[R] Missing interaction effect in binomial GLMM with lmer
Ista Zahn
istazahn at gmail.com
Tue Feb 9 18:19:29 CET 2010
Hi Sam,
Good question. I originally guessed that the "simple effect" (I know
some people on this list don't seem to care for that term, but it's
always made sense to me) coefficients were in the same direction, such
that the effect if Origin at Treat=hot was significantly different
from zero, but not from the effect of Origin at Treat = cold. But a
quick look indicated that is not the case:
contrasts(hatch.frame$Treat) <- contr.treatment(2, base=1)
model1<-lmer(success~Origin*Treat+(1|Female),family=binomial,REML=TRUE,data=hatch.frame)
summary(model1)
Generalized linear mixed model fit by the Laplace approximation
Formula: success ~ Origin * Treat + (1 | Female)
Data: hatch.frame
AIC BIC logLik deviance
95.34 109.3 -42.67 85.34
Random effects:
Groups Name Variance Std.Dev.
Female (Intercept) 0.54993 0.74157
Number of obs: 120, groups: Female, 20
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.609227 1.146844 3.147 0.00165 **
Origin2 -0.004192 1.606214 -0.003 0.99792
Treat2 -5.401703 1.238911 -4.360 1.3e-05 ***
Origin2:Treat2 1.948242 1.697945 1.147 0.25121
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Orign2 Treat2
Origin2 -0.714
Treat2 -0.889 0.635
Orign2:Trt2 0.649 -0.907 -0.730
contrasts(hatch.frame$Treat) <- contr.treatment(2, base=2)
model2<-lmer(success~Origin*Treat+(1|Female),family=binomial,REML=TRUE,data=hatch.frame)
summary(model2)
Generalized linear mixed model fit by the Laplace approximation
Formula: success ~ Origin * Treat + (1 | Female)
Data: hatch.frame
AIC BIC logLik deviance
95.34 109.3 -42.67 85.34
Random effects:
Groups Name Variance Std.Dev.
Female (Intercept) 0.54993 0.74157
Number of obs: 120, groups: Female, 20
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.7925 0.5683 -3.154 0.00161 **
Origin2 1.9441 0.7190 2.704 0.00686 **
Treat1 5.4017 1.2389 4.360 1.3e-05 ***
Origin2:Treat1 -1.9484 1.6979 -1.148 0.25116
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) Orign2 Treat1
Origin2 -0.790
Treat1 -0.385 0.305
Orign2:Trt1 0.281 -0.336 -0.730
So I'm as stumped as you are. How can the effect of Origin at
treat=hot be significantly different from zero, but not significantly
different from -0.004? Clearly there is something here I'm not
understanding. I'm very curious to know the answer.
Best,
Ista
On Tue, Feb 9, 2010 at 12:22 PM, Weber, Sam <Sam.Weber at exeter.ac.uk> wrote:
> Dear all,
>
> I was wondering if anyone could help solve a problem of a missing interaction effect!!
>
> I carried out a 2 x 2 factorial experiment to see if eggs from 2 different locations (Origin = 1 or 2) had different hatching success under 2 different incubation schedules (Treat = 1 or 2). Six eggs were taken from 10 females (random = Female) at each location and split between the treatments, giving 30 eggs from each location in each treatment.
>
> Overall proportions hatching were as follows:
>
> Treat
> 1 2
> Origin
> 1 29/30 5/30
> 2 29/30 16/30
>
>
> I made a binomial response in which hatching was a success and not-hatching was a failure, and analysed as a binomial GLMM. I'm particularly interested in the interaction between the two factors. An expression reproducing the raw data is attached at the end of the post in case it is helpful.
>
> hatch.frame$success<-cbind(hatch.frame$Hatched,hatch.frame$Nothatched)
> model<-lmer(success~Origin*Treat+(1|Female),family=binomial,method="ML",data=hatch.frame)
> model2<-update(model,~.-Origin:Treat)
> anova(model,model2)
>
> Data:
> Models:
> model2: success ~ Origin + Treat + (1 | Female)
> model: success ~ Origin * Treat + (1 | Female)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> model2 4 94.707 105.857 -43.353
> model 5 95.350 109.287 -42.675 1.3572 1 0.244
>
> model3<-update(model2,~.-Origin)
> anova(model2,model3)
>
> Data:
> Models:
> model3: success ~ Treat + (1 | Female)
> model2: success ~ Origin + Treat + (1 | Female)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> model3 3 98.863 107.225 -46.431
> model2 4 94.707 105.857 -43.353 6.1558 1 0.01310 *
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> model4<-update(model2,~.-Treat)
> anova(model2,model4)
>
> Data:
> Models:
> model4: success ~ Origin + (1 | Female)
> model2: success ~ Origin + Treat + (1 | Female)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> model4 3 155.592 163.954 -74.796
> model2 4 94.707 105.857 -43.353 62.885 1 2.191e-15 ***
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> So the model implies that there is a very significant effect of treatment (reduced hatching at treatment 2) with a small effect of origin (improved hatching from origin 2). However the lack of interaction effect implies hatching was better for Origin 2 at both treatments, which looking at the raw values above does not seem to be the case. Identical numbers of eggs hatched from both Origins in Treatment 1, but much more from Origin 2 hatched at Treatment 2.
>
> If you divide the analysis by treatments, Origin only has a significant effect on hatching under Treatment 2 and not with Treatment 1
>
> Hot<-data.frame(hatch.frame[hatch.frame$Treat==2,])
> Cold<-data.frame(hatch.frame[hatch.frame$Treat==1,])
>
> #2
> model<-lmer(success~Origin+(1|Female),family=binomial,method="ML",data=Hot)
> model2<-update(model,~.-Origin)
> anova(model,model2)
> Data: Hot
> Models:
> model2: incubate ~ (1 | Code)
> model: incubate ~ Origin + (1 | Code)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> model2 2 78.633 82.821 -37.316
> model 3 73.697 79.980 -33.848 6.9357 1 0.008449 **
>
>
>
> #1
>
> model<-lmer(success~Origin+(1|Female),family=binomial,method="ML",data=Cold)
> model2<-update(model,~.-Origin)
> anova(model,model2)
>
> Data: Cold
> Models:
> model2: incubate ~ (1 | Code)
> model: incubate ~ Origin + (1 | Code)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> model2 2 21.5086 25.6973 -8.7543
> model 3 23.3472 29.6302 -8.6736 0.1615 1 0.6878
>
>
>
>
>
> So I can't understand where the interaction effect has gone in the full model?! I get the same result in a binomial GLM, without the random effect of Female i.e. a small effect of origin but no interaction with treatment. I'm sure I must be missing something here so I would be very grateful to anyone who can point out my mistakes. I've read previous R Help posts that suggest binomial GLM(M) can create problems when estimated probabilities are close to 0 or 1. In Treatment 1 hatching probability was 0.97 for both Origins, so could this be the source of the problem?
>
>
>
> Thanks for your help
>
>
>
> Sam Weber
>
>
> ----------------------------------------------------------------------------
> University of Exeter
> Centre for Ecology and Conservation
> Tremough Campus
> Penryn
> Cornwall TR10 9EZ
> UK
>
>
>
>
> hatch.frame <-
> structure(list(Female = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L,
> 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L,
> 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L,
> 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 10L,
> 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L,
> 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L,
> 14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L,
> 16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L,
> 18L, 19L, 19L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 20L, 20L
> ), .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10",
> "11", "12", "13", "14", "15", "16", "17", "18", "19", "20"), class = "factor"),
> Origin = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
> 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
> 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
> 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("1", "2"), class = "factor"),
> Treat = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
> 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
> 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
> 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
> 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
> 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
> 2L, 1L, 2L, 1L, 2L), .Label = c("1", "2"), class = "factor"),
> Hatched = c(1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L,
> 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L,
> 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L,
> 1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L,
> 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L,
> 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L,
> 0L, 1L, 1L), Nothatched = c(0, 0, 1, 0, 0, 1, 0, 1, 0, 1,
> 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
> 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
> 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0,
> 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
> 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,
> 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0)), .Names = c("Female",
> "Origin", "Treat", "Hatched", "Nothatched"), row.names = c(NA,
> -120L), class = "data.frame")
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Ista Zahn
Graduate student
University of Rochester
Department of Clinical and Social Psychology
http://yourpsyche.org
More information about the R-help
mailing list