Preserving Dependence Among Virtual-Subject Covariates

Why Dependence Matters in a Virtual Population

A virtual subject is not merely a collection of individually plausible parameter values. Body size, organ size, physiological function, and enzyme abundance may vary together. If a model samples each covariate independently, it can combine a very large body size with unusually small organ volumes, or pair high enzyme abundance with physiological values that would rarely occur alongside it. The individual values may each lie within an acceptable range, yet their combination may not represent the population being simulated.

Preserving dependence means specifying both the distribution of each covariate by itself and how covariates vary together. These are separate tasks. For example, the marginal distribution of liver volume describes liver volume alone; a correlation between liver volume and body size describes one aspect of their joint variation. Correlation does not establish a biological cause, and it does not fully describe every possible form of dependence.

The sampling approach in this tutorial uses a Gaussian copula: correlated standard normal variables supply the dependence structure, and transformations map them to the desired covariate distributions. This is useful when continuous covariates have different shapes, such as approximately normal body size and positively skewed enzyme abundance. It requires explicit choices about the dependence model and its inputs.

Correlated Latent Variables and Marginal Transformations

Suppose there are \(k\) continuous covariates. Let \(\mathbf{Z}=(Z_1,\ldots,Z_k)^\mathsf{T}\) have a multivariate normal distribution with mean vector zero and correlation matrix \(R\). Each \(Z_j\) has a standard normal marginal distribution. A valid \(R\) has ones on its diagonal, is symmetric, and is positive semidefinite; for the Cholesky procedure below, it must be positive definite.

Here \(L\) is a lower-triangular Cholesky factor and \(I\) is the identity matrix. Because \(LL^\mathsf{T}=R\), the generated vector has covariance matrix \(R\). To obtain covariate \(X_j\) with desired cumulative distribution function \(F_j\), transform its latent value using

where \(\Phi\) is the standard normal cumulative distribution function. Each \(U_j\) is uniform on \((0,1)\), so \(F_j^{-1}(U_j)\) has marginal distribution \(F_j\). The transformation preserves the ordering induced by each latent variable, but the resulting Pearson correlations on the original covariate scales will not generally equal the entries of \(R\). Thus \(R\) should be interpreted as the latent normal correlation matrix unless it has been calibrated for the chosen marginals.

For a normal marginal with mean \(\mu_j\) and standard deviation \(\sigma_j\), the transformation simplifies to \(X_j=\mu_j+\sigma_jZ_j\); in that case the Pearson correlation is the corresponding entry of \(R\). For a lognormal marginal, one can transform with \(X_j=\exp(\mu_j+\sigma_j Z_j)\), where \(\mu_j\) and \(\sigma_j\) are the mean and standard deviation of \(\log(X_j)\). The original-scale correlation then differs from the latent correlation.

Calibrating Correlation for Lognormal Covariates

For two lognormal covariates \(X=\exp(\mu_X+\sigma_X Z_X)\) and \(Y=\exp(\mu_Y+\sigma_Y Z_Y)\), assume \((Z_X,Z_Y)\) is bivariate standard normal with correlation \(\rho_Z\). The parameter \(\rho_Z\) is the latent correlation. The original-scale Pearson correlation \(\rho_{XY}\) is

This relationship follows from the moments of the lognormal variables. In particular, \(E(X)=\exp(\mu_X+\sigma_X^2/2)\), and the joint normal moment gives \(E(XY)=\exp(\mu_X+\mu_Y+(\sigma_X^2+\sigma_Y^2)/2+\rho_Z\sigma_X\sigma_Y)\). Subtracting \(E(X)E(Y)\) yields the covariance. Dividing by the product of the marginal standard deviations gives the displayed correlation formula.

Solving for the latent correlation gives

Why Independent Sampling Can Create Implausible Combinations

If each covariate is sampled independently, the joint distribution is the product of its marginal distributions. This construction forces every pairwise correlation to zero when the moments exist, regardless of whether the intended population has positive or negative dependence. It can therefore distort the frequency of combinations at the tails.

For example, suppose body size and an organ-physiology measure each have a defined upper decile. Under independent sampling, the probability that one subject is in both upper deciles is \(0.10\times0.10=0.01\). A model with positive dependence places more joint probability in matching low-low and high-high regions than independence does; with negative dependence, opposite-tail combinations are more common. The appropriate probabilities depend on the specified joint model and cannot be inferred from the marginals alone.

Independent sampling can also be misleading when every draw falls within the marginal ranges. A high body-size value paired with a very low organ volume may be allowed by the separate marginal distributions but uncommon in the target population. Dependence is therefore relevant to simulated subject realism, not merely to summary correlations. The effects on predicted exposure or clearance depend on the PBPK model and should be evaluated in the relevant application.

Construction Checks, Alternatives, and Limitations

Before generating a virtual population, specify which correlations are intended to refer to the original covariate scales and which are latent-normal correlations. Then document each marginal distribution, its units and parameters, the dependence assumptions, and the source or rationale for the target associations. A matrix of pairwise targets must be jointly coherent: valid pairwise values considered one at a time can still form an impossible correlation matrix.

The main alternatives have different strengths. Independent marginal sampling is simple and retains marginal distributions, but it imposes independence. Empirical resampling of complete subject records can preserve observed dependence without imposing a Gaussian copula, but it is limited to the observed covariate combinations and cannot generate combinations beyond the sampled support. Rank-based pairing or sorting can impose an association pattern, but naïve pairing may force an unrealistic dependence, including near-perfect rank agreement. A Gaussian copula is flexible about the marginals and computationally convenient, but it imposes a particular dependence structure and may misrepresent tail dependence.

Diagnostics include confirming that each marginal has the intended units and range, checking the generated correlation matrix against the target on the appropriate scale, and inspecting scatterplots for nonlinear patterns or clusters. A matrix that is nearly singular may generate highly constrained combinations and numerical instability. Extreme transformations can also magnify small tail differences: for example, inverse-CDF transforms applied to probabilities extremely close to zero or one can produce very large or small values. Appropriate numerical safeguards and biologically justified bounds may be needed, but clipping values changes the marginal distribution and should be documented.

The Gaussian copula does not guarantee biological plausibility. Pairwise correlations summarize only linear association on their stated scale; they do not capture all conditional relationships or causal mechanisms. A single correlation matrix can also be inadequate when the population contains distinct subgroups with different covariate relationships. The method represents assumptions about the target population, not facts established merely by producing a simulated cohort.

Common Pitfall

The practical principle is to preserve the intended joint structure, not only plausible marginal ranges. State the assumptions, define the scale of every target association, check that the dependence matrix is feasible, and verify the generated cohort before using it in PBPK simulations.