Assessing interaction effects in GLMMs
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]