# Fitting binomial regression model in R - correct formula, significance testing, and over-dispersion

I'm using generalized linear models to test for the effect of various predictors on some binomial data. My response is a binomial vector of successes and non-successes.

I want to test whether my categorical predictor of interest (P1) is a significant predictor of my response. However, there's another categorical variable (P2) that I want to take account of. Therefore I included it as random effect in my model, which seemed to have not much impact on whether P1 was significant or not:

mod1<-glmer(Response~P1+(1|P2),family=binomial)
mod0<-glmer(Response~1+(1|P2),family=binomial)
anova(mod0,mod1,test="Chisq")


However, reviewers have now asked me for more information on the effect of P2 - they think it will be important and want to know if it's significant or not. If I just test for an effect of P2 as a fixed effect on its own it would be significant, but I think this is just because there is come correlation with P1, not because it is important. Therefore I was thinking I could test whether P2 is a significant predictor within each level of P1, and tried to do this like so:

mod2<-glmer(Response~P2+(1|P1),family=binomial)
mod3<-glmer(Response~1+(1|P1),family=binomial)
anova(mod2,mod3,test="Chisq")


However, I'm worried that the formulas I'm using above are wrong, since I read this guidance on R mixed model formulas:

"Random effects are specified as e|g, where e is an effect and g is a grouping factor"


So should I actually be doing this:

mod4<-glmer(Response~P2|P1,family=binomial)


And if so, how do I test for significance?

As an additional question, my data is heavily overdispersed, but I notice I can't use quasibinomial in lme4. I've read about including an observation level random effect. Using my example above, is this the correct way to do so?

mod5<-glmer(Response~P2+(1|P1)+(1|obs),family=binomial)
mod6<-glmer(Response~1+(1|P1)+(1|obs),family=binomial)
anova(mod5,mod6,test="Chisq")


This doesn't seem to have any obvious affect on the results, so wasn't sure I was doing it right.

Many thanks

There are a variety of issues here:

• hopefully P2 has enough levels (e.g. more than 5 or 6) that it's reasonable to treat it as a random effect; otherwise it's quite likely that the variability associated to P2 will be set to zero (just check VarCorr(mod1)).
• if not, it probably makes sense to treat P2 as fixed, i.e. use
mod2 <- glm(Response~P1+P2,family=binomial)
drop1(mod2,test="Chisq")


or

mod1 <- glm(Response~P1,family=binomial)
anova(mod2,mod1,test="Chisq")


You could use Response~P2|P1, but that actually tests for an interaction - differences in the effect of P2 at different levels of P1 - which is not what you said you wanted. (Also, the random effect is zero-centered, so you should probably add a main (fixed) effect of P2 to the model.)

You could also test for a significant effect of P2 considered as a random effect by comparing glm(Response~P1,...) to glmer(Response~P1+(1|P2),...), but there are some issues about likelihood ratio tests of random effects (e.g. see http://glmm.wikidot.com/faq); the naive test will be conservative.

Note that "nuisance variables" (variables you are not directly interested in) and "random effects" are not synonymous; see e.g. this answer. I wouldn't suggest switching the representation of a variable from fixed to random or vice versa depending on the particular inferences you're trying to make.

Your approach to dealing with overdispersion looks reasonable, although see this PeerJ article for some perspective ... If the observation-level variance was estimated as zero (again see the results of VarCorr(model)), that would explain why the results didn't change ...

• Thanks for the response! Since I posted the question it occurred to me that maybe I should just treat them both as fixed effects, and your answer gives me confidence that this would be the best route to go down.
– rw2
Commented Dec 21, 2015 at 14:23

There was one thing in your description of what you wanted to do that caught my interest, you wrote:

Therefore I was thinking I could test whether P2 is a significant predictor within each level of P1,

That idea sound very much like P1 and P2 are really fixed effects. Perhaps I'm over-interpreting your choice of terms.

Anyway, if the levels of P1 and P2 are interesting in themselves, and you are not only interested in controling for the nuisance of P2, then I suggest, like Ben Bolker already did, that you treat them both as fixed effects. From the glmm-faq:

1. Effects are fixed if they are interesting in themselves or random if there is interest in the underlying population. Searle, Casella and McCulloch [(1992), Section 1.4] explore this distinction in depth. 3. “When a sample exhausts the population, the corresponding variable is fixed; when the sample is a small (i.e., negligible) part of the population the corresponding variable is random” [Green and Tukey (1960)].

If you go this route, then perhaps it would also be illuminating to incorporate an interaction-term in your model.

fm1 <- glm(Response ~ P1 * P2, family = binomial)


Or with some family that accounts for the overdispersion.

• Thanks, I think you're right - both should be fixed effects.
– rw2
Commented Dec 21, 2015 at 14:23