[R-sig-ME] [R] Missing interaction effect in binomial GLMM with lmer
Weber, Sam
Sam.Weber at exeter.ac.uk
Sun Feb 14 18:54:03 CET 2010
Dear all,
I posted the following query to the general R-help list last week, but following the recommendation below I thought it might be worth posting to the mixed-models list as well to see if anyone could give me any pointers. The original post and data are included below. Would be very grateful for any suggestions,
Best regards
Sam Weber
________________________________________
From: Bert Gunter [gunter.berton at gene.com]
Sent: 09 February 2010 18:24
To: 'Ista Zahn'; Weber, Sam
Cc: 'r-help at R-project.org'
Subject: RE: [R] Missing interaction effect in binomial GLMM with lmer
You might do better posting this on the R-sig-mixed-models list.
Bert Gunter
Genentech Nonclinical Statistics
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Ista Zahn
Sent: Tuesday, February 09, 2010 9:19 AM
To: Weber, Sam
Cc: r-help at R-project.org
Subject: Re: [R] Missing interaction effect in binomial GLMM with lmer
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")
More information about the R-sig-mixed-models
mailing list