Search code examples
routputregressionimputationr-mice

Reformat results of pooled imputation results for regression in R

I am using multiple imputation on missing data and then using the pool_mi function to fit the results of the imputation trials to my regression model. However the output is not formatted in an easily interpretable way and I am hoping to get some guidance on how to do so. Below is the example code for what I have done with the function and images of the output vs my desired output.

library(mitools)
data(data.ma05)
dat <- data.ma05

# imputation of the dataset: use six imputations
resp <- dat[, - c(1:2) ]
imp <- mice::mice( resp, method="norm", maxit=3, m=6 )
datlist <- miceadds::mids2datlist( imp )

# linear regression with cluster robust standard errors
mod <- lapply(  datlist, FUN=function(data){
            miceadds::lm.cluster( data=data, formula=denote ~ migrant+ misei,
                    cluster=dat$idclass )
            }  )
# extract parameters and covariance matrix
betas <- lapply( mod, FUN=function(rr){ coef(rr) } )
vars <- lapply( mod, FUN=function(rr){ vcov(rr) } )
# conduct statistical inference
summary( miceadds::pool_mi( qhat=betas, u=vars ) )

Here is the output table I get. enter image description here

But is there a way either using stargazer or some other function/package where I can reformat my results so the p-values are rounded and have stars on the side to indicate significance? The picture below demonstrates more of how I want the output to look. I'm aware that the regression in this image is entirely different function/variables/data, but am including it to provide clarity as to how I would like the output for the p-values to appear. Thank you!

enter image description here

Solution

  • I do not know how to achieve this in stargazer, but it is very easy to do with the modelsummary package for R (disclaimer: I am the author).

    modelsummary supports 100+ models out-of-the-box, but not model objects of class pool_mi, which is what your code produces. Fortunately, it is very easy to add support for new models, as described in detail in the documentation.

    Specifically, we need to define two S3 methods called tidy.CLASSNAME and glance.CLASSNAME. The first method returns a data.frame with one coefficient per row and with column names that follow the standard terminology from the broom package. The second method returns a one-row data.frame with goodness-of-fit or other model characteristics, one per column.

    In your case, these simple methods seem to do the job (obviously, you can customize):

    library(modelsummary)

tidy.pool_mi <- function(x, ...) {
  msg <- capture.output(out <- summary(x, ...))
  out$term <- row.names(out)
  colnames(out) <- c("estimate", "std.error", "statistic", "p.value",
                     "conf.low", "conf.high", "miss", "term")
  return(out)
} 

glance.pool_mi <- function(x, ...) {
  data.frame(nimp = x$m)
}

    When we call these methods on the output of pool_mi, we get this:

    mod_pooled <- miceadds::pool_mi( qhat=betas, u=vars )


glance(mod_pooled)
#>   nimp
#> 1    6

tidy(mod_pooled)
#>               estimate   std.error statistic      p.value    conf.low
#> (Intercept)  2.5875292 0.085702631 30.191946 3.808654e-55  2.41769205
#> migrant      0.5840870 0.087258478  6.693756 1.051145e-10  0.41237872
#> misei       -0.0140385 0.001661665 -8.448452 2.090746e-12 -0.01735046
#>               conf.high   miss        term
#> (Intercept)  2.75736637 22.7 % (Intercept)
#> migrant      0.75579526 13.4 %     migrant
#> misei       -0.01072654 28.2 %       misei

    Finally, as soon as these methods are defined, the modelsummary package should work immediately:

    modelsummary(mod_pooled)
    Model 1
    (Intercept) 2.571
    (0.089)
    migrant 0.571
    (0.086)
    misei -0.014
    (0.002)
    Num.Imp. 6