1
$\begingroup$

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().

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – themeo Nov 2 '14 at 10:23
  • 1
    $\begingroup$ 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. $\endgroup$ – Aaron - Reinstate Monica Nov 3 '14 at 3:08
  • $\begingroup$ 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. $\endgroup$ – themeo Nov 3 '14 at 10:23
  • $\begingroup$ 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. $\endgroup$ – Aaron - Reinstate Monica Nov 3 '14 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.