Preserving Random-Effect Correlations in Simulation

Why Preserve Random-Effect Correlations?

An NLME model can allow a subject’s random effects to vary together. For example, subjects with a larger-than-typical value of one parameter may tend to have a larger or smaller value of another. A simulation that draws each random effect separately retains its marginal distribution but removes this fitted relationship. The resulting simulated subjects can therefore have a different joint distribution from the one specified by the model.

This tutorial focuses on generating draws from the fitted random-effect covariance matrix and checking that the dependence has been preserved. It assumes the random effects are modeled jointly as multivariate normal, as in the standard Gaussian random-effects specification. The earlier material on drawing random effects establishes the basic normal-draw setup; here the key issue is how to draw the components jointly.

Definition: Let \(\boldsymbol{\eta}_i\) be the vector of random effects for subject \(i\). Under a multivariate normal random-effects model, \(\boldsymbol{\eta}_i \sim N(\mathbf{0},\mathbf{\Omega})\), where \(\mathbf{\Omega}\) is the fitted random-effect covariance matrix. Its diagonal entries are the random-effect variances, and its off-diagonal entries are their covariances. The correlation between effects \(j\) and \(k\) is \(\Omega_{jk}/\sqrt{\Omega_{jj}\Omega_{kk}}\).

The covariance matrix must be symmetric and positive semidefinite. In the usual nonsingular case it is positive definite, so it has a Cholesky factor \(\mathbf{L}\), a lower-triangular matrix satisfying \(\mathbf{\Omega}=\mathbf{L}\mathbf{L}^{\mathsf T}\). If \(\mathbf{z}\) has independent standard normal components, then multiplying by this factor creates a random-effect vector with the fitted covariance structure.

\mathbf{z}\sim N(\mathbf{0},\mathbf{I}),\qquad \boldsymbol{\eta}=\mathbf{L}\mathbf{z},\qquad \mathbf{\Omega}=\mathbf{L}\mathbf{L}^{\mathsf T}

Why the Transformation Works

The construction has the required mean and covariance. Since the components of \(\mathbf{z}\) have mean zero, the linear transformation also has mean zero. Its covariance is obtained by multiplying the covariance of \(\mathbf{z}\) on the left by \(\mathbf{L}\) and on the right by its transpose:

E(\boldsymbol{\eta})=\mathbf{L}E(\mathbf{z})=\mathbf{0}, \qquad \operatorname{Cov}(\boldsymbol{\eta}) =\mathbf{L}\operatorname{Cov}(\mathbf{z})\mathbf{L}^{\mathsf T} =\mathbf{L}\mathbf{I}\mathbf{L}^{\mathsf T} =\mathbf{\Omega}.

Because a linear transformation of a multivariate normal vector is multivariate normal, \(\boldsymbol{\eta}\) has the specified multivariate normal distribution. This is the justification for the simulation procedure, not merely a convenient way to create numbers. In contrast, using separate independent draws for every random effect produces a diagonal covariance matrix: all off-diagonal covariances are zero.

Worked Example: Transforming a Three-Effect Draw

Suppose a hypothetical fitted model has the following Cholesky factor for three random effects. All entries are in the units implied by the corresponding random effects.

\mathbf{L}= \begin{pmatrix} 0.2000&0&0\\ 0.0902&0.2000&0\\ 0.0800&0.1386&0.3666 \end{pmatrix}

For one simulated subject, let the independent standard normal draw be \(\mathbf{z}=(1,-1,0.5)^{\mathsf T}\). Multiply by \(\mathbf{L}\), row by row:

\begin{aligned} \eta_1&=0.2000(1)=0.2000,\\ \eta_2&=0.0902(1)+0.2000(-1)=-0.1098,\\ \eta_3&=0.0800(1)+0.1386(-1)+0.3666(0.5)\\ &=0.0800-0.1386+0.1833=0.1247. \end{aligned}

Thus this draw is \(\boldsymbol{\eta}=(0.2000,-0.1098,0.1247)^{\mathsf T}\). The third component depends on all three standard normal inputs in the third row. Omitting the final term or changing its sign would produce a different draw.

The covariance matrix implied by this factor is \(\mathbf{\Omega}=\mathbf{L}\mathbf{L}^{\mathsf T}\). For example, the covariance between effects 1 and 3 is the product of their corresponding rows, \(0.2000(0.0800)=0.0160\). The covariance between effects 2 and 3 is \(0.0902(0.0800)+0.2000(0.1386)=0.034936\). Both are nonzero, so the components are not independent.

From a Covariance Matrix to Simulated Subjects

For each simulated subject, generate a fresh vector of independent standard normal values, multiply by the same fitted Cholesky factor, and retain the resulting random-effect vector as a unit. Reusing one component across subjects or independently reshuffling the components after transformation breaks the intended joint draws. The order of effects in \(\mathbf{z}\), \(\mathbf{L}\), and the model’s random-effect vector must also match.

Assemble the fitted covariance matrix.
Use the fitted variances and covariances in the model’s stated random-effect order. Check symmetry and confirm that the entries have the intended units.

Obtain a factor.
Compute a Cholesky factor \(\mathbf{L}\) such that \(\mathbf{L}\mathbf{L}^{\mathsf T}=\mathbf{\Omega}\), if the matrix is positive definite.

Generate and transform draws.
For each subject, draw independent standard normal components and calculate \(\boldsymbol{\eta}_i=\mathbf{L}\mathbf{z}_i\).

Check the simulated joint distribution.
Across many draws, compare the sample means, variances, covariances, and correlations with their fitted targets, allowing for Monte Carlo variation.

Worked Example: Comparing Correlated and Independent Draws

Suppose a hypothetical two-effect model has standard deviations \(0.3\) and \(0.5\), with fitted correlation \(-0.6\). The covariance is \((-0.6)(0.3)(0.5)=-0.09\), so the covariance matrix is

\mathbf{\Omega}= \begin{pmatrix} 0.09&-0.09\\ -0.09&0.25 \end{pmatrix}, \qquad \mathbf{L}= \begin{pmatrix} 0.3&0\\ -0.3&0.4 \end{pmatrix}.

The factor can be checked by multiplication: its first diagonal entry gives \(0.3^2=0.09\); the off-diagonal entry is \(0.3(-0.3)=-0.09\); and the second diagonal entry is \((-0.3)^2+0.4^2=0.09+0.16=0.25\).

For \(\mathbf{z}=(0.5,-1)^{\mathsf T}\), the correlated draw is

\boldsymbol{\eta} = \begin{pmatrix} 0.3&0\\ -0.3&0.4 \end{pmatrix} \begin{pmatrix}0.5\\-1\end{pmatrix} = \begin{pmatrix}0.15\\-0.55\end{pmatrix}.

If the same marginal standard deviations were instead simulated independently, the factor would be \(\operatorname{diag}(0.3,0.5)\), and the same \(\mathbf{z}\) would give \((0.15,-0.5)^{\mathsf T}\). Both procedures have marginal variances \(0.09\) and \(0.25\), but only the first has covariance \(-0.09\). Across a large simulation, the first method should yield a sample correlation near \(-0.6\), while the independent method should yield a sample correlation near zero. Neither finite-sample correlation is guaranteed to equal its theoretical value exactly.

What Independence Changes

An off-diagonal covariance affects more than a correlation summary. It changes the distribution of combinations of random effects. For example, in the two-effect model above, the variance of their sum is

\operatorname{Var}(\eta_1+\eta_2) =\operatorname{Var}(\eta_1)+\operatorname{Var}(\eta_2) +2\operatorname{Cov}(\eta_1,\eta_2) =0.09+0.25+2(-0.09)=0.16.

Its standard deviation is \(0.4\). Under the incorrect independence assumption, the sum variance would instead be \(0.09+0.25=0.34\), with standard deviation \(\sqrt{0.34}\approx0.5831\). Thus independence would substantially increase the simulated variability of this sum, despite leaving each effect’s marginal variance unchanged.

Worked Example: Checking a Simulated Sum

Consider a second hypothetical fitted covariance matrix:

\mathbf{\Omega}= \begin{pmatrix} 0.16&0.06\\ 0.06&0.09 \end{pmatrix}.

The fitted correlation is \(0.06/(\sqrt{0.16}\sqrt{0.09})=0.06/(0.4\cdot0.3)=0.5\). For the difference of the two random effects, the fitted variance is

\operatorname{Var}(\eta_1-\eta_2) =0.16+0.09-2(0.06)=0.13.

The corresponding standard deviation is \(\sqrt{0.13}\approx0.3606\). If the same marginal variances were simulated independently, the difference variance would be \(0.16+0.09=0.25\), with standard deviation \(0.5\). Therefore, a useful simulation check is not only whether each marginal variance is close to its target, but also whether relationships such as the variance of a sum or difference agree with the full covariance matrix.

Diagnostics, Alternatives, and Limitations

After generating a substantial number of draws, calculate their sample covariance matrix in the same random-effect order as the fitted matrix. The diagonal entries should be near the fitted variances, and the off-diagonal entries should be near the fitted covariances. A quick standardized check is to compare the sample correlations with the fitted correlations. Differences are expected for finite simulations; they should generally shrink as the number of draws increases.

If sample correlations are consistently near zero despite nonzero fitted targets, investigate whether the code uses only the diagonal variances, whether the Cholesky factor is multiplied on the correct side, and whether effect order is consistent. In the convention here, column-vector draws use \(\mathbf{L}\mathbf{z}\), with \(\mathbf{\Omega}=\mathbf{L}\mathbf{L}^{\mathsf T}\). Also check that the factor was not transposed incorrectly and that the same transformed vector was not subsequently split and recombined across subjects.

The main alternative is to simulate random effects independently. That is appropriate only when the model’s covariance structure specifies zero correlations, or when independence is an explicitly chosen approximation. A covariance matrix can also be represented by its standard deviations and correlation matrix; this is equivalent provided the conversion to covariances is done correctly. Preserving a fitted covariance matrix does not establish that its correlations are estimated precisely or that the model is adequate. It only ensures the simulation reproduces the specified covariance structure.

A Cholesky factorization can fail if the estimated matrix is not positive definite, for example because of numerical rounding or a boundary estimate. Do not silently replace problematic off-diagonal entries with zero: that changes the model being simulated. Check the matrix, estimation output, and numerical precision, and use a justified method for a positive-semidefinite covariance matrix when singularity is part of the fitted model. Finally, this construction assumes multivariate normal random effects. Matching a covariance matrix alone does not guarantee the right distribution if that assumption is inappropriate.

Common pitfall: “I drew each random effect with the correct variance, so the simulated subjects follow the fitted model.” This is false when the fitted covariance matrix has nonzero off-diagonal entries. Correct marginal distributions do not imply the correct joint distribution. Use a joint transformation such as \(\mathbf{L}\mathbf{z}\), then verify the simulated covariances as well as the variances.

Check Your Understanding

Try each question, then compare with the brief solutions below.

A fitted covariance matrix has variances \(0.25\) and \(1.00\), and covariance \(0.20\). What is the correlation?
Why does drawing each random effect independently fail to reproduce a fitted nonzero covariance?
For \(\mathbf{L}=\begin{pmatrix}0.4&0\\0.3&0.5\end{pmatrix}\) and \(\mathbf{z}=(1,-2)^{\mathsf T}\), calculate the random-effect draw.
With the covariance matrix in Question 1, calculate the variance of the sum of the two effects. What would it be under independence with the same marginal variances?

Solutions. 1. The standard deviations are \(0.5\) and \(1.0\), so the correlation is \(0.20/(0.5\cdot1.0)=0.4\). 2. Independent draws have zero covariance by construction, so they cannot reproduce a nonzero target covariance. 3. Multiply row by row: \(\eta_1=0.4(1)=0.4\), and \(\eta_2=0.3(1)+0.5(-2)=-0.7\). The draw is \((0.4,-0.7)^{\mathsf T}\). 4. The fitted sum variance is \(0.25+1.00+2(0.20)=1.65\). Under independence it is \(0.25+1.00=1.25\).