[R-sig-ME] post hoc tests following lmer
parn at nt.ntnu.no
parn at nt.ntnu.no
Fri May 18 13:19:04 CEST 2007
I have struggled to find a good way to perform a post hoc test on the
fixed factors in a mixed model using lmer. I found some threads on the
...but I am really curious if anyone have any additional thoughts on
I provide you some dummy data of 200 bird nests.
The male and female attending the nest could be either young or old.
Each chick in a clutch could either have a disease or not.
Some individuals occur in the data set several times and 'identity' of a
bird is included as a random variable.
I would like to test if the proportion of chicks that carries the
disease in a clutch depend on age of the male and female, and the
And (hopefully) be able to perform a posthoc test of which
age-combination of the parents that differs in proportion of sick
# Dummy data
# male age
male.age <- rep(c("Y", "O"), c(100, 100))
# female age
female.age <- rep(c("Y", "O", "Y", "O"), c(50, 50, 50, 50))
# male id
# 100 unique young males
male.id.y <- (1:100)
# 70 unique old males, and 30 of the old captured before
male.id.o <- sample(c(sample(male.id.y, 30, replace=F), 101:170))
male.id <- c(male.id.y, male.id.o)
# female id
female.id.y <- 1:100
female.id.o <- sample(c(sample(female.id.y, 40, replace=F), 101:160))
female.id <- c(female.id.y[1:50], female.id.o[1:50],
# clutch size
clutch <- floor(rnorm(200, mean=9, sd=1.35))
# number of chicks with disease in pairs with of a young male (ym) and a
# young female (yf):
n.epo.ym.yf <- rbinom(50, clutch[1:50], p=0.09)
# and so on (different probabilities of disease in each combination of
male and female age.
n.epo.ym.of <- rbinom(50, clutch[51:100], p=0.25)
n.epo.om.yf <- rbinom(50, clutch[101:150], p=0.075)
n.epo.om.of <- rbinom(50, clutch[151:200], p=0.06)
n.epo <- c(n.epo.ym.yf, n.epo.ym.of, n.epo.om.yf, n.epo.om.of)
# number of healthy offspring
n.wpo <- clutch - n.epo
ep.data <- data.frame(male.id, female.id, male.age, female.age, clutch,
# create response variable
# two-column matrix with the columns giving the numbers of 'successes'
# (n.epo) and 'failures' (n.wpo)
y <- cbind(n.epo, n.wpo)
# plot proportion of epo in the clutch for the
# different age combinations in pairs
p.epo <- n.epo/clutch
pairage <- as.factor(paste(male.age, female.age, sep=""))
plot(p.epo ~ pairage)
# mixed model.
# response: proportion epo in each clutch: y
# fixed factors: male age, female age and the interaction thereof:
# random factors: male.id, female.id
model <- lmer(y ~ male.age*female.age + (1|male.id) + (1|female.id),
But what would be the appropriate way to make a post hoc test of which
age combinations that differs, an lmer-equivalent of
'multicomp', perhaps an mcmcsamp-approach!?
Thanks in advance for any comments!
More information about the R-sig-mixed-models