# Grand-mean centering in GLMM changes estimates for variance (and everything else)?

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

## 1 Answer

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

• 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) – Cuenco Jan 3 '19 at 14:14