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. pigauto workflow classes and markers are caller-controlled metadata; they do not authenticate the fits or their history.

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. For recognized drmTMB and gllvmTMB_multi fits, pigauto still requires convergence code zero and a positive-definite Hessian, even when a custom extractor is supplied.

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. 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. drmTMB and gllvmTMB need to be installed and loadable only when their automatic adapters are used. For either recognized backend class, pooling also requires the fit to report scalar convergence code 0 and pdHess = TRUE; missing or invalid status is rejected before coefficient extraction, including with custom extractors.

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.

See also

multi_impute_analysis(), with_imputations()

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 conf.high
#>  (Intercept)  -0.1680    0.1827   161.6   -0.9192  0.3594  -0.5288    0.1929
#>            x  -0.4498    0.1801 628604.    -2.498 0.01250  -0.8027  -0.09685
#>            z  0.07221    0.2012   26.80    0.3589  0.7225  -0.3407    0.4852
#>       fmi      riv
#>   0.08985  0.08537
#>  0.001264 0.001263
#>    0.2473   0.2394
# }