Tutorials › Optimal Design and Model-Informed Dosing › Correlated Parameter Uncertainty in Emax Models

Dose Optimization Under Uncertainty · Tutorial 217 of 497

Correlated Parameter Uncertainty in Emax Models

Use the full parameter covariance matrix to quantify uncertainty in Emax predictions and target-attaining doses, and recognize when a local approximation needs further checking.

Advanced 11 min read

What You'll Learn

  • Write the Emax response and target-dose functions using baseline, maximum effect, and ED50.
  • Derive the local variance of a predicted response from the full parameter covariance matrix.
  • Calculate how parameter correlations change prediction uncertainty, even when marginal standard errors stay fixed.
  • Propagate covariance into uncertainty for an interior target-attaining dose.
  • Check covariance matrices, parameter scales, and local-approximation limitations.

Why Dependence Among Emax Estimates Matters

An Emax model prediction depends on several estimated parameters at once. If the estimated baseline, maximum effect, and ED50 tend to move together across repeated samples or model fits, uncertainty in a prediction is not determined by their standard errors alone. The direction of that co-movement matters: covariance can make parameter effects reinforce one another or partially cancel.

This tutorial considers uncertainty in the model-predicted mean response and in a dose that attains a specified response target. The focus is the covariance among the three Emax parameter estimates. It is distinct from residual variation in individual outcomes: parameter uncertainty describes uncertainty about the fitted response curve, while residual variation describes scatter around that curve.

Definition: Write the Emax mean response at dose \(D\) as \(f(D;\theta)=B+M D/(K+D)\), where \(B\) is baseline response, \(M\) is the maximum increase above baseline, and \(K\) is ED50, the dose producing half of that maximum increase. The parameter vector is \(\theta=(B,M,K)^\mathsf{T}\). Its estimate is \(\hat{\theta}=(\hat B,\hat M,\hat K)^\mathsf{T}\), with estimated covariance matrix \(\hat{\Sigma}\). Its diagonal entries are the estimated parameter variances; its off-diagonal entries are covariances.

For an increasing Emax curve, the usual parameter conditions are \(M>0\) and \(K>0\). The baseline and maximum response are \(B\) and \(B+M\), respectively; at finite positive doses the mean response lies between them. Parameter units must be tracked: \(B\) and \(M\) have response units, while \(K\) and \(D\) have dose units. The covariance between two estimates has the product of their units.

Propagating Covariance to a Predicted Response

For a fixed dose \(D\), the predicted mean is a function of all three parameters. Its local sensitivity is summarized by the gradient, evaluated at the parameter estimates. A first-order Taylor approximation uses that gradient to translate parameter uncertainty into prediction uncertainty.

$$ \nabla_\theta f(D;\theta) = \begin{pmatrix} 1\\[2pt] \dfrac{D}{K+D}\\[6pt] -\dfrac{MD}{(K+D)^2} \end{pmatrix} \equiv g(D;\theta). $$

The first component says that a one-unit increase in baseline shifts the prediction by one response unit. The second describes sensitivity to the maximum increment, and the third describes sensitivity to ED50. In particular, increasing \(K\) lowers the response at any fixed positive dose when \(M>0\).

Result: If the fitted parameter estimates are approximately distributed around \(\hat{\theta}\) with covariance matrix \(\hat{\Sigma}\), and the response function is sufficiently smooth in the relevant parameter region, the first-order delta-method variance of the estimated mean response is $$ \widehat{\operatorname{Var}}\{f(D;\hat{\theta})\} \approx g(D;\hat{\theta})^\mathsf{T}\hat{\Sigma}\,g(D;\hat{\theta}). $$ This is a local approximation for parameter uncertainty in the mean response, not a prediction variance for a future individual outcome.

To see why the full matrix enters, expand the quadratic form:

$$ \begin{aligned} \widehat{\operatorname{Var}}\{f(D;\hat{\theta})\} \approx{}& g_B^2\hat{\sigma}_{BB} +g_M^2\hat{\sigma}_{MM} +g_K^2\hat{\sigma}_{KK}\\ &+2g_Bg_M\hat{\sigma}_{BM} +2g_Bg_K\hat{\sigma}_{BK} +2g_Mg_K\hat{\sigma}_{MK}. \end{aligned} $$

Here \(g_B,g_M,g_K\) are the gradient components evaluated at \(\hat{\theta}\), and \(\hat{\sigma}_{ij}\) denotes the corresponding covariance entry. The factor of 2 appears because each off-diagonal covariance occurs twice in the matrix product. A covariance term can be positive or negative; its contribution also depends on the signs of the relevant sensitivities.

Worked Response Calculations

Worked Example: A Prediction With Correlated Estimates

Suppose hypothetical estimates are \(\hat B=10\) response units, \(\hat M=20\) response units, and \(\hat K=8\) mg. At \(D=8\) mg, suppose the standard errors are 2 response units, 5 response units, and 4 mg, respectively. Suppose the estimated correlations are \(\rho_{BM}=-0.4\), \(\rho_{BK}=0.2\), and \(\rho_{MK}=0.5\). Calculate the predicted mean and its approximate standard error.

First, the mean response is \(10+20(8)/(8+8)=20\) response units. The gradient at the estimates is \(g=(1,\ 0.5,\ -20(8)/16^2)^\mathsf{T}=(1,\ 0.5,\ -0.625)^\mathsf{T}\).

The diagonal covariance entries are \(4\), \(25\), and \(16\). Converting each correlation to covariance gives \(\hat{\sigma}_{BM}=-0.4(2)(5)=-4\), \(\hat{\sigma}_{BK}=0.2(2)(4)=1.6\), and \(\hat{\sigma}_{MK}=0.5(5)(4)=10\). Thus the approximate variance is \(1^2(4)+0.5^2(25)+(-0.625)^2(16)+2(1)(0.5)(-4)+2(1)(-0.625)(1.6)+2(0.5)(-0.625)(10)\) \(=4+6.25+6.25-4-2-6.25=4.25\) squared response units. The approximate standard error is \(\sqrt{4.25}=2.062\) response units.

If the off-diagonal covariances were incorrectly set to zero, the variance would instead be \(16.5\), with standard error \(\sqrt{16.5}=4.062\) response units. In this hypothetical case, the covariance contributions substantially reduce uncertainty in the mean prediction. This does not mean that correlation always reduces uncertainty.

Worked Example: Same Marginal Errors, Different Correlations

Keep the estimates and standard errors from the previous example, but suppose all three pairwise correlations are \(0.4\). The correlation matrix with all off-diagonal entries \(0.4\) is positive definite, so it is a valid covariance structure for this calculation. Compare the prediction standard error with the one obtained when covariances are ignored.

The covariances are now \(\hat{\sigma}_{BM}=0.4(2)(5)=4\), \(\hat{\sigma}_{BK}=0.4(2)(4)=3.2\), and \(\hat{\sigma}_{MK}=0.4(5)(4)=8\). The diagonal contribution remains \(16.5\). The covariance contributions are \(2(1)(0.5)(4)=4\), \(2(1)(-0.625)(3.2)=-4\), and \(2(0.5)(-0.625)(8)=-5\).

Therefore the variance is \(16.5+4-4-5=11.5\) squared response units, and the standard error is \(\sqrt{11.5}=3.391\) response units. This is smaller than the diagonal-only standard error of \(4.062\), but larger than the \(2.062\) standard error in the first example. The \(B\)-\(M\) covariance increases uncertainty here, while the two covariance terms involving \(K\) decrease it. The result depends on both the covariance signs and the response sensitivities.

Covariance in a Target-Attaining Dose

For a target response \(T\), solve \(T=B+MD/(K+D)\) for \(D\). A finite positive dose exists for an interior target satisfying \(B<T<B+M\), assuming \(M>0\) and \(K>0\). Rearranging gives \((T-B)(K+D)=MD\), so \((T-B)K=[M-(T-B)]D\). Hence the target-attaining dose is

$$ D_T(\theta) = \frac{K(T-B)}{M-(T-B)}. $$

This expression makes clear why target-dose uncertainty can respond differently to covariance than prediction uncertainty. The target dose is a nonlinear function of all three parameters, and its sensitivities can have different signs from those of the predicted response.

Result: Let \(A=T-B\). For a fixed target \(T\) with \(0<A<M\), the gradient of the target-attaining dose is $$ \nabla_\theta D_T = \begin{pmatrix} -\dfrac{KM}{(M-A)^2}\\[6pt] -\dfrac{KA}{(M-A)^2}\\[6pt] \dfrac{A}{M-A} \end{pmatrix}. $$ The first-order variance approximation is \(\widehat{\operatorname{Var}}(D_T)\approx \nabla_\theta D_T(\hat\theta)^\mathsf{T}\hat{\Sigma}\nabla_\theta D_T(\hat\theta)\), provided the target is feasible near the estimates and the local approximation is adequate.

For completeness, write \(D_T=KA/(M-A)\). Its derivative with respect to \(A\) is \(KM/(M-A)^2\). Since \(A=T-B\), its derivative with respect to \(B\) is the negative of this quantity. Differentiating with respect to \(M\) gives \(-KA/(M-A)^2\), and differentiating with respect to \(K\) gives \(A/(M-A)\). These are the stated gradient components. Substitution into the same covariance quadratic form includes all three pairwise covariance contributions.

Worked Example: Uncertainty in the Dose for a Target

Use the hypothetical estimates and covariance values from the first worked example, and set \(T=20\) response units. The target is between the estimated baseline \(10\) and upper asymptote \(30\), so it has a finite positive target dose. Calculate that dose and its approximate standard error.

Here \(A=T-\hat B=20-10=10\) response units, and \(\hat D_T=8(10)/(20-10)=8\) mg. The gradient is \((-8(20)/10^2,\ -8(10)/10^2,\ 10/10)^\mathsf{T}=(-1.6,-0.8,1)^\mathsf{T}\). Its components have units of mg per response unit, mg per response unit, and mg per mg, respectively.

The diagonal variance contribution is \((-1.6)^2(4)+(-0.8)^2(25)+(1)^2(16)=10.24+16+16=42.24\) mg squared. The covariance contributions are \(2(-1.6)(-0.8)(-4)=-10.24\), \(2(-1.6)(1)(1.6)=-5.12\), and \(2(-0.8)(1)(10)=-16\) mg squared. The total approximate variance is \(42.24-10.24-5.12-16=10.88\) mg squared, giving a standard error of \(\sqrt{10.88}=3.298\) mg.

Ignoring covariance would give a standard error of \(\sqrt{42.24}=6.499\) mg. Under a local normal approximation, a 95% interval would be \(8\pm1.96(3.298)\), or approximately 1.536 to 14.46 mg. It is an approximate interval for the target-attaining dose under the fitted model and covariance approximation, not a statement that there is a 95% probability the fixed dose parameter lies in that interval.

Assumptions, Diagnostics, and Alternatives

The delta-method calculation is most defensible when the parameter estimates are reasonably summarized by their covariance matrix, the response or target-dose function is differentiable in the relevant region, and the function is nearly linear over the uncertainty region. For target doses, feasibility also matters: the target should remain between baseline and the upper asymptote for plausible parameter values. A covariance matrix should be symmetric and positive semidefinite; small negative eigenvalues can arise from numerical rounding, but material violations indicate an invalid or corrupted covariance estimate.

Useful checks include verifying that the reported standard errors square to the diagonal variances, converting correlations to covariances with the correct parameter units, and confirming that the result has appropriate units and a nonnegative variance. It is also informative to compare the full quadratic form with its diagonal-only version. A large difference means dependence materially affects the reported uncertainty; it is not a reason to discard the covariance terms.

A local standard error is not a complete description when the function is strongly curved over plausible estimates, or when a substantial share of parameter uncertainty implies an infeasible target. For example, a target near the estimated upper asymptote can make the target-dose function highly sensitive to small changes in \(M\) or \(B\). In such situations, a symmetric interval based on the delta method can be misleading. The earlier material on Propagating Parameter Uncertainty by Monte Carlo provides a way to examine the transformed distribution by drawing parameter vectors jointly, preserving their covariance. Profile likelihood or bootstrap distributions can also be useful alternatives when their assumptions and fitting procedures are appropriate.

1
Specify the estimand.
Decide whether uncertainty is required for the mean response at a fixed dose or for a target-attaining dose. Do not add individual residual variation unless the question concerns a future outcome.
2
Use the joint covariance estimate.
Check parameter order, standard errors, units, and the off-diagonal covariance entries before computing a gradient quadratic form.
3
Evaluate local adequacy.
Check smoothness, feasibility, curvature, and whether the covariance matrix is valid. If the approximation appears doubtful, inspect a joint simulation or another suitable uncertainty method.

A common mistake is to propagate each parameter standard error separately and then add only the three variance contributions. That is equivalent to replacing the estimated covariance matrix with a diagonal matrix, an assumption of zero covariance. It is not justified simply because the model reports marginal standard errors. The correction is to use the full covariance matrix and include each cross-term with its sign. Another mistake is to interpret a smaller covariance-adjusted standard error as proof that the model is more certain in every relevant sense: it only describes local uncertainty for the particular function and dose being evaluated.

Key takeaway: For Emax predictions and target-attaining doses, parameter uncertainty is a joint quantity. The full covariance matrix and the function’s gradient determine whether dependence increases or decreases local uncertainty.

Check Your Understanding

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

  1. At a fixed positive dose, what are the three sensitivities of the Emax mean response to \(B\), \(M\), and \(K\)?
  2. For \(B=10\), \(M=20\), \(K=8\) mg, and target \(T=25\), calculate the target-attaining dose. Is the target feasible?
  3. In the first worked example, what is the prediction variance if all covariances are set to zero, and what standard error does it imply?
  4. Why does a covariance between \(M\) and \(K\) not by itself determine whether prediction variance increases?

Solutions. (1) The sensitivities are \(1\), \(D/(K+D)\), and \(-MD/(K+D)^2\), respectively. The last is negative for \(M,D,K>0\). (2) \(A=T-B=15\), so \(D_T=8(15)/(20-15)=24\) mg. The target is feasible because \(10<25<30\). (3) The diagonal-only variance is \(4+6.25+6.25=16.5\) squared response units, giving standard error \(\sqrt{16.5}=4.062\) response units. (4) Its contribution is \(2g_Mg_K\hat{\sigma}_{MK}\); it depends on the signs of both gradient components as well as the covariance. Thus the covariance’s sign alone does not establish whether it raises or lowers the variance.

Clinical Biostats