Why Estimation Strategy Matters
A PK/PD analysis uses concentration information to describe exposure and response information to describe drug effect. The estimation strategy determines whether those two sources of information are handled in stages or in one joint model. This matters most when concentrations are noisy or sparse, when PK and PD parameters are correlated, or when a dose decision depends on a target exposure. A sequential analysis can be useful and transparent, but treating estimated concentrations as if they were known can make downstream uncertainty look too small. A simultaneous analysis can carry uncertainty through the model, but it relies on a correctly specified joint likelihood and can be harder to diagnose.
The comparison below uses a hypothetical concentration–response dataset and a deliberately simple PK/PD model. The numerical estimates are illustrative rather than results from a real study. The aim is to show what each strategy estimates, how uncertainty is propagated, and why a difference in dose recommendation should be investigated rather than resolved by choosing whichever estimate is more convenient.
A Shared Model and Its Assumptions
Suppose concentration after a dose is described by a one-compartment, first-order decline model, \(C(t)=C_0 e^{-kt}\), where \(C(t)\) is concentration in mg/L at time \(t\) in hours, \(C_0\) is the concentration at time zero, and \(k\) is the elimination rate constant in h\(^{-1}\). For a simple concentration–response relationship, use a Hill model with Hill coefficient fixed at one:
Here \(E(C)\) is response in response units, \(E_0\) is baseline response, \(E_{\max}\) is the maximum drug-related change, and \(EC_{50}\) is the concentration producing half of \(E_{\max}\). In the examples, the response is an increase, \(E_0=0\), and \(E_{\max}=100\) response units. The response model uses the latent concentration curve, not a concentration measurement assumed to be error-free.
A basic joint model assumes that concentration observations have independent, mean-zero errors with constant variance on the concentration scale, and that response observations have independent, mean-zero errors with constant variance on the response scale, conditional on the model parameters. It also assumes the stated structural models adequately describe the sampling times and the response delay: there is no unmodeled hysteresis, time-varying baseline, or systematic bias. If errors are normally distributed, the corresponding residual likelihoods are Gaussian. These assumptions are working approximations; they need scrutiny, especially for concentrations with proportional error or responses with nonconstant variance.
In a sequential plug-in analysis, the PK fit is obtained first, and predicted concentrations are supplied to the PD fit as fixed values. In a simultaneous fit, the joint likelihood contains both the concentration likelihood and the response likelihood. The same PK parameters determine concentrations in both parts, so response data can inform PK parameters when the model links the two. The joint approach does not automatically make a model correct: it can transmit misspecification from one component into the other.
Worked Example: Concentration Fit and PK Uncertainty
Worked Example: Concentration Fit and PK Uncertainty
Consider five hypothetical concentration observations fitted with the exponential model. A nonlinear least-squares fit gives \(\hat{C}_0=15.0\) mg/L and \(\hat{k}=0.280\) h\(^{-1}\). For this illustrative fit, the unrounded concentration residuals, observed minus fitted, are \(-0.231\), \(-0.362\), \(0.307\), \(-0.924\), and \(0.745\) mg/L. The residuals are small relative to concentrations of several mg/L, but the fourth residual is noticeably larger in magnitude than the others.
The residual sum of squares, using the unrounded residuals, is
There are five observations and two fitted PK parameters, so the residual variance estimate is \(s_C^2=SSE/(5-2)\approx 0.562483\;(\mathrm{mg/L})^2\). For nonlinear least squares, the local covariance estimate is \(s_C^2(J^\mathsf{T}J)^{-1}\), where \(J\) is the matrix of derivatives of fitted concentrations with respect to the fitted parameters, evaluated at the estimate. Suppose the \(k,k\) entry of \((J^\mathsf{T}J)^{-1}\) for this hypothetical fit is \(0.00013636\;(\mathrm{h}^{-1})^2/(\mathrm{mg/L})^2\). Then
Thus the estimated elimination rate is \(0.280\) h\(^{-1}\), with an approximate standard error of \(0.008758\) h\(^{-1}\). This standard error is a local, model-based approximation; it depends on the residual and Jacobian assumptions. It is not obtained by dividing the residual sum of squares by an arbitrary constant. For this calculation, using \(1.69\) instead of the more precise \(1.68745\) would slightly inflate the residual variance and the resulting standard error.
A quick diagnostic is to plot residuals against time and fitted concentration. A pattern of residuals that is persistently positive early and negative late could indicate a misspecified decline curve; a fan shape could indicate proportional rather than constant variance. A single large residual merits checking, but it does not by itself justify removing that observation.
What Changes Between the Two Fits?
In the sequential fit, the PK stage produces \(\hat{C}_0\), \(\hat{k}\), and a fitted concentration curve. If the second stage uses that curve as fixed, it behaves as though its values were known exactly. The estimated PD parameters can still be useful, but their reported uncertainty is conditional on the particular PK fit. A two-stage procedure that instead uses measured concentrations directly has a different issue: measurement error in concentration can distort the fitted concentration–response relationship, particularly when that error is substantial.
In a simultaneous fit, the concentration and response contributions are fitted together. In simplified notation, if \(C_i^{obs}\) is a concentration observation at time \(t_i\), and \(E_j^{obs}\) is a response observation at time \(u_j\), a joint likelihood can be written as
Here \(\theta_{PK}\) and \(\theta_{PD}\) denote PK and PD parameters, and the two error distributions are specified by the model assumptions. Maximizing this likelihood estimates the parameters together. The resulting covariance matrix can include covariance between PK and PD estimates; that covariance is relevant when calculating uncertainty in a quantity such as a target dose.
Worked Example: A Hypothetical Concentration–Response Dataset
Worked Example: A Hypothetical Concentration–Response Dataset
Suppose four hypothetical observations, taken at matched sampling occasions, are \(C^{obs}=1.0,\ 2.0,\ 4.0,\ 6.0\) mg/L and \(E^{obs}=32,\ 49,\ 66,\ 75\) response units. The response rises with concentration and does not exceed the assumed \(E_{\max}=100\) response units. This is a plausible pattern for the stated increasing Emax model, though four observations alone cannot establish that the model assumptions hold.
For a sequential analysis, suppose the PK fit supplies predicted concentrations and the PD fit gives \(EC_{50}=2.0\) mg/L. At a predicted concentration of \(4.0\) mg/L, the fitted response is
This is close to the hypothetical observation of 66 response units at 4.0 mg/L. The calculation is conditional on the sequential PK fit: it does not include uncertainty in the fitted concentration of 4.0 mg/L unless that uncertainty is carried into the second-stage analysis.
Now suppose a simultaneous joint fit to these same data, using the latent concentration curve and the assumed error models, gives \(EC_{50}=2.4\) mg/L and a slightly different PK estimate. The fitted response at the same concentration is \(100(4.0)/(2.4+4.0)=62.50\) response units. That fit can still be reasonable across all observations; the difference illustrates that the joint fit balances concentration and response information together rather than holding the PK curve fixed. The estimates are illustrative fit outputs, not estimates uniquely recoverable from the four rounded pairs alone: an actual fit requires the full sampling times, error model, and parameter constraints.
Worked Example: Propagating Uncertainty to a Dose Prediction
Worked Example: Propagating Uncertainty to a Dose Prediction
Suppose a target response is 70 response units, with \(E_0=0\) and \(E_{\max}=100\). Solving the Emax equation for target concentration gives
Assume bioavailability \(F=1\), a dosing interval \(\tau=12\) h, and an average steady-state concentration approximation \(C_{ss,avg}=F D/(CL\tau)\), where \(D\) is dose in mg and \(CL\) is clearance in L/h. With \(CL=Vk\), \(V=10.0\) L, the dose required under this approximation is
For a sequential estimate with \(k=0.280\) h\(^{-1}\) and \(EC_{50}=2.0\) mg/L, the arithmetic is
For an illustrative simultaneous fit with \(k=0.290\) h\(^{-1}\) and \(EC_{50}=2.4\) mg/L,
The point predictions differ by \(194.9-156.8=38.1\) mg, about 24.3% of the sequential prediction. A difference of this size could affect a regimen decision, but it does not establish which fit is preferable; model adequacy, sampling design, and uncertainty must be assessed.
The delta method approximates variance for a smooth function of estimated parameters by \(g^\mathsf{T}\Sigma g\), where \(g\) is the gradient and \(\Sigma\) is their covariance matrix. For \(D=120k\,EC_{50}(70/30)\), the derivatives are \(\partial D/\partial k=D/k\) and \(\partial D/\partial EC_{50}=D/EC_{50}\). For the sequential estimate, suppose \(SE(k)=0.008758\) h\(^{-1}\), \(SE(EC_{50})=0.30\) mg/L, and the covariance is approximated as zero. Then
For comparison, suppose the simultaneous fit gives \(SE(k)=0.010\) h\(^{-1}\), \(SE(EC_{50})=0.25\) mg/L, and \(\mathrm{Cov}(k,EC_{50})=-0.001\) h\(^{-1}\)mg/L. The covariance is feasible because its magnitude is below \(SE(k)SE(EC_{50})=0.0025\) h\(^{-1}\)mg/L. With derivatives \(D/k=672\) mg per h\(^{-1}\) and \(D/EC_{50}=81.2\) mg per mg/L,
These standard errors are illustrative approximations, not dosing intervals. Their assumptions include a locally adequate normal approximation and the stated parameter covariance. In an applied analysis, a likelihood-based interval or bootstrap that refits the PK/PD model and recalculates the dose can better reflect skewness, parameter constraints, and nonlinear uncertainty.
When Estimates Diverge, and How to Diagnose It
Sequential and simultaneous estimates are most likely to differ when the concentration data poorly identify the PK curve, when response data are informative about exposure, or when the model allows PK and PD parameters to compensate for one another. Sparse sampling around the concentration profile, substantial concentration measurement error, a narrow observed concentration range, and strong curvature in the response relationship can all increase that sensitivity. The estimates may also differ if the sequential analysis uses observed concentrations while the simultaneous analysis models latent concentrations with measurement error; these are not equivalent analyses.
A practical comparison should check the following before interpreting a different target dose:
- Residual structure: inspect concentration and response residuals against time, fitted values, and relevant covariates. Trends or changing spread challenge the assumed error or structural model.
- Parameter precision: examine standard errors, confidence intervals, covariance, and whether parameters are weakly identified or near constraints. A large change in a poorly identified parameter may have little evidential meaning.
- Prediction checks: compare both fits on the same concentration and response observations, and examine predictions in the concentration range relevant to the target. Extrapolation beyond observed exposure can magnify differences.
- Uncertainty propagation: compare plug-in dose uncertainty with an approach that carries PK uncertainty through the PD and dose calculations, such as a parametric bootstrap.
- Sensitivity: vary plausible residual error assumptions and model structures. A dose recommendation that changes substantially under small defensible changes is not robust.
A sequential approach is often easier to communicate and can be appropriate when PK parameters are well estimated independently, PK uncertainty is small relative to PD uncertainty, and the PD analysis propagates relevant uncertainty. A simultaneous approach is attractive when the concentration and response data jointly inform shared parameters or when measurement error in concentration must be represented explicitly. Its limitations are greater computational and modeling demands and the possibility that an incorrect error model or structural link biases both components. Neither approach eliminates the need for sound sampling, appropriate model checks, or cautious extrapolation.
Common Pitfall
A second pitfall is to choose the simultaneous model merely because it is “more complete.” Joint estimation is only as credible as its likelihood and structural assumptions. If a joint fit changes the dose prediction substantially, inspect the residuals, parameter correlations, and predictive performance to understand why. A later analysis of PD residual patterns can help assess response-model adequacy, but this comparison already requires checking that the current joint and sequential fits use comparable data and assumptions.
Check Your Understanding
Answer each question, then compare with the short solutions below.
- What does a sequential plug-in PD analysis condition on, and what uncertainty can it omit?
- For \(E_{\max}=100\) response units, \(E_0=0\), and \(EC_{50}=3.0\) mg/L, what concentration corresponds to 60 response units under the stated Emax model?
- Why can response observations affect a PK parameter estimate in a simultaneous fit?
- If a sequential dose calculation uses \(D=V k\tau EC_{50}E/(E_{\max}-E)\), name two sources of uncertainty that should be considered.
- Give one diagnostic sign that the assumed concentration error model may be inadequate.
Solutions. (1) It conditions on the fitted PK parameters or predicted concentrations as fixed and can omit their uncertainty in PD and dose uncertainty. (2) \(C_{target}=3.0(60)/(100-60)=4.50\) mg/L. (3) The response likelihood depends on PK parameters through the concentration curve, so the joint likelihood combines information from both data types. (4) Uncertainty in \(k\) and \(EC_{50}\), and their covariance if nonzero; uncertainty in other estimated inputs may matter too. (5) For example, residual spread that increases with fitted concentration suggests constant additive variance may be inappropriate.