AUC of logistic and ordinal model following multiple imputation using MICE (with R)

Question

I am asking a question concerning the additive predictive benefit of the inclusion of a variable to a logistic and an ordinal model. I am using mice to impute missing covariates and am having difficulty finding ways to calculate the AUC and R squared of the pooled imputed models. Does anyone have any advice?

The summary readout only provides the term, estimate, std.error, statistic, df , p.value

Example code:

imputed_Data <- mice(Cross_sectional, m=10, predictorMatrix=predM, seed=500, method = meth)
Imputedreferecemodel <- with(imputed_Data, glm(Poor ~ age + sex + education + illness + injurycause, family = "binomial", na.action=na.omit) )
summary(pool(Imputedreferecemodel))

Many thanks.

Answer 1

You could use the psfmi package in combination with mice. You could use the function pool_performance to measure performances for logistic regression, according to documentation:

pool_performance Pooling performance measures for logistic and Cox regression models.

I use the nhanes dataset which is standard in mice to show you a reproducible example.

# install.packages("devtools")
# devtools::install_github("mwheymans/psfmi") # for installing package
library(psfmi)
library(mice)

# Make reproducible data with 0 and 1 outcome variable
set.seed(123)
nhanes$hyp <- ifelse(nhanes$hyp==1,0,1)
nhanes$hyp <- as.factor(nhanes$hyp)

# Mice
imp <- mice(nhanes, m=5, maxit=5) 

nhanes_comp <- complete(imp, action = "long", include = FALSE)

pool_lr <- psfmi_lr(data=nhanes_comp, nimp=5, impvar=".imp", 
                    formula=hyp ~ bmi, method="D1")
pool_lr$RR_model
#> $`Step 1 - no variables removed -`
#>          term    estimate std.error   statistic       df   p.value        OR
#> 1 (Intercept) -0.76441322 3.4753113 -0.21995532 16.06120 0.8286773 0.4656071
#> 2         bmi -0.01262911 0.1302484 -0.09696177 15.79361 0.9239765 0.9874503
#>      lower.EXP upper.EXP
#> 1 0.0002947263 735.56349
#> 2 0.7489846190   1.30184

# Check performance
pool_performance(pool_lr, data = nhanes_comp, formula = hyp ~ bmi, 
                 nimp=5, impvar=".imp", 
                 cal.plot=TRUE, plot.indiv="mean", 
                 groups_cal=4, model_type="binomial")
#> Warning: argument plot.indiv is deprecated; please use plot.method instead.

#> $ROC_pooled
#>                     95% Low C-statistic 95% Up
#> C-statistic (logit)  0.2731      0.5207 0.7586
#> 
#> $coef_pooled
#> (Intercept)         bmi 
#> -0.76441322 -0.01262911 
#> 
#> $R2_pooled
#> [1] 0.009631891
#> 
#> $Brier_Scaled_pooled
#> [1] 0.004627443
#> 
#> $nimp
#> [1] 5
#> 
#> $HLtest_pooled
#>        F_value    P(>F) df1      df2
#> [1,] 0.9405937 0.400953   2 31.90878
#> 
#> $model_type
#> [1] "binomial"

^{Created on 2022-12-02 with reprex v2.0.2}

Answer 2

When conducting logistic regression, I believe that it's good practice is to use McFadden's or Tjur's R2, since both of those tend to be used with generalized linear models. mice::pool.r.squared is designed to be only with lm models. A previous StackOverflow user had the same question as you and it seems that the best function for a multiply-imputed glm() model is mfc() from the Github package glmice. The function looks fairly simple and uses McFadden's R2, although the package hasn't been touched for a few years. That previous user wasn't able to get mfc() to work, but it worked for me.

# install.packages("remotes")
# remotes::install_github("noahlorinczcomi/glmice")
library(glmice)
library(mice)
data(nhanes)
nhanes$hyp <- ifelse(nhanes$hyp == 2, 1, 0)
imp <- mice(nhanes, m = 10, seed = 500, printFlag = FALSE)
mod <- with(imp, glm(hyp ~ age + bmi, family = "binomial"))
# summary(pool(mod))
mcf(mod)
#> [1] "34.9656%"

It looks like there are fewer resources on calculating AUC for a multiply-imputed glm(). I did find a vignette from the finalfit package, which calculated area under the curve. Unfortunately, it calculated AUC for each imputation. There might be a way to pool the output, but I'm not sure how (hopefully another SO user might suggest an idea?).

library(finalfit)
mod %>% 
  getfit() %>% 
  purrr::map(~ pROC::roc(.x$y, .x$fitted)$auc)
# not pasting the output because it's a lot

small update

As of 1/23/23, I've noticed that the glmice Github page has been taken down. I'm posting the mcf() function here for reference.

#' Calculates McFadden's Pseudo R-Squared
#'
#' Returns McFadden's pseudo r-squared for logistic regression models performed on 'mice'-imputed data sets.
#' @param model a logit model from which you would like to return McFadden's pseudo r-squared. This can be a model created either with 'glm.mids()' or 'with()'
#' @return mcfs2: McFadden's pseudo r-squared
#' @export

mcf <- function (model) {
  iterations <- model$call1$m
  null_ds <- as.numeric()
  res_ds <- as.numeric()
  for (i in 1:iterations) {
    null_ds[i] <- model$analyses[[i]]$null.deviance
    res_ds[i] <- model$analyses[[i]]$deviance
  }
  ds <- cbind(as.numeric(null_ds), as.numeric(res_ds))
  m_null <- mean(null_ds)
  m_res <- mean(res_ds)
  mcfs <- round(((1 - (m_res / m_null)) * 100), 4)
  mcfs <- paste0(mcfs, "%")
  # end
  return(mcfs)
}