[R-sig-ME] post hoc tests following lmer
Martin Henry H. Stevens
stevenmh at muohio.edu
Fri May 18 14:38:21 CEST 2007
Hi Henrik,
See the languageR package, and a paper by Baayen, Davidson and Bates
(submitted), at Baayen's web site. This may be helpful, as it gives a
couple more bells-and-whistles.
Also, if I understand you correctly, you have two ages each of female
and male? How about creating no-intercept model, generating Bayesian
credible intervals with mcmcsamp and HPDintervals and then comparing
then comparing the four combinations that way? I would be interested
to here the thoughts of others on this.
Hank
On May 18, 2007, at 7:19 AM, parn at nt.ntnu.no wrote:
> 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
>
> --
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list