Search code examples
rgtsummary

Summary of numbers at risk and events for Poisson regression using gtsummary


I have a summary table of IRRs and 95% CI following univariable and multivariable Poisson regression, created with gtsummary that looks a bit like this.

For a logistic regression model it is pretty straightforward to use tbl_summary to create some count data to append to the left of the table. However for a Poisson model I'd like to be able to sum, rather than count, the number of days at risk and the number of events. Each row of the underlying dataset contains a number of days at risk and a number of events, such that the regression model is run like this:

glm(events ~ study_arm + strata_group, 
     offset = log(days_at_risk), 
     family=poisson(link = "log"), 
     data = df)

Is it possible to use gtsummary to create two columns with sums of the number of events and the number of days at risk for each of the rows of the table? (This can then be added to my table using tbl_merge.)

Here is a more complete example of what I'd like to achieve

df = tibble(
  study_arm = c("control", "intervention", "control", "intervention", "control", "intervention", "control", "intervention"),
  events = c(3,4,12,6,0,3,11,9),
  strata_group = c("A", "A", "A", "A", "B", "B", "B", "B"),
  days_at_risk = c(100,100,200,200,300,300,100,100)
)

m=glm(events ~ study_arm + strata_group, 
     offset = log(days_at_risk), 
     family=poisson(link = "log"), 
     data = df)

tbl_regression(m, exponentiate = T)

#this is the summary I wish to be able to generate with tbl_summary so I can merge it with the tbl_regression output
bind_rows(
  df %>% group_by(study_arm) %>% 
  summarise(n_events = sum(events), 
            total_days_at_risk = sum(days_at_risk), 
            rate=n_events/total_days_at_risk) %>% 
  mutate(row_group = "study_arm") %>% rename(characteristic=study_arm),
df %>% group_by(strata_group) %>% 
  summarise(n_events = sum(events), 
            total_days_at_risk = sum(days_at_risk), 
            rate=n_events/total_days_at_risk) %>% 
  mutate(row_group = "strata_group") %>% rename(characteristic=strata_group)
) %>% 
  select(row_group, characteristic, n_events, total_days_at_risk, rate)

Solution

  • Hi and welcome to stackoverflow!

    Here's an example of how to get the table you're looking for.

    1. Use the add_nevent() function to get the sum of the number of observed events.
    2. The sum of the exposure times is already in the table (.$table_body). Add a column header to unhide the exposure coulmn.
    3. Calculate the rate, then assign a column header and a formatting function.

    Happy Programming!

    library(gtsummary)
    library(tidyverse)
    
    df <- tibble(
      study_arm = c("control", "intervention", "control", "intervention", "control", "intervention", "control", "intervention"),
      events = c(3, 4, 12, 6, 0, 3, 11, 9),
      strata_group = c("A", "A", "A", "A", "B", "B", "B", "B"),
      days_at_risk = c(100, 100, 200, 200, 300, 300, 100, 100)
    )
    
    m <- 
      glm(events ~ study_arm + strata_group,
          offset = log(days_at_risk),
          family = poisson(link = "log"),
          data = df
      )
    
    tbl <-
      tbl_regression(m, exponentiate = T) %>%
      # add sum of the number events
      add_nevent(location = "level") %>%
      # add the sum of the exposure times.
      # this column is present in the table by default, but the column is hidden
      # adding the column header unhides the column
      modify_header(exposure ~ "**Exposure**") %>%
      # calculate the rate and add to tbl
      # after the column is added to the table, we need to add
      # a column header and tell gtsummary how to format the new column
      modify_table_body(
        ~.x %>%
          mutate(rate = stat_nevent / exposure, 
                 .after = stat_nevent)
      ) %>%
      modify_header(rate ~ "**Rate**") %>%
      modify_fmt_fun(rate ~ partial(style_percent, symbol = TRUE))
    

    enter image description here

    Created on 2021-07-13 by the reprex package (v2.0.0)