[R] Assessing interaction effects in GLMMs

[Ben Bolker]

Luke Duncan <luke.mangaliso.duncan <at> gmail.com> writes:

> 
> Dear R gurus
> 
> I am running a GLMM that looks at whether chimpanzees spend time in shade
> more than sun (response variable 'y': used cbind() on counts in the sun and
> shade) based on the time of day (Time) and the availability of shade
> (Tertile). I've included some random factors too which are the chimpanzee
> in question (Individual) and where they are in a given area (Zone). 

   How many zones are there? It could be a toss-up between treating
Zone as random vs fixed ...

> There
> are also two continuous predictors (Minimum daily temperature: Min; Maximum
> daily temperature: Max). I have run my GLMM and I know that Time and Min
> are significant predictors of the patterns of shade use while Tertile and
> Max are not. In addition, a Time*Tertile interaction effect is a good
> predictor as well.

  OK, be careful (this is a general point about interpreting models
in R, not specific to GLMMs); the main effects of Time and Tertile
are assessed **at the baseline level of the other** when the interaction
is present in the model.  You have to be very careful interpreting
the meaning of main effects in the presence of interactions.  The
Type III Wald test you use below probably only makes sense if you
have set sum-to-zero contrasts (options(contrasts="sum")).

> I now need to assess how the specific interaction effect conditions differ
> to one another. So, for example, how does shade use differ between 10h00 at
> low shade and 10h00 at high shade? I tried using the package multcomp, but
> that will only allow me to work out the contrasts for the first-order
> effects (Time, Tertile) but won't allow me to do so for the interaction
> effects. Any ideas?

  Are you trying to test a hypothesis, or just to ask about the values?
You can construct the values yourself, or (with the development version
of lme4) use the predict() method, to compute what the probabilities
will be in each of these situations.  By the way, the contrast you
have suggested above (10AM, low vs high shade) is not an interaction;
to test it, I would suggest that you relevel() your time variable
to 10 AM, then look at the estimated effect of Tertile for high vs low.
> 
> My code:
> 
> > m1 <- lmer(y ~ Time*Tertile + (1|Individual) + (1|Zone) + Max +
> Min,family=binomial,REML=F)

> > Anova(m1,type=3,test="Wald")
> Analysis of Deviance Table (Type III tests)
> 
> Response: y
>                Chisq Df Pr(>Chisq)
> (Intercept)   0.9511  1     0.3294
> Time         60.7807  4  1.988e-12 ***
> Tertile       0.3391  1     0.5603
> Max           1.3198  1     0.2506
> Min          77.7736  1  < 2.2e-16 ***
> Time:Tertile 38.9038  4  7.292e-08 ***
> ---

> > summary(m1)
> Generalized linear mixed model fit by the Laplace approximation
> Formula: y ~ Time * Tertile + (1 | Individual) + (1 | Zone) + Max + Min
>   AIC  BIC logLik deviance
>  1168 1224 -569.9     1140
> Random effects:
>  Groups     Name        Variance Std.Dev.
>  Zone       (Intercept) 0.81949  0.90526
>  Individual (Intercept) 0.36417  0.60347
> Number of obs: 412, groups: Zone, 8; Individual, 7
> 

[snip]