[R-sig-ME] post hoc tests following lmer

[parn at nt.ntnu.no]

Dear lmer-users,

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
R-list:
http://www.nabble.com/investigating-interactions-with-mixed-models-tf3272509.html#a9099081
http://www.nabble.com/How-to-use-lmer-function-and-multicomp-package--tf1730665.html#a4702360

...but I am really curious if anyone have any additional thoughts on
this subject.

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
interaction thereof.
And (hopefully) be able to perform a posthoc test of which
age-combination of the parents that differs in proportion of sick
chicks.


# 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],
female.id.y[51:100], female.id.o[51:100])

# 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,
n.epo, n.wpo)


# 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:
# male.age*female.age
# random factors: male.id, female.id


model <- lmer(y ~ male.age*female.age + (1|male.id) + (1|female.id),
family=binomial, data=ep.data)
summary(model)


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!

Best regards,

Henrik

--