library(kstools)

Introduction

Missing data are common in clinical trials, observational studies, registries, and real-world evidence analyses. A standard approach for handling missing values is multiple imputation (MI), where several complete datasets are generated, each containing plausible replacements for missing observations.

After analysing each imputed dataset, the analyst obtains:

• a parameter estimate from each imputation;
• a corresponding standard error.

The final step is to combine these results into a single inference. This process is commonly known as pooling or MI analysis and is described by Rubin (1987).

The mianalyze() function in kstools implements Rubin’s pooling rules and supports the Barnard-Rubin (1999) small-sample degrees-of-freedom adjustment, equivalent to the EDF= option in SAS PROC MIANALYZE. The exported alias ks_mianalyze() provides the same interface.

Multiple Imputation Workflow

The typical workflow is:

1. Create m imputed datasets.
2. Perform the same analysis on each dataset.
3. Extract parameter estimates and standard errors.
4. Pool the results using Rubin’s rules.

Conceptually:

%%{init: {"theme": "dark", "themeVariables": {"lineColor": "#94a3b8", "primaryTextColor": "#f8fafc"}}}%%
flowchart TD
    A([Original Dataset])
    B([Multiple Imputation])
    D1[Imputed Dataset 1<br/>Q₁, U₁]
    D2[Imputed Dataset 2<br/>Q₂, U₂]
    D3[Imputed Dataset 3<br/>Q₃, U₃]
    Ddots[...]
    Dm[Imputed Dataset m<br/>Qₘ, Uₘ]
    R{{Rubin's Rules<br/>Pool Estimates and Variances}}
    F([Pooled Estimate<br/>SE, CI, p-value])
    
    A --> B
    B --> D1
    B --> D2
    B --> D3
    B --> Ddots
    B --> Dm
    D1 --> R
    D2 --> R
    D3 --> R
    Ddots --> R
    Dm --> R
    R --> F
    
    classDef source fill:#1e40af,stroke:#93c5fd,stroke-width:2px,color:#f8fafc;
    classDef process fill:#c2410c,stroke:#fdba74,stroke-width:2px,color:#f8fafc;
    classDef imputed fill:#166534,stroke:#86efac,stroke-width:2px,color:#f8fafc;
    classDef pooling fill:#6d28d9,stroke:#c4b5fd,stroke-width:2px,color:#f8fafc;
    classDef result fill:#15803d,stroke:#86efac,stroke-width:3px,color:#f8fafc;

    class A source;
    class B process;
    class D1,D2,D3,Dm imputed;
    class R pooling;
    class F result;

    linkStyle default stroke:#94a3b8,stroke-width:1.5px

Rubin’s Rules

Suppose Q_i is the estimate from imputation i and U_i is its variance. Let m be the number of imputations.

Pooled Estimate

$$\bar{Q} = \frac{1}{m} \sum_{i=1}^{m}Q_i$$

Within-Imputation Variance

$$\bar{U} = \frac{1}{m} \sum_{i=1}^{m}U_i$$

This quantity reflects uncertainty that would exist even if no data were missing.

Between-Imputation Variance

$$B = \frac{1}{m-1}\sum_{i=1}^{m}(Q_i-\bar{Q})^2$$

Large values indicate substantial uncertainty due to missing data.

Total Variance

$$T = \bar{U} + \left(1+\frac{1}{m}\right) B$$

The pooled standard error is $SE = \sqrt{T}$.

Relative Increase in Variance

$$r = \frac{\left(1+\frac{1}{m}\right) B}{\bar{U}}$$

Values near zero indicate little impact from missing data.

Fraction of Missing Information

For finite Rubin degrees of freedom v_m,

$$\hat\lambda = \frac{r + 2/(v_m + 3)}{r + 1}$$

When v_m is treated as infinite (large-sample case), this simplifies to λ̂ = r/(r + 1). The quantity is often interpreted as the proportion of uncertainty attributable to missingness.

Relative Efficiency

$$RE = \frac{1}{1 + \lambda/m}$$

Values close to 1 indicate that additional imputations would provide little benefit.

Degrees of Freedom

Rubin (1987)

$$v_m = (m-1)\left(1+\frac{\bar{U}}{\left(1+m^{-1}\right)B}\right)^2$$

When m is small and missing information is limited, v_m can become very large.

Problem with Large-Sample Degrees of Freedom

Suppose a model parameter is based on only 40 residual degrees of freedom. Rubin’s formula may yield v_m = 500, which exceeds the information available in the complete-data analysis and can produce confidence intervals that are unrealistically narrow.

Barnard-Rubin Adjustment

Barnard and Rubin (1999) proposed a small-sample correction. Let v₀ be the complete-data degrees of freedom. Define

$$\gamma = \frac{(1 + m^{-1}) B}{T}$$

and

$$\hat{v}_{obs} = (1-\gamma)\, v_0\, \frac{v_0+1}{v_0+3}$$

The adjusted degrees of freedom are

$$v_m^* = \left(\frac{1}{v_m} + \frac{1}{\hat{v}_{obs}}\right)^{-1}$$

This is the method used by SAS PROC MIANALYZE when EDF= is specified, and by mianalyze(edf = ...).

Basic Example

est <- c(0.345, 0.303, 0.298, 0.378)
err <- c(0.056, 0.101, 0.099, 0.034)
out <- mianalyze(est, err)
out

Example: Manual Verification

est <- c(0.345, 0.303, 0.298, 0.378)
err <- c(0.056, 0.101, 0.099, 0.034)

m <- length(est)
Qbar <- mean(est)
Vw <- mean(err^2)
Vb <- var(est)
Vt <- Vw + (1 + 1 / m) * Vb

manual <- c(Estimate = Qbar, Vw = Vw, Vb = Vb, Vt = Vt, SE = sqrt(Vt))
pooled <- c(
  Estimate = out$Estimate, Vw = out$Vw, Vb = out$Vb,
  Vt = out$Vt, SE = out$StdErr
)
rbind(manual = manual, mianalyze = pooled)
#>           Estimate        Vw       Vb       Vt         SE
#> manual       0.331 0.0060735 0.001426 0.007856 0.08863408
#> mianalyze    0.331 0.0060735 0.001426 0.007856 0.08863408

Example: Treatment Difference

est <- c(-2.8, -3.2, -2.5, -3.1, -2.9)
err <- c(1.10, 1.15, 1.05, 1.18, 1.12)
mianalyze(est, err)

Interpretation:

• the pooled estimate is the average treatment difference;
• the pooled standard error combines model and imputation uncertainty;
• confidence intervals incorporate missing-data uncertainty.

Example: Odds Ratio Analysis

Pool results on the log scale, then exponentiate.

log_or <- c(0.18, 0.21, 0.15, 0.19, 0.23)
se_log_or <- c(0.09, 0.10, 0.09, 0.08, 0.11)
res <- mianalyze(log_or, se_log_or)
res

exp(c(OR = res$Estimate, LCL = res$lcl, UCL = res$ucl))
#>        OR       LCL       UCL 
#> 1.2116705 0.9948654 1.4757226

Example: Hazard Ratio Analysis

Cox models typically return log-hazard ratios.

log_hr <- c(-0.22, -0.17, -0.31, -0.28, -0.24)
se <- c(0.14, 0.13, 0.15, 0.12, 0.13)
res <- mianalyze(log_hr, se)
res

exp(c(HR = res$Estimate, LCL = res$lcl, UCL = res$ucl))
#>        HR       LCL       UCL 
#> 0.7834876 0.5861225 1.0473115

Example: Small-Sample Degrees of Freedom

Assume the complete-data analysis would have v₀ = 40:

est <- c(0.345, 0.303, 0.298, 0.378)
err <- c(0.056, 0.101, 0.099, 0.034)

br <- mianalyze(est, err, edf = 40)
rubin <- mianalyze(est, err)

data.frame(
  method = c("Barnard-Rubin (edf=40)", "Rubin large-sample"),
  df = c(br$df, rubin$df),
  lcl = c(br$lcl, rubin$lcl),
  ucl = c(br$ucl, rubin$ucl)
)
#>                   method       df       lcl       ucl
#> 1 Barnard-Rubin (edf=40) 19.57893 0.1458573 0.5161427
#> 2     Rubin large-sample 58.27270 0.1535973 0.5084027

The adjusted degrees of freedom are smaller, producing slightly wider confidence intervals.

Comparison with SAS PROC MIANALYZE

mianalyze() is a univariate analogue of SAS PROC MIANALYZE. The formulas for Q̄, Ū, B, T, r, λ̂, RE, v_m, and the Barnard-Rubin v_m^* match those in the SAS/STAT User’s Guide (Chapter The MIANALYZE Procedure).

Option mapping

Note: SAS documents a default of EDF=1 with no adjustment when EDF= is omitted. In kstools, omitting edf (or using edf = Inf) likewise uses Rubin’s large-sample v_m without Barnard-Rubin adjustment.

Scope

Aligned with SAS univariate MI inference:

• pooled estimate, total variance, t-based CI and test;
• variance diagnostics r, λ̂, RE;
• Barnard-Rubin df via edf.

Not implemented (multivariate SAS features):

• Wald multivariate tests (MULT / MULTIVARIATE);
• TEST linear hypotheses;
• pooling from full covariance matrices alone (without per-parameter SEs).

SAS Example 58.1 (means and standard errors, EDF=30)

The Fitness example in SAS/STAT 12.3 Example 58.1 combines univariate means and standard errors across m = 5 imputations with EDF=30 (sample size n = 31). Reproducing Oxygen from the published input table:

oxygen <- mianalyze(
  est = c(47.0120, 47.2407, 47.4995, 47.1485, 47.0042),
  err = c(0.95984, 0.93540, 1.00766, 0.95439, 0.96528),
  edf = 30
)
oxygen

Comparison with SAS Output 58.1.2 / 58.1.3 (Oxygen). Small differences are expected because the User’s Guide prints input means/SEs with limited precision.

sas_oxygen <- data.frame(
  quantity = c(
    "Estimate", "StdErr", "Between", "Within", "Total",
    "DF", "r", "lambda", "RE", "LCL", "UCL"
  ),
  sas = c(
    47.180993, 0.990266, 0.041478, 0.930853, 0.980626,
    26.298, 0.053471, 0.051977, 0.989712, 45.1466, 49.2154
  ),
  mianalyze = c(
    oxygen$Estimate, oxygen$StdErr, oxygen$Vb, oxygen$Vw, oxygen$Vt,
    oxygen$df, oxygen$r, oxygen$lambda, oxygen$re, oxygen$lcl, oxygen$ucl
  )
)
sas_oxygen$abs_diff <- abs(sas_oxygen$mianalyze - sas_oxygen$sas)
knitr::kable(sas_oxygen, digits = 6)

The same pattern holds for RunTime and RunPulse in that example, and for the Getting Started / Example 58.3 regression coefficients when moments are reconstructed from the published Between/Within variances (with and without EDF=28). Automated checks live in tests/testthat/test-ks_mianalyze-sas.R.

Relative efficiency (SAS Table 58.2)

SAS tabulates RE = 1/(1 + λ/m). For illustration:

re_grid <- expand.grid(
  m = c(3, 5, 10, 20),
  lambda = c(0.1, 0.2, 0.3, 0.5, 0.7)
)
re_grid$RE <- with(re_grid, 1 / (1 + lambda / m))
knitr::kable(
  xtabs(RE ~ m + lambda, data = re_grid),
  digits = 4,
  caption = "Relative efficiency (matches SAS Table 58.2)"
)

API correspondence

SAS:

proc mianalyze data=outuni edf=30;
  modeleffects Oxygen;
  stderr SOxygen;
run;

Equivalent R:

mianalyze(est = oxygen_est, err = oxygen_se, edf = 30)

Typical Clinical Trial Workflow

A common workflow with the mice package:

library(mice)

imp <- mice(adsl, m = 20)

fits <- with(
  imp,
  lm(CHG ~ TRT01A + BASE)
)

est <- sapply(
  fits$analyses,
  function(x) coef(summary(x))["TRT01AActive", "Estimate"]
)

err <- sapply(
  fits$analyses,
  function(x) coef(summary(x))["TRT01AActive", "Std. Error"]
)

mianalyze(est, err)

Interpreting Output

Variance Components

Large V_b relative to V_w indicates substantial missing-data uncertainty.

Missing Information

Inference

Special Cases

Single Imputation

mianalyze(est = 10, err = 2)

Here B = 0, T = Ū, and results reduce to ordinary complete-data inference.

No Between-Imputation Variability

mianalyze(est = rep(5, 5), err = rep(1, 5))

When all imputations produce the same estimate, B = 0 and no additional uncertainty is attributed to missing data.

Summary

• mianalyze() / ks_mianalyze() pool MI point estimates and standard errors with Rubin’s rules: Q̄, T = Ū + (1 + 1/m)B, t-based intervals and tests.
• Diagnostics r, λ̂, and RE quantify the impact of missingness and the efficiency of finite m.
• edf applies the Barnard-Rubin (1999) small-sample df adjustment, matching SAS EDF=.
• Agreement with SAS univariate PROC MIANALYZE was checked against SAS/STAT 12.3 formulas, Table 58.2, Example 58.1 (EDF=30), and Example 58.3 (EDF=28).
• Pool on the analysis scale (e.g. log-OR, log-HR), then transform intervals if needed.
• Prefer edf whenever complete-data residual df are modest so that v_m does not exceed the information in the original sample.

References

Barnard, J. and Rubin, D.B. (1999). Small-sample degrees of freedom with multiple imputation. Biometrika, 86, 948–955.

Rubin, D.B. (1987). Multiple Imputation for Nonresponse in Surveys. Wiley.

SAS Institute Inc. (2013). SAS/STAT 12.3 User’s Guide. Chapter 58, The MIANALYZE Procedure. Cary, NC: SAS Institute Inc.

Van Buuren, S. (2018). Flexible Imputation of Missing Data. Second Edition. Chapman & Hall/CRC.