Worked Example: Two-Stage Versus Joint Emax Estimation

Why Compare the Estimation Strategies?

An Emax estimate can depend on how uncertainty in concentration is handled. In a two-stage analysis, PK data are fitted first, and the resulting concentration estimates are treated as fixed when estimating the PD parameters. A joint analysis instead estimates the PK and PD parameters together from their likelihood contributions. The two procedures can therefore produce different estimates and uncertainty, even when they use the same observations and structural Emax relationship.

The earlier material on sequential versus simultaneous estimation establishes why neither strategy is automatically superior: the result depends on the information in each data stream and on the model assumptions. Here, a compact simulated example makes the comparison explicit. It uses three independent subjects, one PK summary observation per subject, and one PD observation per subject. This is a deliberately simplified model for illustrating the estimation calculation, not a general design recommendation.

Model and Assumptions

For subject \(i\), assume the PK observation and PD observation follow

The PK and PD errors are independent conditional on the parameters, and observations from different subjects are independent. The concentration scales satisfy \(A_i>0\), \(E_{\max}>0\), and \(EC_{50}\) is fixed and known. For the numerical example, \(EC_{50}=2\) mg/L, \(\sigma_C=1\) mg/L, and \(\sigma_Y=0.5\) response units. Treating both residual variances as known keeps the uncertainty comparison focused on the consequence of carrying PK uncertainty forward. In an application, residual variances and other model parameters may also need estimation.

When \(C_i^{\mathrm{obs}}>0\), the two-stage PK estimate under this one-observation normal model is \(\hat A_{i,\mathrm{seq}}=C_i^{\mathrm{obs}}\). If \(C_i^{\mathrm{obs}}\leq 0\), the constrained estimate under the closure \(A_i\geq 0\) is 0; under the stated strict constraint \(A_i>0\), no maximum-likelihood estimate is attained. The second stage uses \(f_i=f(\hat A_{i,\mathrm{seq}})\) and fits \(Y_i=E_{\max}f_i+\epsilon_i\), with independent \(\epsilon_i\sim N(0,\sigma_Y^2)\). This is a zero-baseline Emax model: there is no \(E_0\) term.

This is the unconstrained least-squares estimate and, under the stated normal model with fixed \(f_i\), it is the maximum-likelihood estimate when it is positive, as it is for the supplied data. If it is nonpositive, the strict constraint \(E_{\max}>0\) yields no attained maximum-likelihood estimate: the likelihood approaches its supremum as \(E_{\max}\to0^+\). It does not account for uncertainty in the PK estimates when calculating its standard error.

The Joint Likelihood and Its Score Equations

The joint likelihood retains both data sources. Omitting constants that do not depend on the parameters, its negative log-likelihood is

Minimizing this quantity allows PD observations to inform the fitted \(A_i\) values as well as \(E_{\max}\). That is the defining distinction in this comparison; it does not mean the PK data are discarded or that the PD data directly measure concentration.

A short score derivation provides a useful check on a joint fit. Write \(r_i=E_{\max}f(A_i)-Y_i\), the fitted-minus-observed PD residual. Differentiating the negative log-likelihood with respect to \(E_{\max}\) and \(A_i\) gives

At an interior maximum, both derivatives must be zero. In particular, the second equation shows how the PK estimate can shift: a nonzero PD residual contributes to the score for \(A_i\). The PK and PD score contributions must balance, rather than the PK-only score being zero at the final joint estimate.

Comparing Standard Errors and Interval Coverage

For the two-stage estimate, the conditional model-based standard error is \(\sigma_Y/\sqrt{\sum_i f_i^2}\). It treats the PK-derived fractions as fixed. Using the values above,

A joint standard error accounts for the fact that the \(A_i\) values are also estimated. One standard large-sample calculation uses the expected Fisher information, evaluated at the joint fit. Its \(E_{\max}\)-by-\(E_{\max}\) element is \(\sum_i f_i^2/\sigma_Y^2\). The cross-information for \(E_{\max}\) and \(A_i\) is \(E_{\max}f_if'(A_i)/\sigma_Y^2\), and the \(A_i\)-by-\(A_i\) element is \(1/\sigma_C^2+E_{\max}^2f'(A_i)^2/\sigma_Y^2\). Eliminating the nuisance \(A_i\) parameters gives the profile information

This expression is the Schur complement of the \(A_i\) information block: when nuisance parameters are estimated, their information contribution is removed from the information available for \(E_{\max}\). For the joint fit, the first term is \(5.33778\). The three subtracted terms are approximately \(0.9822\), \(0.8621\), and \(0.3531\), giving \(I_{\mathrm{prof}}\approx3.1404\) and

Using a normal approximation, the two-sided 95% Wald intervals are \(\hat E_{\max}\pm1.96SE\). The two-stage interval is \(9.8458\pm1.96(0.4266)\), or approximately 9.010 to 10.682 response units. The joint interval is \(10.000\pm1.96(0.5643)\), or approximately 8.894 to 11.106 response units. The joint interval is wider here because the PK scales are uncertain and are estimated jointly rather than treated as known.

The data were constructed with a true \(E_{\max}\) of 10 response units. Both intervals contain that value, so each has coverage in this particular simulated replicate. That statement is not an estimate of either method’s long-run coverage. Coverage is the proportion of intervals that contain the true parameter across repeated data sets generated under specified conditions. A single interval can include the truth even when a method is poorly calibrated, or miss it even when its long-run coverage is close to 95%.

Diagnostics, Alternatives, and Limitations

A fitted number alone is not enough to establish that either strategy is working well. For a joint fit, check convergence from multiple reasonable starting values, inspect the score or gradient at the reported solution, and verify that fitted PK and PD means are plausible. Examine standardized residuals against fitted values and relevant design variables. Strong patterns, outliers, or residual spread inconsistent with the assumed variances can indicate model misspecification. A nearly flat likelihood profile or an information matrix close to singular indicates weak identification; in that situation, a Wald standard error may be misleading.

For a two-stage analysis, inspect the PK fits and the PD residuals, but also remember what its reported second-stage standard error omits: uncertainty in the first-stage estimates. A narrow conditional interval is not evidence that the full procedure has equivalent repeated-sampling precision or coverage. Resampling or simulation that repeats both stages can assess the full procedure’s uncertainty, provided the resampling or simulation represents the intended sampling process.

Neither method is universally preferable. Joint fitting can propagate concentration uncertainty, but that benefit depends on a credible PK model and sufficient information to estimate the coupled parameters. Two-stage fitting may be operationally convenient, but treating uncertain exposure estimates as fixed can understate uncertainty and can affect estimates when the exposure-response function is nonlinear. The simple example also omits concentration-time profiles, between-subject random effects, estimated residual variances, and uncertainty in \(EC_{50}\). Its numerical standard errors are illustrative large-sample calculations; with only three subjects, their normal approximation should not be taken as a dependable small-sample interval procedure.

Common Mistake

Wrong reasoning: “The two-stage standard error is smaller, so two-stage estimation is more precise.” The calculation above treats the PK fractions as known; it does not measure uncertainty for the complete two-stage procedure. Conversely, the wider joint interval does not prove that joint estimation is always less precise or that the two-stage estimate is biased. The correction is to identify what uncertainty each reported standard error includes, then compare procedures under the same data-generating assumptions—ideally over repeated simulated data sets when evaluating bias, precision, and coverage.