Varying group coefficients in lme4

Question

All,

I am estimating a multilevel logistic regression with group predictors, but am unclear about some of the advice given by Gelman and Hill (2007) in their book. Therein, they recommend allowing every coefficient to possibly vary, given a large enough N. Does that include group predictors as well? They weren't clear, treating "varying slope" as just another complexity you can incorporate into a mixed effects model in lme4 along with group predictors (see: p. 549 in their book).

For example, I have roughly 50,000 observations with a binary response (plenty large N). Predictors exist at two levels, such that my model looks like:

M1 <- lmer(Y ~ X1 + X2 + X3 + X4 + G1 + G2 + G3 + G4 + (1 | group), family=binomial(link="logit"))

X1:X4 are individual-level predictors and G1:G4 are group-level predictors, thus: a multilevel model. Does their recommendation of treating all coefficients as potentially variable mean including even the group predictors within the random effect, such that:

M2 <- lmer(Y ~ X1 + X2 + X3 + X4 + G1 + G2 + G3 + G4 + (1 + X1 + X2 + X3 + X4 + G1 + G2 + G3 + G4 | group), family=binomial(link="logit"))

I ran M2 and it gave sensible estimates. AIC/BIC suggest much better fit than M1. I'm just unsure if it's appropriate since, unlike individual-level predictors, the group-level predictors are not going to vary in a given group. It will obviously vary across groups, though.

Further, if this is not an incorrect way to approach it, how suspicious should I be if one of the group predictors of interest is statistically insignificant as a stand-alone fixed effect (varying intercept model like M1), but is significant as a fixed effect in a varying slope model like M2?

Thanks for any input and feedback on this topic. I really appreciate it.

Answer 1

First of all, AIC/BIC do not make sense in mixed models. I mean, if you can explain what your $n$ that goes into your BIC is (number of groups? number of observations? something in between? how about level 1 and level 2 variables that obviously have different amount information in them?)... So I wouldn't pay any attention to these.

Second, I am surprised your model with the random effects for group level variables was identified at all. Let us think about an extreme case: a binary group level variable in the a model lmer(Y ~ X + G + (1 + X + G| group). What is it that it describes? That a group has an additional random shift when G==1, i.e., group-level heteroskedasticity. So that appears to be something rather odd to estimate.

So all in all, I would run this as

M2 <- lmer(Y ~ X1 + X2 + X3 + X4 + G1 + G2 + G3 + G4 + (1 + X1 + X2 + X3 + X4 | group), family=binomial(link="logit"))

i.e., only with an individual level covariates having random effects assigned to them.

Answer 2

Its better, I think, that you write your model first, then we can get how to estimate it in R.

Do you wanna the slope of your model to vary by group? And do you think that the second level predictors (ie. predictors at the group level) can be used to predict the varying slopes? If both answers are yes, here how (I think) you can estimate it in R (i'm not 100% sure about this. I just take a look at the book, at chapter 14):

Assuming only one predictor at level one (for the sake of simplicity) and one predictor at level two, and that the second level is given by group:

Model:

$Pr(y_{i}=1) = a + b_{j}*x_{1}$

$b_{j} = \gamma_{0} + \gamma_{1}*g_{1} + e_{j}$

In r:

glmer(y ~ x1 + (x1+0|group) + g1, family=binomial)

As you can see, the difference from your code is that g1 doesn't go within the random effects. You have to expand them to the individual data level to fit the model.

It happens that in cases like this, it's probably better to use a fully Baysian modeling (using WinBugs, Jags or alike), because in this setting I think you can't model properly the second level.

I hope it helps.