--- title: Sending center or not? plain-language-summary: | This be the plain language thing date: last-modified bibliography: references.bib csl: critical-care-medicine.csl number-sections: false --- ```{r setup} # Load dependencies #| output: false #| echo: false #| warning: false library(tidyverse) library(lubridate) library(brms) library(DBI) library(RSQLite) library(flowchart) library(gtsummary) library(ggplot2) library(bayestestR) library(marginaleffects) library(rstan) library(flowchart) library("devtools") library(rcmswe) library(NatParksPalettes) library(showtext) font_add_google("Noto Sans", "Noto Sans") # Downloads from Google Fonts showtext_auto() # Set gtsummart theme theme_gtsummary_compact() # Load data d <- read_delim(file = "~/PhD/nsicu-transfers/data/pre-processed-data/patient-df-2025-05-25 22:12:51.563534.csv", delim = ";")[ , -1] # Prune data d_pruned <- d %>% filter(sir_dsc_time >= "2019-10-01" & sir_dsc_time <= "2024-05-15") %>% # 2019-10-01 -- 2024-05-15 filter(!is.na(TERTIARY_HADM_ID)) %>% # Keep only rows with a tertiary HADM match filter(sir_total_time <= 720) %>% # Keep only rows with a primary ICU stay 12 hours or shorter filter(SIR_PAR_OFFSET_TIGHT %in% (c(-1,0,1))) %>% # Keep only rows with a recent PAR admit filter(road_distance >= 45) %>% # Keep only rows with transfers 45 km or longer filter(DX_GROUP %in% c("ASAH", "ICH", "AIS", "TBI")) %>% #Keep only admits with ASAH, ICH, AIS or TBI filter(min_rank(DX_RANK) == 1, .by = VtfId_LopNr) %>% # Within each unique VtfId_LopNr, keep only the lowest ranking diagnostic PAR admit slice_max(sir_dsc_time, by = LopNr, n = 1, with_ties = FALSE) # This gives the latest date by default (if ties, may give more than one) heavy_users <- d_pruned %>% group_by(formatted_icu_name) %>% summarise( proportion_hems_ift = mean(hems_ift == TRUE, na.rm = TRUE), n = n() ) %>% filter(proportion_hems_ift >= 0.1) %>% select(formatted_icu_name) %>% unique() light_users <- d_pruned %>% group_by(formatted_icu_name) %>% summarise( proportion_hems_ift = mean(hems_ift == TRUE, na.rm = TRUE), n = n() ) %>% filter(proportion_hems_ift < 0.05) %>% select(formatted_icu_name) %>% unique() ``` ```{r} # Load dependencies #| output: false #| echo: false #| warning: false library(tidyverse) library(lubridate) library(brms) library(DBI) library(RSQLite) library(flowchart) library(gtsummary) library(ggplot2) library(bayestestR) library(marginaleffects) library(rstan) library(flowchart) library("devtools") library(rcmswe) library(NatParksPalettes) library(showtext) font_add_google("Noto Sans", "Noto Sans") # Downloads from Google Fonts showtext_auto() # Set gtsummart theme theme_gtsummary_compact() # Load data d <- read_delim(file = "~/PhD/nsicu-transfers/data/pre-processed-data/patient-df-2025-05-25 22:12:51.563534.csv", delim = ";")[ , -1] # Prune data d_pruned <- d %>% filter(sir_dsc_time >= "2019-10-01" & sir_dsc_time <= "2024-05-15") %>% # 2019-10-01 -- 2024-05-15 filter(!is.na(TERTIARY_HADM_ID)) %>% # Keep only rows with a tertiary HADM match filter(sir_total_time <= 720) %>% # Keep only rows with a primary ICU stay 12 hours or shorter filter(SIR_PAR_OFFSET_TIGHT %in% (c(-1,0,1))) %>% # Keep only rows with a recent PAR admit filter(road_distance >= 45) %>% # Keep only rows with transfers 45 km or longer filter(DX_GROUP %in% c("ASAH", "ICH", "AIS", "TBI")) %>% #Keep only admits with ASAH, ICH, AIS or TBI filter(min_rank(DX_RANK) == 1, .by = VtfId_LopNr) %>% # Within each unique VtfId_LopNr, keep only the lowest ranking diagnostic PAR admit slice_max(sir_dsc_time, by = LopNr, n = 1, with_ties = FALSE) # This gives the latest date by default (if ties, may give more than one) heavy_users <- d_pruned %>% group_by(formatted_icu_name) %>% summarise( proportion_hems_ift = mean(hems_ift == TRUE, na.rm = TRUE), n = n() ) %>% filter(proportion_hems_ift >= 0.1) %>% select(formatted_icu_name) %>% unique() light_users <- d_pruned %>% group_by(formatted_icu_name) %>% summarise( proportion_hems_ift = mean(hems_ift == TRUE, na.rm = TRUE), n = n() ) %>% filter(proportion_hems_ift < 0.05) %>% select(formatted_icu_name) %>% unique() icu_lookup <- tribble( ~formatted_icu_name, ~icu_category, "Hudiksvall IVA", "II", "Östersund IVA", "II", "Skövde IVA", "II", "Eskilstuna IVA", "II", "Falun IVA", "II", "Kalix IVA", "I", "Lycksele IVA", "I", "Karlskoga IVA", "II", "Lidköping IVA", "II", "Skellefteå IVA", "II", "Sunderby IVA", "II", "Växjö IVA", "II", "Trollhättan IVA", "II", "Västervik IVA", "II", "Karlskrona IVA", "II", "Alingsås IVA", "I", "Jönköping IVA", "II", "Mora IVA", "I", "Gävle IVA", "II", "Kristianstad IVA", "II", "Helsingborg IVA", "II", "Kalmar IVA", "II", "Nyköping IVA", "II", "Sundsvall IVA", "II", "Gällivare IVA", "I", "Sollefteå IVA", "I", "Västerås IVA", "II", "Borås IVA", "II", "Visby IVA", "II", "Lindesberg IVA", "I", "Karlstad IVA", "II", "Bollnäs IVA", "I", "Norrtälje IVA", "I", "Värnamo IVA", "II", "Örebro IVA", "II", "Halmstad IVA", "II", "Örnsköldsvik IVA", "II", "Piteå IVA", "I", "Norrköping IVA", "II", "Ystad IVA", "II", "Varberg IVA", "II", "Ljungby IVA", "I", "Torsby IVA", "I", "Arvika IVA", "I", "Eksjö IVA", "I", "Södertälje IVA", "II" ) d_pruned <- d_pruned %>% left_join(icu_lookup, by = "formatted_icu_name") ``` ```{r glm-causal} #| output: false #| echo: false #| warning: false m.causal <- d_pruned %>% mutate(z_SAPS_score_minus_age_and_gcs = as.numeric(scale(SAPS_score_minus_age_and_gcs))) %>% mutate( DX_GROUP = forcats::fct_relabel(DX_GROUP, ~ gsub("_", "-", .)) ) %>% mutate(hems_ift = as.integer(hems_ift)) %>% mutate(z_age = as.numeric(scale(age))) %>% brm( formula = MORTALITY_30D ~ hems_ift + (1 + hems_ift | formatted_icu_name) + DX_GROUP + DNR + SAPS_obtunded + z_age + z_SAPS_score_minus_age_and_gcs, data = ., family = bernoulli(), cores = 4) m.causal.exp_vals <- summary(m.causal)$fixed %>% slice(2) %>% select("Estimate", "l-95% CI", "u-95% CI") %>% unlist() %>% exp() ``` ```{r glm-causal-no-cov} #| output: false #| echo: false #| warning: false m.causal.nocov <- d_pruned %>% mutate(z_SAPS_score_minus_age_and_gcs = as.numeric(scale(SAPS_score_minus_age_and_gcs))) %>% mutate( DX_GROUP = forcats::fct_relabel(DX_GROUP, ~ gsub("_", "-", .)) ) %>% mutate(hems_ift = as.integer(hems_ift)) %>% mutate(z_age = as.numeric(scale(age))) %>% brm( formula = MORTALITY_30D ~ hems_ift + (1 + hems_ift || formatted_icu_name) + DX_GROUP + DNR + SAPS_obtunded + z_age + z_SAPS_score_minus_age_and_gcs, data = ., family = bernoulli(), cores = 4 ) m.causal.nocov.exp_vals <- summary(m.causal.nocov)$fixed %>% slice(2) %>% select("Estimate", "l-95% CI", "u-95% CI") %>% unlist() %>% exp() ``` ```{r create-tbl-glm-outcome} #| output: false #| echo: false #| warning: false tbl.glm.causal <- m.causal %>% tbl_regression(exponentiate = TRUE, conf.level = .95, label = c( DX_GROUP ~ "Inferred diagnosis", z_age ~ "Age, std", SAPS_obtunded ~ "SAPS GCS <15", DNR ~ "DNR order", hems_ift ~ "Interfacility transfer by HEMS", z_SAPS_score_minus_age_and_gcs ~ "SAPS score (excl. age/GCS components), std") ) %>% modify_header(estimate = "**aOR**") %>% modify_header(label = "**Predictor**") tbl.glm.causal.nocov <- m.causal.nocov %>% tbl_regression(exponentiate = TRUE, conf.level = .95, label = c( DX_GROUP ~ "Inferred diagnosis", z_age ~ "Age, std", SAPS_obtunded ~ "SAPS GCS <15", DNR ~ "DNR order", hems_ift ~ "Interfacility transfer by HEMS", z_SAPS_score_minus_age_and_gcs ~ "SAPS score (excl. age/GCS components), std") ) %>% modify_header(estimate = "**aOR**") %>% modify_header(label = "**Predictor**") glm_table <- tbl_merge(tbls = list(tbl.glm.causal, tbl.glm.causal.nocov), tab_spanner = c("**Model I**", "**Model I-nocov**")) %>% modify_caption(caption = "**Multivariate Bayesian mixed effect generalized linear models predicting 30-day mortality (I, II)**") ``` ```{r} glm_table ``` ```{r m-causal-contrasts} # calculate marginal effects #| output: false #| echo: false #| warning: false cmp_m.causal <- comparisons(m.causal, variables = "hems_ift", re_formula=NULL) ci_m.causal <- cmp_m.causal %>% tidy() %>% select(estimate) %>% ci(x=., method="HDI") est_contrast <- cmp_m.causal %>% tidy() %>% select(estimate) cmp_m.causal.nocov <- comparisons(m.causal.nocov, variables = "hems_ift", re_formula=NULL) ci_m.causal.nocov <- cmp_m.causal.nocov %>% tidy() %>% select(estimate) %>% ci(x=., method="HDI") est_contrast.nocov <- cmp_m.causal.nocov %>% tidy() %>% select(estimate) ``` ```{r} mean(est_contrast$estimate > 0) ``` ```{r} mean(est_contrast.nocov$estimate > 0) ``` ### Adding receiving center ```{r} d_pruned <- d_pruned %>% mutate( hospital_name_receiving = case_when( hospital_name_receiving == "Akademiska sjukhuset" ~ "Uppsala", hospital_name_receiving == "Karolinska universitetssjukhuset, Solna" ~ "KS", hospital_name_receiving == "Norrlands universitetssjukhus" ~ "Umeå", hospital_name_receiving == "Sahlgrenska universitetssjukhuset" ~ "SU", hospital_name_receiving == "Universitetssjukhuset i Linköping" ~ "Linköping", hospital_name_receiving == "Universitetssjukhuset i Lund" ~ "Lund", hospital_name_receiving == "Universitetssjukhuset i Örebro" ~ "Örebro", TRUE ~ hospital_name_receiving ), # collapse into two meta-groups hospital_group = case_when( hospital_name_receiving %in% c("Örebro", "Lund", "Linköping") ~ "Group A (Örebro/Lund/Linköping)", hospital_name_receiving %in% c("Umeå", "Uppsala", "SU", "KS") ~ "Group B (Umeå/Uppsala/SU/KS)" ), hospital_group = factor(hospital_group, levels = c("Group A (Örebro/Lund/Linköping)", "Group B (Umeå/Uppsala/SU/KS)")) ) ``` ```{r glm-causal-receiving} #| output: false #| echo: false #| warning: false m.causal.receiving <- d_pruned %>% mutate(z_SAPS_score_minus_age_and_gcs = as.numeric(scale(SAPS_score_minus_age_and_gcs))) %>% mutate(hems_ift = as.integer(hems_ift)) %>% mutate( DX_GROUP = forcats::fct_relabel(DX_GROUP, ~ gsub("_", "-", .)) ) %>% mutate(z_age = as.numeric(scale(age))) %>% brm( formula = MORTALITY_30D ~ 1 + hems_ift + hospital_name_receiving + DX_GROUP + DNR + sir_consciousness_level + z_age + z_SAPS_score_minus_age_and_gcs, data = ., family = bernoulli(), cores = 4 ) ``` ```{r} d_pruned %>% filter(hospital_group == "Group B (Umeå/Uppsala/SU/KS)") ``` ```{r glm-causal-receiving-at-group-b} #| output: false #| echo: false #| warning: false m.causal.receiving.b <- d_pruned %>% #filter(hospital_group == "Group B (Umeå/Uppsala/SU/KS)") %>% mutate(z_SAPS_score_minus_age_and_gcs = as.numeric(scale(SAPS_score_minus_age_and_gcs))) %>% mutate(hems_ift = as.integer(hems_ift)) %>% mutate( DX_GROUP = forcats::fct_relabel(DX_GROUP, ~ gsub("_", "-", .)) ) %>% mutate(z_age = as.numeric(scale(age))) %>% brm( formula = MORTALITY_30D ~ 1 + hems_ift + DX_GROUP + DNR + sir_consciousness_level + z_age + z_SAPS_score_minus_age_and_gcs + hospital_name_receiving, data = ., family = bernoulli(), cores = 4 ) ``` ```{r} library(insight) ate_cf <- function(model, ndraws = 4000) { dat <- insight::get_data(model) dat0 <- transform(dat, hems_ift = 0L) dat1 <- transform(dat, hems_ift = 1L) ep0 <- posterior_epred(model, newdata = dat0, ndraws = ndraws, re_formula = NULL) ep1 <- posterior_epred(model, newdata = dat1, ndraws = ndraws, re_formula = NULL) ate <- rowMeans(ep1 - ep0) mean0 <- rowMeans(ep0) mean1 <- rowMeans(ep1) summ <- data.frame( mean = mean(ate), median = median(ate), sd = sd(ate), q025 = quantile(ate, 0.025), q975 = quantile(ate, 0.975), pr_gt0 = mean(ate > 0) ) list(summary = summ, draws = ate, mean0 = mean0, mean1 = mean1) } ``` ```{r} res <- ate_cf(m.causal.receiving.b) res$summary ``` ```{r} d_pruned %>% mutate( hospital_name_receiving = case_when( hospital_name_receiving == "Akademiska sjukhuset" ~ "Uppsala", hospital_name_receiving == "Karolinska universitetssjukhuset, Solna" ~ "KS", hospital_name_receiving == "Norrlands universitetssjukhus" ~ "Umeå", hospital_name_receiving == "Sahlgrenska universitetssjukhuset" ~ "SU", hospital_name_receiving == "Universitetssjukhuset i Linköping" ~ "Linköping", hospital_name_receiving == "Universitetssjukhuset i Lund" ~ "Lund", hospital_name_receiving == "Universitetssjukhuset i Örebro" ~ "Örebro", TRUE ~ hospital_name_receiving ), hospital_name_receiving = factor( hospital_name_receiving, levels = c("Örebro", "Lund", "Linköping", "Umeå", "SU", "Uppsala", "KS") ) ) %>% select( hospital_name_receiving,hems_ift, MORTALITY_30D, age, SAPS_obtunded, SAPS_unconcious, SAPS_total_score ) %>% tbl_summary( by = hospital_name_receiving, missing = "no", label = list( hems_ift ~ "Interfacility transfer by HEMS", MORTALITY_30D ~ "30-day mortality", age ~ "Age (years)", SAPS_obtunded ~ "SAPS GCS <15", SAPS_unconcious ~ "SAPS GCS <9", SAPS_total_score ~ "SAPS Score" ), missing_text = "Missing" ) %>% modify_caption(caption = "**Patient characteristics**") ``` ```{r} d_pruned %>% select( hospital_group, hems_ift, MORTALITY_30D, DX_GROUP, age, SAPS_obtunded, SAPS_unconcious, SAPS_total_score ) %>% tbl_summary( by = hospital_group, missing = "no", label = list( hems_ift ~ "Interfacility transfer by HEMS", MORTALITY_30D ~ "30-day mortality", age ~ "Age (years)", DX_GROUP ~ "Diagnosis", SAPS_obtunded ~ "SAPS GCS <15", SAPS_unconcious ~ "SAPS GCS <9", SAPS_total_score ~ "SAPS Score" ), missing_text = "Missing" ) %>% modify_caption(caption = "**Patient characteristics by receiving hospital group**") %>% add_p() ``` ```{r} d_pruned %>% mutate( hospital_name_receiving = case_when( hospital_name_receiving == "Akademiska sjukhuset" ~ "Uppsala", hospital_name_receiving == "Karolinska universitetssjukhuset, Solna" ~ "KS", hospital_name_receiving == "Norrlands universitetssjukhus" ~ "Umeå", hospital_name_receiving == "Sahlgrenska universitetssjukhuset" ~ "SU", hospital_name_receiving == "Universitetssjukhuset i Linköping" ~ "Linköping", hospital_name_receiving == "Universitetssjukhuset i Lund" ~ "Lund", hospital_name_receiving == "Universitetssjukhuset i Örebro" ~ "Örebro", TRUE ~ hospital_name_receiving ), # collapse into two meta-groups hospital_group = case_when( hospital_name_receiving %in% c("Örebro", "Lund", "Linköping") ~ "Group A (Örebro/Lund/Linköping)", hospital_name_receiving %in% c("Umeå", "Uppsala", "SU", "KS") ~ "Group B (Umeå/Uppsala/SU/KS)" ), hospital_group = factor(hospital_group, levels = c("Group A (Örebro/Lund/Linköping)", "Group B (Umeå/Uppsala/SU/KS)")) ) %>% mutate(hems_ift = case_when( hems_ift == TRUE ~ "HEMS transfer", hems_ift == FALSE ~ "Non-HEMS transfer" )) %>% filter(hospital_group == "Group B (Umeå/Uppsala/SU/KS)") %>% select(hems_ift, age, DX_GROUP, sir_consciousness_level, SAPS_obtunded, any_AMV, SAPS_total_score, SAPS_score_minus_age_and_gcs, MORTALITY_30D) %>% tbl_summary(by=hems_ift, missing="no", label = list(hems_ift ~ "Interfacility transfer by HEMS", age ~ "Age (years)", DX_GROUP ~ "Inferred diagnosis", sir_consciousness_level ~ "SAPS consciousness level", SAPS_obtunded ~ "SAPS GCS <15", any_AMV ~ "Pre-transfer Mechanical Ventilation", SAPS_total_score ~ "SAPS Score", SAPS_score_minus_age_and_gcs ~ "SAPS Score excl. age and GCS components", MORTALITY_30D ~ "30-day mortality" ), missing_text="Missing" ) %>% modify_caption(caption = "**Patient characteristics (Group B)**") %>% add_p() ``` ```{r} d_pruned %>% mutate( icu_category = case_when( icu_category == "I" ~ "I (Små sjukhus)", icu_category == "II" ~ "II (Mellanstora sjukhus") )%>% select( icu_category, hems_ift, MORTALITY_30D, DX_GROUP, age, SAPS_obtunded, SAPS_unconcious, SAPS_total_score ) %>% tbl_summary( by = icu_category, missing = "no", label = list( hems_ift ~ "Interfacility transfer by HEMS", MORTALITY_30D ~ "30-day mortality", age ~ "Age (years)", DX_GROUP ~ "Diagnosis", SAPS_obtunded ~ "SAPS GCS <15", SAPS_unconcious ~ "SAPS GCS <9", SAPS_total_score ~ "SAPS Score" ), missing_text = "Missing" ) %>% modify_caption(caption = "**Patient characteristics by sending hospital type**") %>% add_p() ```