Deriving The Rubin Degrees-Of-Freedom Formula

Why Degrees of Freedom Depend on Imputation Uncertainty

A pooled estimate from multiple imputation has uncertainty from two sources: sampling uncertainty within each completed-data analysis and variation among the results from different imputations. Rubin’s degrees-of-freedom formula reflects how much the estimated total variance itself is uncertain because the between-imputation component is estimated from a finite number of imputations.

The earlier material on large-sample degrees of freedom introduced the formula. Here the focus is its justification: approximate the variance of the estimated total variance, then match that variance to the variance of a scaled chi-squared random variable. This derivation concerns a scalar estimate and the large-sample version of the degrees-of-freedom calculation.

Definition: For \(m\) imputations, let \(\hat{Q}_j\) be the scalar estimate from completed data set \(j\), and let \(U_j\) be its estimated variance, for \(j=1,\ldots,m\). Define the pooled estimate \(\bar{Q}=m^{-1}\sum_{j=1}^{m}\hat{Q}_j\), the pooled within-imputation variance \(\bar{U}=m^{-1}\sum_{j=1}^{m}U_j\), and the between-imputation variance \(B=(m-1)^{-1}\sum_{j=1}^{m}(\hat{Q}_j-\bar{Q})^2\). Rubin’s estimated total variance is \(T=\bar{U}+(1+1/m)B\). The relative increase in variance used here is \(r=(1+1/m)B/\bar{U}\).

The factor \(1+1/m\) accounts for the finite number of imputations in the total-variance estimate. In the degrees-of-freedom derivation, \(\bar{U}\) is treated as fixed, while \(B\) is the estimated component whose sampling variation drives the approximation. This is a simplifying large-sample treatment, not a claim that \(\bar{U}\) is known without error in every analysis.

Assumptions Behind the Moment Matching

The argument uses an approximate model for the imputation-specific estimates. Conditional on the observed data and the quantities being estimated, suppose the \(\hat{Q}_j\) behave approximately like independent draws from a common normal distribution with mean \(Q\) and variance \(\sigma_B^2\). Then their sample variance \(B\) estimates \(\sigma_B^2\). Under this normal model, \((m-1)B/\sigma_B^2\) has a chi-squared distribution with \(m-1\) degrees of freedom.

That distributional result implies \(\operatorname{Var}(B)=2\sigma_B^4/(m-1)\). Since \(\sigma_B^2\) is unknown, the derivation estimates this variance by substituting \(B\) for \(\sigma_B^2\), giving \(\operatorname{Var}(B)\approx 2B^2/(m-1)\). Treating \(\bar{U}\) as fixed, the estimated total variance changes with \(B\) by the multiplier \(1+1/m\). Therefore:

\operatorname{Var}(T)\approx (1+1/m)^2\operatorname{Var}(B) \approx \frac{2(1+1/m)^2B^2}{m-1}.

The assumptions matter. The imputation-specific estimates should be sufficiently well behaved for the normal approximation to be reasonable; the imputations should act like independent draws for purposes of estimating between-imputation variation; and the within-imputation variance should be a meaningful estimate for the scalar quantity. In a MICE analysis, serious chain nonconvergence, strong unexplained skewness in the estimates, or a highly influential imputation can undermine these approximations. The algebra below does not repair such problems.

Deriving the Approximate Degrees of Freedom

A scaled chi-squared random variable can be written as \(cX\), where \(X\sim\chi^2_{\nu}\), \(c\) is a positive scale, and \(\nu\) is its degrees of freedom. Its mean is \(c\nu\), and its variance is \(2c^2\nu\). If it approximates the estimated total variance, its mean is matched to \(T\). Thus \(c\nu=T\), and its variance can be written \(2T^2/\nu\). Matching that variance to the approximation for \(\operatorname{Var}(T)\) gives the degrees of freedom.

Approximate the variance of the total-variance estimate.
Use the chi-squared sampling variance of \(B\), treating \(\bar{U}\) as fixed: \(\operatorname{Var}(T)\approx 2(1+1/m)^2B^2/(m-1)\).

Match the scaled chi-squared variance.
For a scaled chi-squared variable with mean \(T\) and degrees of freedom \(\nu\), the variance is \(2T^2/\nu\).

Solve for the degrees of freedom.
Equate the two variance expressions and substitute \(T=\bar{U}+(1+1/m)B\).

\frac{2T^2}{\nu} \approx \frac{2(1+1/m)^2B^2}{m-1}

\nu \approx (m-1)\left\{\frac{T}{(1+1/m)B}\right\}^2 = (m-1)\left\{1+\frac{\bar{U}}{(1+1/m)B}\right\}^2 = (m-1)\left(1+\frac{1}{r}\right)^2.

Result: Under the large-sample and moment-matching approximations above, Rubin’s degrees of freedom for a pooled scalar estimate are \(\nu_R=(m-1)(1+1/r)^2\), where \(r=(1+1/m)B/\bar{U}\). This is an approximate, generally non-integer degrees of freedom, not a count of observations.

The form of the result is informative. When between-imputation variation is large relative to within-imputation variance, \(r\) is large and \(\nu_R\) is closer to \(m-1\). When \(r\) is small, the factor \((1+1/r)^2\) can make \(\nu_R\) much larger. In that case, the estimated total variance is less dominated by the finite-imputation estimate \(B\), so this particular source of uncertainty is less influential.

Worked Calculations

Worked Example: Moderate Between-Imputation Variation

Hypothetical data. Suppose a scalar analysis uses \(m=5\) imputations and gives \(\bar{U}=4.00\) squared units and \(B=1.00\) squared units. Assume the imputation-specific estimates are approximately normal and independent for this calculation, the chains have converged, and no single imputation is anomalous. These summaries are consistent with using the stated approximation, though summaries alone cannot establish all its assumptions.

First calculate \(r\):

r=(1+1/5)\frac{1.00}{4.00} =1.2(0.25)=0.300.

Then Rubin’s approximate degrees of freedom are:

\nu_R=(5-1)\left(1+\frac{1}{0.300}\right)^2 =4(4.3333)^2 \approx 75.11.

As a check on the moment matching, \(T=4.00+1.2(1.00)=5.20\), and the estimated variance of \(T\) is \(2(1.2)^2(1.00)^2/4=0.720\). The matching degrees of freedom are \(2(5.20)^2/0.720\approx75.11\), agreeing with the formula. The result is an approximate degrees of freedom for subsequent scalar inference, not 75.11 imputations or observations.

Worked Example: Larger Relative Increase in Variance

Hypothetical data. Suppose \(m=10\), \(\bar{U}=9.00\) squared units, and \(B=3.00\) squared units. Assume approximate normality of the imputation-specific estimates, adequate chain mixing and convergence, and no strong outlier among the ten estimates.

The finite-imputation multiplier is \(1+1/10=1.1\), so:

r=1.1\frac{3.00}{9.00}=0.3667.

Substitution gives:

\nu_R=(10-1)\left(1+\frac{1}{0.3667}\right)^2 \approx 9(3.727)^2 \approx 125.0.

Here \(T=9.00+1.1(3.00)=12.30\). The variance approximation is \(2(1.1)^2(3.00)^2/9=2.420\). Moment matching gives \(2(12.30)^2/2.420\approx125.0\). The degrees of freedom exceed \(m-1=9\) because the pooled within-imputation component is substantial relative to the adjusted between-imputation component; the formula is not limited to \(m-1\).

Worked Example: Few Imputations and Strong Between-Imputation Variation

Hypothetical data. Suppose \(m=3\), \(\bar{U}=1.00\) squared units, and \(B=4.00\) squared units. Assume the three estimates are plausible independent draws from an approximately normal distribution, with no convergence or outlier warning. With only three imputations, however, the estimate of between-imputation variance is especially noisy.

The adjusted between-imputation component and relative increase are:

(1+1/3)B=\frac{4}{3}(4.00)=5.333, \qquad r=\frac{5.333}{1.00}=5.333.

Thus:

\nu_R=(3-1)\left(1+\frac{1}{5.333}\right)^2 =2(1.1875)^2 \approx2.820.

For verification, \(T=1.00+5.333=6.333\), and the estimated variance of \(T\) is \((4/3)^2[2(4.00)^2/2]\approx28.444\). Therefore \(2(6.333)^2/28.444\approx2.820\). This small, non-integer result reflects substantial uncertainty in the estimated total variance when the number of imputations is small and \(B\) is large relative to \(\bar{U}\).

Diagnostics, Alternatives, and Limitations

This derivation isolates the uncertainty arising from estimating \(B\). It does not incorporate every possible source of uncertainty in a finite-sample analysis. In particular, the approximation treats \(\bar{U}\) as fixed and uses a large-sample treatment of the within-imputation analysis. The course’s later material on complete-data degrees of freedom and small-sample adjustments addresses additional finite-sample considerations; those adjustments are not part of the formula derived here.

A common alternative is to use a normal reference distribution, effectively treating the degrees of freedom as infinite. That may be a reasonable approximation when the resulting degrees of freedom are very large, but it can understate tail uncertainty when \(\nu_R\) is small. Another tempting shortcut is to set the degrees of freedom equal to \(m-1\), because \(B\) is a sample variance from \(m\) estimates. That ignores the addition of \(\bar{U}\) to the total variance: moment matching is for \(T\), not for \(B\) alone, and the within-imputation component can substantially raise the matched degrees of freedom.

Before relying on the approximation, inspect the imputation process and the imputation-specific estimates. Poor chain mixing, nonconvergence, a highly skewed distribution of estimates, or one extreme estimate can make the normal and chi-squared approximations questionable. Also verify the definitions used by the software: omitting the \(1+1/m\) factor changes \(r\) and therefore changes the result. If \(B\) is exactly zero, the displayed expression involves division by zero; that boundary case requires care rather than mechanical substitution.

Common Mistake

Common mistake: “There are five imputations, so the degrees of freedom must be \(m-1=4\).” The value \(m-1\) is the degrees of freedom in the chi-squared result for the sample variance \(B\). Rubin’s formula matches the uncertainty of the total variance, which also contains \(\bar{U}\). The correct calculation is \((m-1)(1+1/r)^2\), so its value depends on the relative size of the two variance components and can be much greater than \(m-1\).

Check Your Understanding

Use the moment-matching argument and show the relevant calculations.

In one sentence, identify which component’s sampling variation drives this large-sample derivation.
With \(m=5\), \(\bar{U}=4.00\), and \(B=1.00\), calculate \(r\) and \(\nu_R\).
With \(m=10\), \(\bar{U}=9.00\), and \(B=3.00\), calculate \(T\) and the moment-matched degrees of freedom using \(2T^2/\operatorname{Var}(T)\).
Explain why setting \(\nu_R=m-1\) generally does not reproduce the formula.
Name one diagnostic finding that would make the normal or chi-squared approximation less credible.

Solutions. (1) The derivation is driven by the sampling variation of \(B\), with \(\bar{U}\) treated as fixed. (2) \(r=1.2(1.00)/4.00=0.300\), so \(\nu_R=4(1+1/0.300)^2\approx75.11\). (3) \(T=9.00+1.1(3.00)=12.30\); \(\operatorname{Var}(T)\approx2(1.1)^2(3.00)^2/9=2.420\), so \(2T^2/\operatorname{Var}(T)=2(12.30)^2/2.420\approx125.0\). (4) \(m-1\) describes the chi-squared variation of \(B\) alone, whereas the formula matches the uncertainty of \(T=\bar{U}+(1+1/m)B\); the \(\bar{U}\) component changes the matched degrees of freedom. (5) Examples include nonconverged chains, poor chain mixing, a strongly skewed set of imputation-specific estimates, or one extreme estimate.