Created
July 11, 2013 19:23
-
-
Save stevenworthington/5978441 to your computer and use it in GitHub Desktop.
Example of how to create custom contrasts to test hypotheses in lme4 models.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Note: requires loading the "socsub" data frame (not a bundled R dataset) | |
# ------------------------------------------------------------------------------------ | |
# pairwise comparisons including interactions | |
# use lm model to get design matrix | |
model1 <- lm(agro.rec.tot ~ sex*ageclass + loggrpmem, offset = logtimeage, data = socsub) | |
# get contrasts between sexes for each age class at median loggrpmem | |
cc <- contrast(model1, | |
a = list(ageclass = levels(socsub$ageclass), sex = "F", loggrpmem = median(socsub$loggrpmem)), | |
b = list(ageclass = levels(socsub$ageclass), sex = "M", loggrpmem = median(socsub$loggrpmem)) | |
) | |
# function to append new contrasts to existing contrast matrix | |
append.cc <- function(append, dn, cmat = cc$X) { | |
# get attributes of cc matrix | |
att <- attributes(cmat) | |
# add new contrasts to cc matrix (strips attributes) | |
cmat <- rbind(cmat, append) | |
# change dimensions attributes of cc matrix | |
att$dim[1] <- nrow(cmat) | |
att$dimnames[[1]] <- dn | |
# add the attributes back to the raw matrix | |
attributes(cmat) <- att | |
return(cmat) | |
} | |
# create matrix of new contrasts to append to contrast matrix | |
appMat <- matrix(ncol = 9, byrow = TRUE, data = | |
# females | |
c(0, 0, -1, 0, 0, 0, 0, 0, 0, # Inf II - Inf I == 0 | |
0, 0, 0, -1, 0, 0, 0, 0, 0, # Juv I - Inf I == 0 | |
0, 0, 0, 0, -1, 0, 0, 0, 0, # Juv II - Inf I == 0 | |
0, 0, -1, 1, 0, 0, 0, 0, 0, # Juv I - Inf II == 0 | |
0, 0, -1, 0, 1, 0, 0, 0, 0, # Juv II - Inf II == 0 | |
0, 0, 0, -1, 1, 0, 0, 0, 0, # Juv II - Juv I == 0 | |
# males | |
0, 0, -1, 0, 0, 0, -1, 0, 0, # Inf II - Inf I == 0 | |
0, 0, 0, -1, 0, 0, 0, -1, 0, # Juv I - Inf I == 0 | |
0, 0, 0, 0, -1, 0, 0, 0, -1, # Juv II - Inf I == 0 | |
0, 0, -1, 1, 0, 0, -1, 1, 0, # Juv I - Inf II == 0 | |
0, 0, -1, 0, 1, 0, -1, 0, 1, # Juv II - Inf II == 0 | |
0, 0, 0, -1, 1, 0, 0, -1, 1) # Juv II - Juv I == 0 | |
) | |
# create vector of dimnames that reflect contrasts | |
dnames <- c("Inf I: M - F", | |
"Inf II: M - F", | |
"Juv I: M - F", | |
"Juv II: M - F", | |
"F: Inf II - Inf I", | |
"F: Juv I - Inf I", | |
"F: Juv II - Inf I", | |
"F: Juv I - Inf II", | |
"F: Juv II - Inf II", | |
"F: Juv II - Juv I", | |
"M: Inf II - Inf I", | |
"M: Juv I - Inf I", | |
"M: Juv II - Inf I", | |
"M: Juv I - Inf II", | |
"M: Juv II - Inf II", | |
"M: Juv II - Juv I") | |
# append the new contrasts to the matrix | |
cc2$X <- append.cc(append = appMat, dn = dnames) | |
# fit the same model in lme4 | |
model2 <- glmer(agro.rec.tot ~ sex*ageclass + loggrpmem + (1|id), | |
offset = logtimeage, | |
family = poisson, | |
data = socsub | |
) | |
# plug in new appended contrast matrix | |
summary(glht(model2, linfct = cc2$X), test = adjusted(type = "holm")) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hi Steve
Where can I find the socsub dataset? I'd like to go through this example.
Cheers
Llew