You're conflating two slightly different concepts. In short, the output of VarCorr()
is based on the estimated variance-covariance matrix of the unconditional distribution of the random effects, and is not dependent on the observed data $\mathbf{\mathcal{Y}}$. The correlations that you calculated manually are based on the conditional modes of the random effects, which are where the density of $\mathbf{\mathcal{Y}}$ is maximal given the random effects, which means that they are conditional on the observed data (paraphrased from here).
The correlation output of VarCorr()
(what you see in the output of print(model)
) is based on the estimate of $\mathbf{\Sigma}$ for the marginal, or unconditional, distribution of the random effects, given by
$\mathbf{\mathcal{B}} \backsim \mathcal{N}(\mathbf{0},\mathbf{\Sigma})$
Where $\mathbf{\mathcal{B}}$ is the multivariate normal vector of random effects, with mean 0 and q x q variance-covariance matrix $\mathbf{\Sigma}$, and q is the number of columns in the random effects model matrix (see the lme4 paper for details, and you can check in lme4
with getME(model, "q")
).
The conditional modes are where the density of the conditional distribution of $(\mathbf{\mathcal{Y}}|\mathbf{\mathcal{B}} = \mathbf{\mathcal{b}})$ takes the highest value. They are calculated for each level of the grouping factor of the random effect, and so the correlation of these can change depending on what data you've observed.
As for figuring out the manual calculation, see fortunes::fortune(250)
. The relevant code is in nlme:::VarCorr.pdMat
or lme4:::mkVarCorr
(I think the lme4 code is a little easier to read). I think the way lme4
does it is:
Takes $\sigma$ (the residual standard deviation/scale factor), random effect component names, number of terms for each component, and $\theta$ (the random effects parameters as Cholesky factors) as input.
For each random effects term $\mathcal{i}$,
- Takes the cross product of the scale factor and the transpose of $\Lambda_i$ (which is dependent on $\theta$), which gives the variance-covariance matrix
- Calculates the standard deviation as the square root of the diagonal of the vcov matrix
- Calculates the correlation as covariance (non-diagonal terms) divided by the standard deviation
A (not very flexible) version of how mkVarCorr
calculates the correlation:
library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), data = sleepstudy, REML = FALSE)
getRanCorr <- function(mod) {
n_ranefs <- ncol(ranef(mod)[[1]])
theta <- getME(mod, "theta")
sig <- sigma(mod)
Li <- diag(nrow = n_ranefs)
Li[lower.tri(Li, diag = TRUE)] <- theta
val <- tcrossprod(sig * Li)
stddev <- sqrt(diag(val))
corr <- t(val/stddev)/stddev # or cov2cor(val), but I suspect there's a
# reason for this
corr[2,1]
}
getRanCorr(fm1)
##[1] 0.08131969
VarCorr(fm1)
##Groups Name Std.Dev. Corr
##Subject (Intercept) 23.7806
## Days 5.7168 0.081
##Residual 25.5918
And with your model,
mod <- lmer(Y ~ X + (X|ID), data = df, REML = FALSE)
getRanCorr(mod)
##[1] -0.870463
attr(VarCorr(mod)$ID,"correlation")[2,1]
##[1] -0.870463