Presentation-Ready Analytic Result Tables Using gtsummary

library(tidyverse)
library(gtsummary)
library(survival)
data = trial %>% # build in sample dataset in gtsummary package
  mutate(trt = factor(trt))

tbl1 <- data %>%
  select(age, grade, trt) %>% # select outcome and variables you want to include in order
  tbl_summary(
    by = trt, 
    type = all_continuous() ~ "continuous2", # update the continuous variables to summarize on multiple lines, default is median(Q1,Q3)
    statistic = all_continuous() ~ c(
      "{N_nonmiss}",
      "{median} ({p25}, {p75})",
      "{min}, {max}"
    ),
    missing = "ifany", # missing = "no", whether to display a row with the number of missing observations
    missing_text = "(Missing)", # text label for the missing number row
    missing_stat = "{N_miss}/{N_obs} ({p_miss})%)", # stats of missing data
    digits = list(all_categorical() ~ c(0, 1), all_continuous2() ~ c(0, 2, 2)), # digits of stats, can simply be ~1
    percent = "column") %>% # Default is columnwise percentage, could be set as , percent = "row"
  add_p(
    list(all_continuous() ~ "t.test", all_categorical() ~ "chisq.test"),
    pvalue_fun = function(x) style_pvalue(x, digits = 3)
    ) %>% 
  # add header to p.value
  modify_header(p.value ~ "**Parametric p-value**")

## To check the header names
# show_header_names(tbl1)

tbl1 %>% 
  as_kable_extra()

tbl2 <- data  %>% 
  select(age, grade, trt) %>% 
  tbl_summary(by = trt) %>%
  add_p(
    list(all_continuous() ~ "wilcox.test", all_categorical() ~ "fisher.test"),
    pvalue_fun = function(x) style_pvalue(x, digits = 3)
    ) %>% 
  # hide all columns, except p-value
  modify_column_hide(-p.value) %>%
  # add header to p.value
  modify_header(p.value ~ "**Non-Parametric p-value**")

tbl3 <- data %>%
  select(age, grade, marker, trt) %>%
  tbl_uvregression(
    y = trt, # y should be a factor
    method = glm, # specify the regression function
    method.args = list(family = binomial), # specify the outcome type
    # adjust all models for marker level
    formula = "{y} ~ {x} + marker", # remove "marker" to obtain regular univariable models
    exponentiate = TRUE, 
    pvalue_fun = function(x) style_pvalue(x, digits = 3)
  ) %>%
  add_global_p(keep = FALSE) %>% # adds the global p-value for categorical variables, use car::Anova() by default,  # Likelihood ratio test (optional): Anova(fit, test = "LR")
  # hide all columns, except p-value
  modify_column_hide(-p.value) %>%
  # add header to p.value
  modify_header(p.value ~ "**Adjusted GLM p-value**") %>%
  modify_footnote(p.value ~ "Logistic regression adjusted for marker level")

tbl_merge(list(tbl1, tbl2, tbl3)) %>%
  modify_spanning_header(everything() ~ NA) %>% 
  remove_footnote_header(columns = all_stat_cols()) %>% 
  modify_abbreviation(c("Adding abbreviation on top of the existing ones, the function will reorder them by alphabetical"))

# set theme to get MEAN (SD) by default in `tbl_summary()`
# theme_gtsummary_mean_sd()

# function to add pairwise copmarisons to `tbl_summary()`
add_stat_pairwise <- function(data, variable, by, ...) {
  # calculate pairwise p-values
  pw <- pairwise.t.test(data[[variable]], data[[by]], p.adj = "BH") # # c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")

  # convert p-values to list
  index <- 0L
  p.value.list <- list()
  for (i in seq_len(nrow(pw$p.value))) {
    for (j in seq_len(nrow(pw$p.value))) {
      index <- index + 1L
      
      p.value.list[[index]] <- 
        c(pw$p.value[i, j]) %>%
        setNames(glue::glue("**{colnames(pw$p.value)[j]} vs. {rownames(pw$p.value)[i]}**"))
    }
  }
  
  # convert list to data frame
  p.value.list %>% 
    unlist() %>%
    purrr::discard(is.na) %>%
    t() %>%
    as.data.frame() %>%
    # formatting/roundign p-values
    dplyr::mutate(dplyr::across(everything(), style_pvalue))
}

data %>%
  select(grade, age, marker) %>%
  tbl_summary(by = grade, missing = "no", digits = ~1) %>%
  # add pariwaise p-values
  add_stat(everything() ~ add_stat_pairwise)

univariable_tbl = tbl_uvregression(
    data %>% dplyr::select(ttdeath, death, trt, age, stage, grade),
    method = coxph,
    y = Surv(ttdeath, death),
    exponentiate = TRUE,
    label = list(trt = "Treatment", age = "Age", 
                 stage = "Stage", grade = "Grade"), # Can also label in the tbl_function
    pvalue_fun = function(x) style_pvalue(x, digits = 3)
    ) %>% 
  add_global_p(include=c("stage", "grade"), keep=FALSE) %>% 
  modify_column_merge( # Combine the estimate and CI columns
    pattern = "{estimate} [{conf.low}, {conf.high}]",
    rows = !reference_row %in% TRUE & !is.na(conf.low)
  ) %>% 
  modify_table_styling( # Add REF to reference level
    columns = c(estimate),
    rows = reference_row %in% TRUE,
    missing_symbol = "REF"
  ) %>% 
  modify_header(estimate = "**HR [95% CI]**") %>% 
  modify_column_hide("stat_n") # remove the N column
  # modify_caption(caption = "PFS of Cohort A+C") %>% 
  # modify_spanning_header(c(stat_n, estimate, conf.low, p.value) ~ "**Univariable**")

multivariable_tbl <- coxph(
   Surv(ttdeath, death) ~ age + stage, data = data
   ) %>%
  tbl_regression(
      exponentiate=TRUE,
      label = list(age = "Age",
                   stage = "Stage")
   ) %>% 
  add_global_p(include=c("stage"), keep=FALSE) %>% 
  modify_column_merge(
    pattern = "{estimate} [{conf.low}, {conf.high}]",
    rows = !reference_row %in% TRUE & !is.na(conf.low)
  ) %>% 
  modify_table_styling(
    columns = c(estimate),
    rows = reference_row %in% TRUE,
    missing_symbol = "REF"
  ) %>% 
  modify_header(estimate = "**HR [95% CI]**")

tbl_merge(tbls = list(univariable_tbl, multivariable_tbl),
                     tab_spanner = c("**Univariable Models**", "**Multivariable Model**"))

## To check the header names
# show_header_names(pfs_os_tbl)


pfs_os_tbl = tbl_merge(tbls = list(pfs_table, os_table),
                     tab_spanner = FALSE) %>% 
  modify_caption(caption = "PFS and OS table.") %>% 
  modify_abbreviation(c("PFS = Progression Free Survival, OS = Overall Survival")) %>% 
  as_gt() %>% 
  # gt::cols_width(
  #   label	~ px(180),
  #   everything() ~ px(80)
  # ) %>% 
  gt::tab_spanner(
    columns = c(estimate_1_1, conf.low_1_1, p.value_1_1, estimate_2_1, conf.low_2_1, p.value_2_1),
    label = "PFS",
    gather = FALSE
  ) %>% 
  gt::tab_spanner(
    columns = c(estimate_1_2, conf.low_1_2, p.value_1_2, estimate_2_2, conf.low_2_2, p.value_2_2),
    label = "OS",
    gather = FALSE
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(weight = "bold"),
    locations = list(
      gt::cells_column_spanners(matches(c("PFS","OS")))
      )
    ) %>% 
  gt::tab_options(latex.tbl.pos = "ht")

as_gt()
as_kable_extra()
as_flex_table()

t.test(data$age ~ data$trt)
wilcox.test(data$age ~ data$trt)

chisq.test(table(data$grade, data$trt))
fisher.test(table(data$grade, data$trt))

summary(glm(trt ~ age + marker, data = data, family = "binomial"))

fit = glm(trt ~ grade + marker, data = data, family = "binomial")
car::Anova(fit)

summary(coxph(Surv(ttdeath, death) ~ stage, data = data))

sessionInfo()