Why Pool the Contrast, Not Just Its Ingredients?
A regression model may estimate coefficients on a scale that is convenient for fitting but not for interpretation. For example, logistic regression coefficients are on the log-odds scale, while a clinical question may concern a ratio of predicted probabilities. That ratio is a nonlinear function of the fitted coefficients. Across multiply imputed data sets, the contrast can vary both because the coefficients vary and because their estimated covariance varies.
The procedure here is to calculate the same contrast in every imputed data set, estimate its variance there, and then apply Rubin’s rules to those contrast estimates and variances on a chosen analysis scale. If the contrast is defined as a log ratio, pooling takes place on the log scale; exponentiation then expresses the pooled estimate and interval as a ratio. The contrast, its scale, and the target population or covariate values must be held fixed across imputations.
For \(M\) imputations, the pooled estimate is \(\bar q=M^{-1}\sum_{m=1}^M q_m\), the average within-imputation variance is \(\bar U=M^{-1}\sum_{m=1}^M U_m\), and the between-imputation variance is \(B=(M-1)^{-1}\sum_{m=1}^M(q_m-\bar q)^2\). The total variance is \(T=\bar U+(1+1/M)B\), with standard error \(\sqrt{T}\). These are scalar Rubin pooling calculations, applied to the contrast on its chosen scale.
Delta-Method Justification for a Logistic Probability Ratio
Consider the logistic model \(\operatorname{logit}(p)=a+bA\), where \(A\) is a binary exposure and \(p\) is the probability of the outcome at a fixed covariate profile. The predicted probabilities for unexposed and exposed individuals are \(p_0=\operatorname{expit}(a)\) and \(p_1=\operatorname{expit}(a+b)\), where \(\operatorname{expit}(t)=1/(1+\exp(-t))\). A probability ratio is \(p_1/p_0\); define the analysis-scale contrast \(g(a,b)=\log(p_1/p_0)\).
For \(f(t)=\log\{\operatorname{expit}(t)\}\), differentiation gives \(f'(t)=1-\operatorname{expit}(t)\). Thus \(g(a,b)=f(a+b)-f(a)\), so the gradient is \((p_0-p_1,\ 1-p_1)^\mathsf{T}\). The delta method approximates the sampling variance of the smooth function \(g(\hat a,\hat b)\) by multiplying this gradient on both sides of the estimated coefficient covariance matrix. This accounts for uncertainty in both coefficients and their covariance.
The approximation relies on a sufficiently regular fitted model and a sampling distribution for the coefficient estimates that is adequately approximated by a multivariate normal distribution. The covariance matrix must estimate the covariance of the coefficients from that imputation, and the probabilities must be defined at the same covariate profile in every imputation. Pooling also relies on proper imputation and a common estimand across imputations; Rubin’s variance combines within-imputation uncertainty with variation among the completed-data estimates. The delta method is an approximation, and can be poor when the contrast is highly curved over the plausible coefficient values or when estimates are unstable.
Worked Example: Delta-Method Variance for One Imputation
Suppose a hypothetical completed data set gives \(\hat a=-1.00\), \(\hat b=0.70\), and coefficient covariance matrix $$ \mathbf{V}= \begin{pmatrix} 0.040 & 0.010\\ 0.010 & 0.090 \end{pmatrix}. $$ The entries are variances and covariance for \(a,b\), respectively. The contrast compares predicted outcome probabilities at \(A=1\) and \(A=0\), at the specified covariate profile.
First calculate \(p_0=\operatorname{expit}(-1.00)\approx0.268941\) and \(p_1=\operatorname{expit}(-1.00+0.70)=\operatorname{expit}(-0.30)\approx0.425557\). Therefore \(q=\log(p_1/p_0)=\log(0.425557/0.268941)\approx0.4589\). The gradient is $$ \nabla g= \begin{pmatrix} p_0-p_1\\ 1-p_1 \end{pmatrix} = \begin{pmatrix} -0.1566\\ 0.5744 \end{pmatrix}. $$
Using \(U=\nabla g^\mathsf{T}\mathbf{V}\nabla g\), the substituted calculation is $$ U=0.040(-0.1566)^2+2(0.010)(-0.1566)(0.5744)+0.090(0.5744)^2 \approx 0.02888. $$ The delta-method standard error is \(\sqrt{0.02888}=0.1699\) log probability-ratio units. The covariance contribution is negative here because the gradient components have opposite signs and the coefficient covariance is positive. This is an estimated within-imputation variance, not yet a multiple-imputation variance.
Pooling and Back-Transforming
The next example pools the nonlinear log probability ratios from three hypothetical imputations. The same outcome, exposure contrast, covariate profile, regression specification, and log-ratio scale are used in each one. The imputation-specific coefficient estimates and covariance matrices have already been used to calculate the probabilities, contrasts, and delta-method variances shown below.
Worked Example: Pooling a Nonlinear Log Probability Ratio
Suppose the three imputations produce the following quantities. The probabilities are predicted for \(A=0\) and \(A=1\), at the same covariate profile in every imputation.
| Imputation | \(p_0\) | \(p_1\) | \(q_m=\log(p_1/p_0)\) | \(U_m\) |
|---|---|---|---|---|
| 1 | 0.2689 | 0.4256 | 0.4589 | 0.02888 |
| 2 | 0.2891 | 0.4750 | 0.4965 | 0.02424 |
| 3 | 0.2497 | 0.3894 | 0.4441 | 0.03263 |
For example, imputation 2 gives \(q_2=\log(0.4750/0.2891)\approx0.4965\). Its \(U_2\) is calculated from that imputation’s coefficient covariance matrix and the gradient evaluated at that imputation’s coefficients, as in the preceding example. The pooled estimate is $$ \bar q=\frac{0.4589+0.4965+0.4441}{3}=0.4665. $$
The average within-imputation variance is $$ \bar U=\frac{0.02888+0.02424+0.03263}{3}=0.02858. $$ The between-imputation variance is $$ B=\frac{(0.4589-0.4665)^2+(0.4965-0.4665)^2+(0.4441-0.4665)^2}{3-1}\approx0.00072976. $$ Consequently, $$ T=0.02858+\left(1+\frac{1}{3}\right)(0.00072976)=0.02955, \qquad SE=\sqrt{0.02955}=0.1719. $$
For an approximate 95% interval, the usual large-sample normal critical value gives \(0.4665\pm1.96(0.1719)\), or approximately \([0.1296,\ 0.8034]\) on the log probability-ratio scale. Exponentiating the estimate gives \(\exp(0.4665)=1.594\). Exponentiating both interval endpoints gives \([\exp(0.1296),\exp(0.8034)]=[1.138,\ 2.233]\). Thus, at the specified covariate profile, the estimated outcome probability for the exposed group is about 1.59 times that for the unexposed group; the approximate 95% interval for this probability ratio is 1.14 to 2.23. This is an interval for the ratio, not a probability statement about the unknown ratio.
When finite-sample degrees of freedom are needed, the scalar Rubin degrees-of-freedom approximation is \((M-1)\{1+\bar U/[(1+1/M)B]\}^2\), subject to the assumptions and adjustments appropriate to the analysis. In the example, this approximation is about \(1845\), so its \(t\) critical value is close to 1.96. The displayed interval uses the large-sample normal approximation for clarity. In a particular analysis, use the applicable degrees-of-freedom method and its critical value rather than treating 1.96 as universal.
Worked Example: Why the Covariance Term Matters
Return to the first imputation in the preceding example. If the off-diagonal covariance were incorrectly discarded, the variance calculation would retain only the two diagonal terms: $$ 0.040(-0.1566)^2+0.090(0.5744)^2 =0.00098094+0.02969418=0.03067512\approx0.03068. $$
The correct variance includes the two equal cross-products, together contributing \(2(0.010)(-0.1566)(0.5744)=-0.001799\). Therefore \(U=0.03068-0.001799=0.02888\), not 0.03068. A covariance matrix is not merely a list of coefficient variances: its covariances can materially affect the uncertainty of a contrast. The sign and size of this effect depend on both the gradient and the covariance.
Checks, Alternatives, and Limitations
Before reporting a pooled nonlinear contrast, check that the estimand is identical across imputations: the same regression definition, contrast, reference values, and analysis scale must be used. Confirm that each \(U_m\) comes from the full covariance matrix and the gradient evaluated at the matching imputation’s estimates. Inspect the imputation-specific contrast estimates and variances for extreme values or unusual instability. Very large or erratic estimates can indicate sparse data, separation in a logistic model, coding differences, or poor local accuracy of the delta approximation. A small number of imputations can also make the between-imputation component and resulting inference unstable.
The direct alternative is to calculate the contrast on its scientifically interpretable scale in each imputation and pool it there if its sampling distribution is sufficiently well behaved and the variance estimate is valid on that scale. For a positive ratio that is skewed, a log scale is often more suitable: pool log ratios and exponentiate afterward. Pooling fitted regression coefficients first and then applying a nonlinear function is generally not equivalent to calculating the nonlinear contrast within each imputation and pooling those contrasts. The latter is the procedure taught here because it targets the same derived quantity in each completed data set and includes its imputation-specific uncertainty.
The delta method is not exact. If the contrast is near a boundary, strongly skewed, or highly curved in the coefficient region supported by the data, the normal approximation on the chosen scale may be inadequate. The transformation must also be one-to-one over the relevant range for straightforward endpoint back-transformation. For a log ratio, positivity of the underlying probabilities makes the ratio positive, and exponentiation maps the log scale to the ratio scale. Back-transform the pooled estimate and its interval endpoints—not the coefficient estimates before constructing the pooled contrast interval.
Check Your Understanding
Answer each question, then compare with the short solutions below.
- In a logistic model, why is the log of a ratio of predicted probabilities nonlinear in the regression coefficients?
- What quantities enter the delta-method variance of a scalar contrast?
- Three imputation-specific log-scale estimates are 0.20, 0.50, and 0.40. Calculate their pooled point estimate.
- Why must the coefficient covariance be included when calculating a contrast variance?
- If a pooled log-ratio estimate is 0.30, what is the corresponding ratio estimate?
Solutions. (1) Each predicted probability is an expit function of a linear predictor, and the log ratio of two such probabilities is therefore generally a nonlinear function of the coefficients. (2) The gradient of the contrast and the full coefficient covariance matrix enter as \(U_m=\nabla g^\mathsf{T}\mathbf{V}_m\nabla g\). (3) The pooled estimate is \((0.20+0.50+0.40)/3=0.3667\). (4) The contrast generally depends on multiple coefficients, whose estimation errors can covary; the cross-products in the quadratic form contribute to its variance. (5) Exponentiate: \(\exp(0.30)=1.350\), approximately.