Copulas For Non-Gaussian Patient Dependence

Why Separate Marginals From Dependence?

A virtual patient is a joint collection of traits, not a set of unrelated parameter draws. Earlier material on encoding correlated patient traits motivates specifying relationships among traits; copulas provide a way to do this when the traits themselves do not have convenient Gaussian distributions. For example, a clearance parameter may be right-skewed, a biomarker may be bounded, and body weight may have a heavy-tailed distribution. A copula lets each trait retain its specified marginal distribution while a separate model describes how the traits move together.

This separation is useful, but it does not make every dependence structure interchangeable with every set of marginals. The copula determines how probability is distributed across joint combinations, including the relative frequency of jointly high or jointly low values. The marginals determine the scale, units, and one-trait distributions.

Definition: A copula is a joint cumulative distribution function whose univariate marginal distributions are uniform on $[0,1]$. If $X_1,\ldots,X_d$ have continuous marginal cumulative distribution functions $F_1,\ldots,F_d$, their joint distribution can be written as $$ F(x_1,\ldots,x_d)=C\bigl(F_1(x_1),\ldots,F_d(x_d)\bigr), $$ where $C$ is their copula. For continuous marginals, the copula is unique. It represents dependence on the uniform probability scale, separate from the marginal scales.

The probability integral transform explains this representation. If $X_j$ has a continuous distribution $F_j$, then $U_j=F_j(X_j)$ is uniform on $[0,1]$. The joint distribution of $U_1,\ldots,U_d$ is the copula. Conversely, if a vector $U$ has copula $C$, then applying the marginal quantile functions $F_j^{-1}$ produces traits with the requested marginal distributions and copula $C$. Here $F_j^{-1}(u)$ denotes the quantile at probability $u$.

Gaussian Copula Construction

A Gaussian copula starts with a multivariate normal latent vector and transforms each component to a uniform variable. Let $R$ be a symmetric positive-definite correlation matrix, and let $Z\sim N_d(0,R)$. For the standard normal cumulative distribution function $\Phi$, define $U_j=\Phi(Z_j)$, then set $X_j=F_j^{-1}(U_j)$. The resulting $X_j$ has marginal distribution $F_j$. The transformation is increasing, so it preserves the ordering of each trait.

Z\sim N_d(0,R),\qquad U_j=\Phi(Z_j),\qquad X_j=F_j^{-1}(U_j),\quad j=1,\ldots,d.

For simulation, factor $R=LL^\mathsf{T}$, draw independent standard normal components in $\varepsilon$, and calculate $Z=L\varepsilon$. Then apply $\Phi$ componentwise and the selected marginal quantile functions. The matrix must be positive semidefinite to define a valid multivariate normal distribution; a positive-definite matrix permits the usual nonsingular factorization. A matrix assembled from pairwise correlations can fail this condition, so validity must be checked rather than presumed.

Worked Example: Skewed and Heavy-Tailed Marginals

Suppose a hypothetical two-trait virtual population uses a lognormal clearance with median $10$ L/hour and log-scale standard deviation $0.5$, and an exponential biomarker concentration with mean $4$ ng/mL. Use a Gaussian copula with latent correlation $0.6$. To illustrate one simulated pair, take $Z_1=0.5$ and an independent standard-normal draw $E=-0.25$. Since $0.6^2+0.8^2=1$, the second latent component can be formed as $Z_2=0.6Z_1+0.8E$.

First, $Z_2=0.6(0.5)+0.8(-0.25)=0.10$. The corresponding uniforms are $U_1=\Phi(0.5)\approx0.6915$ and $U_2=\Phi(0.10)\approx0.5398$.

For the lognormal marginal, the quantile is $\exp(\ln(10)+0.5\Phi^{-1}(U_1))$. Since $\Phi^{-1}(U_1)=0.5$, this gives $\exp(2.3026+0.25)\approx12.84$ L/hour. For an exponential distribution with mean $4$, the quantile is $-4\ln(1-u)$. Thus the biomarker value is $-4\ln(1-0.5398)\approx3.10$ ng/mL.

This pair is one possible draw, not evidence that the observed traits have Pearson correlation $0.6$. That number is the correlation of the latent normal variables. The nonlinear marginal transformations generally change the Pearson correlation of the traits, while preserving their rank ordering.

What Dependence Does the Gaussian Parameter Describe?

For a Gaussian copula, the off-diagonal entries of $R$ describe latent normal correlations. Because each mapping $Z_j\mapsto\Phi(Z_j)\mapsto F_j^{-1}(\Phi(Z_j))$ is increasing for continuous marginals, the ordering of trait $j$ is the ordering of $Z_j$. Consequently, rank dependence can be related directly to the latent correlation even when the traits are skewed or bounded. Their Pearson correlation, which depends on the marginal scales and shapes, usually cannot be read directly from $R$.

For a bivariate Gaussian copula with latent correlation $\rho$, Kendall’s rank correlation $\tau$ and Spearman’s rank correlation $\rho_S$ are

\tau=\frac{2}{\pi}\arcsin(\rho), \qquad \rho_S=\frac{6}{\pi}\arcsin\!\left(\frac{\rho}{2}\right).

These relationships apply to the continuous Gaussian-copula model. They provide a way to choose a latent correlation when the target is a rank statistic, not a way to force an arbitrary Pearson correlation between non-Gaussian traits.

Worked Example: Calibrating to Kendall’s Tau

Suppose a hypothetical two-trait population is intended to have Kendall’s $\tau=0.50$ under a Gaussian copula. Rearranging the Gaussian-copula relationship gives $\rho=\sin(\pi\tau/2)$. Substitution yields

\rho=\sin\left(\frac{\pi(0.50)}{2}\right) =\sin\left(\frac{\pi}{4}\right) \approx0.7071.

The latent correlation matrix is therefore $R=\begin{pmatrix}1&0.7071\\0.7071&1\end{pmatrix}$. Its eigenvalues are approximately $1.7071$ and $0.2929$, both positive, so it is positive definite. The corresponding Spearman correlation is $\rho_S=(6/\pi)\arcsin(0.7071/2)\approx0.690$.

The calibration targets rank dependence: under the model, the population Kendall correlation is $0.50$. It does not imply that the generated traits have Pearson correlation $0.7071$. In a finite simulated population, the calculated rank correlation will also vary around its population value.

Rank-Based Copulas

A rank-based construction can transfer a dependence pattern observed in paired reference data without assuming that this pattern follows a Gaussian copula. For $n$ paired observations on $d$ traits, rank each trait separately. A rank row records the relative positions of one reference individual on all $d$ traits. To generate a new population with chosen marginals, sort the desired values for each trait and assign them according to the reference rank pattern. The paired ranks, rather than the original measurements, carry the dependence into the generated population.

One common representation uses pseudo-observations $u_{ij}=r_{ij}/(n+1)$, where $r_{ij}$ is the rank of individual $i$ on trait $j$. Their empirical joint distribution estimates the copula over the observed ranks. When transferring ranks to newly generated marginal values, a convenient finite-sample convention is to use sorted marginal quantiles at probabilities $(r-0.5)/n$, or to use a supplied set of $n$ sorted candidate values. These are finite-sample choices; they do not define a unique population copula beyond the observed range.

Worked Example: Transferring a Rank Pattern

Suppose five hypothetical reference patients have rank pairs for traits A and B of $(1,2),(2,1),(3,4),(4,3),(5,5)$. A new candidate set has sorted A values $[2,3,5,9,20]$ mg/L and sorted B values $[1,2,4,8,16]$ hours. Assign each B value according to its rank in the paired pattern for the corresponding A rank.

The generated pairs are $(2,2),(3,1),(5,8),(9,4),(20,16)$, with A in mg/L and B in hours. Each marginal contains exactly the five requested values. The rank pairs retain the reference ordering pattern.

For these five pairs, the B ranks in A-rank order are $2,1,4,3,5$. Spearman’s sample rank correlation is $1-6\sum d_i^2/[n(n^2-1)]$, where $d_i$ is the difference between the paired ranks. Here $\sum d_i^2=1+1+1+1+0=4$, so $1-6(4)/[5(25-1)]=0.80$. This statistic describes this small rank pattern; five observations provide only coarse information about dependence, particularly in the tails.

Justification and Practical Choices

The marginal-preservation claim follows directly from the construction. If $Z_j$ is standard normal, then for $0<u<1$,

P(U_j\leq u) =P(\Phi(Z_j)\leq u) =P(Z_j\leq\Phi^{-1}(u)) =\Phi(\Phi^{-1}(u)) =u.

Thus $U_j$ is uniform. For a continuous marginal $F_j$, the quantile transform gives $P(F_j^{-1}(U_j)\leq x)=P(U_j\leq F_j(x))=F_j(x)$. Hence each generated trait has the specified marginal. The joint dependence is determined by the joint distribution of the uniforms, namely the chosen copula.

The Gaussian approach is convenient when a smooth parametric dependence model is appropriate and a valid correlation matrix can be specified. The rank-based approach is useful when paired data provide a dependence pattern worth retaining, without imposing a Gaussian copula shape. A Gaussian copula gives a compact model that can be simulated beyond a finite observed rank table, but imposes a particular dependence structure. An empirical rank construction retains observed orderings more directly, but its resolution and extrapolation are limited by the reference sample.

In more than two dimensions, matching a collection of pairwise rank correlations does not by itself ensure that a valid joint dependence model exists. For a Gaussian copula, the calibrated latent matrix must be positive semidefinite. A rank-based sample provides a joint rank pattern when its rows are retained together, but a small sample may not represent rare combinations reliably.

Diagnostics, Failure Modes, and Limitations

Check the marginals. Compare generated marginal quantiles and ranges with the specified distributions. A correct copula does not repair a mistaken marginal model.
Check dependence on the intended scale. Report whether a target is Pearson, Spearman, or Kendall correlation. Compare generated rank summaries with a rank target; do not compare latent $\rho$ directly with trait Pearson correlation.
Check matrix validity. A non-positive-semidefinite latent correlation matrix cannot define the intended multivariate normal draw. Review the matrix eigenvalues or use a documented valid-matrix procedure rather than silently proceeding.
Inspect tail combinations. For a Gaussian copula with $|\rho|<1$, asymptotic tail dependence is zero: extreme co-occurrence becomes relatively unlikely as the threshold approaches the marginal extremes. Heavy-tailed marginals alone do not remove this limitation of the copula.
Account for ties and discreteness. Continuous strictly increasing transformations preserve ranks. Discrete marginals create ties, so rank dependence may depend on the tie-handling convention and the copula is not uniquely identified by the margins.
Recognize empirical resolution. An empirical rank pattern has only the observed number of rank levels. It cannot reveal unobserved tail behavior, and a small reference sample can make its apparent dependence unstable.

Copulas describe dependence, not biological validity. A mathematically valid joint draw may still be physiologically implausible. Likewise, preserving marginal distributions and rank dependence does not establish that either was estimated precisely. Those questions require separate evaluation of the population targets and the biological restrictions relevant to the model.

Common pitfall: “The Gaussian copula has $\rho=0.7$, so the generated traits must have Pearson correlation $0.7$.” This is generally false after nonlinear marginal transformations. The value $0.7$ is latent normal correlation; calibrate and assess the dependence measure actually required, and verify it on generated traits.

Check Your Understanding

Answer each question, then compare with the brief solutions.

If $Z_j$ is standard normal and $U_j=\Phi(Z_j)$, what is the distribution of $U_j$, and why?
A bivariate Gaussian copula is calibrated to Kendall’s $\tau=0.30$. Calculate its latent correlation $\rho$.
Why does a heavy-tailed marginal not, by itself, guarantee strong joint tail dependence under a Gaussian copula?
What is one advantage and one limitation of transferring dependence with an empirical rank pattern?

Solutions. (1) $U_j$ is uniform on $[0,1]$, since $P(U_j\leq u)=P(Z_j\leq\Phi^{-1}(u))=u$. (2) $\rho=\sin(\pi\tau/2)=\sin(0.15\pi)\approx0.4540$. This is the latent correlation, not necessarily the traits’ Pearson correlation. (3) The copula governs joint tail behavior separately from marginal tail heaviness; for a Gaussian copula with $|\rho|<1$, asymptotic tail dependence is zero. (4) It can retain an observed multivariate ordering pattern without assuming a Gaussian copula, but its resolution and tail information are limited by the reference sample size.