1
$\begingroup$

I know that when running a linear mixed effects model, centering around the grand mean should change the estimates for the coefficients, but not the estimate for the variance.

For example, I have some data in which subjects respond to some yes/no question over 6 time periods (time is perfectly balanced, so group-mean would basically be the same). I am interested in the effect of time over their answer. I run a logistic mixed effect model with and without centering of the variable time.

Without centering:

> m1 = glmer(Yes ~ time + (1+time|id), data = dd, family = "binomial", 
+            control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5)))
> summary(m1)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Yes ~ time + (1 + time | id)
   Data: dd
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

     AIC      BIC   logLik deviance df.resid 
  1516.8   1542.7   -753.4   1506.8     1303 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8536 -0.5075 -0.2720  0.5576  2.0201 

Random effects:
 Groups Name        Variance Std.Dev. Corr 
 id     (Intercept) 4.98084  2.2318        
        time        0.07198  0.2683   -0.37
Number of obs: 1308, groups:  id, 218

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.82498    0.23031  -3.582 0.000341 ***
time         0.10282    0.04683   2.196 0.028125 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
time -0.713

Now with centering:

> dd$c_time = dd$time - mean(dd$time)
> m2 = glmer(Yes ~ c_time + (1+c_time|id), data = dd, family = "binomial", 
+            control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5)))
> summary(m2)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Yes ~ c_time + (1 + c_time | id)
   Data: dd
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

     AIC      BIC   logLik deviance df.resid 
  1516.8   1542.7   -753.4   1506.8     1303 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8536 -0.5076 -0.2719  0.5576  2.0201 

Random effects:
 Groups Name        Variance Std.Dev. Corr
 id     (Intercept) 4.32746  2.0803       
        c_time      0.07198  0.2683   0.06
Number of obs: 1308, groups:  id, 218

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.46512    0.16151  -2.880  0.00398 **
c_time       0.10286    0.04683   2.196  0.02806 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
c_time -0.002

And indeed the estimated variance for the random slope did not change (0.07198 in both). And neither did the estimate for the fixed effect. Good.

But, for my data it makes little sense that the random intercept and slope will correlate. So I also ran the same models with independent random components. Without centering:

> m1 = glmer(Yes ~ time + (1|id) + (0+time|id), data = dd, family = "binomial", 
+            control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5)))
> summary(m1)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Yes ~ time + (1 | id) + (0 + time | id)
   Data: dd
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

     AIC      BIC   logLik deviance df.resid 
  1515.6   1536.3   -753.8   1507.6     1304 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8887 -0.5087 -0.2901  0.5653  2.0188 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 3.65259  1.9112  
 id.1   time        0.05097  0.2258  
Number of obs: 1308, groups:  id, 218

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.77117    0.20871  -3.695  0.00022 ***
time         0.09083    0.04439   2.046  0.04077 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
time -0.648

And with centering:

> m2 = glmer(Yes ~ c_time + (1|id) + (0+c_time|id), data = dd, family = "binomial", 
+            control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5)))
> summary(m2)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Yes ~ c_time + (1 | id) + (0 + c_time | id)
   Data: dd
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

     AIC      BIC   logLik deviance df.resid 
  1514.8   1535.5   -753.4   1506.8     1304 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8525 -0.5093 -0.2686  0.5554  2.0206 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 4.32549  2.0798  
 id.1   c_time      0.07164  0.2677  
Number of obs: 1308, groups:  id, 218

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.46647    0.16131  -2.892  0.00383 **
c_time       0.10497    0.04581   2.292  0.02193 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
       (Intr)
c_time -0.016

And I was surprised to see that the estimates for the variance of time now changed. Moreover, the estimate for the fixed effect also changed. Since I know centering is a complex issue is random coefficient models, I am now worried about the interpretation of these results. Can someone:

  1. Explain intuitively and/or mathematically why there's a change in variance estimates when the model is specified with a diagonal covariance matrix, but not when it is specified with a general matrix?
  2. Help me interpret the meaning of the coefficients in the centered model without correlation between the random components? Particularly, why the change in the estimation of the fixed effect? I thought that when grand-mean centering it should be similar interpretation as when using the raw scores...

Thanks for any help

$\endgroup$
1
$\begingroup$

I would guess that in the second case one of the models has not appropriately converged. You could try fitting with different optimizers and check if you observe the same behavior.

On the other hand, note that changing the scale of your time variable may help the algorithm converge more appropriately. For example, if the estimated variance for the random slopes is very close to zero, this may be because it is truly zero or because in the timescale you have selected there is very little variability between subjects. In the latter, you can get a more stable estimate of this variance by rescaling (e.g., dividing or multiplying your time variable with a suitable constant).

$\endgroup$
  • $\begingroup$ Convergence does not appear to be the problem. I just checked with multiple other optimizers (basically tried all options in here, and those that converged all agreed with the above results. I dont follow your second suggestion. First, I am not sure the estimated variance is that close to zero (do you think it is?). Even still, what do you mean by rescaling the time? Literally multiplying time by a constant (say, 10), then running the same models? (I did, it only decreased var estimates by factor of 100) $\endgroup$ – Cuenco Jan 3 at 14:14

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.