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. To pool an analysis, first use with_imputations() to fit one eligible downstream model to each completed dataset, then pass those fits here. Start with the analysis-aware multi_impute_analysis() workflow or, for eligible continuous-trait data, multi_impute() with draws_method = "posterior". Check those help topics for input restrictions and inferential limits. This pooling function 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
)

# S3 method for class 'pigauto_pooled'
print(x, digits = 4, ...)

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. Failed try-error and pigauto_mi_error elements are removed with a warning; at least two retained fits are required. Inspect these failures before interpreting results: dropping analyses can change the imputation sample.
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. Called once on the first retained fit after failures are removed; it must return one finite positive number assumed appropriate for the common analysis. Supplying it selects the Barnard & Rubin (1999) correction for finite complete-data degrees of freedom. 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.
x: A pooled coefficient table.
digits: Number of significant digits used for the displayed table.
...: Ignored. Included for compatibility with the base print() generic.

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.

The m attribute records the retained fit count used in all pooling calculations; it can be smaller than the supplied list.

The print method returns x invisibly.

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$. The implementation uses the following Rubin-style pooling equations: $$\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).$$ Here nu is the pooled degrees of freedom returned in df: the classical value when df_fun = NULL, or the finite-complete-data adjustment below when df_fun is supplied. When df_fun returns finite complete-data df nu_com, pigauto applies a finite-df adjustment motivated by Barnard and Rubin (1999). The equations below specify pigauto's implementation. It uses $\lambda = (1 + 1/M) B / T$, the proportion of total variance from between-imputation uncertainty, which differs from the finite-df fmi above. It 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. For a coefficient with W = 0 and B > 0, all total variance is between imputations: the pooled coefficient estimates differ across imputed datasets, while each fit reports zero variance. This implementation returns df = M - 1 and fmi = 1, including when df_fun returns a finite positive complete-data df. The observed-data df term is zero in this limit, so the harmonic adjustment is skipped. When both variance components are zero (W = B = 0), it returns df = Inf and fmi = 0.

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. Posterior-draw combination is a separate procedure, not a blanket recommendation: use it only with an eligible imputation method and a compatible downstream model, check convergence for every imputation-specific fit, and give each imputed dataset equal representation in the combined draws. For brms, use brms::brm_multiple(formula, data = mi$datasets), or brms::combine_models() on a list of existing fits. For MCMCglmm, concatenate equally represented Sol samples only after checking each fit. Combined draws do not establish calibrated inference for pigauto. See the multiple-imputation guide for workflow limits and the current brms diagnostics guidance. No Bayesian downstream workflow is supported by the analysis-aware backend.

The references below provide general multiple-imputation and ecological context. The exact book passages for these equations and the ecological review passages have not been checked against the primary full texts; these citations do not validate pigauto's imputation workflows.

Printing

Prints coefficients, retained-fit count, confidence level and an FMI-based imputation-count suggestion. That suggestion is a stability heuristic, not a validity check. The ... argument is currently unused.

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_
# Workflow example only: m = 2L is not an analysis recommendation.
# Choose enough imputations to make Monte Carlo error small relative to
# the standard errors you report.
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
# }