PK-PD Random-Effect Correlation and Joint Estimation

Why Clearance–Sensitivity Dependence Matters

In a PK/PD model, individuals can differ both in pharmacokinetic parameters, such as clearance, and in pharmacodynamic parameters, such as sensitivity to concentration. These differences need not be independent. If individuals with higher clearance also tend to have lower sensitivity, for example, their reduced exposure and weaker response per unit exposure may reinforce one another. If the association has the opposite sign, the two differences may partly offset one another. Either pattern can affect the distribution of exposure-response outcomes and how population parameters are estimated.

The focus here is a correlation between individual random effects, not a correlation between residual errors and not merely a correlation between two estimated parameters. The distinction matters: random-effect correlation describes a population-level distribution of latent individual differences. Joint estimation can use PK and PD observations together to estimate that distribution, but the data must contain information about both effects and their association.

Definition: Let $\eta_{CL,i}$ be individual $i$'s log-scale random effect on clearance and $\eta_{S,i}$ the log-scale random effect on PD sensitivity. A standard bivariate normal model is $$ \begin{pmatrix}\eta_{CL,i}\\ \eta_{S,i}\end{pmatrix} \sim N\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix} \omega_{CL}^{2} & \rho\,\omega_{CL}\omega_{S}\\ \rho\,\omega_{CL}\omega_{S} & \omega_{S}^{2} \end{pmatrix} \right], $$ where $\omega_{CL}$ and $\omega_S$ are the random-effect standard deviations and $\rho$, between $-1$ and $1$, is their correlation. For example, $CL_i=CL_{\mathrm{typ}}\exp(\eta_{CL,i})$ and $S_i=S_{\mathrm{typ}}\exp(\eta_{S,i})$. The correlation is defined on the log-random-effect scale.

A positive $\rho$ means that above-average clearance random effects tend to occur with above-average sensitivity random effects; a negative $\rho$ means that above-average clearance effects tend to occur with below-average sensitivity effects. Because clearance and sensitivity are different parameters, their numerical scales cannot be compared directly. The standard deviations put the covariance on a scale appropriate to each effect.

How Dependence Enters an Exposure–Response Relationship

Consider linear PK with total exposure $AUC_i=\mathrm{Dose}/CL_i$, and a simple individual exposure-linked effect $Q_i=S_i AUC_i$. This illustration isolates the role of the two random effects; it is not a claim that every PD endpoint is linear in exposure. Substituting the log-normal parameterizations gives

\log Q_i =\log\left(\frac{S_{\mathrm{typ}}\mathrm{Dose}}{CL_{\mathrm{typ}}}\right) +\eta_{S,i}-\eta_{CL,i}.

Thus, the variation in log effect depends on the variance of a difference, not the sum of two independent variances. The random effects have mean zero, so the mean of $\eta_S-\eta_{CL}$ is zero, while its variance includes their covariance.

\operatorname{Var}(\eta_S-\eta_{CL}) =\omega_S^2+\omega_{CL}^2-2\rho\omega_S\omega_{CL}.

The minus sign is important. Positive correlation can reduce variation in $Q_i$: people with higher sensitivity also tend to have higher clearance, which lowers exposure. Negative correlation can increase variation: higher clearance tends to coincide with lower sensitivity, so both effects reduce response for some individuals and increase it for others relative to the population center. This statement concerns the specified $S_i AUC_i$ quantity; the consequence for a nonlinear endpoint depends on its dose-response function and the part of that function sampled by the data.

Worked Examples

Worked Example: Interpreting a Fitted Correlation

Suppose a hypothetical joint model estimates $\omega_{CL}=0.30$, $\omega_S=0.40$, and $\rho=-0.60$. Both random effects are on the log scale. Find the covariance and interpret the association.

The covariance is $$ \operatorname{Cov}(\eta_{CL},\eta_S) =\rho\omega_{CL}\omega_S =(-0.60)(0.30)(0.40)=-0.0720. $$

The negative sign indicates that individuals with higher-than-typical log clearance effects tend, on average, to have lower-than-typical log sensitivity effects. In a bivariate normal model, the conditional mean is $$ E(\eta_S\mid \eta_{CL}=x) =\rho\frac{\omega_S}{\omega_{CL}}x. $$

At $\eta_{CL}=0.30$, one clearance standard deviation above zero, the conditional mean sensitivity effect is $$ (-0.60)\frac{0.40}{0.30}(0.30)=-0.240. $$ A clearance effect of $0.30$ corresponds to a clearance ratio $\exp(0.30)=1.350$ relative to typical clearance. The conditional mean sensitivity effect of $-0.240$ corresponds to a sensitivity ratio $\exp(-0.240)=0.7866$. These are conditional population summaries, not deterministic predictions for every individual at that clearance value.

Worked Example: Covariance Changes the Spread of Exposure-Linked Effect

Use a hypothetical model with $\omega_{CL}=0.30$, $\omega_S=0.40$, and $\rho=-0.60$. Let $Q_i=S_i AUC_i$ under fixed dose and linear PK. Compare the multiplicative spread in $Q_i$ with the spread that would follow if the same random-effect standard deviations were independent.

The variance of $\log Q_i$, apart from its fixed typical-value term, is $$ 0.30^2+0.40^2-2(-0.60)(0.30)(0.40) =0.09+0.16+0.144=0.394. $$

If the effects were independent, the variance would instead be $0.09+0.16=0.250$. Under the stated log-normal model, the expected multiplicative factor relative to the typical-value term is $\exp(v/2)$, where $v$ is the variance of $\log Q_i$. With the negative correlation, this factor is $$ \exp(0.394/2)=\exp(0.197)=1.218. $$ Under independence it is $$ \exp(0.250/2)=\exp(0.125)=1.133. $$

The ratio is $1.218/1.133=1.075$, so the negative correlation increases this expected factor by about $7.5\%$ in this simplified model. This calculation illustrates a distributional consequence, not a universal percentage bias in an estimated PD parameter. For a saturable or otherwise nonlinear PD model, the effect must be evaluated through that model rather than by applying this linear-exposure result directly.

The conditional-mean expression used above follows directly from the bivariate normal model. For jointly normal $X$ and $Y$, the conditional mean of $Y$ given $X=x$ is $E(Y)+\operatorname{Cov}(X,Y)\operatorname{Var}(X)^{-1}(x-E(X))$. Here both means are zero, $\operatorname{Var}(\eta_{CL})=\omega_{CL}^2$, and $\operatorname{Cov}(\eta_{CL},\eta_S)=\rho\omega_{CL}\omega_S$. Therefore, $$ E(\eta_S\mid\eta_{CL}=x) =\frac{\rho\omega_{CL}\omega_S}{\omega_{CL}^2}x =\rho\frac{\omega_S}{\omega_{CL}}x. $$ This result describes the mean of the conditional distribution; individual sensitivity effects still vary around it.

Worked Example: A Sample Association Is Not a Joint-Model Estimate

Consider a hypothetical set of five known latent-effect pairs, included only to show the arithmetic of a sample correlation. Let the clearance effects be $(-0.20,-0.10,0,0.10,0.20)$ and the sensitivity effects be $(0.15,0.02,-0.04,-0.08,-0.05)$. Both sample means are zero. Calculate their sample correlation.

The sum of cross-products is $$ \sum_{i=1}^{5}(\eta_{CL,i}-\bar{\eta}_{CL})(\eta_{S,i}-\bar{\eta}_S) =(-0.20)(0.15)+(-0.10)(0.02)+(0)(-0.04)+(0.10)(-0.08)+(0.20)(-0.05)=-0.0500. $$ The sums of squared deviations are $0.100$ for clearance effects and $0.0334$ for sensitivity effects. Since the common divisor $n-1$ cancels in a sample correlation, $$ r=\frac{-0.0500}{\sqrt{(0.100)(0.0334)}}=-0.865. $$

This set shows a strong negative association among these five pairs. It does not establish that a population random-effect correlation is $-0.865$: five observations provide little precision, and in real analyses the latent effects are not observed. Individual empirical Bayes estimates are typically shrunk toward zero and have unequal uncertainty, so their raw sample correlation is not a substitute for estimating $\rho$ in the joint likelihood. The calculation is a descriptive diagnostic at most.

Joint Estimation and What the Data Must Contribute

A joint model represents the PK and PD observations conditional on the same individual random effects and then accounts for the population distribution of those effects. For individual $i$, let $y_{PK,i}$ and $y_{PD,i}$ denote that individual's observed PK and PD data. A schematic likelihood contribution is

L_i(\theta)=\int p(y_{PK,i}\mid\eta_{CL,i},\theta_{PK})\, p(y_{PD,i}\mid\eta_{CL,i},\eta_{S,i},\theta_{PD})\, \phi_2\left( \begin{pmatrix}\eta_{CL,i}\\\eta_{S,i}\end{pmatrix};\,0,\,\Omega \right) \,d\eta_{CL,i}\,d\eta_{S,i},

where $\theta$ collects fixed-effect and residual parameters, $\Omega$ is the random-effect covariance matrix, and $\phi_2$ is the bivariate normal density. The PD likelihood depends on clearance because clearance helps determine the concentration or exposure driving the response. The PK data inform clearance, and the PD data inform sensitivity and the response relationship. Their combined contribution can inform the covariance between the effects.

Joint estimation does not guarantee that $\rho$ is well identified. The data need to distinguish variation in clearance from variation in sensitivity, and need enough information across individuals to learn whether those variations co-occur. If the PK observations weakly inform an individual's clearance, the PD data may be consistent with several combinations of exposure and sensitivity. An estimated covariance can then be unstable, highly uncertain, or sensitive to model assumptions. This is a limitation of information in the data, not something solved merely by choosing a joint fitting algorithm.

When the fitted $\rho$ is near zero, that alone does not prove independence. The uncertainty in the estimate matters, as does whether the model can distinguish the two effects. Conversely, a large estimated magnitude is not automatically credible: inspect uncertainty, boundary behavior, stability across reasonable starting values, and whether the model reproduces relevant PK and PD patterns. A covariance matrix must also be valid; in the two-effect case this requires $-1\leq\rho\leq1$ and nonnegative variances.

Sequential Fits, Diagnostics, and Limitations

A sequential analysis that estimates PK first and then treats individual PK estimates or predicted exposures as fixed in the PD stage can fail to carry forward clearance uncertainty and clearance–sensitivity dependence. The resulting PD fit may attribute response differences to sensitivity that partly reflect uncertain exposure. The earlier material on shrinkage-induced bias in sequential PD fits addresses one related mechanism: individual PK estimates pulled toward the population mean can distort the exposure information supplied to a later PD fit. A sequential approach is not necessarily biased in every setting, but it cannot estimate the joint covariance if that dependence is excluded from its model.

A joint analysis can represent the covariance directly, but its estimate is conditional on the PK model, PD model, random-effect distribution, and residual-error assumptions. In particular, correlation among latent random effects must not be confused with correlation among PK and PD residuals. A residual correlation, if scientifically and structurally appropriate, concerns unexplained observation-level variation after the model effects are accounted for; it is a different feature and requires its own justification.

Useful checks include verifying that the estimated covariance matrix is admissible, examining uncertainty in $\rho$, and checking whether simulated or model-based predictions reproduce the joint distribution of PK exposure and PD response across individuals. Extreme estimates near $-1$ or $1$, large changes under small modeling choices, or a substantial difference between the fitted correlation and the pattern visible in individual data summaries warrant scrutiny. Individual random-effect estimates are noisy and shrunken, however, so a plot of those estimates is not a definitive test.

The main alternatives answer different questions. A model that fixes $\rho=0$ is simpler and may be reasonable when independence is scientifically defensible or the data cannot support a free covariance, but it does not represent dependence. A sequential plug-in analysis is computationally straightforward, but ordinarily loses some PK uncertainty and does not estimate this cross-domain covariance. A joint model with free $\rho$ represents dependence and propagates information through the shared model, at the cost of an additional parameter and greater demands on the data. No option is uniformly best: the choice should reflect the estimand, scientific plausibility, data information, and the consequences of misspecifying dependence.

Common pitfall: “The fitted correlation is negative, so higher clearance causes lower sensitivity.” This reverses the meaning of a random-effect association. A negative $\rho$ describes a population tendency for higher clearance effects to accompany lower sensitivity effects under the model; it does not establish causation. Also avoid concluding that $\rho=0$ because a small sample of shrunken individual estimates looks uncorrelated. Evaluate the joint-model estimate and its uncertainty, and state the association in terms of the modeled random effects.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

If $\omega_{CL}=0.25$, $\omega_S=0.30$, and $\rho=0.40$, what is the covariance of the two random effects, and what does its sign mean?
For $Q=S\,AUC$, with $AUC$ inversely proportional to clearance, explain whether a positive clearance–sensitivity correlation increases or decreases $\operatorname{Var}(\log Q)$, relative to independence.
Why is the sample correlation of empirical Bayes estimates not generally an adequate estimate of the population random-effect correlation?
Name one way joint estimation can be more informative than a plug-in sequential fit and one reason its correlation estimate may nevertheless be unreliable.

Solutions. (1) The covariance is $(0.40)(0.25)(0.30)=0.0300$, in the product of the two log-effect scales. Its positive sign means above-average clearance effects tend to accompany above-average sensitivity effects. (2) Since $\operatorname{Var}(\log Q)=\omega_S^2+\omega_{CL}^2-2\rho\omega_{CL}\omega_S$, positive correlation decreases this variance relative to independence. (3) Empirical Bayes estimates are uncertain and usually shrunk toward zero, so their sample correlation need not represent the latent population correlation. (4) Joint estimation can use PK and PD observations together while representing covariance and propagating uncertainty; its estimate can still be unreliable when the data do not separately inform individual clearance, sensitivity, and their association.