How to add a q value (adjusted p-value) to a modelsummary table after pooling the results of a multinomial model over multiple imputed datasets

Question

I am using modelsummary to display the results of several multinomial models, each pooled over 5 datasets using the mice::pool function. It works great, but I want to add the q-value / adjusted p-value for false discovery rate.

I understand I need to create a tidy_custom.mipo function to add this statistic but I can't get it to work.

Below is the code to get the 'pool_univariate' list of mipo objects, which I then pass to modelsummary. It works great, I just want to add the q-value statistic.

Any idea how to do that?

Thanks a lot!

# list of exposures 
exposures <- c(
    Cs(exposure1,exposure2,exposure3)

## model function
models <- function(x) {
    lapply(imputed_data, function(y)
        multinom(as.formula(
            paste0(
                "outcome ~ ",
                x
            )
        ), data = y, model = TRUE)
    )
}

## run models
models_univariate <- as.list(seq(1,length(exposures))) 
models_univariate <- pblapply(exposures, models)

## pool 

pool_univariate <- as.list(seq(1,length(exposures))) 

# run pool
for(j in seq_along(exposures)) {
    pool_univariate[[j]] <- pool(models_univariate[[j]])
}

Answer 1

It is difficult to answer this question without a minimal working example. Here I give a simpler example than the original, for the linear regression context.

First, load the package and estimate a regression model:

library(modelsummary)
mod <- lm(mpg ~ hp + drat + vs + am, data = mtcars)

Second, since we want to summarize a model of class lm, we define a new method called tidy_custom.lm. This function takes a statistical model as input, and returns a data frame that conforms to the broom package specification, with one column called term and other columns containing matching statistics. In the current example, the data frame will include three new statistics (q.value, bonferroni and holm). These values are computed using R’s p.adjust function, which adjusts p values for multiple comparison:

tidy_custom.lm <- function(x, ...) {
  out <- broom::tidy(x)
  out$q.value <- p.adjust(out$p.value, n = 10, method = "fdr")
  out$bonferroni <- p.adjust(out$p.value, n = 10, method = "bonferroni")
  out$holm <- p.adjust(out$p.value, n = 10, method = "holm")
  return(out)
}

Now, we can call modelsummary with our lm model, and request the statistics:

modelsummary(mod, statistic = "q.value")

We can also compare different p values and label them nicely using glue strings:

modelsummary(mod,
  statistic = c(
    "p = {p.value}",
    "q = {q.value}",
    "p (Bonferroni) = {bonferroni}",
    "p (Holm) = {holm}"
  )
)