amices / mice

Multivariate Imputation by Chained Equations
https://amices.org/mice/
GNU General Public License v2.0
433 stars 107 forks source link

add flexible pooling for the total variance of (Q-Qbar) #509

Closed gerkovink closed 1 year ago

gerkovink commented 1 year ago

I had need for a custom total variance calculation, other than $T = \bar{U} + B + B/m$. Have updated pool() with argument custom.t to accept a custom string for $T$ that takes the following usage:

library(mice)
library(magrittr)
set.seed(123)
imp <- mice(nhanes, printFlag = FALSE)
fit <- with(imp, lm(age ~ bmi))

# regular calculation
pool(fit)  %T>% 
  print %>% summary()
#> Class: mipo    m = 5 
#>          term m    estimate        ubar            b           t dfcom       df
#> 1 (Intercept) 5  3.91587204 1.154102714 0.3341192913 1.555045863    23 12.48683
#> 2         bmi 5 -0.08122935 0.001629126 0.0003996515 0.002108708    23 13.53213
#>         riv    lambda       fmi
#> 1 0.3474068 0.2578336 0.3536785
#> 2 0.2943799 0.2274293 0.3208922
#>          term    estimate  std.error statistic       df    p.value
#> 1 (Intercept)  3.91587204 1.24701478  3.140197 12.48683 0.00816489
#> 2         bmi -0.08122935 0.04592067 -1.768906 13.53213 0.09943129

# custom calculation
pool(fit, rule = "reiter2003", custom.t = ".data$b + .data$b/m") %T>% 
  print %>% summary()
#> Class: mipo    m = 5 
#>          term m    estimate        ubar            b            t dfcom
#> 1 (Intercept) 5  3.91587204 1.154102714 0.3341192913 0.4009431496    23
#> 2         bmi 5 -0.08122935 0.001629126 0.0003996515 0.0004795819    23
#>         df       riv lambda fmi
#> 1 1335.291 0.3474068     NA  NA
#> 2 1828.730 0.2943799     NA  NA
#>          term    estimate  std.error statistic       df      p.value
#> 1 (Intercept)  3.91587204 0.63320072  6.184251 1335.291 8.278363e-10
#> 2         bmi -0.08122935 0.02189936 -3.709212 1828.730 2.141111e-04

# another custom calculation
outsideoffunction <- 14
pool(fit, rule = "rubin1987", custom.t = "outsideoffunction")  %T>% 
  print %>% summary()
#> Class: mipo    m = 5 
#>          term m    estimate        ubar            b  t dfcom       df
#> 1 (Intercept) 5  3.91587204 1.154102714 0.3341192913 14    23 20.53591
#> 2         bmi 5 -0.08122935 0.001629126 0.0003996515 14    23 21.22865
#>         riv       lambda       fmi
#> 1 0.3474068 2.863880e-02 0.3209004
#> 2 0.2943799 3.425585e-05 0.2912026
#>          term    estimate std.error   statistic       df   p.value
#> 1 (Intercept)  3.91587204  3.741657  1.04656083 20.53591 0.3074648
#> 2         bmi -0.08122935  3.741657 -0.02170946 21.22865 0.9828825

Created on 2022-11-03 with reprex v2.0.2