Solved – Random effect on the intercept in models with categorical predictor

categorical datamixed model

I'm doing a mixed model analysis of a model with a categorical predictor (using glmer, binomial model). To my surprise, the variance I obtain for the intercept does not change when I change the reference level of this categorical predictor. This goes against my understanding that as the intercept is the log-odd for the reference level of the categorical predictor, then the variance for the random intercept should change when changing the reference level.

What am I missing here ? Many thanks for your assistance.

Update: @charles Many thanks for your reply. I still don't get it. The following code, in which I simulate data similar to those that I'm analyzing, might help pursuing the discussion.

require(lme4)
set.seed(123)

# Create 6 groups of 100 observations each
group <- gl(6,100)
# Create a categorical predictor
category <- as.factor(round(runif(600,1,5))) # unbalanced design per group
# Create a binary response variable
used <- runif(600)
 used[1:100] <- used[1:100]+runif(100,0.1,0.3) # to increase the number of 1s in the first category, just to create a bit of variability
 used[101:200] <- used[101:200]-runif(100,0.1,0.3) # to decrease the number of 1s in the second category, just to create a bit of variability 
used <- round(used)
used[used>1] <- 1; used[used<0] <- 0

# Fit two models with different reference categories
mod.ref1 <- glmer(used~category+(1|group),family="binomial")
category <- relevel(category,"3") # change reference category
mod.ref3 <- glmer(used~category+(1|group),family="binomial")

# Compare coefficients and random effects for both models
coef(mod.ref1); ranef(mod.ref1)
coef(mod.ref3); ranef(mod.ref3)

And the results are:

coef(mod.ref1)
$group
  (Intercept)   category2  category3  category4   category5
1  0.72458362 -0.02654181 -0.2808615 0.02515119 -0.09999195
2 -0.57971784 -0.02654181 -0.2808615 0.02515119 -0.09999195
3  0.14790658 -0.02654181 -0.2808615 0.02515119 -0.09999195
4  0.14356195 -0.02654181 -0.2808615 0.02515119 -0.09999195
5 -0.02681879 -0.02654181 -0.2808615 0.02515119 -0.09999195
6  0.11810340 -0.02654181 -0.2808615 0.02515119 -0.09999195

ranef(mod.ref1)
$group
  (Intercept)
1  0.63660977
2 -0.66769169
3  0.05993273
4  0.05558810
5 -0.11479264
6  0.03012955

coef(mod.ref3)
$group
  (Intercept) category1 category2 category4 category5
1   0.4436543 0.2808285 0.2543196 0.3060796 0.1809765
2  -0.8605623 0.2808285 0.2543196 0.3060796 0.1809765
3  -0.1329837 0.2808285 0.2543196 0.3060796 0.1809765
4  -0.1373277 0.2808285 0.2543196 0.3060796 0.1809765
5  -0.3076979 0.2808285 0.2543196 0.3060796 0.1809765
6  -0.1627864 0.2808285 0.2543196 0.3060796 0.1809765

ranef(mod.ref3)
$group
   (Intercept)
1  0.63656425
2 -0.66765237
3  0.05992621
4  0.05558221
5 -0.11478796
6  0.03012352

and the variance of the random intercepts is 

Random effects:
 Groups Name        Variance Std.Dev.
 group  (Intercept) 0.18119  0.42567 
Number of obs: 600, groups: group, 6

for both models.

So, (1) variance for the random intercepts is the same in both models, (2) random intercepts appear to be the same in both models, (3) intercepts vary between models as they refer to the reference category. I don't understand how (2) and (3) can be true simultaneously, and I'm wondering if I'm misunderstanding the meaning of random intercepts when categorical predictors come into play. Could anyone put down the equation underlying this model ?

Best Answer

Brief: You can't say the random intercepts are the same in both models. They aren't estimated. If they were they would be intercepts rather than random intercepts. On the other hand the variance may be the same, since this is estimated.

Longer: I think I need to know Latex and have more time to give a reasonable answer. But hopefully someone can build on this.

The univariate random effects model with a single random intercept is the following, where $b_j$ is the random intercept and $\epsilon$ is the error term:

$$ y_{ij}=\beta_{0}+\beta_{1}x1_{ij}+b_{j}+\epsilon_{ij}. $$

If $x1$ is a dummy variable and if $x2=0$, the intercept is: $$ \beta_{0} + b_j $$

If $x2=1$, the intercept is: $$ \beta_{0} + \beta_{1} + b_j $$ Note that changing the coding has no effect on $b_j$, the random intercept.

The random intercept can be considered a latent variable that is not estimated along with the fixed parameters $\beta_1, \beta_2, \dots$, but whose variance is estimated together with the variance of the error term.

Practically speaking, this is why it is such a pain to get predictions for individuals out of mixed-effects models. What is a BLUP or EBLUP?

Should i understand that including a random intercept will allow me to control for variability between groups in ONLY the reference category of the predictor, and that I need to include random 'slopes' to control for variability between groups in ALL categories of the predictor?

Random intercept should be sufficient I think.

  • Random intercept model: the overall level of the response is allowed to vary between clusters after controlling for covariates (it does not only apply to the reference group).
  • Random slope/coefficient model: often used with longitudinal data. In addition to the above, this allows the effects of the covariates to vary between clusters.
Related Question