Tbl_regression from the gtsummary package for negative binomial regressions

To assess if there is an association between certain groups of patients (patient_group; categorical variable) and a disease (disease_outcome; count variable) I am running negative binomial regression models (due to overdispersion). To check for confounding by other variables I am running 3 models with increasing amounts of covariates.

To display the IRRs and CIs i want to use the tbl_regression function from the package gtsummary (I am using the latest version 1.3.7.9022). However, calling the function returns the IRR and the corresponding 95% CIs non-exponentiated, even though I put exponentiate=TRUE:

# Load packages
library(haven)
library(magrittr)
library(MASS)
library(dplyr)
install.packages("gtsummary")
remotes::install_github("ddsjoberg/gtsummary")
library(gtsummary)

# Load example data.
dat <- read_dta("https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta")
dat <- within(dat, {
  prog <- factor(prog, levels = 1:3, labels = c("General", "Academic", "Vocational"))
  id <- factor(id)
})


# Run negative binomial regression and pipe in the tbl_regression function
Model 1 <-  
  glm.nb(data=dat, formula=daysabs ~ prog) %>%
  tbl_regression(exponentiate=TRUE) 

Model 1

This returns the summary table, but the regression coefficients have not been exponentiated. Is there a way to get gtsummary to return exponentiated coefficients and CIs?

Thanks!

I was just doing some poking around to see what is going on. The tbl_regression() function uses broom::tidy() in the background. Support for negbin models was just added 7 days ago, but for some reason an exponentiate= argument was not added for this type of model.

I am going to request that it be added. In the meantime, this code should get you up and going with negbin models.

library(gtsummary)
library(tidyverse)

# add a custom tidying function
my_negbin_tidy <- function(x, exponentiate = FALSE, ...) {
  df_tidy <- broom::tidy(x, ...)
  
  # exponentiate coef and CI if requested
  if (exponentiate) {
    df_tidy <-
      df_tidy %>%
      mutate_at(vars(any_of(c("estimate", "conf.low", "conf.high"))), exp)
  }
  
  df_tidy
}

# build model
mod <- MASS::glm.nb(response ~ age, gtsummary::trial) 

# summarize model results
tbl <- 
  tbl_regression(
    mod, 
    exponentiate = TRUE,
    tidy_fun = my_negbin_tidy
  ) 

^{Created on 2021-04-12 by the reprex package (v2.0.0)}