# Coding of categorical random effects in R: int vs factor

I have a problem with coding of a 2-level categorical predictor variable in R, and subsequently using it as a random slope in lmer().

I can keep the factor as numeric, coded using the treatment coding:

> unique (b$multi) [1] 0 1  Running lmer() using a dataset coded in this way yields: > l1 = glmer(OK ~ multi + (0 + multi|item) + (1|subject)+ (1|item), family="binomial", data=b) > summary(l1) Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod] Family: binomial ( logit ) Formula: OK ~ multi + (0 + multi | item) + (1 | subject) + (1 | item) Data: b AIC BIC logLik deviance df.resid 4806.5 4838.9 -2398.3 4796.5 4792 Scaled residuals: Min 1Q Median 3Q Max -7.8294 -0.5560 -0.1548 0.5623 14.3342 Random effects: Groups Name Variance Std.Dev. subject (Intercept) 1.84379 1.3579 item (Intercept) 2.40306 1.5502 item.1 multi 0.04145 0.2036 Number of obs: 4797, groups: subject, 123; item, 39 [...]  Above there is only one random slope related to multi. However, something very different happens when I convert the variable into a factor: > b$multi = as.factor(b$multi) > levels (b$multi)
[1] "0" "1"


When I fit a model using multi as a random slope variable:

l2 = glmer(OK ~ multi + (0+multi|item) + (1|subject)+ (1|item), family="binomial", data=b) Warning message: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 1 negative eigenvalues

... the model fails to converge and I get a very different random effects structure:

> summary(l2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: OK ~ multi + (0 + multi | item) + (1 | subject) + (1 | item)
Data: b

AIC      BIC   logLik deviance df.resid
4807.8   4853.1  -2396.9   4793.8     4790

Scaled residuals:
Min      1Q  Median      3Q     Max
-8.3636 -0.5608 -0.1540  0.5627 15.2515

Random effects:
Groups  Name        Variance Std.Dev. Corr
subject (Intercept) 1.8375   1.3555
item    (Intercept) 0.9659   0.9828
item.1  multi0      1.5973   1.2638
multi1      1.0224   1.0111   1.00
Number of obs: 4797, groups:  subject, 123; item, 39
[...]


The number of parameters in the model clearly change (reflected by the change in AIC, etc.), and I get two random slopes.

My question is which way of coding the categorical variable is better? Intuition tells me that it is the first one, but I have seen recommendations for both ways of coding in various tutorials and classes about running GLMMs in R and this is why it baffles me. Both types of the predictor variable work identically in ordinary regression using lm().

The difference is because you're separating the intercept and the slope in the random effect. That's an odd thing to do; the usual way to fit this model would be

OK ~ multi + (multi | item) + (1 | subject)


with multi being a factor.

What happens is that in the first model you get what you expect; the 0+multi|item term gives one parameter and the 1|item term gives one parameter, but in the second model the 0 + multi | item term results in two parameters, which are simply the estimate for each condition. If you take the 1|item term out of that model you should get a result that is equivalent to both your first model and the one I give above, except for differences in parameterization.

Note also the correlation of exactly one in your second model; this is a clue that you've overparameterized it and that one of those parameters is not necessary.

• Thanks, stating the model as (multi|item) or (1+multi|item) indeed made the two versions of coding the variable work in the same way. The reason I specified the model as (0+multi|item) + (1|item) is that it removes the correlation term between the random intercept and the random slope - a parameter that increases model complexity but does not help in better controlling type I error rate. So my take-home message is that if I want to keep the correlation parameter out, I should recode the variable in the first version, i.e. have the contrasts defined as integers, not as factors. – themeo Nov 2 '14 at 10:23
• Well... I'm a little suspect of intentionally not including the correlation term. I think you need it, but I haven't fully thought it through for this model. So don't take it out simply because it makes the model more complex (it's not that complex anyway); take it out instead only if that makes the model match the design of your study better. – Aaron left Stack Overflow Nov 3 '14 at 3:08
• That is at least what these guys found out: Random effects structure for confirmatory hypothesis testing: Keep it maximal . Through Monte-Carlo modeling they showed that models with correlation term removed retain high power and near-nominal p-values of the full model. My model BTW is much more complex that that presented above, I just stripped it to demonstrate the problem clearly. – themeo Nov 3 '14 at 10:23
• Looks like an interesting article, thanks. I just read the abstract so far, and will simply note that they emphasize that the model should be "justified by the design," so in that respect, we completely agree. – Aaron left Stack Overflow Nov 3 '14 at 15:32