Why Covariates Must Sometimes Be Sampled Together
A virtual population is defined not only by the separate distributions of its covariates, but also by how those covariates occur together. For example, body weight and renal function may each have plausible marginal distributions, while their joint pattern determines whether simulated patients resemble the intended population. Sampling each variable separately can pair values that were uncommon or absent in the source population.
That mismatch matters when model predictions depend on more than one covariate, especially through nonlinear terms or interactions. Even if independent sampling preserves the mean and standard deviation of every covariate, it removes their dependence. The resulting average prediction can therefore differ from the target-population average.
The target population should guide the choice of sampling method. If representative records are available, resampling complete records can preserve the dependence present in those data. If a suitable parametric joint model is available, it can generate new combinations while representing the estimated dependence. Neither method is automatically adequate: each depends on the quality of the target data or model.
Empirical Joint Sampling: Keep Records Intact
Suppose a data set contains one row per individual and columns for several continuous covariates. Empirical joint sampling draws rows, with replacement, from that data set. Each selected row contributes its entire covariate vector. The method preserves the observed combinations and, in expectation, the empirical marginal distributions and dependence.
Drawing each column separately is not equivalent. It combines values from different rows and destroys the original pairing. The following small hypothetical data set makes that difference concrete.
Worked Example: Resampling Paired Covariate Records
A hypothetical source population has four covariate records. Weight is measured in kilograms and renal function in milliliters per minute.
| Record | Weight (kg) | Renal function (mL/min) |
|---|---|---|
| A | 50 | 120 |
| B | 60 | 105 |
| C | 70 | 90 |
| D | 80 | 75 |
To generate a virtual population of four by empirical joint resampling, draw four row labels with replacement, such as B, D, B, and C. The resulting covariate pairs are (60 kg, 105 mL/min), (80 kg, 75 mL/min), (60 kg, 105 mL/min), and (70 kg, 90 mL/min). Each pair remains intact because it came from one observed record.
The four unique source records have mean weight \((50+60+70+80)/4=65\) kg and mean renal function \((120+105+90+75)/4=97.5\) mL/min. In this particular resample, the means are \((60+80+60+70)/4=67.5\) kg and \((105+75+105+90)/4=93.75\) mL/min. A finite resample need not reproduce the source means exactly; repeated resampling reproduces the empirical distribution in expectation.
Here the source records show an exact negative linear pattern: each 10 kg increase in weight accompanies a 15 mL/min decrease in renal function. Sampling the two columns independently could instead pair the 80 kg record’s value with 120 mL/min. That combination is not among the observed records and conflicts with the pattern in this small source data set.
Empirical resampling is most defensible when records represent the intended population and the data set covers its relevant range. It cannot generate covariate values outside the observed records, and its empirical distribution may be too coarse when the source sample is small. Those limitations are important when the simulation target differs from the data used to construct the empirical distribution.
Parametric Joint Sampling of Continuous Covariates
A parametric approach specifies a distribution for the entire covariate vector. For two continuous covariates, a common choice is a bivariate normal distribution. Its parameters are the two means, two standard deviations, and their correlation. The multivariate normal model is an assumption, not a consequence of knowing those summaries: it specifies a particular joint shape and dependence structure.
To draw from this model, first generate independent standard normal variables \(Z_1,Z_2\), each with mean 0 and variance 1. Transform them using a matrix \(L\) satisfying \(\boldsymbol{\Sigma}=LL^\mathsf{T}\), then add the mean vector. This is the Cholesky construction. It creates correlated covariates from independent random inputs; the final covariates are not independent.
The covariance identity follows because adding a constant does not change covariance and \(\operatorname{Cov}(\mathbf{Z})=I\). Thus \(\operatorname{Cov}(L\mathbf{Z})=L\operatorname{Cov}(\mathbf{Z})L^\mathsf{T}=LIL^\mathsf{T}=\boldsymbol{\Sigma}\). This verifies that the transformation has the specified covariance, rather than merely the desired marginal means.
Worked Example: Generating a Correlated Pair
Assume, hypothetically, that weight \(X_1\) has mean 70 kg and standard deviation 12 kg, while renal function \(X_2\) has mean 90 mL/min and standard deviation 20 mL/min. Suppose the target correlation is \(\rho=0.60\), and assume a bivariate normal model is reasonable for these covariates over the target population.
The covariance is \(0.60(12)(20)=144\) kg·mL/min. The covariance matrix and one valid Cholesky factor are \[ \boldsymbol{\Sigma}= \begin{pmatrix} 144&144\\144&400 \end{pmatrix}, \qquad L= \begin{pmatrix} 12&0\\12&16 \end{pmatrix}. \] Indeed, \(LL^\mathsf{T}\) has diagonal entries \(144\) and \(12^2+16^2=400\), with off-diagonal entry \(12(12)=144\).
For an illustrative draw with \(Z_1=1\) and \(Z_2=-0.5\), the generated values are \[ X_1=70+12(1)=82\text{ kg}, \qquad X_2=90+12(1)+16(-0.5)=94\text{ mL/min}. \] The draw is one possible pair, not a claim that these covariate values are common. Across repeated draws, the specified means, variances, and covariance describe the model distribution.
For more than two continuous covariates, the same construction uses a vector of independent standard normal variables and a Cholesky factor of the full covariance matrix. A valid factorization requires a valid covariance matrix: diagonal entries must be nonnegative, the matrix must be symmetric, and its eigenvalues must be nonnegative. A matrix that fails these conditions cannot represent a covariance structure as specified.
Why Dependence Changes Predictions
A direct way to see the prediction consequence is to consider a model containing a product of two covariates. The expected product depends on their covariance, not only on their separate means. By the definition of covariance,
Rearranging gives \(E(X_1X_2)=\mu_1\mu_2+\operatorname{Cov}(X_1,X_2)\). Under independent sampling, the covariance is zero, so the expected product becomes \(\mu_1\mu_2\). Therefore, replacing a joint distribution with independent marginal draws can change the expected value of a prediction that includes a product term.
Worked Example: Covariance and Average Model Prediction
Consider a hypothetical prediction model, with outcome measured in an arbitrary concentration unit: \[ Y=100+0.5W+2R+0.04WR. \] Here \(W\) is weight deviation in kilograms from a reference value and \(R\) is renal-function deviation in liters per hour from a reference value. Assume \(E(W)=E(R)=0\), standard deviations 10 kg and 4 L/hour, and correlation 0.50. These are assumptions for this example, not data from a real study.
The covariance is \(\operatorname{Cov}(W,R)=0.50(10)(4)=20\) kg·L/hour. Taking expectations term by term gives \[ E(Y)=100+0.5E(W)+2E(R)+0.04E(WR). \] Because both means are zero, and \(E(WR)=\operatorname{Cov}(W,R)=20\), the joint-sampling expectation is \[ E(Y)=100+0+0+0.04(20)=100.8. \]
If \(W\) and \(R\) are instead sampled independently while retaining these same marginal distributions, their covariance is zero. The resulting expectation is \(100+0.04(0)=100.0\). Thus independent sampling lowers the mean prediction by 0.8 outcome units in this example. This difference is caused by the lost dependence interacting with the product term; it does not imply that every model or every dependence pattern produces bias in the same direction.
More generally, for a linear prediction \(Y=\beta_0+\beta_1X_1+\beta_2X_2\), its mean is \(\beta_0+\beta_1\mu_1+\beta_2\mu_2\). If only the marginal means are preserved, independence does not change that particular expectation. But it can still change the distribution of predictions, and nonlinear functions or interaction terms can change the mean as well. Assess the model structure and the target estimand rather than assuming dependence is harmless.
Choosing a Joint Sampling Approach
| Approach | What it preserves or assumes | Key limitation |
|---|---|---|
| Empirical row resampling | Observed covariate records and their pairings | Cannot generate values beyond observed records; quality depends on source data |
| Parametric joint model | Dependence and marginal shapes specified by a joint distribution | Results depend on the adequacy of the distributional assumptions and estimates |
| Independent marginal sampling | Each specified marginal distribution separately | Imposes zero dependence and may create implausible combinations or alter predictions |
The multivariate normal model is useful when its joint shape is plausible, but matching means and covariances alone does not establish that plausibility. Skewness, heavy tails, nonlinear dependence, or unusual subpopulations can make a normal joint model inadequate. If the dependence structure is not well represented by the chosen parametric family, a simulation may have realistic-looking marginal distributions and still produce unrealistic pairs.
A useful review compares the generated joint sample with the target data, not just one-variable summaries. Inspect scatterplots, pairwise correlations, ranges, and whether combinations in the tails are credible. For empirical sampling, check whether the source records cover the target population and whether rare but important combinations are represented. For a parametric model, check that its covariance matrix is valid and that the assumed shape is plausible. Marked differences in pairwise patterns or frequent implausible combinations are warning signs, not problems to conceal by reporting marginal summaries alone.
When the dependence structure is uncertain, it can be useful to compare predictions under more than one defensible joint model, including a clearly identified independent-sampling scenario as a sensitivity analysis. The independent case is not a neutral default: it is a particular assumption that may be unrealistic. The comparison should be framed as sensitivity to covariate dependence, not as proof that one scenario is correct.
Common Mistakes and Practical Limits
- Resampling columns independently and calling it empirical joint sampling. Separate column draws break the observed pairing. Draw complete records when the objective is to preserve empirical joint structure.
- Assuming correct marginal distributions guarantee a correct virtual population. Marginals say nothing about how values co-occur. Examine joint patterns and the prediction consequences of dependence.
- Using a correlation matrix that is not mathematically valid. Pairwise correlations that each lie between \(-1\) and \(1\) do not by themselves guarantee a valid multivariate covariance matrix. Check positive semidefiniteness before generating draws.
- Assuming correlation proves causation or a biological mechanism. Correlation describes statistical dependence in the target distribution. It does not establish why covariates are associated.
- Treating a finite simulation as an exact copy of the target. Random draws vary from run to run; a particular sample need not match target summaries exactly. A substantial mismatch can also indicate a coding error or an unsuitable target model.
Joint sampling is limited by what is known about the target population. Empirical records may be sparse or unrepresentative; parametric models may be misspecified; and a plausible dependence pattern does not ensure every simulated covariate profile is clinically meaningful. The simulation design should make these assumptions explicit and examine whether conclusions are sensitive to them.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- A source data set contains complete records of weight and renal function. What should be sampled to preserve their empirical dependence?
- A bivariate normal model has standard deviations 5 and 8 and correlation \(-0.25\). What is the covariance, including its units if the covariates are measured in kg and L/hour?
- For a model \(Y=10+3X_1-2X_2+0.5X_1X_2\), both covariate means are zero and their covariance is 4. What is \(E(Y)\) under joint sampling? What is it under independent sampling with the same marginal means?
- Name two checks that can reveal a problem with a proposed parametric joint sample.
Solutions. (1) Sample complete rows, so both covariate values come from the same record. (2) The covariance is \((-0.25)(5)(8)=-10\) kg·L/hour. (3) Under joint sampling, \(E(X_1X_2)=4\), so \(E(Y)=10+0+0+0.5(4)=12\). Under independent sampling, the covariance and expected product are zero, giving \(E(Y)=10\). (4) For example, check that the covariance matrix is positive semidefinite and compare generated scatterplots or pairwise correlations with the target data.