Can I exclude a fixed factor from a glmer to test the main effect of different continuous variable?

Question

I'm unsure how to proceed to report my effects.

I have a glmer to run a logistic regression:

glmres.RU.conf.full = glmer(correct ~ choiceConf.scaled * condition + (1 + choiceConf.scaled|subject)  + (1|cityPair), data = myDataForglm.RU, family = binomial)

I'm interested in two things: (1) the interaction between choiceConf.scaled * condition and (2) the main effect of choiceConf.scaled.

Testing my interaction effect doing a likelihood ratio test is straightforward:

glmres.RU.conf.sum  = glmer(correct ~ choiceConf.scaled + condition + (1 + choiceConf.scaled|subject)  + (1|cityPair), data = myDataForglm.RU, family = binomial)
anova(glmres.RU.conf.full, glmres.RU.conf.sum)

                    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
glmres.RU.conf.sum   7 4467.4 4513.5 -2226.7   4453.4                         
glmres.RU.conf.full  8 4469.3 4522.0 -2226.7   4453.3 0.0606      1     0.8056

I am unsure about what the correct null model to test the main effect of choiceConf.scaled is. I can think of two approaches. First, I could compare a model including the two fixed effects vs. one including only the fixed factor I'm not interested in.

glmres.RU.conf.sum  = glmer(correct ~ choiceConf.scaled + condition + (1 + choiceConf.scaled|subject)  + (1|cityPair), data = myDataForglm.RU, family = binomial)
glmres.RU.null      = glmer(correct ~ condition                     + (1 |subject)                     + (1|cityPair), data = myDataForglm.RU, family = binomial)

               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
glmres.RU.null      4 4471.8 4498.2 -2231.9   4463.8                           
glmres.RU.conf.sum  7 4467.4 4513.5 -2226.7   4453.4 10.443      3    0.01515 *

Note that the AIC and BIC go in different directions. For the BIC, the null model is better (presumably because it's penalizing the 3 Df difference too strongly). So I tried to find a way to compare models that don't differ so much in complexity. I could do this by never considering the effect of condition in my model at all (I'm not interested in it for now, after all).

glmres.confOnly.full     = glmer(correct ~ choiceConf.scaled            + (1 + choiceConf.scaled|subject)  + (1|cityPair), data = myDataForglm.RU, family = binomial)
glmres.confOnly.null     = glmer(correct ~ 1                            + (1 |subject)                     + (1|cityPair), data = myDataForglm.RU, family = binomial)


anova(glmres.confOnly.full, glmres.confOnly.null)

                 Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
glmres.confOnly.null  3 5178.4 5198.2 -2586.2   5172.4                            
glmres.confOnly.full  6 5082.8 5122.4 -2535.4   5070.8 101.6      3  < 2.2e-16 ***

But this sounds incorrect to me. Is it ok to exclude factors (that are known to have a significant effect) from a model? Given that the BIC values are so different between the first and the second approach, I suspect it's not correct.

Any other comments welcome.

Answer 1

This really isn't specific to GLMMs, or even to GLMs. It's a basic issue with model testing in any kind of models descended from linear models (i.e. [G]LM[M]s ...) Ideally, you wouldn't use AIC/BIC to "test the effect of a parameter" - while the answers often go in similar directions, they are not designed for hypothesis testing.

Let's take a look at an example with a strong interaction:

beta <- c(0,1,0,-2)
set.seed(101)
dd <- data.frame(x=rnorm(100),
                 y=factor(sample(c("a","b"),size=100,replace=TRUE)))
dd$z <- rnorm(100,mean=model.matrix(~x*y,dd) %*% beta,sd=0.2)
with(dd,plot(z~x,col=as.numeric(y)))
legend("bottomright",legend=paste0("y=",c("a","b")),col=1:2,pch=1)

The question is: what do you mean when you say "test the main effect"? Whether the main effect of $x$ is significant, and in which direction, depends entirely on what value of $y$ you decide on. If you use R's defaults, you are implicitly testing the effect of x at the baseline level of y ("a" in this case)

## helper function
pp <- function(...) {
    cc <- coef(summary(lm(z~x*y,dd,...)))
    printCoefmat(cc,digits=3)
}
pp()    
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.0249     0.0236    1.05     0.29    
## x             1.0242     0.0255   40.23   <2e-16 ***
## yb           -0.0338     0.0373   -0.91     0.37    
## x:yb         -2.0045     0.0401  -49.99   <2e-16 ***

So in this case, x is strongly significant (slope=1.02 with SE of 0.0255).

If we change to "SAS-style" contrasts where the last level ($y$="b") is the baseline, we get a strongly negative effect (-0.98 with SE of SE 0.03).

pp(contrasts=list(y=contr.SAS))
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.00889    0.02892   -0.31     0.76    
## x           -0.98022    0.03098  -31.64   <2e-16 ***
## y1           0.03380    0.03734    0.91     0.37    
## x:y1         2.00445    0.04010   49.99   <2e-16 ***

Or, if we use sum-to-zero contrasts, we get the average of the slopes when y=a and y=b - about zero (slope = 0.02, SE = 0.02)

pp(contrasts=list(y=contr.sum))
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.00801    0.01867    0.43     0.67    
## x            0.02201    0.02005    1.10     0.28    
## y1           0.01690    0.01867    0.91     0.37    
## x:y1         1.00223    0.02005   49.99   <2e-16 ***

The question here is whether "halfway between conditions" means anything at all; if my levels were (e.g.) "male" and "female", is looking at the halfway point sensible?

Another reasonable choice would be to test the combination of the main effect and its interaction (i.e. the total effect of x) by comparing models with and without x:

anova(lm(z~x*y,dd),lm(z~y,dd))
## Analysis of Variance Table
## 
## Model 1: z ~ x * y
## Model 2: z ~ y
##   Res.Df    RSS Df Sum of Sq    F    Pr(>F)    
## 1     96  3.208                                
## 2     98 90.748 -2   -87.541 1310 < 2.2e-16 ***

The bottom line is that when you want to test main effects in the presence of interactions you have to decide what you mean. This is why some old-school statisticians start to foam at the mouth and mutter about the principle of marginality when this topic comes up ...