Meta-Analysis: Cox-PH Sensitivity Analysis

This notebook evaluates the performance of Cox proportional hazard models constructed using three censor cutoff periods (30, 60, and 90 days) and three methods of adjusting for pre-admission health conditions (i.e., pre-existing comorbidity burden): (1) the inclusion of the 29 individual covariates (i.e., health conditions) comprising the Elixhauser Comorbidity Index, (2) the Elixhauser summary score, and (3) the top 10 principal components computed from logistic principal component analysis (LPCA).

We also perform additional meta-analyses in order to evaluate the estimates of the individual covariates (e.g., age, comoridities, etc) on each outcome.

Of note, this analysis is performed only on adult patients due to the pediatric cohort’s low incidence of both poor health outcomes and neurological diagnoses during COVID-19 hospitalization.

library(tidyverse)
library(metafor)
library(meta)
library(cowplot)
library(DT)
library(ggbreak) 
library(RColorBrewer)
library(viridis)
devtools::install_github("rdboyes/forester")
devtools::install_github("thomasp85/patchwork")
source("R/forester_custom.R")
source("R/analyze_survival.R")

Import Data from each Healthcare system

# read in files from results folder 
# this folder contains all of the local healthcare system level analyses
rdas <- list.files(
  path = "results",
  pattern = ".rda",
  full.names = TRUE
)
for (rda in rdas) {
  load(rda)
}

rm(rdas, rda)

# create a list of participating healthcare systems from our study tracking spreadsheet
site_google_url <- "https://docs.google.com/spreadsheets/d/1epcYNd_0jCUMktOHf8mz5v651zy1JALD6PgzobrGWDY/edit?usp=sharing"

# load site parameters
site_params <- googlesheets4::read_sheet(site_google_url, sheet = 1)
site_avails <- googlesheets4::read_sheet(site_google_url, sheet = 2)

# filter the list of sites who ran the analysis
sorted_sites <- site_avails %>%
  filter(!is.na(date_v4_received)) %>%
  pull(siteid) %>%
  paste("results", sep = "_")

# list sites without race
sites_wo_race <- site_params %>%
  filter(!include_race) %>%
  pull(siteid)

# combine all rda files with 'results' in name
results <- mget(ls(pattern = "results"))

## read in pt counts
pt_counts_df <- read.csv('tables/site_pt_counts.csv')

Identify sites to include in the analysis

We will only include a site if they have >= 3 neuro patients

sites_adult <- pt_counts_df %>% 
  filter(population == "Adult") %>% 
  mutate(neuro_sum = n_var_Central + n_var_Peripheral) %>% 
  filter(!neuro_sum < 3) %>% 
  distinct(site) %>% 
  mutate(site = gsub("_results", "", site))

sites_adult_results <- sites_adult %>% mutate(site = paste0(site, "_results")) %>% as.vector()

adult_results <- results[grepl(paste(sites_adult_results$site, collapse = "|"), names(results))]

Define parameters

outcomes <- c("time_first_discharge_reg_elix", "deceased_reg_elix")

comorb_adj <- c("lpca", "score", "ind")

censor_cut <- c("30", "60", "90")

Get Cox model results - Adults

cox_results_adults <- list()

for (outcome_i in outcomes) {
  for (comorb_adj_i in comorb_adj) {
    for (censor_cut_i in censor_cut) {
      
      cox_results_adults[[outcome_i]][[comorb_adj_i]][[censor_cut_i]] <-
        adult_results %>%
        lapply(get_cox_row, population = "adults", comorb_method = comorb_adj_i, censor_cutoff = censor_cut_i, outcome = outcome_i) %>%
        bind_rows() %>% 
        mutate(outcome = outcome_i,
               comorb_method = comorb_adj_i,
               censor_cutoff = censor_cut_i)
    }
  }
}

Evaluate concordance

cox_concordance_results_adults <- list()

for (outcome_i in outcomes) {
  for (comorb_adj_i in comorb_adj) {
    for (censor_cut_i in censor_cut) {
      
      cox_concordance_results_adults[[outcome_i]][[comorb_adj_i]][[censor_cut_i]] <-
        adult_results %>%
        lapply(get_summary_stats, population = "adults", comorb_method = comorb_adj_i, censor_cutoff = censor_cut_i, outcome = outcome_i, cox_stat = "concordance") %>%
        bind_rows() 
    }
  }
}

Bind all concordance results together

bind_results <-  function(results_list) {
  
  lpca_discharge <- results_list[["time_first_discharge_reg_elix"]][["lpca"]] %>% bind_rows()

  score_discharge <- results_list[["time_first_discharge_reg_elix"]][["score"]] %>% bind_rows()
  
  ind_discharge <- results_list[["time_first_discharge_reg_elix"]][["ind"]] %>% bind_rows()
  
  lpca_mortality <- results_list[["deceased_reg_elix"]][["lpca"]] %>% bind_rows()
  
  score_mortality <- results_list[["deceased_reg_elix"]][["score"]] %>% bind_rows()
  
  ind_mortality <- results_list[["deceased_reg_elix"]][["ind"]] %>% bind_rows()
  
  combined_results <- bind_rows(lpca_discharge, score_discharge, ind_discharge,
                                   lpca_mortality, score_mortality, ind_mortality)
  
  return(combined_results)
  
  
}

cox_results <- bind_results(results_list = cox_results_adults)
concordance_results <- bind_results(results_list = cox_concordance_results_adults)

Concordance

c_plot <- ggplot(concordance_results %>% 
         mutate(comorb_method = factor(comorb_method,
                                       levels = c("ind", "score", "lpca")) %>% 
                  fct_recode(
                    `Individual\nCovariates` = "ind",
                    `Summary\nScore` = "score",
                    `LPCA` = "lpca"
                    ),
                outcome = as.factor(outcome) %>% 
                  fct_recode(
                    Mortality = "deceased_reg_elix",
                    Discharge = "time_first_discharge_reg_elix"),
                censor_cutoff = as.factor(censor_cutoff) %>% 
                  fct_recode(
                    `30 Days` = "30",
                    `60 Days` = "60", 
                    `90 Days` = "90")), 
       aes(x = comorb_method, y = C, group = comorb_method, fill = comorb_method)) + 
  geom_violin(trim = FALSE) + 
  facet_grid(outcome ~ censor_cutoff, scales = "free") + 
  scale_y_continuous("Concordance", position="left") +   
  xlab("") +
  theme_bw() +
  theme(legend.position = "none") +
  scale_fill_brewer(palette="Dark2") +
  geom_boxplot(width=0.3, fill="white", outlier.shape=NA) +
  geom_jitter(size = 0.4, alpha = 0.5, color = "black") +
  theme(strip.text.x = element_text(size = 15, face = "bold"),
        strip.text.y = element_text(size = 15, face = "bold"),
        axis.title.y = element_text(size = 15, face = "bold"),
        axis.text.x = element_text(size = 13, face = "bold", color = "black"),
        axis.text.y = element_text(size = 12, face = "bold"))

c_plot

ggsave("figures/Comorbidity_adjustment.png", c_plot, width = 12, height = 7, dpi = 300)

c_results <- concordance_results %>% 
  group_by(outcome, censor_cutoff, comorb_method) %>% 
  summarize(min = min(C),
            q1 = quantile(C, 0.25),
            median = median(C),
            mean = mean(C),
            q3 = quantile(C, 0.75),
            max = max(C)) %>% 
    select(outcome, censor_cutoff, comorb_method, median, mean)

## `summarise()` has grouped output by 'outcome', 'censor_cutoff'. You can
## override using the `.groups` argument.

write.csv(c_results, "tables/comorbidity_adj_c_results.csv", row.names=FALSE)


# summary stats
datatable(concordance_results %>% 
  group_by(outcome, censor_cutoff, comorb_method) %>% 
  summarize(min = min(C),
            q1 = quantile(C, 0.25),
            median = median(C),
            mean = mean(C),
            q3 = quantile(C, 0.75),
            max = max(C)),
  filter="top") 

## `summarise()` has grouped output by 'outcome', 'censor_cutoff'. You can
## override using the `.groups` argument.

Evaluate comorb_method with highest OR for CNS and PNS

# summary stats
datatable(cox_results %>% 
            filter(variable == "neuro_postCentral" | variable == "neuro_postPeripheral") %>% 
            group_by(outcome, variable, censor_cutoff, comorb_method) %>% 
            mutate(exp.coef. = round(exp.coef.)) %>% 
            summarize(min = min(exp.coef., na.rm = TRUE),
                      q1 = quantile(exp.coef., 0.25, na.rm = TRUE),
                      median = median(exp.coef., na.rm = TRUE),
                      mean = mean(exp.coef. , na.rm = TRUE),
                      q3 = quantile(exp.coef., 0.75, na.rm = TRUE),
                      max = max(exp.coef., na.rm = TRUE)),
          filter="top")

## `summarise()` has grouped output by 'outcome', 'variable', 'censor_cutoff'. You
## can override using the `.groups` argument.

Random effects meta-analysis - CNS

set sm = "HR" when estimate is logHR (which is coef in our model): https://cran.r-project.org/web/packages/meta/meta.pdf - (search logHR in cran to see example)

*Note: We excluded sites with < 3 neuro patients in a category

## adults 
meta_results_adults_cns <- list()

for (outcome_i in outcomes) {
  for (comorb_adj_i in comorb_adj) {
    for (censor_cut_i in censor_cut) {
      
      meta_results_adults_cns[[outcome_i]][[comorb_adj_i]][[censor_cut_i]] <-
        cox_results %>% 
        filter(outcome == outcome_i,
               comorb_method == comorb_adj_i,
               censor_cutoff == censor_cut_i
               ) %>% 
        bind_rows() %>%
        data.frame() %>%
        filter(variable == "neuro_postCentral") %>%
        metagen(
          TE = coef,
          seTE = se.coef.,
          data = .,
          sm = "HR", # hazard ratios
          comb.random = TRUE,
          comb.fixed = FALSE,
          method.tau = "DL", # default tau method
          hakn = FALSE,
          prediction = TRUE,
          studlab = site
        ) 
    }
  }
}

Random effects meta-analysis - PNS

## adults
meta_results_adults_pns <- list()

for (outcome_i in outcomes) {
  for (comorb_adj_i in comorb_adj) {
    for (censor_cut_i in censor_cut) {
      
      meta_results_adults_pns[[outcome_i]][[comorb_adj_i]][[censor_cut_i]] <-
        cox_results %>% 
        filter(outcome == outcome_i,
               comorb_method == comorb_adj_i,
               censor_cutoff == censor_cut_i
               ) %>% 
        bind_rows() %>%
        data.frame() %>%
        filter(variable == "neuro_postPeripheral") %>%
        metagen(
          TE = coef,
          seTE = se.coef.,
          data = .,
          sm = "HR", # hazard ratios
          comb.random = TRUE,
          comb.fixed = FALSE,
          method.tau = "DL", # default tau method
          hakn = FALSE,
          prediction = TRUE,
          studlab = site
        ) 
    }
  }
}

Format meta results

format_meta <- function(meta_results) {

  ma_combine <- list()

  for (outcome_i in outcomes) {
    for (comorb_adj_i in comorb_adj) {
      for (censor_cut_i in censor_cut) {
        ma_combine[[outcome_i]][[comorb_adj_i]][[censor_cut_i]] <- data.frame(
        TE = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["TE"]],
        #seTE =meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["seTE"]],
        studlab = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["studlab"]],
        upper = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["upper"]],
        lower = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["lower"]],
        TE.random = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["TE.random"]],
        lower.random = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["lower.random"]],
        upper.random = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["upper.random"]],
        pval.random = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["pval.random"]],
        Weight.random = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["w.random"]],
        lower.predict = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["lower.predict"]],
        upper.predict = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["upper.predict"]],
        I2 = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["I2"]],
        lower.I2 = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["lower.I2"]],
        upper.I2 = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["upper.I2"]],
        H = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["H"]],
        lower.H = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["lower.H"]],
        upper.H = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["upper.H"]],
        tau2 = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["tau2"]],
        lower.tau2 = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["lower.tau2"]],
        upper.tau2 = meta_results[[outcome_i]][[comorb_adj_i]][[censor_cut_i]][["upper.tau2"]],
        analysis = paste0(outcome_i,"_", comorb_adj_i, "_", censor_cut_i)
        ) %>% 
          mutate(outcome = outcome_i,
                 comorb_method = comorb_adj_i,
                 censor_cutoff = censor_cut_i) %>% 
          select(studlab, outcome, comorb_method, censor_cutoff, everything())
    }
  }
    }

  return(ma_combine)
  }

meta_results_adults_cns_list = format_meta(meta_results_adults_cns)
meta_results_adults_pns_list = format_meta(meta_results_adults_pns)

meta_results_cns <- bind_results(results_list = meta_results_adults_cns_list)
meta_results_pns <- bind_results(results_list = meta_results_adults_cns_list)

meta_results_cns_table <- meta_results_cns %>% 
  distinct(TE.random, outcome, comorb_method, censor_cutoff) %>% 
  group_by(outcome, censor_cutoff, comorb_method) %>% 
  mutate(TE.random = as.numeric(TE.random))

write.csv(meta_results_cns_table, "tables/comorbidity_adj_cns_metanalysis_TE.csv", row.names=FALSE)


# summary stats
datatable(meta_results_cns %>% 
            distinct(TE.random, outcome, comorb_method, censor_cutoff) %>% 
            group_by(outcome, censor_cutoff, comorb_method) %>% 
            mutate(TE.random = as.numeric(TE.random)),
          filter="top")

meta_results_pns_table <- meta_results_pns %>% 
  distinct(TE.random, outcome, comorb_method, censor_cutoff) %>% 
  group_by(outcome, censor_cutoff, comorb_method) %>% 
  mutate(TE.random = as.numeric(TE.random))

write.csv(meta_results_pns_table, "tables/comorbidity_adj_pns_metanalysis_TE.csv", row.names=FALSE)


# summary stats
datatable(meta_results_pns %>% 
            distinct(TE.random, outcome, comorb_method, censor_cutoff) %>% 
            group_by(outcome, censor_cutoff, comorb_method) %>% 
            mutate(TE.random = as.numeric(TE.random)),
          filter="top")

Hazard-Ratio

pd <- position_dodge(width = 0.6)

meta_results_cns$neuro <- "CNS"
meta_results_pns$neuro <- "PNS"

meta_mortality_results <- rbind(meta_results_cns, meta_results_pns)

meta_comorb_comparison <- ggplot(meta_mortality_results %>%
                                  mutate(comorb_method = factor(comorb_method,
                                                                levels = c("lpca", "score", "ind")) %>% 
                                          fct_recode(
                                            `LPCA` = "lpca",
                                            `Summary\nScore` = "score",
                                            `Individual\nCovariates` = "ind"
                                             ),
                                         outcome = as.factor(outcome) %>% 
                                          fct_recode(
                                            Mortality = "deceased_reg_elix",
                                            Discharge = "time_first_discharge_reg_elix"),
                                        censor_cutoff = as.factor(censor_cutoff) %>% 
                                          fct_recode(
                                            `30 Days` = "30",
                                            `60 Days` = "60", 
                                            `90 Days` = "90"), 
                                        TE.random = exp(TE.random),
                                        lower.random = exp(lower.random),
                                        upper.random = exp(upper.random)) %>%  
                                   rename(`Days since admission` = "censor_cutoff") %>% 
                                   distinct(TE.random, lower.random, upper.random, outcome, comorb_method, `Days since admission`, neuro),
                                 aes(x=TE.random, y=comorb_method, colour=`Days since admission`, group=`Days since admission`)) +
  scale_color_viridis(discrete=TRUE, direction = 1) +
  geom_point(aes(x=TE.random), shape=15, size=3, position = pd) +
  geom_linerange(aes(xmin=lower.random, xmax=upper.random), position = pd) +
  labs(x = "Hazard Ratio", y = "") +  
  geom_vline(xintercept = 1, linetype="dashed")  + 
  facet_grid(outcome ~ neuro) + 
  theme_bw() +
  theme(strip.text.x = element_text(size = 15, face = "bold"),
        strip.text.y = element_text(size = 15, face = "bold"),
        axis.title.x = element_text(size = 15, face = "bold"),
        axis.text.x = element_text(size = 13, face = "bold", color = "black"),
        axis.text.y = element_text(size = 12, face = "bold"),
       legend.title = element_text(size = 15, face = "bold"),
       legend.text = element_text(size = 13, face = "bold"))

meta_comorb_comparison

ggsave("figures/Comorbidity_adjustment_HR.png", meta_comorb_comparison, width = 12, height = 7, dpi = 300)

Evaluate individual covariates

variables <- cox_results %>% 
  filter(comorb_method == "ind") %>% 
  distinct(variable)

variables <- unique(variables$variable)

meta_results_all <- list()

for(variable_i in variables) {
  for (outcome_i in outcomes) {
        
    meta_results_all[[outcome_i]][[variable_i]] <-
      cox_results %>% 
      filter(outcome == outcome_i,
             comorb_method == "ind",
             censor_cutoff == "90",
             variable == variable_i) %>% 
      bind_rows() %>%
      data.frame() %>%
      metagen(
        TE = coef,
        seTE = se.coef.,
        data = .,
        sm = "HR", # hazard ratios
        comb.random = TRUE,
        comb.fixed = FALSE,
        method.tau = "DL", # default tau method
        hakn = FALSE,
        prediction = TRUE,
        studlab = site
      ) 
  }
}

format_meta_ind <- function(meta_results) {

ma_ind_combine <- list()

for(variable_i in variables) {
  for (outcome_i in outcomes) {
    ma_ind_combine[[outcome_i]][[variable_i]] <- data.frame(
      TE = meta_results[[outcome_i]][[variable_i]][["TE"]],
      #seTE =meta_results[[outcome_i]][[variable_i]][["seTE"]],
      studlab = meta_results[[outcome_i]][[variable_i]][["studlab"]],
      upper = meta_results[[outcome_i]][[variable_i]][["upper"]],
      lower = meta_results[[outcome_i]][[variable_i]][["lower"]],
      TE.random = meta_results[[outcome_i]][[variable_i]][["TE.random"]],
      lower.random = meta_results[[outcome_i]][[variable_i]][["lower.random"]],
      upper.random = meta_results[[outcome_i]][[variable_i]][["upper.random"]],
      pval.random = meta_results[[outcome_i]][[variable_i]][["pval.random"]],
      Weight.random = meta_results[[outcome_i]][[variable_i]][["w.random"]],
      lower.predict = meta_results[[outcome_i]][[variable_i]][["lower.predict"]],
      upper.predict = meta_results[[outcome_i]][[variable_i]][["upper.predict"]],
      I2 = meta_results[[outcome_i]][[variable_i]][["I2"]],
      lower.I2 = meta_results[[outcome_i]][[variable_i]][["lower.I2"]],
      upper.I2 = meta_results[[outcome_i]][[variable_i]][["upper.I2"]],
      H = meta_results[[outcome_i]][[variable_i]][["H"]],
      lower.H = meta_results[[outcome_i]][[variable_i]][["lower.H"]],
      upper.H = meta_results[[outcome_i]][[variable_i]][["upper.H"]],
      tau2 = meta_results[[outcome_i]][[variable_i]][["tau2"]],
      lower.tau2 = meta_results[[outcome_i]][[variable_i]][["lower.tau2"]],
      upper.tau2 = meta_results[[outcome_i]][[variable_i]][["upper.tau2"]]
      ) %>% 
        mutate(outcome = outcome_i,
               variable = variable_i) %>% 
        select(studlab, outcome, variable, everything())
  }
  }
  return(ma_ind_combine)
  }

meta_ind_results_all <- format_meta_ind(meta_results_all)

meta_ind_discharge <- bind_rows(meta_ind_results_all[['time_first_discharge_reg_elix']]) %>% 
  distinct(variable, TE.random, lower.random, upper.random, pval.random) %>% 
  mutate(TE.random = exp(TE.random),
         lower.random = exp(lower.random),
         upper.random = exp(upper.random)) %>% 
  filter(!variable ==   "DMcx") %>% 
  mutate(`Estimate` = round(TE.random, 2),
         lower.random = round(lower.random, 2),
         upper.random = round(upper.random, 2),
    "Hazard Ratio" = paste(Estimate, "(", lower.random, ",", upper.random, ")"),
    p.value_format = ifelse(pval.random > 0.01, round(pval.random, 2), round(pval.random, 3)),
    "P-value" = ifelse(pval.random < 0.001, "< .001*", as.character(p.value_format)),
    "P-value" = ifelse(pval.random >= 0.001 & round(pval.random,2) < 0.05, paste(`P-value`, "*"), `P-value`)
  ) %>%
  select(-p.value_format) %>% 
  rename(Covariate = "variable")
         

meta_ind_deceased <- bind_rows(meta_ind_results_all[['deceased_reg_elix']]) %>% 
  distinct(variable, TE.random, lower.random, upper.random, pval.random) %>% 
  mutate(TE.random = exp(TE.random),
         lower.random = exp(lower.random),
         upper.random = exp(upper.random)) %>% 
  filter(!variable ==   "DMcx") %>% 
  mutate(Estimate = round(TE.random, 2),
         lower.random = round(lower.random, 2),
         upper.random = round(upper.random, 2),
    "Hazard Ratio" = paste(Estimate, "(", lower.random, ",", upper.random, ")"),
    p.value_format = ifelse(pval.random > 0.01, round(pval.random, 2), round(pval.random, 3)),
    "P-value" = ifelse(pval.random < 0.001, "< .001*", as.character(p.value_format)),
    "P-value" = ifelse(pval.random >= 0.001 & round(pval.random,2) < 0.05, paste(`P-value`, "*"), `P-value`)
  ) %>%
  select(-p.value_format) %>% 
  rename(Covariate = "variable")

Evaluating Outcomes by Age

meta_ind_deceased_age <- meta_ind_deceased %>% 
  filter(grepl(x = Covariate, pattern = "(age|neuro)")) %>% 
  mutate(
         Covariate = as.factor(Covariate),
         Covariate = factor(Covariate,
                            levels = rev(c("neuro_postPeripheral",
                                       "neuro_postCentral",
                                       "age_group26to49",
                                       "age_group50to69",
                                       "age_group70to79",
                                       "age_group80plus")),
                            labels = rev(c("PNS", "CNS", "26-49", "50-69", "70-79", "80+"))
         ),
         Outcome = "Mortality"
  )

meta_ind_discharge_age <- meta_ind_discharge %>% 
  filter(grepl(x = Covariate, pattern = "(age|neuro)")) %>% 
  mutate(
         Covariate = as.factor(Covariate),
         Covariate = factor(Covariate,
                            levels = rev(c("neuro_postPeripheral",
                                       "neuro_postCentral",
                                       "age_group26to49",
                                       "age_group50to69",
                                       "age_group70to79",
                                       "age_group80plus")),
                            labels = rev(c("PNS", "CNS", "26-49", "50-69", "70-79", "80+"))
         ),
         Outcome = "Discharge"
  )

meta_ind_age <- rbind(meta_ind_deceased_age, meta_ind_discharge_age)


color_palette <- c(PNS = "slateblue", CNS ="tomato", `26-49` = "#fde725",
                   `50-69` = "#35b779", `70-79` = "#31688e", `80+` = "#440154")

age_hr <- ggplot(meta_ind_age, aes(x=TE.random, y=Covariate, colour=Covariate, group = Outcome)) +
  scale_color_manual("", values = color_palette) +
  geom_point(aes(x=TE.random), shape=15, size=3, position = pd) +
  geom_linerange(aes(xmin=lower.random, xmax=upper.random), position = pd) +
  labs(x = "Hazard Ratio", y = "") + 
  geom_vline(xintercept = 1, linetype="dashed")  + 
  facet_grid(~Outcome, scales = "free_x") + 
  theme_bw() +
  theme(strip.text.x = element_text(size = 15, face = "bold"),
        strip.text.y = element_text(size = 15, face = "bold"),
        axis.title.x = element_text(size = 15, face = "bold"),
        axis.text.x = element_text(size = 13, face = "bold", color = "black"),
        axis.text.y = element_text(size = 12, face = "bold"),
        legend.position = "none") 

ggsave("figures/Hazard_ratio_age.png", age_hr, width = 12, height = 6, dpi = 300)

Evaluate all covariates

Create a scale function to identify the min and max CI for each site, excluding sites with values of ‘Inf’ or our invalid sites

scale_plot <- function(combined_results_df) {
  scale <- combined_results_df %>%
    filter(
      !lower.random == "Inf",
      !upper.random == "Inf"
    ) %>%
    mutate(
      min_est = min(lower.random, na.rm = TRUE) - 1,
      max_est = max(upper.random, na.rm = TRUE) + 1,
      min_est = if_else(min_est + 1 >= -1, min(lower.random, na.rm = TRUE)-0.1, min_est),
      min_est = if_else(min_est > 1, 0.5, min_est),
      max_est = if_else(max_est - 1 < 1, 1.1, max_est)) %>%
    distinct(min_est, max_est)

  return(scale)
}

Mortality Results

## define forest plot shapes
# shape #16 is normal circle; #18 is diamond
shapes <- rep(16, times = nrow(meta_ind_deceased))
#shapes[length(shapes)] <- 18
sizes <- rep(3.25, times = nrow(meta_ind_deceased))
#sizes[length(sizes)] <- 5
scale <- scale_plot(meta_ind_deceased)

forester(
  left_side_data = meta_ind_deceased %>% select(Covariate),
  estimate = meta_ind_deceased$Estimate,
  ci_low = meta_ind_deceased$lower.random,
  ci_high = meta_ind_deceased$upper.random,
  right_side_data = meta_ind_deceased[, c("Hazard Ratio", "P-value")],
  display = TRUE,
  font_family = "arial",
  arrows = TRUE,
  arrow_labels = c("Decreased Risk", "Increased Risk"),
  null_line_at = 1,
  xlim = c(scale$min_est,25),
  point_sizes = sizes,
  point_shapes = shapes, 
  x_scale_linear = FALSE,
  add_plot_width = 5,
  ggplot_width = 100,
  set_png_height = 12.5,
  file_path = here::here(paste0("figures/meta-analysis-coefficients-ind-90-mortality.png")))

Discharge Results

## define forest plot shapes
# shape #16 is normal circle; #18 is diamond
shapes <- rep(16, times = nrow(meta_ind_discharge))
#shapes[length(shapes)] <- 18
sizes <- rep(3.25, times = nrow(meta_ind_discharge))
#sizes[length(sizes)] <- 5
scale <- scale_plot(meta_ind_discharge)

forester(
  left_side_data = meta_ind_discharge %>% select(Covariate),
  estimate = meta_ind_discharge$Estimate,
  ci_low = meta_ind_discharge$lower.random,
  ci_high = meta_ind_discharge$upper.random,
  right_side_data = meta_ind_discharge[, c("Hazard Ratio", "P-value")],
  display = TRUE,
  font_family = "arial",
  arrows = TRUE,
  arrow_labels = c("Decreased Risk", "Increased Risk"),
  null_line_at = 1,
  xlim = c(scale$min_est, 3.5),
  point_sizes = sizes,
  point_shapes = shapes, 
  x_scale_linear = FALSE,
  add_plot_width = 5,
  ggplot_width = 100,
  set_png_height = 12.5,
  file_path = here::here(paste0("figures/meta-analysis-coefficients-ind-90-discharge.png")))