Why the Order of Calculation Matters
Suppose a regression model is fitted in each of several multiply imputed data sets, and the scientific quantity of interest is a nonlinear function of its coefficients. Two possible workflows may appear reasonable: first pool the coefficients and then calculate the function, or calculate the function in each imputation and then pool those results. In general, they answer different numerical questions. The distinction matters for the point estimate and for its uncertainty.
For a linear contrast, pooling coefficients first is compatible with calculating the contrast in each imputation; this is the subject of the earlier material on linear contrasts of pooled coefficients. Nonlinear functions do not have that general property. A sound analysis therefore begins by defining the estimand: is the target the nonlinear function of a common parameter, or a quantity defined on a transformed scale? Then the estimates and their within-imputation variances must be formed for that target consistently.
For scalar \(Q_j\), the pooled estimate is \(\bar{Q}=m^{-1}\sum_{j=1}^{m}Q_j\). The average within-imputation variance is \(\bar{U}=m^{-1}\sum_{j=1}^{m}U_j\), and the between-imputation variance is \(B=(m-1)^{-1}\sum_{j=1}^{m}(Q_j-\bar{Q})^2\). The total variance is \(T=\bar{U}+(1+1/m)B\), with pooled standard error \(\sqrt{T}\). The factor \(1+1/m\) accounts for finite Monte Carlo sampling of the imputations in the pooled variance.
These rules require a common estimand across imputations, proper imputations that represent uncertainty about missing values, and a suitable complete-data analysis in each imputation. In particular, the completed-data estimates \(Q_j\) should be approximately unbiased for the same \(q\), and their estimated variances \(U_j\) should reflect the sampling uncertainty of those estimates. For confidence intervals and tests, the usual Rubin small-sample \(t\) reference distribution also relies on an adequate approximation to normality for the complete-data estimate of \(q\). The procedure cannot repair a poorly specified imputation model, an incompatible analysis model, or an unsuitable variance estimate.
Nonlinear Functions Do Not Generally Commute With Pooling
Write \(\bar{\boldsymbol{\theta}}=m^{-1}\sum_j\hat{\boldsymbol{\theta}}_j\). Pooling the coefficient estimates and then transforming gives \(g(\bar{\boldsymbol{\theta}})\). Pooling the imputation-specific nonlinear estimates gives \(m^{-1}\sum_j g(\hat{\boldsymbol{\theta}}_j)\). These expressions are generally unequal because nonlinear functions do not preserve averages. For a convex function, Jensen’s inequality gives \(g(\bar{\boldsymbol{\theta}})\leq m^{-1}\sum_jg(\hat{\boldsymbol{\theta}}_j)\), when the function and the estimates lie in a domain where convexity applies. That inequality is not universal for non-convex functions; the general point is that equality cannot be assumed.
The same distinction applies to uncertainty. Once \(Q_j\) has been defined, \(U_j\) must estimate the complete-data variance of that quantity in imputation \(j\). It is not generally correct to pool the original coefficient variances and then transform the result by applying the same function used for the point estimate. The completed-data analysis must supply a suitable estimate of \(U_j\), using a valid variance calculation for the nonlinear estimand. This tutorial does not prescribe one universal variance method: what is appropriate depends on the model, the function, and the estimand.
Worked Example: Pooling a Squared Estimate
Consider four hypothetical imputations. A positive parameter \(\theta\) has estimates \(0.20,\ 0.30,\ 0.40,\ 0.50\). Suppose the target is \(q=\theta^2\), and the completed-data analyses provide valid estimates of the within-imputation variances for this target: \(0.0016,\ 0.0025,\ 0.0036,\ 0.0049\). The estimates are assumed to come from sufficiently large completed-data samples for the squared estimates to have an adequate approximate normal sampling distribution.
First calculate the target estimate in each imputation: \(Q_j=\hat{\theta}_j^2\), giving \(0.0400,\ 0.0900,\ 0.1600,\ 0.2500\). Then
The average within-imputation variance is \(\bar{U}=(0.0016+0.0025+0.0036+0.0049)/4=0.00315\). The mean of the four \(Q_j\) values is \(0.1350\), so
Thus \(T=0.00315+(1+1/4)(0.008300)=0.013525\), and the pooled standard error is \(\sqrt{0.013525}=0.1163\), in squared units of \(\theta\). By contrast, pooling \(\theta\) first gives \(\bar{\theta}=0.3500\), and squaring gives \(0.1225\), not \(0.1350\). The difference is not an arithmetic error: the two workflows average different quantities.
The variance values in this example are stated as completed-data variance estimates for \(q\); they are not obtained by simply squaring or otherwise transforming the variances of \(\hat{\theta}\). In an applied analysis, the validity of those \(U_j\) values must be justified for the chosen nonlinear estimand.
When Transformation After Pooling Is Valid
There are important exceptions to the general warning. First, an affine transformation \(g(\theta)=a+b\theta\), for fixed constants \(a\) and \(b\), commutes exactly with averaging: \(m^{-1}\sum_j(a+b\hat{\theta}_j)=a+b\bar{\theta}\). Its within-imputation variance is \(b^2U_j\). As a result, applying this transformation before or after Rubin pooling gives the same pooled point estimate and total variance, provided the same complete-data estimates and variances are used.
Worked Example: An Affine Contrast Gives the Same Result Either Way
Suppose three hypothetical imputations produce estimates \(-0.40,\ 0.10,\ 0.50\) for \(\theta\), with within-imputation variances \(0.040,\ 0.050,\ 0.030\). Let \(q=5+2\theta\). Pooling \(\theta\) first gives \(\bar{\theta}=(-0.40+0.10+0.50)/3=0.06667\), so the transformed estimate is \(5+2(0.06667)=5.133\).
Transforming separately gives \(Q_j=4.20,\ 5.20,\ 6.00\), whose mean is also \(5.133\). For \(\theta\), \(\bar{U}=0.0400\). Its between-imputation variance is
Therefore \(T=0.0400+(1+1/3)(0.2033)=0.3111\). On the \(q\) scale, the within-imputation variances are \(4U_j\), so their mean is \(0.1600\), and the between-imputation variance is \(4B=0.8133\). Pooling directly gives \(0.1600+(1+1/3)(0.8133)=1.244\), which equals \(4T=4(0.3111)=1.244\), apart from rounding. Both the point estimate and the variance transform consistently.
A second useful situation is reporting a monotone transformation of a pooled estimate and its confidence interval. If a confidence interval has adequate coverage for a scalar parameter \(\theta\), then applying a continuous, one-to-one, monotone function \(g\) to both endpoints gives an interval for \(g(\theta)\) with the same coverage. For a decreasing function, the transformed endpoints must be reordered. This supports reporting a transformed pooled estimate with a transformed interval when the target is explicitly \(g(\theta)\). It does not make the transformed point estimate equal to the average of \(g(\hat{\theta}_j)\), nor does it justify transforming a standard error as if the function were affine.
Worked Example: Transforming a Pooled Estimate and Interval
Suppose four hypothetical imputations estimate a parameter \(\lambda\) on an analysis scale where the target is \(\tau=\exp(\lambda)\). The estimates are \(0.10,\ 0.20,\ 0.00,\ 0.30\), and each within-imputation variance is \(0.0400\). Here the intended reporting procedure is to pool \(\lambda\), then exponentiate its estimate and confidence limits.
The pooled estimate is \(\bar{\lambda}=(0.10+0.20+0.00+0.30)/4=0.1500\), and \(\bar{U}=0.0400\). The between-imputation variance is
The total variance is \(T=0.0400+1.25(0.01667)=0.06083\), so \(SE=\sqrt{0.06083}=0.2466\). For illustration, use the usual Rubin degrees of freedom \(\nu=(m-1)(1+1/r)^2\), where \(r=(1+1/m)B/\bar{U}\). Here \(r=1.25(0.01667)/0.0400=0.5208\), giving \(\nu=3(1+1/0.5208)^2=25.58\). The two-sided 95% interval on the \(\lambda\) scale is approximately \(0.1500\pm2.059(0.2466)=[-0.3578,\ 0.6578]\).
Exponentiating gives the estimate \(\exp(0.1500)=1.162\) and interval \([\exp(-0.3578),\exp(0.6578)]=[0.6992,\ 1.930]\), in the units of \(\tau\). This interval inherits the coverage of the interval on the \(\lambda\) scale, subject to the adequacy of that interval. The average of the four separately exponentiated estimates is about \(1.169\), not \(1.162\); transforming the pooled estimate and pooling transformed estimates are distinct point-estimation procedures.
Nonlinear Contrasts of Several Coefficients
For a vector of coefficients, the same logic applies with a function \(g(\boldsymbol{\theta})\). In imputation \(j\), compute \(Q_j=g(\hat{\boldsymbol{\theta}}_j)\), then obtain a within-imputation variance \(U_j\) for that scalar contrast. Pool these \(Q_j\) and \(U_j\). The covariance between coefficients may be essential to \(U_j\), even though the final estimand is scalar. Using the separate coefficient variances while ignoring their covariance can produce an incorrect variance for the contrast.
A practical sequence is to define the scientific quantity and its scale before looking at the pooled output, calculate its estimate within each completed-data analysis, and retain the corresponding variance estimate. Then apply the scalar pooling rules and report the scale, estimate, standard error or interval, and interpretation. This approach makes it possible to check that every imputation is estimating the same quantity.
Write the nonlinear function \(q=g(\boldsymbol{\theta})\) and specify its units and parameterization.
Calculate \(Q_j\) and a valid complete-data variance estimate \(U_j\) for that same target.
Use \(\bar{Q}\), \(\bar{U}\), \(B\), and \(T\); verify that the estimates and variances are finite and the estimand is consistent across imputations.
If transforming a pooled estimate or interval, state that choice and ensure it matches the defined target and inferential procedure.
Diagnostics, Alternatives, and Limitations
Before pooling, inspect the imputation-specific estimates and variances. A few extreme \(Q_j\) values can dominate \(B\), particularly for functions that grow rapidly or behave irregularly near a boundary. Non-finite values, estimates outside a function’s domain, or wildly varying \(U_j\) values are warnings to investigate the model fit, convergence, data support, and variance calculation. A smooth-looking pooled answer does not make these problems disappear.
| Approach | What it estimates | When it is appropriate |
|---|---|---|
| Calculate \(g(\hat{\boldsymbol{\theta}}_j)\) in each imputation and pool | The common nonlinear target, with between- and within-imputation uncertainty represented on its scale | Usually preferred when the target itself is the nonlinear contrast and valid \(U_j\) values are available |
| Pool coefficients, then calculate \(g(\bar{\boldsymbol{\theta}})\) | The nonlinear function evaluated at the pooled coefficient estimate | Exact agreement for affine functions; may be appropriate for a clearly defined transformed-scale reporting procedure, but is not generally the pooled transformed estimate |
| Pool on a scale with a valid interval, then transform the estimate and interval | A transformed parameter, with inference inherited from the original-scale interval | Useful for a one-to-one monotone transformation when the interval on the original scale is valid and the target is explicitly \(g(\theta)\) |
The first approach depends on having a defensible complete-data variance for the nonlinear quantity. If that variance is unavailable or unreliable, simply transforming a pooled variance is not a general substitute. A different inferential method may be needed, or the analysis may need to be framed on a scale for which valid estimates and intervals are available. The introductory material on Rubin’s rules supplies the pooling framework; the central issue here is choosing the right quantity to enter that framework.
Common Mistakes
- Assuming every function commutes with averaging. The equality holds for affine functions, not generally for nonlinear ones. Check the algebra or calculate both quantities in a small example.
- Transforming the estimate but not its uncertainty coherently. A nonlinear point transformation does not justify applying that function to a variance or standard error. Use a valid within-imputation variance for the target, or a justified transformed-interval procedure.
- Ignoring covariance in a multivariable contrast. A nonlinear function of several coefficients can depend on their joint sampling variability. Use a variance calculation that accounts for the required covariance information.
- Changing the estimand without saying so. The mean of imputation-specific transformed estimates and the transformed pooled estimate are not generally the same. State which target and scale are being reported.
A useful rule is: pool the quantity whose estimate and uncertainty you intend to interpret. Transforming after pooling is defensible in special cases, notably affine transformations and monotone transformations of a valid pooled interval. Outside those cases, calculate the nonlinear quantity within each imputation and pool its estimates using appropriate within-imputation variances.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Three imputation-specific estimates of \(\theta\) are \(1,\ 2,\ 3\). For \(g(\theta)=\theta^2\), calculate the mean of the transformed estimates and the transformed pooled estimate.
- For \(q=4-3\theta\), does transforming before pooling give the same pooled estimate as transforming after pooling? What happens to each within-imputation variance?
- What must be supplied in each imputation before Rubin pooling a nonlinear contrast of two regression coefficients?
- A valid 95% confidence interval for \(\theta\) is \([2,5]\), and the target is \(\exp(\theta)\). Give the transformed interval and state a condition needed for this operation to preserve coverage.
Solutions. (1) The transformed estimates are \(1,\ 4,\ 9\), whose mean is \(14/3=4.667\). The pooled estimate of \(\theta\) is \(2\), and its square is \(4\); the values differ. (2) Yes. For an affine function, \(4-3\bar{\theta}\) equals the mean of \(4-3\hat{\theta}_j\). Each within-imputation variance is multiplied by \((-3)^2=9\); the constant adds no variance. (3) Calculate the nonlinear contrast estimate \(Q_j\) and a valid within-imputation variance \(U_j\) for that same scalar contrast, accounting for relevant coefficient covariance. (4) The interval is \([\exp(2),\exp(5)]\), approximately \([7.389,\ 148.4]\). Coverage is preserved when the transformation is continuous and one-to-one and monotone over the interval, and the original interval is valid for the parameter.