Why the Joint Distribution Matters
A simulated pharmacometric cohort is more than a list of typical covariate values and random-effect variances. The covariates and interindividual random effects together determine which individuals are represented and how their model parameters vary. If body size and renal function are associated in the target population, simulating them independently can create too many unusual combinations. If clearance and volume random effects are correlated in the model but are generated independently, simulated concentration profiles will have the wrong between-individual structure.
The goal is to specify and generate the joint distribution required by the data-generating model. That includes the marginal distribution of each quantity, such as a lognormal distribution for weight or a normal distribution for a random effect, and the dependence among quantities. The simulation should reproduce those assumptions—not impose a convenient correlation structure merely because it is easy to generate.
For a vector of normally distributed random effects, a common specification is \(\eta_i \sim N(0,\Omega)\), where \(\eta_i\) is the random-effect vector for individual \(i\), and \(\Omega\) is its covariance matrix. The diagonal entries of \(\Omega\) are variances; the off-diagonal entries are covariances. A correlation is the covariance divided by the product of the corresponding standard deviations. Covariates need not be normal, and their distribution and relationship to the random effects must be specified separately.
Generating Correlated Normal Variables
Suppose \(Z\) is a vector of independent standard normal variables, so that \(Z \sim N(0,I)\). If a lower-triangular matrix \(L\) satisfies \(\Omega=LL^\mathsf{T}\), then \(\eta=LZ\) has mean zero and covariance matrix \(\Omega\). This provides a direct way to generate correlated normal random effects. The same construction can be used for standardized normal covariates before they are transformed to the required scales.
The equality \(E(ZZ^\mathsf{T})=I\) follows because each component has variance 1 and distinct components are independent, hence have covariance 0. A valid covariance matrix must be symmetric and positive semidefinite; for a nonsingular Cholesky factorization, it must be positive definite. The factorization fails if the requested covariances are incompatible with the variances.
Worked Example: Generating Two Correlated Normal Covariates
Consider a hypothetical population in which weight has mean 70 kg and standard deviation 12 kg, while estimated glomerular filtration rate (eGFR) has mean 90 mL/min/1.73 m² and standard deviation 20 mL/min/1.73 m². For illustration, assume both are normally distributed with correlation 0.60. Generate two illustrative covariate pairs from independent standard normal draws.
Let \(Z_1,Z_2\) be independent standard normal draws. Construct standardized variables \(U=Z_1\) and \(V=0.60Z_1+\sqrt{1-0.60^2}Z_2\). Since \(\sqrt{1-0.60^2}=0.80\), \(V=0.60Z_1+0.80Z_2\). Both variables have variance 1, and their covariance is \(0.60\).
For the first illustrative draw, take \(Z_1=1.20\) and \(Z_2=-0.50\). Then \(U=1.20\), and \(V=0.60(1.20)+0.80(-0.50)=0.72-0.40=0.32\). Transforming to the stated scales gives weight \(70+12(1.20)=84.4\) kg and eGFR \(90+20(0.32)=96.4\) mL/min/1.73 m².
For a second draw, take \(Z_1=-0.40\) and \(Z_2=1.10\). Then \(U=-0.40\), and \(V=0.60(-0.40)+0.80(1.10)=-0.24+0.88=0.64\). The resulting weight is \(70+12(-0.40)=65.2\) kg, and eGFR is \(90+20(0.64)=102.8\) mL/min/1.73 m². These two pairs illustrate the transformation; they are not expected to have a sample correlation of exactly 0.60. The specified correlation describes the population distribution, not the exact correlation in a small simulated sample.
Generating Correlated Interindividual Random Effects
Random effects are often modeled jointly as multivariate normal, even when the covariates are generated from a different distribution. For example, a model might specify random effects on log clearance and log volume with a nonzero covariance. The correlation belongs to the random-effect distribution; it does not by itself state that either random effect is correlated with a covariate.
Worked Example: Cholesky Generation of Two Random Effects
Suppose a hypothetical pharmacometric model specifies two mean-zero normal random effects with covariance matrix
The standard deviations are \(\sqrt{0.0625}=0.25\) and \(\sqrt{0.0400}=0.20\). The correlation is \(0.025/(0.25\times0.20)=0.50\). A lower-triangular factor is
To verify the second diagonal element, the bottom-right entry of \(LL^\mathsf{T}\) must be 0.0400. Thus \(0.10^2+(\sqrt{0.0300})^2=0.0100+0.0300=0.0400\). The off-diagonal entry is \(0.25(0.10)=0.025\), and the top-left entry is \(0.25^2=0.0625\), so this factor reproduces \(\Omega\).
For independent standard normal draws \(Z_1=1.00\) and \(Z_2=-0.50\), calculate \(\eta=LZ\):
The generated pair is one possible realization, not a claim that the effects must move together for every individual. The positive covariance describes a population tendency: across repeated draws, above-average values of one random effect tend to occur with above-average values of the other.
Non-Normal Covariates and Dependence
A normal transformation is not suitable for every covariate. Weight may be right-skewed, age may be bounded, and genotype or disease status may be categorical. A useful strategy for some settings is to generate correlated latent normal variables and transform them to the desired marginal distributions. The transformation preserves the chosen marginal distribution, but the resulting correlation on the observed scale is generally not equal to the latent normal correlation.
For example, a standard normal variable \(Z\) can be converted to a uniform variable by applying its cumulative distribution function, \(\Phi(Z)\). Applying the inverse cumulative distribution function of a target distribution to \(\Phi(Z)\) produces a draw with that target marginal distribution. This construction is often called a Gaussian copula approach. Its suitability depends on whether the resulting dependence, after transformation, matches the intended joint distribution; matching a latent correlation alone is not enough.
Worked Example: A Binary Covariate Associated With a Random Effect
Consider a hypothetical binary covariate \(C\), with \(C=1\) representing a characteristic present in half the population. Suppose the model also has a mean-zero random effect \(\eta=0.20Z_2\). To induce dependence, draw standard normal variables \(Z_1,Z_2\) with correlation 0.50, and set \(C=1\) when \(Z_1>0\), otherwise \(C=0\). This generates \(P(C=1)=0.50\), because a standard normal variable is positive half the time.
The association can be checked from the joint distribution. For jointly standard normal \(Z_1,Z_2\) with correlation 0.50, \(E[Z_2\mid Z_1]=0.50Z_1\). Therefore,
For a standard normal variable, \(E\{Z_1 I(Z_1>0)\}=\phi(0)=1/\sqrt{2\pi}\approx0.3989\). Thus \(E(C\eta)\approx0.20(0.50)(0.3989)=0.03989\). Since \(E(\eta)=0\), \(\operatorname{Cov}(C,\eta)=0.03989\). Also, \(\operatorname{SD}(C)=\sqrt{0.5(1-0.5)}=0.50\), and \(\operatorname{SD}(\eta)=0.20\). The observed-scale correlation is therefore \(0.03989/(0.50\times0.20)\approx0.3989\), not 0.50. This example shows why the dependence must be checked after transforming latent variables.
Specifying the Dependence You Intend to Simulate
Before generating individuals, distinguish three assumptions that can otherwise be conflated:
- Dependence among covariates: for example, an association between weight and renal function.
- Dependence among random effects: represented by off-diagonal elements in \(\Omega\), when a joint random-effect model is specified.
- Dependence between covariates and random effects: a separate feature of the joint population model. It is not implied by a correlation among random effects.
A common modeling assumption is that random effects have mean zero and covariance \(\Omega\) conditional on the covariates, with no additional covariate–random-effect dependence. That is an assumption to encode, not a mathematical default that must be true in every population. If a dependence is intended, specify its conditional or joint distribution explicitly. For example, a covariate-dependent random-effect mean may be represented as \(E(\eta\mid X)=m(X)\), where \(X\) denotes covariates and \(m(X)\) is specified by the model. Generating \(\eta\) independently of \(X\) would not reproduce that assumption.
For design comparisons, the simulated population should represent the same target population in each candidate design. Otherwise, a difference attributed to a design may instead reflect different covariate or random-effect distributions. This requires a deliberate population-generation specification for every scenario; it does not require pretending that small simulated samples will have identical empirical summaries.
Diagnostics and Failure Modes
After simulation, inspect the realized cohort against the specified population model. Useful checks include means, standard deviations, quantiles, category proportions, covariate ranges, and the intended correlations or other dependence measures. For a large sample from a continuous distribution, empirical summaries should generally be close to their targets, but they will not match exactly. For rare categories or small samples, noticeable deviations may occur by chance.
Several failures have recognizable causes:
- Wrong variance entered as a standard deviation. If a variance of 0.04 is used as a standard deviation, simulated variability is too small. Check units and square roots explicitly.
- Invalid covariance matrix. A symmetric-looking matrix can still be impossible. Check positive definiteness or positive semidefiniteness before factorization; a failed Cholesky decomposition is a warning, not a reason to silently alter the inputs.
- Dependence lost during transformation. A latent normal correlation does not generally equal the correlation after converting variables to binary, skewed, or bounded scales. Check dependence on the observed scales relevant to the model.
- Implausible values from an unsuitable marginal distribution. A normal model for a strictly positive variable can generate negative values. Select a suitable marginal distribution or transformation, and verify its implications for the dependence structure.
- Truncation changes the joint distribution. Simply discarding or clipping out-of-range values can change means, variances, and correlations. If truncation is part of the population definition, specify and simulate the truncated distribution rather than treating it as a harmless correction.
A diagnostic discrepancy is not automatically evidence of a coding error: finite samples vary. It is a reason to compare observed summaries with the target and assess whether the discrepancy is consistent with the sample size and the assumed distribution. Conversely, agreement in means alone does not establish that the joint distribution is correct; dependence and tails may still be wrong.
Choosing a Generation Approach
| Approach | Useful when | Main check or limitation |
|---|---|---|
| Multivariate normal with Cholesky factor | Variables are modeled jointly as normal with a specified covariance matrix | Confirm matrix validity, units, and that normal marginals are suitable |
| Latent normal variables with marginal transformations | Target marginals are non-normal and a copula dependence model is appropriate | Check dependence after transformation; latent correlation is not observed-scale correlation |
| Conditional generation | A model specifies a conditional distribution, such as a random effect given covariates | Verify the conditional model and the resulting joint population distribution |
| Empirical resampling of covariate records | An observed covariate sample is an appropriate representation of the target population | Resampling retains only the information and population coverage in the source sample |
No single approach is universally best. A convenient multivariate normal construction is inappropriate if the target variables are strongly skewed or categorical; an empirical resample is inappropriate if its source population does not represent the intended setting. The approach should follow the data-generating assumptions and the simulation question.
Common Mistake
A frequent error is to generate each covariate and random effect from the correct marginal distribution independently, then conclude that the simulated cohort is correct because every mean and standard deviation looks right. That reasoning ignores the joint distribution. If two quantities are dependent in the specified population, independent draws will remove that dependence and change the frequency of combinations that drive model predictions. The correction is to state the joint or conditional distribution, generate accordingly, and check both marginal summaries and dependence in the resulting cohort.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A random-effect covariance matrix has variances 0.09 and 0.16 and covariance 0.06. What is the correlation, and is the matrix positive definite?
- Why does a latent normal correlation of 0.50 not guarantee an observed-scale correlation of 0.50 after turning one variable into a binary indicator?
- A covariate has the correct mean and standard deviation in a simulated cohort, but its correlation with another covariate is near zero even though the target correlation is 0.70. Is the population simulation adequate? Explain.
- Give one reason that clipping negative values from a normally generated positive covariate can produce an unintended population distribution.
Solutions. (1) The standard deviations are 0.30 and 0.40, so the correlation is \(0.06/(0.30\times0.40)=0.50\). The determinant is \(0.09(0.16)-0.06^2=0.0144-0.0036=0.0108>0\), and both diagonal entries are positive, so this symmetric two-by-two matrix is positive definite. (2) Thresholding changes the distribution and its dependence structure; the correlation must be evaluated after the transformation on the binary and continuous scales. (3) No, not if 0.70 is part of the target joint distribution. Correct marginal summaries do not compensate for lost dependence. Check that the generation method includes the association and that the observed-scale correlation is assessed with adequate sample size. (4) Clipping creates a point mass at zero and changes the mean, variance, and generally the dependence with other variables; it does not sample from the intended truncated distribution.