Joint Likelihood Factorization for PK/PD Data

Why factor a joint PK/PD likelihood?

In a joint PK/PD analysis, concentration observations help estimate the pharmacokinetic (PK) model, while response observations inform the pharmacodynamic (PD) model. The two sources of information are connected: the PD model often uses the concentration trajectory predicted by the PK model, and individual characteristics may affect both PK and PD. A likelihood must represent this connection without counting observations twice or treating uncertain individual effects as known.

The key construction is conditional. Given an individual's latent random effects and model parameters, write a density for the concentration data and a density for the response data. Under conditional independence of their observation errors, multiply those densities. Then integrate over the individual's random-effects distribution. The resulting marginal likelihood is a joint model for the observed data.

Definition: For individual $i$, let $Y_i^C=(Y_{i1}^C,\ldots,Y_{in_i^C}^C)$ denote the observed concentrations and $Y_i^R=(Y_{i1}^R,\ldots,Y_{in_i^R}^R)$ the observed responses. Let $b_i=(b_i^{PK},b_i^{PD})$ be the individual's latent random effects, with joint density $g(b_i;\Omega)$. The conditional joint likelihood contribution is the density $p(Y_i^C,Y_i^R\mid b_i;\theta)$. The marginal individual likelihood contribution is $$ L_i(\theta,\Omega)=\int p(Y_i^C,Y_i^R\mid b_i;\theta)\,g(b_i;\Omega)\,db_i. $$ Here $\theta$ collects fixed PK, PD, and residual-error parameters; $\Omega$ describes the random-effects distribution.

Building the conditional contributions

Write the predicted concentration at the $j$th PK sampling time as $m_{ij}^C(\theta_{PK},b_i^{PK})$. This prediction comes from the PK model, dose history, and relevant covariates. Write the predicted response at the $k$th PD observation time as $m_{ik}^R(\theta_{PD},C_i(\cdot),b_i^{PD})$, where $C_i(\cdot)$ is the concentration trajectory implied by the PK model for that individual. The PD prediction may depend on a concentration at a particular time or on a summary of the trajectory, such as an interval exposure measure.

A common illustrative observation model assumes normally distributed residual errors, conditionally independent across observations and between the PK and PD measurements:

Y_{ij}^C\mid b_i \sim N\left(m_{ij}^C,\sigma_C^2\right), \qquad Y_{ik}^R\mid b_i \sim N\left(m_{ik}^R,\sigma_R^2\right).

These are modeling assumptions, not universal properties of PK/PD data. For example, concentrations may be modeled on a log scale, and repeated measurements may require correlated residual errors. The chosen densities must match the scale and dependence structure of the observations. Under the displayed normal model, the individual conditional likelihood factors as

p(Y_i^C,Y_i^R\mid b_i;\theta) = \left\{\prod_{j=1}^{n_i^C} f_C(Y_{ij}^C\mid b_i;\theta)\right\} \left\{\prod_{k=1}^{n_i^R} f_R(Y_{ik}^R\mid b_i;\theta)\right\}.

The PK factor describes concentration observations around the PK predictions. The PD factor describes responses around predictions that use the PK-implied exposure. Although these factors multiply conditionally, they are not generally independent after the random effects are integrated out.

Assumptions for this product form: Given $b_i$ and the model parameters, the specified concentration and response observations have the stated densities, and their residual errors are conditionally independent across the two data types and across observations. The same $b_i$ values are used wherever the model says they affect PK or PD. If residual errors are correlated, replace the product with the appropriate joint density; do not impose a product merely for convenience.

Why integration is part of the likelihood

The random effects are not observed data. For a marginal likelihood based on the measured concentrations and responses, all possible values of $b_i$ must be accounted for, weighted by their joint population distribution. For $N$ independent individuals, the likelihood is therefore

L(\theta,\Omega) = \prod_{i=1}^{N} \int \left[ \prod_{j=1}^{n_i^C} f_C(Y_{ij}^C\mid b_i;\theta) \right] \left[ \prod_{k=1}^{n_i^R} f_R(Y_{ik}^R\mid b_i;\theta) \right] g(b_i;\Omega)\,db_i.

The log likelihood is the sum of the logarithms of these individual integrals. In general, the integral of a product is not the product of the separate integrals. Consequently, one should not usually integrate the PK and PD factors separately and then multiply the resulting marginal contributions. Shared or correlated random effects connect the two factors through the integration.

Justification by conditional probability: For fixed $b_i$, conditional independence gives $p(Y_i^C,Y_i^R\mid b_i)=p(Y_i^C\mid b_i)p(Y_i^R\mid b_i)$. Averaging this conditional joint density over the random-effects distribution gives $$ p(Y_i^C,Y_i^R)=\int p(Y_i^C\mid b_i)p(Y_i^R\mid b_i)g(b_i)\,db_i. $$ The two data types can be dependent marginally because both conditional densities vary with the same $b_i$. Thus conditional factorization does not imply marginal factorization.

Worked Example: Multiplying the Conditional PK and PD Densities

Consider one hypothetical individual with two concentration observations, 8 and 4 mg/L, and one response observation, 6 response units. At a fixed value of the random effects, suppose the PK predictions are 7 and 5 mg/L, the PD prediction is 5 response units, and the residual standard deviations are $\sigma_C=1$ mg/L and $\sigma_R=2$ response units. Assume the normal observation model and conditional independence stated above.

For a normal density with mean $\mu$ and standard deviation $\sigma$, the density at $y$ is $(\sigma\sqrt{2\pi})^{-1}\exp\{-(y-\mu)^2/(2\sigma^2)\}$. The two concentration residuals are $8-7=1$ mg/L and $4-5=-1$ mg/L. Each concentration density is

\frac{1}{1\sqrt{2\pi}}\exp\left(-\frac{1^2}{2(1^2)}\right) =0.24197\;(\text{mg/L})^{-1}.

The response residual is $6-5=1$ response unit. Its density is

\frac{1}{2\sqrt{2\pi}}\exp\left(-\frac{1^2}{2(2^2)}\right) =0.17603\;(\text{response unit})^{-1}.

Thus the conditional joint likelihood contribution at this fixed random-effect value is

0.24197\times 0.24197\times 0.17603 =0.01031\;(\text{mg/L})^{-2}(\text{response unit})^{-1}.

This is a conditional density contribution, not yet the marginal contribution for the individual. The latter requires averaging this product over the random-effects distribution.

Worked Example: Integrating a Shared Random Effect

Suppose, for illustration, an individual has a shared latent effect with two possible states, A and B, each having probability 0.5. The conditional PK likelihood contribution is 0.08 in state A and 0.20 in state B; the conditional PD likelihood contribution is 0.30 in state A and 0.10 in state B. These values are hypothetical density contributions, with the same measurement units understood across the states. Assume conditional independence given the state.

For state A, the conditional joint contribution is $0.08\times0.30=0.024$. For state B, it is $0.20\times0.10=0.020$. Integrating over the two-state random-effects distribution gives

$$ L_i=0.5(0.024)+0.5(0.020)=0.022. $$

If the marginal PK and PD contributions were calculated separately, they would be $0.5(0.08)+0.5(0.20)=0.14$ and $0.5(0.30)+0.5(0.10)=0.20$. Multiplying those gives $0.14\times0.20=0.028$, not 0.022. That multiplication incorrectly combines PK information from one possible random-effect state with PD information from another. The shared state must be integrated jointly.

Conditioning on measured concentrations is a different representation

The probability chain rule also gives $p(Y_i^C,Y_i^R)=p(Y_i^C)p(Y_i^R\mid Y_i^C)$. This identity is exact. But the conditional response distribution in that factorization is not generally the PD density evaluated at measured concentrations while ignoring PK uncertainty. In a joint random-effects model, PK observations update what is plausible about $b_i$; the correct conditional response density accounts for that information. For example,

p(Y_i^R\mid Y_i^C) = \int p(Y_i^R\mid b_i,Y_i^C)\,p(b_i\mid Y_i^C)\,db_i.

Under the conditional observation assumptions used above, the response density given $b_i$ is specified through the PK-implied trajectory, and the distribution $p(b_i\mid Y_i^C)$ reflects the observed concentrations. Simply inserting noisy measured concentrations into a PD equation is a different modeling choice. It can be appropriate in a design where exposure is treated as observed without error, but it is not equivalent to the latent-trajectory joint model.

Worked Example: Effect of Changing Concentration Units

Suppose one hypothetical concentration observation is 8 mg/L, with a normal model having mean 7 mg/L and standard deviation 1 mg/L. Its density is $0.24197\;(\text{mg/L})^{-1}$. Express the same observation in micrograms per liter: the observation is 8000 µg/L, the mean is 7000 µg/L, and the standard deviation is 1000 µg/L.

The density in the new units is

\frac{1}{1000\sqrt{2\pi}} \exp\left(-\frac{(8000-7000)^2}{2(1000)^2}\right) =0.00024197\;(\text{µg/L})^{-1}.

The numerical density is one thousandth as large because the measurement unit is one thousandth as large. The log density shifts by $-\log(1000)=-6.9078$. More generally, if a fixed one-to-one unit conversion transforms each continuous observation, the density changes by the corresponding Jacobian factor. Applying the same conversion consistently adds a parameter-independent constant to the log likelihood. Therefore, parameter estimates and likelihood comparisons are unaffected, even though the numerical log likelihood depends on the units used for the continuous observations.

Practical checks and limitations

A correctly written factorization is only as credible as its observation models and random-effects structure. Check whether residual patterns suggest a variance that changes with concentration or response magnitude, whether repeated residuals show temporal dependence, and whether the PK and PD predictions systematically miss particular time regions. Compare conditional residual checks with checks based on predictions that account for random effects; the two views answer different questions. If residual dependence is present, a covariance model or another joint residual density may be required.

Also check whether the assumed random-effects distribution and the PK-to-PD exposure link are scientifically defensible. A shared latent effect can induce dependence, but it does not automatically represent every relevant source of association. Sparse sampling, weakly informed random-effects dimensions, or a misspecified PK model can make the likelihood difficult to estimate or interpret. The integral may require numerical approximation in realistic nonlinear mixed-effects models, and an optimizer's convergence does not by itself establish that the model is adequate.

The joint likelihood is distinct from a sequential procedure that first estimates PK parameters and then treats fitted exposure as fixed in a PD fit. The joint construction retains both observation contributions in one marginal likelihood. That distinction defines the model being fitted; by itself it does not guarantee better precision, better fit, or reliable inference. Those properties depend on the data, assumptions, and identifiability.

Common mistakes

Multiplying marginal PK and PD likelihoods. This generally loses the dependence induced by shared random effects. Integrate the product of conditional contributions over the joint random-effects distribution.
Treating the PD factor as a function of error-free measured concentration by default. In a latent-trajectory model, PD predictions use the PK-implied exposure, and the PK observations contribute information about individual effects.
Assuming that conditional independence means marginal independence. Conditional independence holds at fixed random effects; integration can make PK and PD observations dependent.
Calling the numerical log likelihood unit-free. Continuous-data densities have units, and numerical log likelihoods change with the units of measurement. A consistent fixed conversion adds a parameter-independent constant.

Key takeaway: Form the PK and PD density contributions conditional on the same individual random effects, multiply them only when conditional independence is justified, and integrate their product over the joint random-effects distribution to obtain the marginal likelihood.

Check Your Understanding

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

What is the difference between a conditional joint likelihood contribution and a marginal individual likelihood contribution?
Why can the product of the separately integrated PK and PD contributions differ from the joint marginal contribution?
If PK and PD residual errors are correlated conditional on random effects, what must change in the likelihood construction?
A concentration is reported in a unit that is 100 times smaller than the original unit. What happens to its density value, and what happens to parameter estimates if the conversion is applied consistently?
Why is inserting measured concentrations directly into the PD model not generally equivalent to using a latent PK trajectory?

Solutions. (1) The conditional contribution is $p(Y_i^C,Y_i^R\mid b_i)$ at specified random effects; the marginal contribution averages that density over $g(b_i;\Omega)$. (2) In general, $\int p(Y_i^C\mid b_i)p(Y_i^R\mid b_i)g(b_i)\,db_i$ is not equal to $\{\int p(Y_i^C\mid b_i)g(b_i)\,db_i\}\{\int p(Y_i^R\mid b_i)g(b_i)\,db_i\}$, because the same random effects enter both factors. (3) Use a joint conditional residual density that represents the covariance; do not multiply independent-error densities. (4) The density is multiplied by 100 for that observation, with a corresponding additive constant in its log density. A fixed consistent conversion does not change parameter estimates. (5) The measured concentrations include observation error, whereas the latent trajectory is implied by the PK model and individual effects; substituting measurements changes the conditioning and generally ignores uncertainty and the PK model structure.