Why Averaging P-Values Is Tempting—and Misleading
After fitting the same model to each of \(m\) completed data sets, an analyst obtains \(m\) estimates and \(m\) p-values for a hypothesis such as \(H_0:q=q_0\). Taking the arithmetic mean of those p-values is easy, but it does not generally produce a p-value for the multiply imputed analysis. It fails to account correctly for uncertainty between imputations, and it discards information about the estimates’ direction and scale.
A p-value is calculated from a test statistic and its reference distribution under a null hypothesis. It is not an estimate of the effect, nor is it a quantity that Rubin’s rules pool by taking an arithmetic mean. For a scalar estimand, the usual replacement is to pool the estimates and their variances, then form one test using the pooled standard error and an appropriate degrees-of-freedom reference distribution.
The term \(\bar{U}\) represents estimated sampling variance within completed data sets. The term \(B\) represents variation in the estimates across imputations; its multiplier includes a finite-\(m\) adjustment. The pooled test therefore reflects both sources of uncertainty, unlike a calculation that averages test results made separately in each completed data set.
Scalar Rubin Inference
For testing \(H_0:q=q_0\), the scalar Rubin test statistic is \(t=(\bar{q}-q_0)/\sqrt{T}\). In the usual large-complete-data-sample setting, it is compared with a \(t\) distribution whose degrees of freedom are estimated from the relative contribution of between-imputation uncertainty. When \(B>0\), the conventional Rubin degrees of freedom are
For a two-sided test, the pooled p-value is \(2P(t_{\nu}\geq |t|)\), where \(t_{\nu}\) denotes a random variable with a \(t\) distribution and \(\nu\) degrees of freedom. The reference distribution and any finite-complete-data degrees-of-freedom adjustment should match the analysis setting. The formula above is the conventional large-sample complete-data version; it should not be treated as a universal substitute for a finite-sample adjustment when that adjustment is needed.
A useful way to see why total variance has these components is to consider the uncertainty in an estimate when the missing values are not known. Even within a particular completed data set, the estimate has variance \(U_j\). Across plausible completed data sets, the estimate itself also varies. Averaging the within-data-set variance estimates gives \(\bar{U}\); the observed spread of the \(m\) estimates gives \(B\). Because only \(m\) imputations approximate that spread, the between-imputation contribution to total variance is inflated by \(1+1/m\). Thus inference uses \(T\), rather than just a typical \(U_j\), as the variance of the pooled estimate.
This calculation assumes that every completed data set estimates the same scalar quantity on the same scale, that each \(U_j\) estimates the corresponding within-data-set variance, and that the imputations and analysis are adequate for the intended missing-data analysis. The usual \(t\)-based inference also relies on a suitable approximate sampling distribution for the pooled estimate. If these conditions fail—for example, if the estimand changes between imputations—the arithmetic pooling formulas cannot repair the mismatch.
Worked Example: A Mean of P-Values Hides Between-Imputation Uncertainty
Suppose five hypothetical completed data sets give estimates of a regression contrast, in outcome units, of \(0.20,\ 0.25,\ 0.15,\ 0.30,\ 0.10\). Each has estimated variance \(U_j=0.0100\), so each within-imputation standard error is \(0.100\). Test \(H_0:q=0\), using normal-reference p-values for the completed-data tests for simplicity.
The five individual test statistics are \(2.00,\ 2.50,\ 1.50,\ 3.00,\) and \(1.00\). Their two-sided normal-reference p-values are approximately \(0.0455,\ 0.0124,\ 0.1336,\ 0.0027,\) and \(0.3173\). Their arithmetic mean is
That number is not the pooled p-value. For the pooled inference, \(\bar{q}=(0.20+0.25+0.15+0.30+0.10)/5=0.200\), and \(\bar{U}=0.0100\). The squared deviations from \(0.200\) sum to \(0+0.0025+0.0025+0.0100+0.0100=0.0250\), giving \(B=0.0250/(5-1)=0.00625\). Therefore,
The pooled test statistic is \(t=(0.200-0)/0.1323=1.512\). The degrees of freedom are \(\nu=4[1+0.0100/(1.2\times0.00625)]^2=4(2.333)^2=21.78\). The two-sided \(t\)-reference p-value is approximately \(0.145\). Thus scalar Rubin inference gives a p-value of about \(0.145\), not the arithmetic mean \(0.1023\). The difference arises in part because the completed-data estimates vary, adding between-imputation uncertainty to the pooled standard error.
Why the Average Can Give the Wrong Impression
An average of p-values is especially difficult to interpret when completed-data estimates differ in sign. A two-sided p-value depends on the magnitude of a test statistic, not its sign. As a result, estimates pointing in opposite directions can each produce small p-values even when their pooled estimate is near zero and the variation across imputations is substantial.
Worked Example: Opposite Signs Can Produce Small Individual P-Values
Consider four hypothetical estimates of the same contrast, in outcome units: \(0.40,\ -0.40,\ 0.40,\ -0.40\). Let every within-imputation variance be \(0.0400\), so the within-imputation standard error is \(0.200\). Under \(H_0:q=0\), each completed-data test has \(|z|=0.40/0.20=2.00\), giving a two-sided normal-reference p-value of approximately \(0.0455\). The mean of these four p-values is also \(0.0455\).
For pooled inference, \(\bar{q}=0\), \(\bar{U}=0.0400\), and \(B=[4(0.40^2)]/(4-1)=0.6400/3=0.2133\). Hence
The pooled test statistic is \(t=(0-0)/0.5538=0\), so its two-sided p-value is \(1.000\). The pooled estimate is zero, while the individual p-values ignore the opposing signs when converted to two-sided evidence. The unusually large spread of estimates is also a warning to examine the imputation results and analysis assumptions rather than summarize the individual p-values.
The point is not that a pooled p-value must always be larger or smaller than the mean of the individual p-values. It can be either: the two calculations are different quantities, and there is no general rule connecting them. An average can happen to be numerically close to a pooled p-value in a particular analysis, but that agreement does not make the average a calibrated test.
Using a Pooled Test in Practice
For a scalar hypothesis, the relevant inputs are the completed-data estimates and their estimated variances—not just the completed-data p-values. A disciplined calculation keeps the estimate, variance, reference distribution, and hypothesis together:
Confirm that every completed-data analysis estimates the same \(q\), on the same scale, and tests the same null value \(q_0\).
Calculate \(\bar{q}\), \(\bar{U}\), \(B\), and \(T=\bar{U}+(1+1/m)B\).
Calculate \(t=(\bar{q}-q_0)/\sqrt{T}\), determine an appropriate degrees-of-freedom reference, and obtain the p-value from that distribution.
Give the pooled estimate, standard error or interval, degrees of freedom, and p-value. Do not describe the p-value as the probability that \(H_0\) is true.
If the scientific question concerns several parameters jointly, separate scalar Rubin tests do not generally answer the joint question. Use an appropriate pooled multi-parameter procedure, such as the D1 test for a joint regression hypothesis, rather than averaging individual p-values or treating several scalar tests as one joint test. The choice depends on the hypothesis being tested, not on which summary is easiest to calculate.
Worked Example: A Scalar Test of a Pooled Slope
Suppose five hypothetical imputations give a regression slope estimate, in outcome units per one-unit increase in a predictor, of \(1.10,\ 0.90,\ 1.00,\ 1.20,\ 0.80\). Each estimated variance is \(0.1600\). Test \(H_0:q=0\) using the conventional scalar Rubin calculation.
The pooled estimate is \(\bar{q}=(1.10+0.90+1.00+1.20+0.80)/5=1.00\), and \(\bar{U}=0.1600\). The squared deviations are \(0.0100,\ 0.0100,\ 0,\ 0.0400,\ 0.0400\), so \(B=0.1000/4=0.02500\). Thus,
The degrees of freedom are \(\nu=4[1+0.1600/(1.2\times0.02500)]^2=4(6.333)^2=160.4\). The pooled statistic is \(t=(1.00-0)/0.4359=2.294\), with two-sided p-value approximately \(0.023\). At the 5% significance level, this result is inconsistent with \(H_0:q=0\) under the stated model and reference distribution. It does not establish that the null hypothesis has a 2.3% probability of being true.
Diagnostics, Alternatives, and Limitations
A practical diagnostic is to inspect the estimates \(\hat{q}_j\) and their variances \(U_j\), not just a list of p-values. Look for changes in sign, unusually large between-imputation variation, very different within-imputation standard errors, or apparent shifts in the estimand or model across completed data sets. Such patterns may indicate instability, sparse information, or a mismatch in the analysis; averaging p-values conceals them.
The main alternative for a scalar estimand is Rubin inference based on pooled estimates and variances. For a joint hypothesis, use a pooled joint test suited to that hypothesis. Methods that combine p-values are designed for settings with their own assumptions about the tests being combined. Imputation-specific tests are not independent studies: they analyze multiple completed versions of the same incomplete data, so a generic p-value combination rule does not automatically provide valid multiply imputed inference.
Rubin’s scalar procedure also has limits. It requires estimates of the same target and appropriate variance estimates in each completed data set; it does not cure a poorly specified analysis or inadequate imputations. With small samples, unusual sampling distributions, or finite complete-data degrees of freedom, the conventional large-sample degrees-of-freedom formula may be inadequate, and an appropriate adjusted procedure is needed. The pooled p-value is an inferential summary under its assumptions, not a diagnostic for whether the imputations themselves are satisfactory.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Why is the arithmetic mean of completed-data p-values not generally a pooled p-value?
- For \(m=5\), \(\bar{U}=0.0200\), and \(B=0.00500\), calculate \(T\).
- In a two-sided test, four completed-data estimates have opposite signs but equal absolute values and equal within-imputation standard errors. What information does averaging their p-values lose?
- When the scientific hypothesis concerns several regression parameters jointly, should the analyst average their scalar p-values? Name an appropriate type of alternative.
Solutions. (1) The mean does not pool the estimates and variance components or generally have a calibrated reference distribution; it also discards signs and does not correctly account for between-imputation uncertainty. (2) \(T=0.0200+(1+1/5)(0.00500)=0.0200+0.00600=0.02600\). (3) It loses the direction of the estimates and hides whether they cancel in the pooled estimate; it also fails to reflect their across-imputation spread. (4) No. Use a pooled joint test appropriate to the multi-parameter hypothesis, such as the D1 test for a joint regression hypothesis.