Pool downstream model fits across multiple imputations (Rubin's rules)

Combine regression coefficients from M model fits – one per imputed dataset – into a single pooled table using Rubin's rules. The pooled standard errors combine within-imputation sampling variance and between-imputation variance. Correct Rubin arithmetic does not make an incompatible imputation model inferentially valid. For pigauto inference, fits should come from the documented multi_impute_analysis() workflow, or, for continuous traits, from multi_impute() with draws_method = "posterior" (see with_imputations() for the analyses those draws cover), via with_imputations(). The backend is experimental and supports fixed-effect coefficients only. Variance components, correlations, random-effect predictions, BLUPs/conditional modes, latent loadings, and other structured parameters are unsupported.

Usage

pool_mi(
  fits,
  conf.level = 0.95,
  coef_fun = NULL,
  vcov_fun = NULL,
  df_fun = NULL,
  tidy_fun = NULL
)

Arguments

fits: A list of model fits of length M >= 2. Automatic fixed-effect adapters are provided for stats::lm, stats::glm, nlme::gls, nlme::lme, lme4::merMod, glmmTMB::glmmTMB, drmTMB, and gllvmTMB_multi fits. Other classes implementing coef() and vcov() also work. The output of with_imputations() is accepted directly. Bayesian fits (brms, MCMCglmm) are rejected; see Details. Bare user-supplied fit lists are accepted with an explicit warning because their imputation provenance cannot be verified.
conf.level: Confidence level for the pooled interval (default 0.95).
coef_fun: Optional function extracting a named numeric fixed-effect vector from one fit. NULL uses the automatic class adapter. Custom coefficient and covariance extractors can be supplied independently; callers must ensure custom extractors select fixed effects only.
vcov_fun: Optional function extracting the fixed-effect covariance matrix. NULL uses the automatic class adapter. Base matrices and Matrix objects are accepted; the selected covariance block must be square, finite, symmetric, and have non-negative diagonal entries.
df_fun: Optional function returning the complete-data residual degrees of freedom nu_com from one fit. When supplied, pooled degrees of freedom use the Barnard & Rubin (1999) small-sample correction, which is less biased for short series. When NULL (the default) the classical Rubin (1987) formula is used.
tidy_fun: Optional function returning a data.frame with unique term, numeric estimate, and non-negative finite std.error columns. This is an alternative to coef_fun and vcov_fun, not a supplement; combining them is an error. Callers must ensure it selects fixed effects only.

Value

A data.frame with one row per coefficient and columns:

term: Coefficient name.
estimate: Pooled point estimate (mean across fits).
std.error: Pooled standard error sqrt(T) where T = W + (1 + 1/M) * B.
df: Pooled degrees of freedom (Barnard-Rubin if df_fun supplied, else classical Rubin).
statistic: estimate / std.error.
p.value: Two-sided p-value from a t distribution on df.
conf.low, conf.high: Pooled conf.level interval.
fmi: Fraction of missing information.
riv: Relative increase in variance due to non-response.

Details

Let $\hat\theta_i$ be the coefficient vector from fit i and $U_i = \mathrm{vcov}(\mathrm{fit}_i)$, for $i = 1, \ldots, M$. Rubin's rules (Rubin 1987) give $$\bar\theta = M^{-1} \sum_i \hat\theta_i$$ $$W = M^{-1} \sum_i \mathrm{diag}(U_i)$$ $$B = (M-1)^{-1} \sum_i (\hat\theta_i - \bar\theta)^2$$ $$T = W + (1 + 1/M) B$$ with pooled standard error $\sqrt{T}$. The relative increase in variance is $r = (1 + 1/M) B / W$, the classical pooled df is $\nu_{\text{old}} = (M-1)(1 + 1/r)^2$, and the fraction of missing information is $$\mathrm{fmi} = (r + 2/(\nu + 3)) / (r + 1).$$ When df_fun returns finite complete-data df nu_com, the Barnard-Rubin (1999) correction combines $\nu_{\text{obs}} = ((\nu_{\text{com}}+1)/(\nu_{\text{com}}+3)) \nu_{\text{com}} (1 - \lambda)$ with nu_old via $\nu_{\text{BR}} = 1/(1/\nu_{\text{old}} + 1/\nu_{\text{obs}})$. With no between-imputation variation (B = 0), the classical limit is df = Inf, riv = 0, and fmi = 0. If finite complete-data df are supplied, Barnard–Rubin instead retains finite observed-data df and its small-sample FMI adjustment. A completely deterministic quantity (B = W = 0) always has zero FMI and infinite df.

Only fixed-effect coefficients are pooled in version 0.10.0. Random-effect variances and correlations, BLUPs/conditional modes, latent loadings, and other structured parameters require parameter-specific transformations and are not supported by the automatic pool_mi() adapters. Custom extractors are an expert escape hatch and cannot be inspected by pigauto; using them to select unsupported structured parameters is outside the documented scope. The glmmTMB adapter selects conditional fixed effects only. The drmTMB adapter includes named distributional fixed-effect blocks such as regression coefficients for mu and sigma; those are fixed coefficients, not random-effect variance components.

Supported classes. lm, glm, gls, lme, merMod (lme4), glmmTMB (conditional component), drmTMB, and gllvmTMB_multi are handled automatically. Other classes are pooled through their coef() / vcov() methods, or through coef_fun / vcov_fun / tidy_fun.

Bayesian fits (brms, MCMCglmm) are rejected because Rubin's rules are not the right tool for posterior samples: they are pooled by concatenating posterior draws across imputations. For brms, use brms::brm_multiple(formula, data = mi$datasets), or brms::combine_models() on a list of existing fits. For MCMCglmm, concatenate the Sol samples across fits. No Bayesian downstream workflow is supported by the analysis-aware backend.

References

Rubin DB (1987). Multiple Imputation for Nonresponse in Surveys. Wiley.

Barnard J, Rubin DB (1999). "Small-sample degrees of freedom with multiple imputation." Biometrika 86(4): 948-955.

Nakagawa S, Freckleton RP (2008). "Missing inaction: the dangers of ignoring missing data." Trends in Ecology & Evolution 23(11): 592-596.

Nakagawa S, Freckleton RP (2011). "Model averaging, missing data and multiple imputation: a case study for behavioural ecology." Behavioral Ecology and Sociobiology 65(1): 103-116.

Examples

# \donttest{
dat <- data.frame(y = stats::rnorm(30L), x = stats::rnorm(30L),
                  z = stats::rnorm(30L))
dat$x[seq(3L, 30L, by = 5L)] <- NA_real_
mi <- multi_impute_analysis(
  data = dat, formula = y ~ x + z, missing = "x",
  model = "lm", m = 2L
)
fits <- with_imputations(mi, function(d) lm(y ~ x + z, data = d))
pool_mi(fits)
#> Pooled estimates from 2 multiply-imputed fits (Rubin's rules)
#> Confidence level: 95%
#> 
#>         term estimate std.error          df statistic p.value conf.low
#>  (Intercept)  0.04342    0.1846 4212192328.    0.2352  0.8141  -0.3184
#>            x  -0.1296    0.1959       363.1   -0.6616  0.5086  -0.5150
#>            z -0.07846    0.1675      96755.   -0.4684  0.6395  -0.4067
#>  conf.high        fmi        riv
#>     0.4053 0.00001541 0.00001541
#>     0.2557    0.05765    0.05538
#>     0.2498   0.003235   0.003225
# }