Back to Article
Sending center or not?
Download Source

Sending center or not?

Published

September 27, 2025

In [1]:
# Load dependencies 
library(tidyverse)
Warning: package 'tidyr' was built under R version 4.2.3
Warning: package 'readr' was built under R version 4.2.3
Warning: package 'dplyr' was built under R version 4.2.3
Warning: package 'stringr' was built under R version 4.2.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(brms)
Warning: package 'brms' was built under R version 4.2.3
Loading required package: Rcpp
Loading 'brms' package (version 2.21.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: 'brms'

The following object is masked from 'package:stats':

    ar
library(DBI)
library(RSQLite)
library(flowchart)
library(gtsummary)
library(ggplot2)
library(bayestestR)
Warning: package 'bayestestR' was built under R version 4.2.3
library(marginaleffects)
Warning: package 'marginaleffects' was built under R version 4.2.3
library(rstan)
Warning: package 'rstan' was built under R version 4.2.3
Loading required package: StanHeaders

rstan version 2.32.6 (Stan version 2.32.2)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
change `threads_per_chain` option:
rstan_options(threads_per_chain = 1)


Attaching package: 'rstan'

The following object is masked from 'package:tidyr':

    extract
library(flowchart)
library("devtools")
Loading required package: usethis
library(rcmswe)
library(NatParksPalettes)
library(showtext)
Warning: package 'showtext' was built under R version 4.2.3
Loading required package: sysfonts
Warning: package 'sysfonts' was built under R version 4.2.3
Loading required package: showtextdb
font_add_google("Noto Sans", "Noto Sans")  # Downloads from Google Fonts
showtext_auto()

# Set gtsummart theme
theme_gtsummary_compact()
Setting theme `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]
New names:
Rows: 109107 Columns: 127
── Column specification
──────────────────────────────────────────────────────── Delimiter: ";" chr
(18): DX_GROUP, hospital_name_receiving, sir_icu_name, sir_hospital_typ... dbl
(91): ...1, VtfId_LopNr, TERTIARY_HADM_ID, LopNr, SIR_PAR_OFFSET, SIR_P... lgl
(8): daylight, hems_minima_sending, hems_minima_window, daylight_recei... dttm
(9): sir_adm_time, sir_dsc_time, sir_adm_time_UTC, sir_dsc_time_UTC, w... date
(1): DODSDAT_ROUND_UP
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...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()
In [2]:
# Load dependencies 
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()
Setting theme `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]
New names:
Rows: 109107 Columns: 127
── Column specification
──────────────────────────────────────────────────────── Delimiter: ";" chr
(18): DX_GROUP, hospital_name_receiving, sir_icu_name, sir_hospital_typ... dbl
(91): ...1, VtfId_LopNr, TERTIARY_HADM_ID, LopNr, SIR_PAR_OFFSET, SIR_P... lgl
(8): daylight, hems_minima_sending, hems_minima_window, daylight_recei... dttm
(9): sir_adm_time, sir_dsc_time, sir_adm_time_UTC, sir_dsc_time_UTC, w... date
(1): DODSDAT_ROUND_UP
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...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")
In [3]:

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)
Warning: Rows containing NAs were excluded from the model.
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Start sampling

m.causal.exp_vals <- summary(m.causal)$fixed %>%
  slice(2) %>%
  select("Estimate", "l-95% CI", "u-95% CI") %>%
  unlist() %>%
  exp()
In [4]:

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
  )
Warning: Rows containing NAs were excluded from the model.
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Start sampling

m.causal.nocov.exp_vals <- summary(m.causal.nocov)$fixed %>%
  slice(2) %>%
  select("Estimate", "l-95% CI", "u-95% CI") %>%
  unlist() %>%
  exp()
In [5]:
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**")
Warning in tidy.brmsfit(x, ..., effects = "fixed"): some parameter names
contain underscores: term naming may be unreliable!
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**")
Warning in tidy.brmsfit(x, ..., effects = "fixed"): some parameter names
contain underscores: term naming may be unreliable!
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)**")
In [6]:
glm_table
Multivariate Bayesian mixed effect generalized linear models predicting 30-day mortality (I, II)
Predictor Model I Model I-nocov
aOR 95% CI1 aOR 95% CI1
Interfacility transfer by HEMS 0.80 0.51, 1.23 0.83 0.55, 1.24
Inferred diagnosis



    AIS
    ASAH 0.43 0.28, 0.66 0.43 0.28, 0.65
    ICH 0.30 0.19, 0.48 0.30 0.19, 0.48
    TBI 0.32 0.21, 0.47 0.32 0.21, 0.48
DNR order 3.70 1.62, 8.63 3.68 1.63, 8.33
SAPS GCS <15 7.39 4.25, 13.8 7.35 4.26, 13.8
Age, std 1.69 1.44, 2.00 1.69 1.43, 2.01
SAPS score (excl. age/GCS components), std 1.67 1.44, 1.93 1.67 1.44, 1.93
1 CI = Credible Interval
In [7]:
# calculate marginal effects
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)
In [8]:
mean(est_contrast$estimate > 0)
[1] 0.09441341
In [9]:
 mean(est_contrast.nocov$estimate > 0)
[1] 0.05307263

Adding receiving center

In [10]:
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)"))
  ) 
In [11]:
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
  )
Warning: Rows containing NAs were excluded from the model.
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Start sampling
In [12]:
d_pruned %>% filter(hospital_group == "Group B (Umeå/Uppsala/SU/KS)")
# A tibble: 1,311 × 128
   VtfId_LopNr TERTIARY_HADM_ID LopNr SIR_PAR_OFFSET SIR_PAR_OFFSET_TIGHT
         <dbl>            <dbl> <dbl>          <dbl>                <dbl>
 1      166403            10585  7030              0                    0
 2      166458            47230 32121              0                    0
 3      166966            90697 61656              0                    0
 4      166962            39386 27001              0                    0
 5      166675            14277  9502              0                    0
 6      167961            29577 20189              0                    0
 7      168033            13167  8654              0                    0
 8      166959            45235 30619              0                    0
 9      166610            48412 33040             -1                   -1
10      166672            57120 38752              0                    0
# ℹ 1,301 more rows
# ℹ 123 more variables: OREBRO_INTERNAL <dbl>, DX_GROUP <chr>, DX_RANK <dbl>,
#   HADM_ID <dbl>, age <dbl>, hospital_name_receiving <chr>, sex_female <dbl>,
#   sir_icu_name <chr>, sir_hospital_type <chr>, sir_adm_time <dttm>,
#   sir_dsc_time <dttm>, sir_total_time <dbl>, admission_height <dbl>,
#   admission_weight <dbl>, BMI <dbl>, DNR <dbl>, icu_discharge_daytime <dbl>,
#   icu_discharge_nighttime <dbl>, icu_discharge_afterhours <dbl>, …
In [13]:
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
  )
Warning: Rows containing NAs were excluded from the model.
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1
Start sampling
In [14]:
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)
}
In [15]:
res <- ate_cf(m.causal.receiving.b)
res$summary
              mean        median         sd        q025       q975  pr_gt0
2.5% -0.0006701878 -0.0007074971 0.02054095 -0.04012373 0.04117289 0.48625
In [16]:
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**")
Patient characteristics
Characteristic Örebro, N = 511 Lund, N = 1841 Linköping, N = 2661 Umeå, N = 4361 SU, N = 1961 Uppsala, N = 5531 KS, N = 1261
Interfacility transfer by HEMS 0 (0%) 1 (0.5%) 3 (1.1%) 89 (20%) 60 (31%) 221 (40%) 57 (45%)
30-day mortality 16 (31%) 40 (22%) 47 (18%) 100 (23%) 34 (17%) 55 (9.9%) 18 (14%)
Age (years) 65 (53, 75) 60 (46, 70) 64 (52, 74) 65 (54, 74) 63 (52, 72) 61 (50, 71) 69 (59, 76)
SAPS GCS <15 42 (82%) 154 (87%) 203 (77%) 339 (79%) 150 (77%) 382 (69%) 64 (51%)
SAPS GCS <9 19 (37%) 91 (51%) 91 (35%) 193 (45%) 93 (48%) 160 (29%) 29 (23%)
SAPS Score 54 (45, 68) 54 (45, 67) 55 (46, 62) 57 (48, 67) 55 (47, 63) 51 (44, 58) 52 (45, 57)
1 n (%); Median (IQR)
In [17]:
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()
Patient characteristics by receiving hospital group
Characteristic Group A (Örebro/Lund/Linköping), N = 5011 Group B (Umeå/Uppsala/SU/KS), N = 1,3111 p-value2
Interfacility transfer by HEMS 4 (0.8%) 427 (33%) <0.001
30-day mortality 103 (21%) 207 (16%) 0.016
Diagnosis

0.10
    AIS 56 (11%) 199 (15%)
    ASAH 131 (26%) 362 (28%)
    ICH 107 (21%) 248 (19%)
    TBI 207 (41%) 502 (38%)
Age (years) 63 (50, 73) 63 (52, 73) 0.3
SAPS GCS <15 399 (81%) 935 (72%) <0.001
SAPS GCS <9 201 (41%) 475 (37%) 0.089
SAPS Score 54 (45, 64) 53 (45, 62) 0.2
1 n (%); Median (IQR)
2 Pearson’s Chi-squared test; Wilcoxon rank sum test
In [18]:
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()
Patient characteristics (Group B)
Characteristic HEMS transfer, N = 4271 Non-HEMS transfer, N = 8841 p-value2
Age (years) 61 (48, 72) 64 (53, 73) 0.009
Inferred diagnosis

0.5
    AIS 70 (16%) 129 (15%)
    ASAH 108 (25%) 254 (29%)
    ICH 87 (20%) 161 (18%)
    TBI 162 (38%) 340 (38%)
SAPS consciousness level

0.4
    I (GCS ≥13) 203 (48%) 429 (49%)
    II (GCS 7-12) 125 (29%) 229 (26%)
    III (GCS 6) 26 (6.1%) 77 (8.8%)
    IV (GCS 5) 22 (5.2%) 46 (5.3%)
    V (GCS ≤4) 48 (11%) 94 (11%)
SAPS GCS <15 302 (71%) 633 (72%) 0.7
Pre-transfer Mechanical Ventilation 234 (55%) 486 (55%) >0.9
SAPS Score 53 (46, 60) 53 (45, 63) 0.6
SAPS Score excl. age and GCS components 41 (37, 46) 41 (36, 48) >0.9
30-day mortality 60 (14%) 147 (17%) 0.2
1 Median (IQR); n (%)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test
In [19]:
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()
Patient characteristics by sending hospital type
Characteristic I (Små sjukhus), N = 3331 II (Mellanstora sjukhus, N = 1,4791 p-value2
Interfacility transfer by HEMS 83 (25%) 348 (24%) 0.6
30-day mortality 61 (18%) 249 (17%) 0.5
Diagnosis

<0.001
    AIS 82 (25%) 173 (12%)
    ASAH 69 (21%) 424 (29%)
    ICH 62 (19%) 293 (20%)
    TBI 120 (36%) 589 (40%)
Age (years) 65 (54, 73) 63 (51, 73) 0.066
SAPS GCS <15 238 (73%) 1,096 (75%) 0.5
SAPS GCS <9 108 (33%) 568 (39%) 0.056
SAPS Score 52 (45, 62) 54 (46, 62) 0.3
1 n (%); Median (IQR)
2 Pearson’s Chi-squared test; Wilcoxon rank sum test