library(mice)
#> 
#> Attaching package: 'mice'
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> The following objects are masked from 'package:base':
#> 
#>     cbind, rbind
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
# example with manual and automatic pooling
imp <- mice(nhanes, 
            maxit = 2, 
            m = 2, 
            print = FALSE, 
            seed = 123)
fit <- with(data = imp, lm(bmi ~ age))

# manual pooling
summary(fit$analyses[[1]]) # the first analysis
#> 
#> Call:
#> lm(formula = bmi ~ age)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -6.9316 -3.5058  0.1684  2.3200  7.9684 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  28.4575     1.8159  15.671 9.11e-14 ***
#> age          -1.1258     0.9365  -1.202    0.242    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.811 on 23 degrees of freedom
#> Multiple R-squared:  0.05912,    Adjusted R-squared:  0.01822 
#> F-statistic: 1.445 on 1 and 23 DF,  p-value: 0.2415
summary(fit$analyses[[2]]) # the second analysis
#> 
#> Call:
#> lm(formula = bmi ~ age)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -9.5297 -1.5094  0.1703  2.2906  5.3703 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   33.290      1.827  18.225 3.63e-15 ***
#> age           -3.360      0.942  -3.567  0.00164 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 3.833 on 23 degrees of freedom
#> Multiple R-squared:  0.3562, Adjusted R-squared:  0.3282 
#> F-statistic: 12.72 on 1 and 23 DF,  p-value: 0.001637
pool.scalar(Q = c(-1.1258, -3.360), # estimates
            U = c(0.9365^2, .942^2), # variances (SE squared)
            n = 25, 
            k = 2) %>% unlist()
#>          m      qhat1      qhat2         u1         u2       qbar       ubar 
#>  2.0000000 -1.1258000 -3.3600000  0.8770322  0.8873640 -2.2429000  0.8821981 
#>          b          t         df          r        fmi 
#>  2.4958248  4.6259354  1.1087230  4.2436468  0.9021233

# automatic pooling using broom
# differences are due to rounding by summary()
pool(fit)
#> Class: mipo    m = 2 
#>          term m  estimate      ubar         b        t dfcom        df      riv
#> 1 (Intercept) 2 30.873671 3.3169948 11.675886 20.83082    23 0.9973228 5.280029
#> 2         age 2 -2.242995 0.8821795  2.496047  4.62625    23 1.1086620 4.244115
#>      lambda       fmi
#> 1 0.8407651 0.9204359
#> 2 0.8093100 0.9021334