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