Uncertainty for Treatment Contrasts Under Selection Models

Why treatment contrasts need joint uncertainty

A treatment contrast under a selection model is usually a comparison of treatment-specific outcome means in the full target population, not merely a comparison among participants whose outcomes were observed. A selection model describes both the outcome distribution and the probability that an outcome is observed. The response-model parameters may therefore affect how well the outcome parameters are estimated, even when those response parameters do not appear directly in the treatment contrast.

This distinction matters for uncertainty calculations. Treating response parameters as known, or calculating a contrast standard error from only one part of the fitted likelihood, can understate uncertainty. The appropriate starting point is the covariance matrix for the jointly estimated outcome and response parameters. A Wald standard error then follows by propagating that covariance through the treatment-contrast function.

Joint likelihood and the treatment contrast

Let $Y_i$ be a possibly missing outcome for subject $i$, $A_i$ the treatment assignment, $X_i$ observed covariates, and $R_i=1$ if $Y_i$ is observed and $R_i=0$ otherwise. Write the outcome model as $f(y\mid A_i,X_i;\beta)$ and the selection model as $\pi(y,A_i,X_i;\gamma)=P(R_i=1\mid y,A_i,X_i)$. The selection model allows response probability to depend on the outcome, which is the feature that can represent missing not at random (MNAR) response.

Assume subjects are independent, the covariates and treatment assignments are observed, and the specified outcome and response models are adequate for the analysis. For an observed outcome, the observed-data likelihood contribution is the joint probability of the outcome and response. For a missing outcome, the unknown outcome is integrated out. Thus the contributions are:

L_i(\beta,\gamma)= \begin{cases} f(Y_i\mid A_i,X_i;\beta)\,\pi(Y_i,A_i,X_i;\gamma), & R_i=1,\\[4pt] \displaystyle\int f(y\mid A_i,X_i;\beta)\{1-\pi(y,A_i,X_i;\gamma)\}\,dy, & R_i=0. \end{cases}

Let $\theta=(\beta^\mathsf{T},\gamma^\mathsf{T})^\mathsf{T}$ collect all parameters, and let $\hat{\theta}$ be the joint maximum-likelihood estimate. Define the target contrast as $\Delta=g(\theta)$. For example, if $m_a(\beta)$ is the full-data mean outcome under treatment $a$, then $\Delta=m_1(\beta)-m_0(\beta)$. The target covariate distribution, if standardization is used, must be specified as part of the estimand.

Definition: A joint-likelihood Wald standard error for a scalar contrast $\Delta=g(\theta)$ uses the covariance matrix of the jointly estimated parameters. If $\widehat V$ estimates that covariance matrix and $\nabla g$ is the gradient of $g$, evaluated at $\hat{\theta}$, then $$ \widehat{SE}(\hat{\Delta})= \sqrt{\nabla g(\hat{\theta})^\mathsf{T}\widehat V\,\nabla g(\hat{\theta})}. $$ For a confidence level of $100(1-\alpha)\%$, the corresponding Wald interval is $\hat{\Delta}\pm z_{1-\alpha/2}\widehat{SE}(\hat{\Delta})$, where $z_{1-\alpha/2}$ is the standard normal quantile.

From the joint Hessian to a contrast interval

Let $\ell(\theta)$ be the observed-data log likelihood and let $J(\hat{\theta})=-\nabla^2\ell(\hat{\theta})$ be the observed information matrix, evaluated at the joint maximum. Under regularity conditions—including an interior maximum, local identifiability, and a sufficiently accurate likelihood and Hessian—the usual local covariance estimate is $\widehat V=J(\hat{\theta})^{-1}$. The likelihood must be on the total-sample scale: if software reports a Hessian for the average log likelihood, it must be multiplied by the sample size before it is inverted for this covariance calculation.

The delta-method formula follows by a first-order expansion. Near $\hat{\theta}$, the contrast changes approximately by its gradient multiplied by the parameter change. Substituting the covariance of that parameter change gives the variance expression:

g(\theta)\approx g(\hat{\theta})+ \nabla g(\hat{\theta})^\mathsf{T}(\theta-\hat{\theta}), \qquad \widehat{\operatorname{Var}}\{g(\hat{\theta})\} = \nabla g(\hat{\theta})^\mathsf{T}\widehat V\nabla g(\hat{\theta}).

For a linear contrast, the gradient is constant and this calculation is exact for the covariance approximation. For a nonlinear contrast, such as a difference of exponentiated parameters, it is a first-order approximation and works best when the likelihood is locally close to quadratic and the contrast function is not changing too sharply over the relevant parameter range.

A key point is that the full inverse information matrix—not the inverse of only the outcome-parameter block—is used. If the contrast is a function of outcome parameters alone, its gradient entries for response parameters are zero. Nevertheless, response-model uncertainty can still affect its standard error because the outcome-parameter covariance block of the full inverse information reflects joint estimation and parameter dependence.

Worked calculations

Worked Example: A Difference in Arm Means

Suppose a hypothetical study estimates treatment-arm means $\mu_0$ and $\mu_1$, in mmHg, together with a response parameter $\gamma$. The fitted treatment contrast is $\hat{\Delta}=\hat{\mu}_1-\hat{\mu}_0=2.40$ mmHg. Suppose the joint fit gives this covariance matrix, with parameters ordered as $(\mu_0,\mu_1,\gamma)$:

\widehat V= \begin{pmatrix} 0.16&0.04&0.06\\ 0.04&0.25&0.10\\ 0.06&0.10&0.36 \end{pmatrix}.

The gradient for $\Delta=\mu_1-\mu_0$ is $(-1,1,0)^\mathsf{T}$. Therefore,

\widehat{\operatorname{Var}}(\hat{\Delta}) =(-1,1,0)\widehat V(-1,1,0)^\mathsf{T} =0.16+0.25-2(0.04)=0.33, \qquad \widehat{SE}(\hat{\Delta})=\sqrt{0.33}=0.5745\text{ mmHg}.

A 95% Wald interval is $2.40\pm1.96(0.5745)$, or approximately $(1.274,3.526)$ mmHg. Under the model and the Wald approximation, the interval estimates the full-data mean difference between treatment arms. Its endpoints are in mmHg. The response parameter has zero direct gradient here, but the covariance used is from the joint fit and includes the effects of jointly estimating response and outcome parameters.

As a simple check on the calculation, omitting the covariance between the two arm means would give variance $0.16+0.25=0.41$, not $0.33$. The covariance between estimates is part of a contrast variance and cannot generally be discarded.

Worked Example: A Nonlinear Contrast of Positive Means

Consider a hypothetical positive-outcome model parameterized so the full-data mean in arm $a$ is $m_a=\exp(\beta_a)$. Let the response parameter be $\gamma$, and suppose the joint estimates are $\hat{\beta}_0=\log(10)$, $\hat{\beta}_1=\log(12)$, and $\hat{\gamma}=0.20$. The joint covariance matrix, ordered as $(\beta_0,\beta_1,\gamma)$, is:

\widehat V= \begin{pmatrix} 0.010&0.004&0.003\\ 0.004&0.014&0.006\\ 0.003&0.006&0.040 \end{pmatrix}.

The mean difference is $g(\theta)=\exp(\beta_1)-\exp(\beta_0)$, so the fitted contrast is $12-10=2.00$ outcome units. Its gradient is $\nabla g=(-\exp(\beta_0),\exp(\beta_1),0)^\mathsf{T}$, which at the estimates is $(-10,12,0)^\mathsf{T}$. The delta-method variance is:

\widehat{\operatorname{Var}}(\hat{\Delta}) =(-10)^2(0.010)+12^2(0.014)+2(-10)(12)(0.004) =1.000+2.016-0.960=2.056.

Thus $\widehat{SE}(\hat{\Delta})=\sqrt{2.056}=1.434$ outcome units. The 95% Wald interval is $2.00\pm1.96(1.434)$, or approximately $(-0.810,4.810)$ outcome units. It is wide relative to the estimated difference, indicating substantial uncertainty. As in the first example, $\gamma$ has zero direct derivative in this estimand, but its joint estimation is reflected in the covariance for the outcome parameters. The interval also relies on a local linear approximation to the nonlinear contrast.

Worked Example: Seeing the Cost of Fixing the Response Parameter

Suppose a hypothetical fitted selection model has two parameters: a treatment contrast $\beta$, measured in score points, and a response parameter $\gamma$. The observed information at the joint maximum is:

J= \begin{pmatrix} 25&10\\ 10&16 \end{pmatrix}, \qquad \hat{\beta}=1.10.

The determinant is $25(16)-10(10)=300$, so the joint covariance estimate is:

J^{-1} =\frac{1}{300} \begin{pmatrix} 16&-10\\ -10&25 \end{pmatrix}.

The variance estimate for $\hat{\beta}$ is the upper-left entry, $16/300=0.05333$, giving $SE=\sqrt{0.05333}=0.2309$ score points. The 95% Wald interval is $1.10\pm1.96(0.2309)$, or approximately $(0.647,1.553)$ score points.

If $\gamma$ were fixed rather than estimated, the curvature for $\beta$ alone would give variance $1/25=0.04000$ and standard error $0.2000$. That is not the joint-fit standard error: it omits uncertainty from estimating the response parameter. The example illustrates why the full information matrix is needed, even when the reported contrast is itself $\beta$.

Why nuisance-parameter uncertainty enters

The role of the response parameters can be seen from a block partition of the observed information. Write its blocks as $J_{\beta\beta}$, $J_{\beta\gamma}$, $J_{\gamma\beta}$, and $J_{\gamma\gamma}$, with $\beta$ denoting outcome parameters and $\gamma$ response parameters. When the relevant inverses exist, the outcome-parameter block of the joint covariance is:

\widehat V_{\beta\beta} = \left(J_{\beta\beta} -J_{\beta\gamma}J_{\gamma\gamma}^{-1}J_{\gamma\beta}\right)^{-1}.

The subtracted term accounts for the association between the outcome and response parameters in the likelihood. If one instead fixes $\gamma$, the covariance is based on $J_{\beta\beta}^{-1}$. These are generally different. In a well-behaved regular problem, allowing nuisance parameters to be estimated typically reduces the effective information about the parameters of interest; the block formula makes the source of that adjustment explicit.

The gradient must also match the estimand. If the target is a marginal mean contrast defined solely by outcome-model parameters, its response-parameter gradient is zero. If a specified estimand is a function of both outcome and response parameters, the gradient must include its response-parameter derivatives as well. Do not add a separate “response variance” term by hand: the gradient formula using the full covariance already combines direct dependence and joint estimation.

Diagnostics, alternatives, and limitations

A Wald interval is a local approximation, not a universal guarantee. Before reporting one, check that the joint maximum is interior and that the observed information is numerically stable. A nearly singular Hessian, very large standard errors, strong sensitivity to starting values, or substantial changes under modest numerical-integration refinements can signal weak identification or an unreliable curvature estimate. In those circumstances, the resulting standard error and interval may be misleading even if the arithmetic is correct.

Also check whether the interval is plausible for the contrast. A symmetric Wald interval can extend beyond the possible range of an outcome or produce a negative lower bound for a nonnegative mean. For a nonlinear contrast, a gradient-based interval may be particularly inaccurate if the likelihood is skewed, the estimate is near a boundary, or a likelihood ridge makes the local quadratic approximation poor. An interval that is unexpectedly narrow is not evidence by itself that the selection uncertainty is negligible; verify that the joint covariance—not a conditional or partial covariance—was used.

The main alternatives answer related but not identical inferential needs. A profile-likelihood interval evaluates likelihood support over fixed values of the contrast while optimizing nuisance parameters, and can better represent asymmetry or non-quadratic curvature. A parametric bootstrap uses simulated data under the fitted joint model to approximate the sampling distribution of the contrast. It can reflect nonlinear estimation behavior, but inherits dependence on the fitted model and may be computationally demanding. These methods complement a Wald calculation; they do not remove the need to assess model fit, identifiability, and numerical accuracy.

Common pitfall: “The contrast only involves the outcome parameters, so I can invert the outcome block of the Hessian and ignore the response parameters.” This treats the response parameters as fixed. For joint estimation, invert the full observed information and then use the relevant entries of that joint covariance in the contrast variance. Zero gradient entries for response parameters do not make their joint estimation irrelevant.

Check Your Understanding

Answer each question using the joint-likelihood and contrast calculations above. Solutions follow the questions.

For $\Delta=\mu_1-\mu_0$, why does the gradient have a zero entry for a response parameter, and why can the response model still affect the standard error?
If the two estimated arm means have variances $0.09$ and $0.16$, and covariance $0.03$, calculate the variance and standard error of their difference.
A fitted scalar contrast is $0.80$ units with standard error $0.25$ units. Calculate its 95% Wald interval.
Name two conditions or diagnostics that should make an analyst cautious about a joint-likelihood Wald interval.
When might a profile-likelihood or parametric-bootstrap interval be preferable to a Wald interval?

Solutions. (1) The contrast function does not directly include the response parameter, so its derivative with respect to that parameter is zero. But the joint covariance of the outcome estimates comes from the inverse full information matrix and reflects their joint estimation with the response parameters. (2) The variance is $0.09+0.16-2(0.03)=0.19$, so the standard error is $\sqrt{0.19}=0.4359$ units. (3) The interval is $0.80\pm1.96(0.25)=0.80\pm0.490$, or $(0.310,1.290)$ units. (4) Examples include a nearly singular Hessian, a boundary estimate, instability across starting values, or sensitivity to numerical integration. (5) They may be preferable when the likelihood is asymmetric or poorly approximated by a local quadratic, or when the contrast is nonlinear; each still relies on the specified selection model and a reliable fit.