[R-sig-ME] [R] Missing interaction effect in binomial GLMM with lmer
Ben Bolker
bolker at ufl.edu
Mon Feb 15 04:15:47 CET 2010
I don't think this is necessarily a mixed model question (as you
pointed out, one gets about the same effect in the plain non-mixed GLM).
I also don't think it's the Hauck-Donner effect (i.e., a problem with
estimated probabilities being too close to 0 or 1 so that Wald tests are
unreliable) -- in any case, this issue applies only to Wald tests (e.g.
Z- and t-statistics reported by summary()), not to likelihood ratio
tests (e.g. results of anova()). Is it not just a reasonable conclusion
that there is not a significant difference between (the difference
between 29/30 and 16/30) and (the difference between 29/30 and 5/30)?
The likelihood ratio test is known to be a little unreliable for small
sample sizes (I can't work out the approximate "denominator df" for this
case without thinking too hard), but a quick parametric bootstrap
suggests that the p-value of 0.24 is not too far off ...
Weber, Sam wrote:
> 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")
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / people.biology.ufl.edu/bolker
GPG key: people.biology.ufl.edu/bolker/benbolker-publickey.asc
More information about the R-sig-mixed-models
mailing list