# Running ANOVAS in a generalized linear mixed effect model

My experiment (named Exp.2 hereafter) is a 2x2x2 design. In Exp. 2, the data set has three fixed effects predictor variables (threat, distance and modulation) and the dependent variable is a binomial output named Code. The data has two random effects variables named Subject and Group. I'm using the lme4 package to run a binomial linear mixed effects model.

To make my upcoming questions as clear as possible, I'll give a bit of background. I also ran a binomial linear mixed effects model for Exp1, which was a 2x3 design. The full model is:

glmfull <- glmer(Code ~ Threat * Distance + (1|Subject) + (1|Group),
data= exdata, family = "binomial")


From this I saw that whilst Threat and Distance were significant the interaction wasn't. I then built 3 more models:

glm.a <- glmer(Code ~ Threat + (1|Subject) + (1|Group),
data = exdata, family= "binomial")
glm.b <- glmer(Code ~ Distance + (1|Subject) + (1|Group),
data = exdata, family= "binomial")
glm.a.b <- glmer(Code ~ Threat:Distance + (1|Subject) + (1|Group),
data = exdata, family= "binomial")


Then, I ran these anovas:

anova(glm.a,glmfull)
anova(glm.b,glmfull)
anova(glm.a.b,glmfull)


From a, I found that Distance was significant: ($\chi^2(3) = 74.9, p < .001$) From *b**, Threat was significant: ($\chi^2(3) = 25.3, p < .001$) From a.b, the interaction was not significant: p > 0.5

The intercepts I reported, for example to show that difference between the two levels of Distance ($\beta$ = 1.65, SE = 0.2, p < 0.001) were taken from the output I got when running the model of best fit, which was:

 glmer(Code ~ Threat + Distance + (1|Subject) + (1|Group),
data = exdata, family= "binomial")


QUESTION 1: I just want to check it's definitely correct to report $\beta$ values from the model of best fit, rather than the 'full' model?

Moving on to Exp.2 -- essentially, I've worked out that the model of best fit is:

glmbest <- glmer(Code ~ Threat + Distance + Modulation +
Threat:Modulation + (1|Subject) + (1|Group),
data= exdata2, family = "binomial")


Running the full model (Code ~ Threat * Distance * Modulation etc) showed that the three-way interaction wasn't significant, nor was the interaction between Modulation and Distance. The interaction between Threat and Distance was slightly significant (p = 0.013) but I decided to remove it, as it improved model fit.

QUESTION 2 really is this: which models do I include in the ANOVA to get the correct chi-squared statistic to report?

For example, to get the $\chi^2$ for Threat do I run a model where Threat is removed, against the 'full' model, shown in 1[], or against the model of best fit as in 2[]? Or is it neither of the two examples below and something totally different?

## 1.

glm1 <- glmer(Code ~ Distance * Modulation + (1|Subject) + (1|Group),
data=exdata2, family = "binomial")
glmfull1 <- glmer(Code ~ Threat* Distance * Modulation + (1|Subject) + (1|Group),
data=exdata2, family = "binomial")
anova(glm1,glmfull1)


## 2.

glm2 <- glmer(Code ~ Distance+Modulation+Threat:Modulation + (1|Subject) + (1|Group),
data=exdata2, family = "binomial")
anova(glm2,glmbest)


There are so many variations, and I'm not sure which is correct. I hope all this makes sense.

Oh and just a note: Although I know interactions should usually be removed if they aren't accounted for in the hypothesis, the interaction I've kept in Exp.2 (Threat:Modulation) is an interaction I was predicting.

• You might be interested in the lmerTest package. – Roland Sep 19 '16 at 14:11

There are very cogent arguments against stepwise model reduction (see Frank Harrell's Regression Modeling Strategies book); in particular, if you do conditional inference after selecting a model (i.e., you pretend that the selected model was the one you meant to use from the beginning, and that there wasn't a model selection process involved), then all your conclusions will be overconfident (confidence intervals narrower and $p$ values smaller than they should be).

However, it may be OK/not so bad to do stepwise reduction at least to the level of removing non-significant interaction terms (this is done in many otherwise excellent treatments). So in answer to question 1

it's definitely correct to report β values from the model of best fit, rather than the 'full' model?

I would say, "yes; if you've removed some interactions, report the coefficients from the models without those interactions" (reporting main effects from models that also include those effects in an interaction is at best tricky and at worst wrong).

The interaction between Threat and Distance was slightly significant (p = 0.013) but I decided to remove it, as it improved model fit.

?? (1) in my field (ecology and evolution), $p=0.013$ would generally be considered "moderately strong" ... (2) how are you characterizing model fit? Adding interactions (or terms of any sort) will always improve model fit, at least slightly, if you measure goodness-of-fit by some unpenalized measure such as (unadjusted) $R^2$ or likelihood ...

For example, to get the $\chi^2$ for Threat do I run a model where Threat is removed, against the 'full' model, shown in 1[], or against the model of best fit as in 2[]?

To test a hypothesis about a term you should always compare two models that differ only in the term you're interested in testing, so narrowly speaking 2[] is the more correct choice. However, you need to be SUPER careful if/when testing main effects in the presence of interaction. If you're going to violate marginality in this way, you need to be aware that the meaning of your main effect depends on the way you have parameterized the interacting variable. In particular, what you are doing when comparing Distance+Modulation+Threat:Modulation to Distance+Modulation+Threat+Threat:Modulation is testing the effect of Threat when Modulation is equal to its baseline value (i.e. when Modulation=0 if continuous, or when Modulation is equal to the baseline factor value if categorical and using default treatment contrasts). See Schielzeth 2010 "Simple means to improve the interpretability of regression coefficients" ...