---
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()

```



