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 scalevar.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.
