Sampling Correlated Random Effects for Virtual Subjects

Why the Joint Distribution Matters

When a population PK model assigns interindividual variability to more than one parameter, a virtual subject is defined by a joint set of random effects, not by a collection of unrelated draws. For example, clearance and volume random effects may be negatively correlated. Drawing them independently preserves each parameter's marginal spread but removes the modeled tendency for subjects with a higher clearance effect to have a lower volume effect. That can change the simulated distribution of concentrations and exposure.

The earlier material on simulating individual profiles describes how sampled individual parameters are used to generate concentration-time profiles. Here the focus is narrower: how to sample the interindividual random effects that determine those parameters, while preserving their estimated covariance.

The distribution is specified on the random-effects scale used by the model. For instance, in a log-normal parameter model, a typical individual clearance may be written \(CL_i=CL_{\mathrm{typ}}\exp(\eta_{i,CL})\). The covariance in \(\Omega\) is then on the \(\eta\) scale; it is not the covariance of clearance values measured in physical units. A simulated \(\boldsymbol{\eta}_i\) is transformed through the model's parameter equations before the profile is calculated.

The usual model has mean-zero random effects, and independent subjects receive independent draws from the same joint distribution. The random effects within each subject are generally not independent unless the corresponding covariances in \(\Omega\) are zero. In this tutorial, \(\Omega\) is treated as fixed at its estimated value. Sampling uncertainty in that estimate is a separate issue.

Generating Correlated Draws

A standard construction begins with independent standard normal variables and applies a matrix transformation. Suppose \(\Omega\) is positive definite and has Cholesky factor \(L\), so that \(\Omega=LL^\mathsf{T}\), with \(L\) lower triangular. Draw a vector \(\mathbf{z}_i\) whose components are independent \(N(0,1)\) variables, then set \(\boldsymbol{\eta}_i=L\mathbf{z}_i\).

This construction gives the desired mean and covariance. Since \(E(\mathbf{z}_i)=\mathbf{0}\), linearity of expectation gives \(E(\boldsymbol{\eta}_i)=L\mathbf{0}=\mathbf{0}\). Using the covariance of a linear transformation,

Thus the factor does not merely create effects with the right marginal variances: it creates the full covariance matrix. The same construction can be written using a standard normal vector \(\mathbf{z}_i\) and any matrix \(A\) satisfying \(AA^\mathsf{T}=\Omega\); the Cholesky factor is a convenient, widely used choice when the matrix is positive definite.

Implementation and Checks

For each virtual subject, generate a fresh vector of independent standard normal values, multiply by the same factor \(L\), and use the resulting random-effects vector in the model. The following sequence makes the key operations explicit.

The empirical values will not exactly equal their targets in a finite simulation. They should be close enough to be plausible given the simulation size; a discrepancy that remains large when the number of subjects increases suggests an implementation problem. It is useful to inspect the covariance matrix and the correlation matrix, because a covariance's magnitude depends on the scales of its two effects. Also check whether the simulated marginal distributions and parameter transformations behave sensibly—for example, whether a log-scale transformation produces plausible positive parameter values.

Valid Matrices and Edge Cases

A covariance matrix must be symmetric and positive semidefinite: for every vector \(\mathbf{a}\), \(\mathbf{a}^\mathsf{T}\Omega\mathbf{a}\geq0\). A positive-definite matrix gives a conventional Cholesky factor with strictly positive diagonal entries. A positive-semidefinite matrix can be singular, so a standard Cholesky routine may fail even though the matrix is mathematically valid. This can arise, for example, when two effects are perfectly correlated or when numerical estimation produces a zero-variance direction.

A singular covariance matrix can be handled by a factorization suited to positive-semidefinite matrices, such as an eigenvalue decomposition. If \(\Omega=Q\Lambda Q^\mathsf{T}\), with nonnegative eigenvalues on the diagonal of \(\Lambda\), then \(A=Q\Lambda^{1/2}\) satisfies \(AA^\mathsf{T}=\Omega\); use \(\boldsymbol{\eta}_i=A\mathbf{z}_i\). A small negative eigenvalue caused only by numerical rounding may require careful numerical treatment, but a materially negative eigenvalue indicates an invalid matrix, not a harmless factorization inconvenience. Do not silently replace a problematic covariance matrix with a diagonal one: that changes the model.

Alternatives, Diagnostics, and Limitations

Independent normal draws are appropriate only when the intended covariance matrix is diagonal. They are easy to implement, but using them for a non-diagonal estimate discards parameter compensation encoded by the fitted model. A second approach is to use a multivariate-normal random-number routine that accepts a covariance matrix directly. This can be convenient, but the user must still supply the correct parameter order and confirm that the routine interprets the input as a covariance matrix rather than a standard-deviation or correlation matrix. The factor-and-transform construction makes the dependence structure visible and is straightforward to audit.

The multivariate normal model is an assumption, not a guarantee that every estimated random-effects distribution is adequately represented by a normal distribution. Earlier work on heavy-tailed effects and subpopulation mixtures highlights possible departures from a single normal population. If those features are supported by the model or data, sampling from a single multivariate normal distribution may not represent the intended population. The simulation must reflect the selected model rather than imposing normality simply because a Cholesky procedure is convenient.

Common warning signs include an asymmetric or incorrectly ordered matrix, a factor that does not reconstruct \(\Omega\), generated correlations near zero when the target is nonzero, or simulated variances that systematically disagree with the diagonal entries. A frequent mistake is to put standard deviations on the diagonal of \(\Omega\) instead of variances. Another is to enter a correlation matrix where the software expects a covariance matrix. Check the units and scale of every entry, and verify \(LL^\mathsf{T}=\Omega\) numerically before generating a cohort.