How to transform variance-covariance matrix into `theta` parameters in `lme4`

Question

I am trying to understand more about the parametrization of random effects in lme4 models. For example, from this model:

library(lme4)
data(Machines, package="MEMSS")
m <- lmer(score ~ Machine + (Machine|Worker), Machines)

I can extract the variance-covariance matrix for the random effects using VarCorr(m)$Worker and I obtain the following

            (Intercept)  MachineB  MachineC
(Intercept)   16.639159 11.602394 -5.497988
MachineB      11.602394 34.554029  6.431676
MachineC      -5.497988  6.431676 13.616956

From this answer I understand that internally the random are parametrized as the vector $\theta$, which should be the columwise unpacking of the lower triangular Cholesky factor. So I naively thought that I could transform the matrix above into the vector theta doing something like the following:

VarCov2Theta <- function(X){
  X <- t(chol(X)) 
  X_unpacked <- X[lower.tri(X, diag=T)]
  return(X_unpacked)
}

However the values are slightly different:

> VarCov2Theta(VarCorr(m)$Worker)
[1]  4.079113  2.844343 -1.347839  5.144292  1.995492  2.796122

> unname(getME(m, "theta"))
[1]  4.242177  2.958047 -1.401720  5.349938  2.075262  2.907898

Clearly I am missing something; I am also unable to recover the variance-covariance matrix from the factors

> getME(m, "ST")$Worker %*%  t(getME(m, "ST")$Worker)
          [,1]      [,2]      [,3]
[1,] 17.996065  2.958047 -1.401720
[2,]  2.958047 29.108053  1.844859
[3,] -1.401720  1.844859  8.715522

What am I missing here?

Could anyone shed some light on this and show me how to calculate the variance-covariance matrix from the $\theta$ and vice-versa?

Answer 1

This is due to the formulation of the likelihood function that is optimized. Lme4 uses REML by default, which profiles out the fixed effect and the residual term. This reformulation reduces the number of parameters that need to be estimated and therefore increases computational efficiency.

By the way, you also need to do the same thing to get random effect from $\theta$.

test_matrix<-matrix(0,3,3)

test_matrix[!upper.tri(test_matrix)]<-(getME(m, "theta"))


diag(test_matrix %*% t(test_matrix))* sigma(m)^2
[1] 16.63916 34.55403 13.61696

Please read this original paper for more technical explanation. You can find the formula you are looking for at (58).

Answer 2

Looking at the code of the function mkVarCorr it seems that I was missing scaling by residual variance

VarCov2Theta <- function(X, sigma){
  X <- t(chol(X/sigma^2)) 
  X_unpacked <- X[lower.tri(X, diag=T)]
  return(X_unpacked)
}

I am now able to transform the variance-covariance into the theta parametrization of lme4

> VarCov2Theta(W,sigma(m))
[1]  4.242177  2.958047 -1.401720  5.349938  2.075262  2.907898

> unname(getME(m, "theta"))
[1]  4.242177  2.958047 -1.401720  5.349938  2.075262  2.907898

I am still not clear why the residual variance appear there, so I will leave the question open for the moment, in case someone can provide some clarification.