# Reparameterization of the variance-covariance matrix (apVar) of the random-effect parameters estimated by lme

The question is about computing the variance of the random-effect parameters estimated when fitting a linear mixed-effect model when the parameterization of the random-effect parameters changes. This link describes how to use the delta rule to compute the variance of the estimated random effects. I would like to extend it to compute the variance of the covariance parameters but I am not obtaining the expected result. The post is long because it includes all the derivations as well as a concrete example with R.

Let $$\Sigma_S$$ be the variance-covariance matrix of two random effects (e.g. the random intercept and random slope of a linear mixed-effect model): $$\begin{bmatrix} b_{0i} \\\\ b_{1i} \end{bmatrix} \sim \mathcal{N}(0, \Sigma_S)$$ where $$\Sigma_S = \begin{bmatrix} \sigma_{11} & \sigma_{12} \\\\ \sigma_{12} & \sigma_{22} \end{bmatrix}$$ Let $$S_{ij}$$ be the estimated values of the random effect parameters $$\sigma_{ij}$$. The $$S_{ij}$$ are random variables and we can define the 3 by 3 variance-covariance matrix $$\Psi_S$$ of estimated values $$S_{ij}$$ of the random-effect parameters: $$\Psi_S = \text{Var} \left( \begin{bmatrix} S_{ii} \\\\ S_{jj} \\\\ S_{ij} \end{bmatrix} \right) = \begin{bmatrix} \text{Var}(S_{ii}) & \text{Cov}(S_{ii},S_{jj}) & \text{Cov}(S_{ii},S_{ij}) \\\\ \text{Cov}(S_{ii}, S_{jj})& \text{Var}(S_{jj}) & \text{Cov}(S_{jj},S_{ij}) \\\\ \text{Cov}(S_{ii}, S_{ij})& \text{Cov}(S_{jj},S_{ij}) & \text{Var}(S_{ij}) \end{bmatrix}$$ The diagonal elements of this matrix $$\text{Var}(S_{ii})$$, $$\text{Var}(S_{jj})$$ and $$\text{Var}(S_{ij})$$ are the variances of the estimated values of the random-effect parameters $$\sigma_{ij}$$ in matrix $$\Sigma_S$$.

Now, let's assume that a different parameterization of the variance and covariance parameters $$\sigma_{ij}$$ is used (e.g., the so-called "natural" parameterization in Pinheiro, p. 93):
\begin{aligned} \eta_i &= \log(\sqrt{\sigma_{ii}}) \\\\ \eta_{ij} &= \text{logit} \left( \rho_{ij} \right) \end{aligned} where $$\rho_{ij} = \frac{\sigma_{ij}}{\sqrt{\sigma_{ii} \sigma_{jj}}}$$ and $$\text{logit}(x) = \log \left( \frac{x+1}{x-1} \right)$$. Let $$T_i$$ and $$T_{ij}$$ be the estimated values of $$\eta_i$$ and $$\eta_{ij}$$ obtained when fitting the mixed-effect model using this parameterization.

The inverse transformation can be used to get the parameters of $$\Sigma_S$$ from their natural parameterization \begin{aligned} \sigma_{ii} &= g_1(\eta_{ii}) = \exp(\eta_i)^2 \\\\ \sigma_{ij} &= g_2(\eta_{ij}, \eta_i, \eta_j) = \text{logit}^{-1}( \eta_{ij} ) \exp(\eta_i) \exp(\eta_j) \end{aligned} \tag{1}\label{1} where $$\rho_{ij}=\text{logit}^{-1}(\eta_{ij}) = \frac{\exp(\eta_{ij})-1}{\exp(\eta_{ij})+1}$$ is the correlation of coefficient between the random effects and $$\exp(\eta_i) = \sqrt{\sigma_{ii}}$$ is the standard devation of the random effects.

Now, assume that we have estimated the random effect parameters $$T_i$$, $$T_j$$ and $$T_{ij}$$ in the so-called natural scale and that we have computed the corresponding variance-covariance matrix: $$\Psi_T = \text{Var} \left( \begin{bmatrix} T_i \\\\ T_j \\\\ T_{ij} \end{bmatrix} \right) = \begin{bmatrix} \text{Var}(T_{ii}) & \text{Cov}(T_{ii},T_{jj}) & \text{Cov}(T_{ii},T_{ij}) \\\\ \text{Cov}(T_{ii}, T_{jj})& \text{Var}(T_{jj}) & \text{Cov}(T_{jj},T_{ij}) \\\\ \text{Cov}(T_{ii}, T_{ij})& \text{Cov}(T_{jj},T_{ij}) & \text{Var}(T_{ij}) \end{bmatrix}$$

The general question is how to transform $$\Psi_T$$ into $$\Psi_S$$. To start,let's focus on the diagonal terms, i.e the variances $$\text{Var}(T_i)$$, $$\text{Var}(T_j)$$ and $$\text{Var}(T_{ij})$$.

It is quite simple to use the delta-rule to transform the first two terms, i.e. the variance of the variance of the estimated random-effect parameters $$\text{Var}(T_i)$$ into $$\text{Var}(S_{ij})$$ (see also here). The delta rule is based on first-order Taylor approximation of the transformation $$Y = g(X)$$ $$Y \approx g(\mu_X) + g'(X)( X - \mu_X)$$ where $$X$$ and $$Y$$ are two random variables and $$\mu_x = E[X]$$. It is easy to see that $$\text{Var}(Y) \approx g'(\mu_X)^2 \text{Var}(X)$$ since $$\mu_X$$ and $$g(\mu_X)$$ are constants.

In this case, the derivative of the transformation $$S_{ii} = g_1(T_i)$$ is $$g_1'(T_i) = \frac{d}{dT_i} \exp(T_i)^2 = 2 \exp(T_i)^2$$ (see eq. $$\eqref{1}$$). Therefore, the transformation from $$\text{Var}(T_j)$$ to $$\text{Var}(S_{ij})$$ is: $$\text{Var}(S_{ii}) \approx \left(2 \exp(T_i)^2 \right)^2 \text{Var}(T_i) \tag{2}\label{2}$$

The delta rule can also be used to transform the variance of correlation/covariance term $$\text{Var}(T_{ij})$$ but, in this case, the transformation $$g$$ has three arguments (see eq. $$\eqref{1}$$).

The first-order Taylor approximation of $$Y=g(X_1, X_2, X_3)$$ is $$Y \approx g(\mu_{X_1}, \mu_{X_2}, \mu_{X_3}) + \sum_i \frac{dg}{dX_i}(\mu_{X_1}, \mu_{X_2}, \mu_{X_3}) (X_i-\mu_{X_i}) = \nabla_g(\mu_X)^T (X-\mu_X)$$ The variance of $$Y$$ is $$\text{Var}(Y) \approx \sum_i \left( \frac{dg}{dX_i}(\mu_{X_1}, \mu_{X_2}, \mu_{X_3}) \right) ^2 \text{Var}(X_i) = \left( \nabla_g(\mu_X) ^2 \right)^T \text{Var}(X)$$ where $$\text{Var}(X) = \left[ \text{Var}(X_1), \text{Var}(X_2), \text{Var}(X_3) \right]$$.

In this case, the gradient is $$\nabla_{g} = \begin{bmatrix} \frac{\delta g_2}{\delta T_{ij}} \\\\ \frac{\delta g_2}{\delta T_i} \\\\ \frac{\delta g_2}{\delta T_j} \end{bmatrix} = \begin{bmatrix} \left[ \frac{d}{d T_{ij}}\text{logit}^{-1}(T_{ij}) \right] \exp(T_i) \exp(T_j) \\\\ \text{logit}^{-1}(T_{ij}) \left[ \frac{d}{d T_i} \exp(T_i) \right] \exp{T_j} \\\\ \text{logit}^{-1}(T_{ij}) \exp(T_i) \left[ \frac{d}{d T_j} \exp(T_j) \right] \end{bmatrix}= \begin{bmatrix} \frac{2 \exp(T_{ij})}{(\exp(T_{ij})+1)^2} \exp(T_i) \exp(T_j) \\\\ S_{ij} \\\\ S_{ij} \end{bmatrix}$$ since $$\frac{d}{dx}\exp(x)=\exp(x)$$ and $$S_{ij} = R_{ij} \sqrt{S_{ii} S_{jj}}$$ where $$R_{ij} = \text{logit}^{-1}(T_{ij})$$ is the sample correlation of $$\rho_{ij}$$, and $$\sqrt{S_{ii}} = \exp(T_i)$$ is the sample standard deviation $$\sqrt{\sigma_{ii}}$$.

Therefore, the variance of the sample covariance $$S_{ij}$$ is $$\text{Var}(S_{ij}) = \left( \frac{2 \exp(T_{ij})}{(\exp(T_{ij})+1)^2} \exp{T_i} \exp{T_j} \right)^2 \text{Var}(T_{ij}) + S_{ij}^2 \text{Var}(T_i) + S_{ij}^2 \text{Var}(T_j) \tag{3}\label{3}$$

I include an example to be concrete and check the formulae

library(nlme)
library(lmeInfo)

# data set
data("Orthodont",package="MEMSS")
colnames(Orthodont) <- tolower(colnames(Orthodont))
contrasts(Orthodont$sex) <- "contr.sum" # Lme model fit.lme <- nlme::lme(fixed = distance ~ age*sex, random=~age|subject, data=Orthodont, method="REML")  The function VarCorr and the function lmeInfo::extract_varcomp returns the estimated values of the random effects $$T_{ij}$$ (matrix $$\Sigma_S$$): > VarCorr(fit.lme) subject = pdLogChol(age) Variance StdDev Corr (Intercept) 5.78643480 2.4055009 (Intr) age 0.03252449 0.1803455 -0.668 Residual 1.71620366 1.3100396 > extract_varcomp(fit.lme)$Tau
$subject subject.var((Intercept)) subject.cov(age,(Intercept)) subject.var(age) 5.78643480 -0.28962733 0.03252449  Internally, lme use the natural parameterization of the variance-covariance matrix. The slot apVar of the fitted model contains both the estimates $$T_i$$ and $$T_{ij}$$ of the random effect parameters (in the attribute Pars) and their covariance matrix $$\Psi_T$$ (the matrix aPVar). > (aV <- fit.lme$apVar)
reStruct.subject1 reStruct.subject2 reStruct.subject3       lSigma
reStruct.subject1        0.19534695        0.20547534       -0.47083968 -0.017153480
reStruct.subject2        0.20547534        0.32692918       -0.57907986 -0.024253069
reStruct.subject3       -0.47083968       -0.57907986        1.35401393  0.044474133
lSigma                  -0.01715348       -0.02425307        0.04447413  0.009241268
attr(,"Pars")
reStruct.subject1 reStruct.subject2 reStruct.subject3            lSigma
0.8777582        -1.7128810        -1.6128707         0.2700573
attr(,"natural")
[1] TRUE


The function ranef_par implements the transformation from the natural parameterization to the standard parameterization of the random effect parameters (eq. $$\eqref{1}$$) as well as the transformations between the variances of their estimates (eq. $$\eqref{2}$$ and $$\eqref{3}$$).

ranef_params <- function(lmeObject) {
aV <- lmeObject$apVar # number of random effect parameters including covariance terms k <- ncol(aV)-1 # number of random effects (k = q(q+1)/2) q <- 0.5*(sqrt(8*k+1)-1) # # standard deviation terms # mu.eta.sd <- attr(aV,"Pars")[1:q] var.eta.sd <- diag(aV)[1:q] # inverse transformation mu.sigma.ii <- exp(mu.eta.sd)^2 g.exp2 <- function(x) 2*exp(x)^2 # derivative var.sigma.ii <- g.exp2(mu.eta.sd)^2 * var.eta.sd # # correlation terms # mu.eta.r <- attr(aV,"Pars")[(q+1):k] var.eta.r <- diag(aV)[(q+1):k] # inverse transformation invlogit <- function(x) (exp(x)-1)/(exp(x)+1) mu.sigma.ij <- invlogit(mu.eta.r)*prod(exp(mu.eta.sd)) g.invlogit <- function(x) 2*exp(x)/(exp(x)+1)^2 grad.g <- rbind(g.invlogit(mu.eta.r)*prod(exp(mu.eta.sd)), mu.sigma.ij, mu.sigma.ij) var.sigma.ij <- t(grad.g^2) %*% c(var.eta.r, var.eta.sd) # res <- cbind( eta = c(mu.eta.sd, mu.eta.r), var.eta = c(var.eta.sd, var.eta.r), sigma = c(mu.sigma.ii, mu.sigma.ij), var.sigma = c(var.sigma.ii, var.sigma.ij)) aux <- outer(1:q,1:q,paste0) rownames(res) <- c(paste0("re.",diag(aux)),paste0("re.",aux[upper.tri(aux)])) attr(res,"rho.ij") <- unname(invlogit(mu.eta.r)) attr(res,"drho.ij_deta.ij") <- unname(g.invlogit(mu.eta.r)) attr(res,"grad.rhoij") <- as.vector(grad.g) res }  > ranef_params(fit.lme) eta var.eta sigma var.sigma re.11 0.8777582 0.1953469 5.78643480 26.163072722 re.22 -1.7128810 0.3269292 0.03252449 0.001383358 re.12 -1.6128707 1.3540139 -0.28962733 0.063383293 attr(,"rho.ij") [1] -0.6676191 attr(,"drho.ij_deta.ij") [1] 0.2771423 attr(,"grad") [1] 0.1202302 -0.2896273 -0.2896273  The columns corresponds to • eta: estimated values of the random-effect parameters in the natural scale • var.eta: variances of the random effect parameters in the natural scale (diagonal of $$\Psi_T$$) • sigma: estimates values of random-effect parameters (elements of matrix $$\Sigma_S$$) • var.sigma: variances of the random effect parameters in the natural scale (diagonal of $$\Psi_S$$) The two first rows refer to the variances of the random effects while the third row refers to their covariance (note: the function is limited to models with only two random effects). The attributes give some intermediary results in the computation of the sample correlation of $$\rho_{ij}$$ and variance of the sample covariance $$\text{Var}(S_{ij})$$. The function lmerInfo::varcomp_vcovcomputes the variance-covariance matrix $$\Psi_S$$ of the sample random-effect parameters and can be used to check the results (I excluded the fourth dimension that corresponds to the residual error and reordered the terms): > # variance-covariance matrix > re.vv <- varcomp_vcov(fit.lme, type="expected") > # variance terms (diagonal elements) > diag(re.vv)[c(1,3,2)] Tau.subject.var((Intercept)) Tau.subject.var(age) Tau.subject.cov(age,(Intercept)) 26.370535502 0.001392966 0.172430864 > # full matrix > round(re.vv[c(1,3,2),c(1,3,2)], 4) Tau.subject.var((Intercept)) Tau.subject.var(age) Tau.subject.cov(age,(Intercept)) Tau.subject.var((Intercept)) 26.3705 0.1561 -2.0160 Tau.subject.var(age) 0.1561 0.0014 -0.0147 Tau.subject.cov(age,(Intercept)) -2.0160 -0.0147 0.1724  The diagonal show the variances of the random effects ($$\text{Var}(S_{ii})$$ and $$\text{Var}(S_{jj})$$) and the of the covariance $$\text{Var}(S_{ij})$$. Comparing the results show a discrepancy with the values returned by the function ranef_params for the variance of the covriance term (0.1724 vs 0.0634). I would appreciate help the error or reason for this discrepancy. I also wonder how to obtain the off-diagonal terms of $$\Psi_S$$ from $$\Psi_T$$. Update After the answer, it is easy the correct the line in the function above function that computes $$\text{Var}(S_{ij})$$: > var.sigma.ij <- t(grad.g) %*% aV[3:1,3:1] %*% grad.g > var.sigma.ij reStruct.subject3 reStruct.subject3 0.170976  The rows and column of the matrix aV need to be reordered to be consistent with the gradient definition. ## 1 Answer Your derivation of $$\text{Var}(S_{ij})$$ is incorrect because it treats $$T_i$$, $$T_j$$, and $$T_{ij}$$ as uncorrelated. The correct delta method formula would use $$\text{Var}(S_{ij}) \approx \left(\nabla_g(\mathbf{T})\right)' \boldsymbol\Psi_T \left(\nabla_g(\mathbf{T})\right)$$ This blog post provides more background on the multivariate delta method approximation. You can use the formulas there to get the off-diagonal terms also. One last note: The fit.lme$apVar matrix is a numerical approximation (using numerical derivative approximations) to the observed information matrix, whereas varcomp_vcov() calculates the expected information matrix using analytic derivatives. Thus, even with the correct formulas, you will not get an exact match between the fit.lme\$apVar matrix and the (suitably transformed) lmeInfo::varcomp_vcov() matrix.

• Thank you for the answer and the link toward blog post. I get now the correct result. I have update the question to include the change in the code. Commented Nov 18, 2023 at 22:15