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


# data set
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.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
reStruct.subject1 reStruct.subject2 reStruct.subject3            lSigma 
        0.8777582        -1.7128810        -1.6128707         0.2700573 
[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)
> 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
[1] -0.6676191
[1] 0.2771423
[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$.


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          0.170976

The rows and column of the matrix aV need to be reordered to be consistent with the gradient definition.


1 Answer 1


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.

  • $\begingroup$ 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. $\endgroup$
    – gavril
    Nov 18, 2023 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.