# ==============================================================================
# TITLE: R Code for ENIGMA Project Analysis
# AUTHOR: Giacomo Biganzoli & Co-authors
# DATE: 2026
#
# DATA AVAILABILITY NOTE:
# In compliance with ethical and GDPR requirements, the original Excel file
# containing sensitive patient data is not included in this repository.
# This script automatically generates a simulated/mock dataset (N=1500) inline
# to ensure full reproducibility of the code and models without external files.
# ==============================================================================
# ==============================================================================
# 0. ENVIRONMENT SETUP AND LIBRARIES
# ==============================================================================
# Grouping of all necessary dependencies for the analysis
library(dplyr)
library(tidyr)
library(ggplot2)
library(survival)
library(mstate)
library(mice) # Multiple imputation
library(readxl)
library(gtsummary)
library(ggsurvfit)
library(ggsci) # Color palette for ggplot
library(nnet) # Multinomial models
library(splines)
library(tableone)
library(mitools)
library(survey)
library(gt)
library(tibble)
library(gdata) # For the keep() function
# Automatic creation of the output folder if it doesn't exist
if (!dir.exists("output")) dir.create("output")
if (!dir.exists("data")) dir.create("data")
# ==============================================================================
# 1. INLINE MOCK DATA GENERATION (Replaces external Excel loading)
# ==============================================================================
library(writexl)
# Fix seed for total reproducibility
set.seed(2026)
# Generate 1500 initial patients (matches original ~1587 cases before exclusions)
n <- 1500
# 1. IDENTIFIERS & CENTERS
# Mixing the considered centers with a few fake ones to test the exclusion filter
centers_considered <- c(1911, 1913, 1912, 1927, 1934, 1932)
centers_all <- c(centers_considered, 9999, 8888)
center_id <- sample(centers_all, n, replace = TRUE, prob = c(rep(0.15, 6), 0.05, 0.05))
patient_id <- sprintf("%04d", 1:n)
record_id <- paste(center_id, patient_id, sep = "-")
# 2. BASELINE DEMOGRAPHICS & CLINICALS
index_eta <- round(runif(n, 40, 85))
demo_sesso <- sample(c("M", "F"), n, replace = TRUE, prob = c(0.7, 0.3))
index_eziologia <- sample(c("AMI", "HF-CS", "Other", "TEP"), n, replace = TRUE, prob = c(0.40, 0.40, 0.15, 0.05))
Oh_score_scai <- sample(c("B", "C", "D", "E", "A", NA), n, replace = TRUE, prob = c(0.2, 0.35, 0.25, 0.1, 0.05, 0.05))
Oh_lab_gfr <- round(rnorm(n, mean = 55, sd = 20), 1)
Oh_par_vit_ega_latt <- round(runif(n, 1, 12), 1)
demo_diabete <- sample(c("Yes", "No"), n, replace = TRUE, prob = c(0.3, 0.7))
# Enigma baseline generic drugs (must be character)
enigma_arb <- sample(c("Yes", "No", NA), n, replace = TRUE)
enigma_mra <- sample(c("Yes", "No", NA), n, replace = TRUE)
enigma_beta_block <- sample(c("Yes", "No", NA), n, replace = TRUE)
enigma_sglt2 <- sample(c("Yes", "No", NA), n, replace = TRUE)
# 3. TIMELINES (Admission and Devices)
# Generating admission dates strictly before 2025-04-01 to pass the filter
index_data_ricovero <- as.Date("2023-01-01") + sample(0:600, n, replace = TRUE)
# Acute Mechanical Support (Assigned sequentially to ensure logical timelines)
tratt_iabp_data <- as.Date(ifelse(runif(n) < 0.25, index_data_ricovero + sample(0:2, n, replace = TRUE), NA), origin = "1970-01-01")
tratt_impella_data <- as.Date(ifelse(runif(n) < 0.15, index_data_ricovero + sample(0:2, n, replace = TRUE), NA), origin = "1970-01-01")
tratt_ecmo_data <- as.Date(ifelse(runif(n) < 0.10, index_data_ricovero + sample(0:2, n, replace = TRUE), NA), origin = "1970-01-01")
tratt_trapianto_data <- as.Date(ifelse(runif(n) < 0.02, index_data_ricovero + sample(1:10, n, replace = TRUE), NA), origin = "1970-01-01")
# 4. IN-HOSPITAL OUTCOMES (Discharge or Death)
# Assuming ~30% in-hospital mortality
died_in_hosp <- runif(n) < 0.30
los <- sample(3:40, n, replace = TRUE) # Length of Stay
fu_dimiss_osp_data <- as.Date(ifelse(!died_in_hosp, index_data_ricovero + los, NA), origin = "1970-01-01")
valid_dimiss <- c("Home", "Home with service", "Hospice (Home)", "Hospice (Inpatient)",
"Inpatient Rehab", "Outpatient Rehab", "Hospital / Inpatient facility")
fu_dimiss_osp_to <- ifelse(!is.na(fu_dimiss_osp_data), sample(valid_dimiss, n, replace = TRUE), NA)
# 5. DISCHARGE VARIABLES & GDMT (Present only in alive discharged)
disch_fe <- ifelse(!is.na(fu_dimiss_osp_data), round(rnorm(n, 38, 12)), NA)
disch_tapse <- ifelse(!is.na(fu_dimiss_osp_data), round(rnorm(n, 16, 4)), NA)
disch_betabl <- ifelse(!is.na(fu_dimiss_osp_data), sample(c("Yes", "No"), n, replace = TRUE, prob=c(0.6, 0.4)), NA)
disch_acei <- ifelse(!is.na(fu_dimiss_osp_data), sample(c("Yes", "No"), n, replace = TRUE, prob=c(0.4, 0.6)), NA)
disch_sart <- ifelse(!is.na(fu_dimiss_osp_data), sample(c("Yes", "No"), n, replace = TRUE, prob=c(0.3, 0.7)), NA)
disch_entresto <- ifelse(!is.na(fu_dimiss_osp_data), sample(c("Yes", "No"), n, replace = TRUE, prob=c(0.2, 0.8)), NA)
disch_risp_pot <- ifelse(!is.na(fu_dimiss_osp_data), sample(c("Yes", "No"), n, replace = TRUE, prob=c(0.5, 0.5)), NA)
disch_sglt_2 <- ifelse(!is.na(fu_dimiss_osp_data), sample(c("Yes", "No"), n, replace = TRUE, prob=c(0.5, 0.5)), NA)
# 6. FOLLOW-UP & LONG-TERM MORTALITY
death_days_in_hosp <- los + sample(-2:0, n, replace = TRUE)
death_days_post <- los + sample(10:350, n, replace = TRUE)
fu_died_date <- as.Date(rep(NA, n))
# Dates for in-hospital death
fu_died_date[died_in_hosp] <- index_data_ricovero[died_in_hosp] + death_days_in_hosp[died_in_hosp]
# Dates for post-discharge death (~15% of the cohort)
died_post <- !died_in_hosp & (runif(n) < 0.15)
fu_died_date[died_post] <- index_data_ricovero[died_post] + death_days_post[died_post]
# Date of last follow-up for survivors
enigma_date_fu <- as.Date(ifelse(is.na(fu_died_date) & !is.na(fu_dimiss_osp_data),
fu_dimiss_osp_data + sample(30:365, n, replace = TRUE), NA), origin = "1970-01-01")
# 7. CRITICAL FORMATTING FOR ENIGMA SCRIPT
# The main script explicitly calls: as.Date(x, format = '%d/%m/%Y')
enigma_date_lvad_raw <- as.Date(ifelse(runif(n) < 0.05, index_data_ricovero + sample(5:20, n, replace = TRUE), NA), origin = "1970-01-01")
enigma_date_lvad <- format(enigma_date_lvad_raw, "%d/%m/%Y")
enigma_date_trapianto_raw <- as.Date(ifelse(runif(n) < 0.05, fu_dimiss_osp_data + sample(10:100, n, replace = TRUE), NA), origin = "1970-01-01")
enigma_date_trapianto <- format(enigma_date_trapianto_raw, "%d/%m/%Y")
enigma_death_date <- format(fu_died_date, "%d/%m/%Y")
# 8. ANCILLARY / DUMMY VARIABLES
tratt_trapianto <- "No"
fu_lvad <- "No"
enigma_trapianto_sino <- "No"
enigma_lvad_sino <- "No"
# 9. ASSEMBLE DATAFRAME directly into 'dt'
dt <- data.frame(
record_id, index_eta, demo_sesso, index_eziologia, Oh_score_scai, Oh_lab_gfr, Oh_par_vit_ega_latt, demo_diabete,
index_data_ricovero, tratt_iabp_data, tratt_impella_data, tratt_ecmo_data, tratt_trapianto_data,
fu_dimiss_osp_data, fu_dimiss_osp_to, disch_fe, disch_tapse,
disch_betabl, disch_acei, disch_sart, disch_entresto, disch_risp_pot, disch_sglt_2,
enigma_arb, enigma_mra, enigma_beta_block, enigma_sglt2,
enigma_date_lvad, enigma_date_trapianto, fu_died_date, enigma_date_fu, enigma_death_date,
tratt_trapianto, fu_lvad, enigma_trapianto_sino, enigma_lvad_sino
)
# Optional: Save a copy to disk for inspection
if (requireNamespace("writexl", quietly = TRUE)) {
writexl::write_xlsx(dt, "data/estrazione_registry_anonymized.xlsx")
cat("✅ Inline Mock Dataset Successfully Generated (N=1500) and saved to disk.\n")
} else {
cat("✅ Inline Mock Dataset Successfully Generated (N=1500) in memory.\n")
}
# ==============================================================================
# 2. CLEANING AND INITIAL PREPARATION
# ==============================================================================
# ID creation and filtering
dt <- dt %>% tidyr::separate(record_id, into = c('center', 'id'), sep = '-')
dt$idd <- paste(dt$center, dt$id, sep = '_')
dt_or <- dt
dt <- dt %>% mutate(across(c(enigma_arb, enigma_mra, enigma_beta_block, enigma_sglt2), as.character))
ids_na_SCAI <- with(dt %>% filter(is.na(`Oh_score_scai`)), paste(center, id, sep = '_'))
ids_etiology_TEP <- with(dt %>% filter(index_eziologia == 'TEP'), paste(center, id, sep = '_'))
ids_SCAI_A <- with(dt %>% filter(`Oh_score_scai` == 'A'), idd)
ids_centers_excluded <- with(dt %>% filter(!center %in% centers_considered), idd)
ids_date_beyond <- with(dt %>% filter(!is.na(index_data_ricovero) & as.Date(index_data_ricovero) > as.Date('2025-04-01')), idd)
ids_date_NA <- with(dt %>% filter(is.na(index_data_ricovero)), idd)
totake <- dt$idd[!dt$idd %in% unique(c(ids_SCAI_A, ids_na_SCAI, ids_date_beyond, ids_date_NA, ids_centers_excluded, ids_etiology_TEP))]
dt <- dt %>% filter(idd %in% totake)
# Check exclusions
length(dt_or$idd %in% ids_date_beyond)
dt_or %>% filter(!idd %in% ids_centers_excluded) %>% filter(!idd %in% ids_date_beyond) %>%
filter(!idd %in% ids_date_NA) %>% filter(!idd %in% ids_na_SCAI) %>%
filter(!idd %in% ids_SCAI_A) %>%
filter(!idd %in% ids_etiology_TEP)
diff(c(1587, 1077, 864, 841, 795, 762, 744))
# ------------------------------------------------------------------------------
# PREDICTORS MAPPING (Baseline) + AUXILIARY VARIABLES (Discharge for MICE)
# ------------------------------------------------------------------------------
dt <- dt %>%
mutate(
# Analysis Predictors (Baseline and Acute Phase)
var_eta = as.numeric(index_eta),
var_sesso = as.factor(demo_sesso),
var_eziologia = as.factor(index_eziologia),
var_scai_base = as.factor(`Oh_score_scai`),
var_egfr_base = as.numeric(Oh_lab_gfr),
var_lact_base = as.numeric(Oh_par_vit_ega_latt),
var_center = center,
# Auxiliary Variables for MICE (Measured at discharge)
var_fe_disch = as.numeric(disch_fe),
var_tapse_disch = as.numeric(disch_tapse),
# Calculation of OMT Score (0-4) at discharge (NAs will be imputed by MICE)
var_omt_score = as.numeric(disch_betabl == "Yes") +
as.numeric(disch_acei == "Yes" | disch_sart == "Yes" | disch_entresto == "Yes") +
as.numeric(disch_risp_pot == "Yes") +
as.numeric(disch_sglt_2 == "Yes"),
var_diabete = demo_diabete,
) %>%
mutate(
var_fe_disch = ifelse(var_fe_disch == 0, NA, var_fe_disch),
var_tapse_disch = ifelse(var_tapse_disch ==0, NA, var_tapse_disch),
)
t0 <- dt %>%
dplyr::select(
center, id, idd, index_data_ricovero,
fu_dimiss_osp_data, fu_dimiss_osp_to,
tratt_ecmo_data, tratt_iabp_data, tratt_impella_data,
tratt_trapianto_data, enigma_date_trapianto, enigma_date_lvad,
fu_died_date, enigma_date_fu, enigma_death_date,
tratt_trapianto, fu_lvad, enigma_trapianto_sino, enigma_lvad_sino,
# Inclusion of all modeled variables
var_eta, var_sesso, var_eziologia, var_scai_base, var_egfr_base,
var_lact_base, var_fe_disch, var_tapse_disch, var_omt_score, var_diabete
) %>%
mutate(across(
c(enigma_date_lvad, enigma_date_trapianto, enigma_death_date),
function(x) as.Date(x, format = '%d/%m/%Y')
)) %>%
# Protect clinical variables from Date conversion
mutate(across(c(everything(), -center, -id, -idd, -fu_dimiss_osp_to, -fu_lvad,
-enigma_trapianto_sino, -enigma_lvad_sino, -tratt_trapianto,
-var_eta, -var_sesso, -var_eziologia, -var_scai_base,
-var_egfr_base, -var_lact_base, -var_diabete,
-var_fe_disch, -var_tapse_disch, -var_omt_score), as.Date)) %>%
arrange(index_data_ricovero) %>%
mutate(
date_lvad = enigma_date_lvad,
last_fu = case_when(
!is.na(fu_died_date) | !is.na(enigma_death_date) ~ as.Date(pmin(fu_died_date, enigma_death_date, na.rm = TRUE)),
is.na(fu_died_date) & is.na(enigma_death_date) ~ as.Date(pmax(fu_dimiss_osp_data, enigma_date_fu, na.rm = TRUE))
),
death_date = case_when(
!is.na(fu_died_date) & !is.na(enigma_death_date) ~ fu_died_date,
TRUE ~ pmin(fu_died_date, enigma_death_date, na.rm = TRUE)
)
)
# ==============================================================================
# 3. CREATION OF TIME-DEPENDENT VARIABLES AND LOGICAL ALIGNMENT (df_clean_base)
# ==============================================================================
calc_days <- function(date_event, date_start) {
as.numeric(difftime(as.Date(date_event), as.Date(date_start), units = "days"))
}
valid_dimiss <- c("Home", "Home with service", "Hospice (Home)",
"Hospice (Inpatient)", "Inpatient Rehab",
"Outpatient Rehab", "Hospital / Inpatient facility")
grp_home_hospice <- c("Home with service", "Hospice (Home)", "Hospice (Inpatient)")
df_clean_base <- t0 %>%
mutate(
raw_t_iabp = calc_days(tratt_iabp_data, index_data_ricovero),
raw_t_impella = calc_days(tratt_impella_data, index_data_ricovero),
raw_t_ecmo = calc_days(tratt_ecmo_data, index_data_ricovero),
raw_t_dimiss = calc_days(fu_dimiss_osp_data, index_data_ricovero),
raw_t_death = calc_days(death_date, index_data_ricovero),
t_fu = calc_days(last_fu, index_data_ricovero),
raw_t_trap_pre_check = calc_days(tratt_trapianto_data, index_data_ricovero),
raw_t_trap_post_check = calc_days(enigma_date_trapianto, index_data_ricovero),
raw_t_lvad = calc_days(date_lvad, index_data_ricovero)
) %>%
mutate(across(c(starts_with("raw_t_"), t_fu), ~ as.numeric(ifelse(.x < 0, 0, .x)))) %>%
rowwise() %>%
mutate(
mcs_times = list(c(raw_t_iabp, raw_t_impella, raw_t_ecmo)),
mcs_valid = list(mcs_times[!is.na(mcs_times)]),
t_first_mcs = if(length(mcs_valid) > 0) min(mcs_valid) else NA_real_,
n_first_day = if(length(mcs_valid) > 0) sum(mcs_valid == t_first_mcs) else 0,
raw_t_multi_mcs = case_when(
n_first_day >= 2 ~ t_first_mcs,
length(mcs_valid) >= 2 ~ sort(mcs_valid)[2],
TRUE ~ NA_real_
),
status_multi_mcs = ifelse(!is.na(raw_t_multi_mcs), 1, 0),
status_iabp = ifelse(n_first_day == 1 & !is.na(raw_t_iabp) & raw_t_iabp == t_first_mcs, 1, 0),
status_impella = ifelse(n_first_day == 1 & !is.na(raw_t_impella) & raw_t_impella == t_first_mcs, 1, 0),
status_ecmo = ifelse(n_first_day == 1 & !is.na(raw_t_ecmo) & raw_t_ecmo == t_first_mcs, 1, 0)
) %>%
ungroup() %>%
mutate(
status_dimiss_base = ifelse(!is.na(raw_t_dimiss), 1, 0),
status_death = ifelse(!is.na(raw_t_death), 1, 0)
) %>%
# ----------------------------------------------------------------------------
# LOGICAL SANITY CHECK ON OBSERVATIONS AND DEATHS (TAB1 AND MSTATE ALIGNMENT)
# ----------------------------------------------------------------------------
mutate(
# 1. If a patient is discharged alive, death CANNOT precede discharge.
raw_t_death = ifelse(status_dimiss_base == 1 & status_death == 1 & raw_t_death < raw_t_dimiss, raw_t_dimiss, raw_t_death),
# 2. Fill registry gaps: Non-discharged patients without death date.
# Forced as in-hospital deaths at the day of their last follow-up.
status_death = ifelse(status_dimiss_base == 0 & is.na(raw_t_death), 1, status_death),
raw_t_death = ifelse(status_dimiss_base == 0 & is.na(raw_t_death), t_fu, raw_t_death)
) %>%
# ----------------------------------------------------------------------------
mutate(
# New covariate: Post-discharge time (until follow-up or death)
t_fu_post_dimiss = case_when(
status_dimiss_base == 1 & !is.na(t_fu) & !is.na(raw_t_dimiss) ~ pmax(0, t_fu - raw_t_dimiss),
TRUE ~ NA_real_
),
# Target Variable for MICE Imputation
fu_dimiss_osp_to_imp = case_when(
status_dimiss_base == 1 & !(fu_dimiss_osp_to %in% valid_dimiss) ~ NA_character_,
status_dimiss_base == 1 & is.na(fu_dimiss_osp_to) ~ NA_character_,
TRUE ~ fu_dimiss_osp_to
),
fu_dimiss_osp_to_imp = factor(fu_dimiss_osp_to_imp, levels = valid_dimiss),
# Total MCS duration
has_tmcs = ifelse(status_iabp == 1 | status_impella == 1 | status_ecmo == 1 | status_multi_mcs == 1, 1, 0),
durata_tmcs = case_when(
has_tmcs == 1 ~ pmax(0, raw_t_dimiss - t_first_mcs, na.rm = TRUE),
TRUE ~ 0
),
# HRT Temporal Dynamics (Pre/Post)
raw_t_lvad_pre = case_when(!is.na(raw_t_lvad) & (is.na(raw_t_dimiss) | raw_t_lvad <= raw_t_dimiss) ~ raw_t_lvad, TRUE ~ NA_real_),
raw_t_trap_pre = case_when(!is.na(raw_t_trap_pre_check) & (is.na(raw_t_dimiss) | raw_t_trap_pre_check <= raw_t_dimiss) ~ raw_t_trap_pre_check, TRUE ~ NA_real_),
raw_t_lvad_post = case_when(!is.na(raw_t_lvad) & !is.na(raw_t_dimiss) & raw_t_lvad > raw_t_dimiss ~ raw_t_lvad, TRUE ~ NA_real_),
raw_t_trap_post = case_when(!is.na(raw_t_trap_post_check) & status_dimiss_base == 1 & raw_t_trap_post_check >= raw_t_dimiss ~ raw_t_trap_post_check, TRUE ~ NA_real_),
raw_t_hrt_pre = pmin(raw_t_lvad_pre, raw_t_trap_pre, na.rm = TRUE),
raw_t_hrt_post = pmin(raw_t_lvad_post, raw_t_trap_post, na.rm = TRUE),
status_hrt_pre = ifelse(!is.na(raw_t_hrt_pre), 1, 0),
status_hrt_post = ifelse(!is.na(raw_t_hrt_post), 1, 0)
)
# Nelson-Aalen addition to properly handle survival in imputation
df_clean_base$cum_hazard <- nelsonaalen(df_clean_base, raw_t_death, status_death)
df_clean_base <- as.data.frame(df_clean_base)
# ==============================================================================
# 4. MULTIPLE IMPUTATION
# ==============================================================================
df_imputation <- df_clean_base %>%
select(
idd,
fu_dimiss_osp_to_imp, # Target
# Baseline Covariates
var_eta, var_sesso, var_eziologia, var_scai_base,
var_egfr_base, var_lact_base, var_diabete,
# tMCS Variables: Single devices separated and duration for accurate imputation
status_iabp, status_impella, status_ecmo, status_multi_mcs, durata_tmcs,
# HRT status before discharge
status_hrt_pre,
# Auxiliary Temporal Variables
t_fu_post_dimiss,
# AUXILIARY variables at discharge
var_fe_disch, var_tapse_disch, var_omt_score,
raw_t_dimiss, status_death, cum_hazard
)
ini <- mice(df_imputation, maxit = 0)
pred_mat <- ini$predictorMatrix
pred_mat[, "idd"] <- 0 # ID should not predict anything
# --- WHERE MATRIX: BLOCK IMPUTATION FOR IN-HOSPITAL DECEASED PATIENTS ---
where_mat <- is.na(df_imputation)
idx_no_dimiss <- df_clean_base$status_dimiss_base == 0
# Discharge and post-discharge follow-up variables are not imputed if there is no discharge
colonne_da_bloccare <- c("fu_dimiss_osp_to_imp", "var_fe_disch", "var_tapse_disch", "var_omt_score", "t_fu_post_dimiss")
where_mat[idx_no_dimiss, colonne_da_bloccare] <- FALSE
# ------------------------------------------------------------------
imputed_data <- mice(
df_imputation,
m = 10,
predictorMatrix = pred_mat,
where = where_mat,
seed = 2026,
printFlag = FALSE
)
# ==============================================================================
# 5. IMPUTATION DIAGNOSTICS (Plot for Supplementary Material)
# ==============================================================================
var_name <- "fu_dimiss_osp_to_imp"
dist_obs <- df_imputation %>%
filter(!is.na(.data[[var_name]])) %>%
count(.data[[var_name]]) %>%
mutate(
Prop = n / sum(n),
Type = "1_Observed (Known data)"
) %>%
rename(Categoria = all_of(var_name))
dist_imp <- complete(imputed_data, action = "long") %>%
filter(!is.na(.data[[var_name]])) %>%
count(.imp, .data[[var_name]]) %>%
group_by(.imp) %>%
mutate(Prop = n / sum(n)) %>%
group_by(Categoria = .data[[var_name]]) %>%
summarise(Prop = mean(Prop), .groups = "drop") %>%
mutate(Type = "2_Imputed (Pooled M=10)")
plot_data_cat <- bind_rows(dist_obs, dist_imp) %>%
mutate(
Categoria = as.character(Categoria),
Categoria = case_when(
Categoria == "Home with service" ~ "Home (Service)",
Categoria == "Hospice (Home)" ~ "Hospice (Home)",
Categoria == "Hospice (Inpatient)" ~ "Hospice (Inpat)",
Categoria == "Hospital / Inpatient facility" ~ "Hospital Transfer",
Categoria == "Inpatient Rehab" ~ "Inpat Rehab",
Categoria == "Outpatient Rehab" ~ "Outpat Rehab",
TRUE ~ Categoria
)
)
plot_dimissione <- ggplot(plot_data_cat, aes(x = Categoria, y = Prop, fill = Type)) +
geom_bar(stat = "identity", position = position_dodge(width = 0.8), width = 0.7, color = "black") +
scale_y_continuous() +
scale_fill_manual(values = c("1_Observed (Known data)" = "grey70",
"2_Imputed (Pooled M=10)" = "tomato")) +
labs(
title = "Imputation Diagnostics: Discharge Setting",
subtitle = "Distributions computed solely on patients surviving to hospital discharge",
x = "Discharge Setting",
y = "Proportion (%)",
fill = "Dataset"
) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1, size = 11, face = "bold"),
panel.grid.major.x = element_blank()
)
print(plot_dimissione)
# ==============================================================================
# 6. MSTATE EXPANSION FOR EACH IMPUTED DATASET
# ==============================================================================
stati_clinici <- c(
"Init", "IABP", "Impella", "ECMO", "Multi_MCS", "HRT_Pre",
"Home", "Home_Serv_Hospice", "Inpat_Rehab", "Outpat_Rehab", "Hospital",
"HRT_Post", "Death"
)
idx_multi <- 5; idx_hrt_pre <- 6; idx_dimiss <- 7:11; idx_hrt_post <- 12; idx_morte <- 13
tmat <- transMat(x = list(
c(2, 3, 4, idx_multi, idx_hrt_pre, idx_dimiss, idx_morte), # 1. Init
c(idx_multi, idx_hrt_pre, idx_dimiss, idx_morte), # 2. IABP
c(idx_multi, idx_hrt_pre, idx_dimiss, idx_morte), # 3. Impella
c(idx_multi, idx_hrt_pre, idx_dimiss, idx_morte), # 4. ECMO
c(idx_hrt_pre, idx_dimiss, idx_morte), # 5. Multi_MCS
c(idx_dimiss, idx_morte), # 6. HRT_Pre
c(idx_hrt_post, idx_morte), # 7. Home
c(idx_hrt_post, idx_morte), # 8. Home_Serv_Hospice
c(idx_hrt_post, idx_morte), # 9. Inpat_Rehab
c(idx_hrt_post, idx_morte), # 10. Outpat_Rehab
c(idx_hrt_post, idx_morte), # 11. Hospital
c(idx_morte), # 12. HRT_Post
c() # 13. Death
), names = stati_clinici)
col_tempi <- c(NA, "time_iabp", "time_impella", "time_ecmo", "time_multi_mcs", "time_hrt_pre",
rep("time_dimiss_all", 5), "time_hrt_post", "time_death")
col_status <- c(NA, "status_iabp", "status_impella", "status_ecmo", "status_multi_mcs", "status_hrt_pre",
"st_home", "st_home_srv_hosp", "st_inpat_rehab", "st_outpat_rehab",
"st_hospital", "status_hrt_post", "status_death")
col_to_force <- unique(c(col_tempi[-1], col_status[-1]))
ms_datasets_imputed <- list()
for (i in 1:imputed_data$m) {
df_comp <- complete(imputed_data, action = i)
df_i <- df_clean_base
# Safe Global Overwrite with imputed data
cols_imputed <- intersect(names(df_comp), names(df_i))
df_i[cols_imputed] <- df_comp[cols_imputed]
df_i <- df_i %>%
mutate(
st_home = ifelse(status_dimiss_base == 1 & fu_dimiss_osp_to_imp %in% "Home", 1, 0),
st_home_srv_hosp = ifelse(status_dimiss_base == 1 & fu_dimiss_osp_to_imp %in% grp_home_hospice, 1, 0),
st_inpat_rehab = ifelse(status_dimiss_base == 1 & fu_dimiss_osp_to_imp %in% "Inpatient Rehab", 1, 0),
st_outpat_rehab = ifelse(status_dimiss_base == 1 & fu_dimiss_osp_to_imp %in% "Outpatient Rehab", 1, 0),
st_hospital = ifelse(status_dimiss_base == 1 & fu_dimiss_osp_to_imp %in% "Hospital / Inpatient facility", 1, 0),
max_obs_time = pmax(raw_t_iabp, raw_t_impella, raw_t_ecmo, raw_t_multi_mcs,
raw_t_hrt_pre, raw_t_dimiss,
raw_t_hrt_post, raw_t_death, t_fu, na.rm = TRUE),
max_obs_time = ifelse(is.na(max_obs_time) | is.infinite(max_obs_time), 0, max_obs_time),
t_censor_global = max_obs_time + 0.008,
time_iabp = ifelse(status_iabp == 1, raw_t_iabp + 0.001, t_censor_global),
time_impella = ifelse(status_impella == 1, raw_t_impella + 0.001, t_censor_global),
time_ecmo = ifelse(status_ecmo == 1, raw_t_ecmo + 0.001, t_censor_global),
time_multi_mcs = ifelse(status_multi_mcs == 1, raw_t_multi_mcs + 0.0015, t_censor_global),
time_hrt_pre = ifelse(status_hrt_pre == 1, raw_t_hrt_pre + 0.002, t_censor_global),
# SIMULTANEOUS EVENTS MANAGEMENT (Tie-breaking for Mstate)
# Discharge wins temporal tie over death (guaranteeing transition to settings)
time_dimiss_all = ifelse(status_dimiss_base == 1, raw_t_dimiss + 0.003, t_censor_global),
time_hrt_post = ifelse(status_hrt_post == 1, pmax(raw_t_hrt_post, raw_t_dimiss, na.rm=TRUE) + 0.004, t_censor_global),
time_death = ifelse(status_death == 1, pmax(raw_t_death, raw_t_dimiss, raw_t_hrt_post, na.rm=TRUE) + 0.005, t_censor_global)
)
for(col in col_to_force) { df_i[[col]] <- as.numeric(unlist(df_i[[col]])) }
ms_data_i <- msprep(time = col_tempi, status = col_status, data = df_i, trans = tmat, id = "idd")
# --- CARRYING COVARIATES INTO LONG FORMAT ---
ms_data_i$var_eta <- df_i$var_eta[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_sesso <- df_i$var_sesso[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_eziologia <- df_i$var_eziologia[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_scai_base <- df_i$var_scai_base[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_egfr_base <- df_i$var_egfr_base[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_lact_base <- df_i$var_lact_base[match(ms_data_i$idd, df_i$idd)]
ms_data_i$status_hrt_pre <- df_i$status_hrt_pre[match(ms_data_i$idd, df_i$idd)]
ms_data_i$status_iabp <- df_i$status_iabp[match(ms_data_i$idd, df_i$idd)]
ms_data_i$status_impella <- df_i$status_impella[match(ms_data_i$idd, df_i$idd)]
ms_data_i$status_ecmo <- df_i$status_ecmo[match(ms_data_i$idd, df_i$idd)]
ms_data_i$status_multi_mcs<- df_i$status_multi_mcs[match(ms_data_i$idd, df_i$idd)]
ms_data_i$durata_tmcs <- df_i$durata_tmcs[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_fe_disch <- df_i$var_fe_disch[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_tapse_disch <- df_i$var_tapse_disch[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_omt_score <- df_i$var_omt_score[match(ms_data_i$idd, df_i$idd)]
ms_data_i$var_diabete <- df_i$var_diabete[match(ms_data_i$idd, df_i$idd)]
ms_data_i$center <- df_i$center[match(ms_data_i$idd, df_i$idd)]
ms_datasets_imputed[[i]] <- ms_data_i
}
# ==============================================================================
# 7. EVENTS POOLING AND PLOT OF THE AVERAGED TRANSITION MATRIX
# ==============================================================================
base_ev <- events(ms_datasets_imputed[[1]])
pool_freq <- base_ev$Frequencies * 0
pool_prop <- base_ev$Proportions * 0
for (i in 1:imputed_data$m) {
ev <- events(ms_datasets_imputed[[i]])
pool_freq <- pool_freq + ev$Frequencies
pool_prop <- pool_prop + ev$Proportions
}
pool_freq <- pool_freq / imputed_data$m
pool_prop <- pool_prop / imputed_data$m
df_freq <- as.data.frame(as.table(pool_freq))
df_prop <- as.data.frame(as.table(pool_prop))
df_plot <- df_freq %>%
inner_join(df_prop %>% rename(prop = Freq), by = c('from', 'to')) %>%
filter(from != "total entering", to != "no event") %>%
filter(Freq > 0)
plot_transizioni <- ggplot(df_plot, aes(x = to, y = from, fill = prop)) +
geom_tile() +
scale_fill_gradientn(colors = terrain.colors(10)) +
geom_text(aes(label = round(Freq)), col = 'black') +
theme_minimal() +
theme(legend.position = 'none',
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
print(plot_transizioni)
# ==============================================================================
# 8. TABLE 1: POST-HOC AND ALIGNED WITH MSTATE (SPLIT & STACKED)
# ==============================================================================
table1_data <- df_clean_base %>%
dplyr::select(
center, var_eta, var_sesso, var_eziologia, var_scai_base, var_egfr_base, var_lact_base, var_diabete,
status_dimiss_base, fu_dimiss_osp_to,
var_fe_disch, var_tapse_disch, var_omt_score
) %>%
mutate(
# Binary status transformation
across(starts_with("status_"), ~ factor(.x, levels = c(0, 1), labels = c("No", "Yes"))),
fu_dimiss_osp_to = as.character(fu_dimiss_osp_to),
# Since status_dimiss_base = 0 indicates those who did not leave alive and the
# Mstate logic funnels them into Death, we force it directly to "Died in hospital".
fu_dimiss_osp_to = case_when(
status_dimiss_base == 'No' ~ 'Died in hospital',
TRUE ~ fu_dimiss_osp_to
),
fu_dimiss_osp_to = as.factor(fu_dimiss_osp_to)
)
t1_base <- table1_data %>%
dplyr::select(-var_fe_disch, -var_tapse_disch, -var_omt_score) %>%
tbl_summary(
statistic = list(
all_continuous() ~ "{median} ({IQR})",
all_categorical() ~ "{n} ({p}%)"
),
digits = all_continuous() ~ 1,
missing = "ifany",
missing_text = "Missing",
missing_stat = "{N_miss} ({p_miss}%)",
label = list(
var_eta ~ "Age (Years)",
var_sesso ~ "Sex",
var_eziologia ~ "Etiology",
var_scai_base ~ "SCAI Shock Stage",
var_egfr_base ~ "Baseline eGFR (ml/min)",
var_lact_base ~ "Baseline Lactate (mmol/L)",
var_diabete ~ "Type 2 Diabetes",
fu_dimiss_osp_to ~ "Discharge Setting"
)
)
t1_disch <- table1_data %>%
filter(status_dimiss_base == 'Yes') %>%
dplyr::select(var_fe_disch, var_tapse_disch, var_omt_score) %>%
tbl_summary(
statistic = list(
all_continuous() ~ "{median} ({IQR})",
all_categorical() ~ "{n} ({p}%)"
),
digits = all_continuous() ~ 1,
missing = "ifany",
missing_text = "Missing",
missing_stat = "{N_miss} ({p_miss}%)",
label = list(
var_fe_disch ~ "Discharge LVEF (%)",
var_tapse_disch ~ "Discharge TAPSE (mm)",
var_omt_score ~ "Discharge OMT Score (0-4)"
)
)
table1_final <- tbl_stack(
list(t1_base, t1_disch),
group_header = c("Overall Cohort (Admission & Clinical Course)",
"Discharge Parameters (Among Hospital Survivors)")
) %>%
modify_header(label = "**Clinical Characteristic**") %>%
bold_labels()
print(table1_final)
# --- OVERALL SURVIVAL PLOT ---
survfit(Surv(t_fu, status_death) ~ 1, data = df_clean_base) -> fit
summary(fit, t = c(30, 60, 180, 365)) -> s
plot_km <- ggsurvfit(survfit(Surv(t_fu, status_death) ~ 1, data = df_clean_base))+
add_confidence_interval()+
geom_segment(inherit.aes = F, data = data.frame(time = s$time, surv = s$surv), aes(x = 0, xend = time, y = surv, yend= surv), linetype = 'dotted')+
geom_text(inherit.aes = F, data = data.frame(time = s$time, surv = s$surv), aes(x = -6, y = surv, label = round(surv,2) ), size = 2)+
coord_cartesian(xlim = c(0,366))+
ggpubr::theme_classic2()+
labs(title = "Overall Cohort Survival", x = 'Time from index event (days)', y = 'Survival Probability')+
add_risktable()+
scale_x_continuous(breaks = c(0, 30, 60, 180, 365))+
add_censor_mark()
print(plot_km)
# ==============================================================================
# 9. MULTI-STATE ANALYSIS AND POOLING OF BASELINE HAZARDS (MICE-COMPATIBLE)
# ==============================================================================
list_msfit <- list()
for (i in 1:imputed_data$m) {
ms_data_i <- ms_datasets_imputed[[i]]
cox_mod_i <- coxph(Surv(Tstart, Tstop, status) ~ strata(trans),
data = ms_data_i, method = "breslow")
list_msfit[[i]] <- msfit(cox_mod_i, trans = tmat)
}
pooled_hazards <- list_msfit[[1]]
pooled_hazards$Haz$Haz <- 0
pooled_hazards$varHaz$varHaz <- 0
for (i in 1:imputed_data$m) {
pooled_hazards$Haz$Haz <- pooled_hazards$Haz$Haz + list_msfit[[i]]$Haz$Haz
pooled_hazards$varHaz$varHaz <- pooled_hazards$varHaz$varHaz + list_msfit[[i]]$varHaz$varHaz
}
pooled_hazards$Haz$Haz <- pooled_hazards$Haz$Haz / imputed_data$m
pooled_hazards$varHaz$varHaz <- pooled_hazards$varHaz$varHaz / imputed_data$m
prob_transizioni <- probtrans(pooled_hazards, predt = 0)
# ==============================================================================
# DATA TRANSFORMATION INTO LONG FORMAT FOR GGPLOT
# ==============================================================================
pt_df <- prob_transizioni[[1]]
pt_long <- pt_df %>%
dplyr::select(time, starts_with("pstate")) %>%
pivot_longer(cols = starts_with("pstate"), names_to = "state_col", values_to = "prob") %>%
mutate(
state_num = as.numeric(gsub("pstate", "", state_col)),
state_name = factor(stati_clinici[state_num], levels = stati_clinici)
) %>%
mutate(macro = case_when(
state_name %in% c('IABP', 'Impella', 'Multi_MCS', 'ECMO') ~ 'MCS',
state_name %in% c('Inpat_Rehab', 'Outpat_Rehab', 'Home', 'Hospital',
'Home_Serv_Hospice') ~ 'Post-Discharge',
TRUE ~ 'base'
))
# Plot 1: State Occupation Probabilities (Faceted by Macro)
plot_occupazione <- ggplot(pt_long %>% filter(time <= 365),
aes(x = time, y = prob, col = state_name)) +
facet_wrap(~macro, scales = 'free_y', ncol = 2) +
geom_line(lwd = 1.2) +
scale_x_continuous(limits = c(0, 365), expand = c(0, 0)) +
labs(
title = "State Occupation Probabilities",
x = "Days from admission",
y = "Probability",
col = "Clinical states"
) +
ggpubr::theme_pubclean() +
theme(
legend.position = "bottom",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) +
coord_cartesian(xlim = c(0, 180)) +
ggsci::scale_color_simpsons(labels = c('Still in Hospital', 'IABP', 'Impella',
'ECMO', 'Mulit-MCS', 'HRT-pre discharge',
'Discharged at Home', 'Palliative care / Hospice',
'Inpatient rehab', 'Outpatient Rehab', 'Discharged in Hospital',
'HRT post discharge', 'Death'))
print(plot_occupazione)
# ==============================================================================
# 10. MODELS: IN-HOSPITAL & POST-DISCHARGE (CONVOLUTION)
# ==============================================================================
list_multinom <- list()
for (i in 1:imputed_data$m) {
ms_data_i <- ms_datasets_imputed[[i]]
data_intervals <- ms_data_i %>%
filter(from %in% 1:6) %>%
group_by(idd, Tstart, Tstop) %>%
summarise(
event_occurred = max(status),
event_to = if(any(status == 1)) to[status == 1] else 0,
var_eta = first(var_eta),
var_sesso = first(var_sesso),
var_eziologia = first(var_eziologia),
var_scai_base = first(var_scai_base),
var_egfr_base = first(var_egfr_base),
var_lact_base = first(var_lact_base),
var_diabete = first(var_diabete),
var_center = first(center),
current_state = first(from),
.groups = "drop"
) %>%
mutate(
outcome_class = case_when(
event_occurred == 0 ~ "0_None",
event_to == 7 ~ "1_Home",
event_to %in% c(9, 10) ~ "2_Rehab",
event_to %in% c(11) ~ "3_Hospital",
event_to %in% 8 ~ '4_Hospice/Palliative Care',
event_to == 13 ~ "5_Death"
),
var_eta_10 = var_eta / 25,
var_egfr_10 = var_egfr_base / 20,
var_lact_base = var_lact_base / 5,
var_eziologia = as.character(var_eziologia),
var_eziologia = case_when(var_eziologia == 'AMI'~'AMI', var_eziologia == 'HF-CS'~'HF-CS', TRUE~'Other'),
tdc_mcs_type = case_when(
current_state %in% c(1, 6) ~ "1_None",
current_state == 2 ~ "2_IABP",
current_state == 3 ~ "3_Impella",
current_state == 4 ~ "4_ECMO",
current_state == 5 ~ "5_Multi_MCS"
),
tdc_mcs_type = factor(tdc_mcs_type, levels = c("1_None", "2_IABP", "3_Impella", "4_ECMO", "5_Multi_MCS")),
tdc_hrt_pre = factor(ifelse(current_state == 6, "Yes", "No"), levels = c("No", "Yes"))
)
data_intervals$dummy_event <- ifelse(data_intervals$event_occurred == 1, 1, 0)
max_giorni <- ceiling(max(data_intervals$Tstop, na.rm = TRUE))
data_daily <- survSplit(Surv(Tstart, Tstop, dummy_event) ~ ., data = data_intervals, cut = seq(1, max_giorni, by = 1)) %>%
mutate(
daily_outcome = ifelse(dummy_event == 1, as.character(outcome_class), "0_None"),
daily_outcome = factor(daily_outcome, levels = c("0_None", "1_Home", "2_Rehab", "3_Hospital", "4_Hospice/Palliative Care", '5_Death')),
day = Tstop
)
list_multinom[[i]] <- multinom(
daily_outcome ~ ns(day, df = 3) + var_eta_10 + var_sesso + var_eziologia + var_scai_base + var_egfr_10 +
var_lact_base + tdc_mcs_type + tdc_hrt_pre + var_center,
data = data_daily, trace = FALSE, maxit = 500, decay = 0.05
)
}
list_glm_post = list()
list_data_post <- list()
for (i in 1:imputed_data$m) {
ms_data_i <- ms_datasets_imputed[[i]]
storico_pazienti <- ms_data_i %>%
filter(from %in% 1:5) %>%
group_by(idd) %>%
summarise(max_mcs_state = suppressWarnings(max(from, na.rm = TRUE)), .groups = 'drop') %>%
mutate(
storico_mcs = case_when(
max_mcs_state == 1 ~ "1_None",
max_mcs_state == 2 ~ "2_IABP",
max_mcs_state == 3 ~ "3_Impella",
max_mcs_state == 4 ~ "4_ECMO",
max_mcs_state == 5 ~ "5_Multi_MCS",
TRUE ~ "1_None"
),
storico_mcs = factor(storico_mcs, levels = c("1_None", "2_IABP", "3_Impella", "4_ECMO", "5_Multi_MCS"))
)
data_post <- ms_data_i %>%
filter(from %in% 7:11 & to == 13) %>%
left_join(storico_pazienti, by = "idd") %>%
mutate(
var_center = center,
var_eta_10 = var_eta / 25,
var_egfr_10 = var_egfr_base / 20,
var_lact_base = var_lact_base / 5,
var_fe_disch = case_when(var_fe_disch < 40~'<40', var_fe_disch >=40~'≥40'),
var_eziologia = as.character(var_eziologia),
var_eziologia = case_when(var_eziologia == 'AMI' ~ 'AMI', var_eziologia == 'HF-CS' ~ 'HF-CS', TRUE ~ 'Other'),
Tstop_reset = Tstop - Tstart,
Tstart_reset = 0,
dummy_event = status,
setting_disch = case_when(
from == 7 ~ "1_Home",
from %in% c(9, 10) ~ "2_Rehab",
from %in% c(11) ~ "3_Hospital",
from %in% 8 ~ "4_Hospice/Palliative Care"
),
setting_disch = factor(setting_disch, levels = c("1_Home", "2_Rehab", "3_Hospital", "4_Hospice/Palliative Care")),
setting_disch = relevel(setting_disch, ref = '1_Home'),
storico_mcs = if_else(is.na(storico_mcs), factor("1_None", levels = levels(storico_mcs)), storico_mcs)
)
list_data_post[[i]] <- data_post
max_giorni_post <- ceiling(max(data_post$Tstop_reset, na.rm = TRUE))
data_post_daily <- survSplit(Surv(Tstart_reset, Tstop_reset, dummy_event) ~ .,
data = data_post, cut = seq(1, max_giorni_post, by = 1)) %>%
mutate(day = Tstop_reset)
list_glm_post[[i]] <- glm(
dummy_event ~ ns(day, df = 3) +
var_eta_10 + var_sesso + var_eziologia + var_scai_base +
var_egfr_10 + var_lact_base + status_hrt_pre +
setting_disch + var_fe_disch + var_omt_score + var_diabete +
storico_mcs+var_center,
data = data_post_daily,
family = binomial(link = "logit")
)
}
# ==============================================================================
# 11. FOREST PLOT AND PLOTS SAVING (USING /OUTPUT FOLDER)
# ==============================================================================
pulisci_termini <- function(term) {
case_when(
term == "var_eta_10" ~ "Age (+25 years)",
term == "var_sessoM" ~ "Gender: Male (vs Female)",
term == "var_eziologiaHF-CS" ~ "Etiology: HF-CS (vs AMI)",
term == "var_eziologiaOther" ~ "Etiology: Others (vs AMI)",
term == "var_scai_baseB" ~ "SCAI stage: B (vs B)",
term == "var_scai_baseC" ~ "SCAI stage: C (vs B)",
term == "var_scai_baseD" ~ "SCAI stage: D (vs B)",
term == "var_scai_baseE" ~ "SCAI stage: E (vs B)",
term == "var_egfr_10" ~ "Baseline eGFR (+20 ml/min)",
term == "var_lact_base" ~ "Baseline Lactate (+5 mmol/L)",
term == "tdc_mcs_type2_IABP" ~ "MCS (Acute): IABP (vs None)",
term == "tdc_mcs_type3_Impella" ~ "MCS (Acute): Impella (vs None)",
term == "tdc_mcs_type4_ECMO" ~ "MCS (Acuto): ECMO (vs None)",
term == "tdc_mcs_type5_Multi_MCS" ~ "MCS (Acute): Multi-MCS (vs None)",
term == "storico_mcs2_IABP" ~ "MCS (Max): IABP (vs None)",
term == "storico_mcs3_Impella" ~ "MCS (Max): Impella (vs None)",
term == "storico_mcs4_ECMO" ~ "MCS (Max): ECMO (vs None)",
term == "storico_mcs5_Multi_MCS" ~ "MCS (Max): Multi-MCS (vs None)",
term == "setting_disch2_Rehab" ~ " Discharfed Rehab (vs Home)",
term == "setting_disch3_Hospital" ~ "Hosopital transfer (vs Home)",
term == "setting_disch4_Hospice/Palliative Care" ~ "Discharged in Hospice (vs Home)",
term == "var_fe_disch<40" ~ "LVEF at discharge: <40% (vs ≥40%)",
term == "var_omt_score" ~ "GDMT at discharge (+1)",
term == "var_diabeteYes" ~ "Type 2 Diabetes: Yes (vs No)",
term %in% c("status_hrt_pre", "status_hrt_pre1", "status_hrt_preYes") ~ "HRT Pre-discharge (vs None)",
term == "tdc_hrt_preYes" ~ "Acute Phase: HRT Pre-discharge (vs None)",
TRUE ~ term
)
}
# ==============================================================================
# 11.1 IN-HOSPITAL FOREST PLOT (Multinomial)
# ==============================================================================
cat("Pooling in-hospital model...\n")
M <- length(list_multinom)
# EXACT MODEL NAMES
esiti <- c("1_Home", "2_Rehab", "3_Hospital", "4_Hospice/Palliative Care", "5_Death")
lista_coefs <- lapply(list_multinom, function(m) summary(m)$coefficients)
lista_vars <- lapply(list_multinom, function(m) (summary(m)$standard.errors)^2)
termini <- colnames(lista_coefs[[1]])
df_intra <- data.frame()
for (e in esiti) {
pooled_beta <- numeric(length(termini))
pooled_se <- numeric(length(termini))
for (j in 1:length(termini)) {
betas_m <- sapply(1:M, function(i) lista_coefs[[i]][e, j])
vars_m <- sapply(1:M, function(i) lista_vars[[i]][e, j])
Q_bar <- mean(betas_m); U_bar <- mean(vars_m); B <- var(betas_m)
T_var <- U_bar + B + (B / M)
pooled_beta[j] <- Q_bar; pooled_se[j] <- sqrt(T_var)
}
df_intra <- bind_rows(df_intra, data.frame(
Esito = e, Termine = termini, Odds.Ratio = exp(pooled_beta),
CI_inf = exp(pooled_beta - 1.96 * pooled_se), CI_sup = exp(pooled_beta + 1.96 * pooled_se)
))
}
df_intra <- df_intra %>%
filter(!grepl("ns\\(|Intercept|var_center", Termine)) %>%
mutate(
Termine_Clean = pulisci_termini(Termine),
Termine_Clean = factor(Termine_Clean, levels = rev(unique(Termine_Clean))),
Esito = factor(Esito, levels = c("5_Death", "4_Hospice/Palliative Care", "3_Hospital", "2_Rehab", "1_Home"),
labels = c("Acute Phase Mortality", "Hospice", "Hospital transfer", "Rehabilitation", "Home"))
)
plot_forest_intra <- ggplot(df_intra, aes(xmin = CI_inf, xmax = CI_sup, x = Odds.Ratio, y = Termine_Clean)) +
geom_pointrange(size = 0.1, color = "black") +
geom_vline(xintercept = 1, linetype = "dashed", color = "red") +
scale_x_log10(breaks = c(0.1, 0.5, 1, 2, 5, 10, 20)) +
facet_wrap(~ Esito, ncol = 2) +
labs(x = "Odds Ratio (Logarithmic Scale)", y = "") +
theme_bw(base_size = 10) +
theme(axis.text.y = element_text(face = "bold"),
strip.background = element_rect(fill = "#e0e0e0", color = "black"),
strip.text = element_text(face = "bold", color = "black", size = 11),
panel.grid.minor = element_blank())
png('output/fig_intra.png', width = 30, height = 30, unit = 'cm', res = 1000)
print(plot_forest_intra)
dev.off()
# Reshape data and apply the two-sided significance check
df_wide <- df_intra %>%
mutate(
# TWO-SIDED CHECK: Significant if the entire 95% CI is strictly above 1 OR strictly below 1
is_significant = (CI_inf > 1.00) | (CI_sup < 1.00),
# Format all numeric columns strictly to 2 decimal places
across(c(Odds.Ratio, CI_inf, CI_sup), ~ sprintf("%.2f", .x)),
# Construct the classic clinical string format
base_est = paste0(Odds.Ratio, " (", CI_inf, "–", CI_sup, ")"),
# Dynamically append asterisk if either condition is met
OR_CI = if_else(is_significant, paste0(base_est, "*"), base_est)
) %>%
select(Termine_Clean, Esito, OR_CI) %>%
pivot_wider(
names_from = Esito,
values_from = OR_CI
)
# Build the horizontal presentation table
gt_horizontal <- df_wide %>%
gt() %>%
cols_label(
Termine_Clean = "Predictor / Variable"
) %>%
tab_header(
title = md("**Multivariable Predictors Across Clinical Outcomes**"),
subtitle = "Values are expressed as Odds Ratio (95% Confidence Interval)"
) %>%
tab_footnote(
footnote = "* Statistically significant: 95% Confidence Interval does not overlap 1.00 (p < 0.05).",
locations = cells_title(groups = "subtitle")
) %>%
cols_align(
align = "left",
columns = Termine_Clean
) %>%
cols_align(
align = "center",
columns = -Termine_Clean
) %>%
tab_options(
heading.title.font.size = px(16),
heading.subtitle.font.size = px(12),
table.border.top.color = "black",
table.border.bottom.color = "black",
table_body.hlines.color = "#eaeaea",
column_labels.font.weight = "bold"
)
# Print the table
print(gt_horizontal)
# ==============================================================================
# 11.2 POST-DISCHARGE FOREST PLOT
# ==============================================================================
cat("Pooling post-discharge model...\n")
pooled_post <- pool(as.mira(list_glm_post))
df_post <- summary(pooled_post, exponentiate = TRUE, conf.int = TRUE) %>%
rename(Odds.Ratio = estimate, CI_inf = `2.5 %`, CI_sup = `97.5 %`) %>%
filter(!grepl("ns\\(|Intercept|var_center", term)) %>%
mutate(
Termine_Clean = pulisci_termini(term),
Termine_Clean = factor(Termine_Clean, levels = rev(unique(Termine_Clean)))
)
plot_forest_post <- ggplot(df_post, aes(xmin = CI_inf, xmax = CI_sup, x = Odds.Ratio, y = Termine_Clean)) +
geom_pointrange(color = "black", size = 0.2) +
geom_vline(xintercept = 1, linetype = "dashed", color = "red") +
scale_x_log10(breaks = c(0.1, 0.5, 1, 2, 5, 10, 20)) +
labs(x = "Odds Ratio (Logarithmic Scale)", y = "") +
theme_minimal(base_size = 13) +
theme(axis.text.y = element_text(face = "bold"), panel.grid.minor = element_blank())
png('output/fig_post.png', width = 30, height = 30, unit = 'cm', res = 1000)
print(plot_forest_post)
dev.off()
# Reshape data and apply the two-sided significance check
df_wide <- df_post %>%
mutate(
# TWO-SIDED CHECK: Significant if the entire 95% CI is strictly above 1 OR strictly below 1
is_significant = (CI_inf > 1.00) | (CI_sup < 1.00),
# Format all numeric columns strictly to 2 decimal places
across(c(Odds.Ratio, CI_inf, CI_sup), ~ sprintf("%.2f", .x)),
# Construct the classic clinical string format
base_est = paste0(Odds.Ratio, " (", CI_inf, "–", CI_sup, ")"),
# Dynamically append asterisk if either condition is met
OR_CI = if_else(is_significant, paste0(base_est, "*"), base_est)
) %>%
select(Termine_Clean, OR_CI)
# Build the horizontal presentation table
gt_horizontal_post <- df_wide %>%
gt() %>%
cols_label(
Termine_Clean = "Predictor / Variable"
) %>%
tab_header(
title = md("**Multivariable Predictors of Post Discharge Mortality**"),
subtitle = "Values are expressed as Odds Ratio (95% Confidence Interval)"
) %>%
tab_footnote(
footnote = "* Statistically significant: p < 0.05.",
locations = cells_title(groups = "subtitle")
) %>%
cols_align(
align = "left",
columns = Termine_Clean
) %>%
cols_align(
align = "center",
columns = -Termine_Clean
) %>%
tab_options(
heading.title.font.size = px(16),
heading.subtitle.font.size = px(12),
table.border.top.color = "black",
table.border.bottom.color = "black",
table_body.hlines.color = "#eaeaea",
column_labels.font.weight = "bold"
)
# Print the table
print(gt_horizontal_post)
# ==============================================================================
# 11.3 STATE OCCUPATION PROBABILITIES FOR 4 SPECIFIC PROFILES (Clinical Covariates)
# ==============================================================================
cat("Calculating S.O.P. for the 4 profiles...\n")
max_giorni_totali <- 365
liv_sesso <- list_multinom[[1]]$xlevels$var_sesso
liv_ezio <- list_multinom[[1]]$xlevels$var_eziologia
liv_scai <- list_multinom[[1]]$xlevels$var_scai_base
liv_mcs <- list_multinom[[1]]$xlevels$tdc_mcs_type
liv_hrt <- list_multinom[[1]]$xlevels$tdc_hrt_pre
liv_center <- list_multinom[[1]]$xlevels$var_center
centro_default <- liv_center[1]
# Define constant background levels (Reference 65yo Patient)
val_eta <- 65 / 25
val_scai <- "C"
val_egfr <- 60 / 20
val_lact <- 4.3 / 5
val_mcs <- "1_None"
val_hrt <- "No"
val_diab <- "No"
val_omt <- 2
val_fe <- "≥40"
crea_profilo <- function(nome, sesso, ezio) {
data.frame(
Profilo = nome,
day = 1:max_giorni_totali,
var_eta_10 = val_eta,
var_sesso = factor(sesso, levels = liv_sesso),
var_eziologia = factor(ezio, levels = liv_ezio),
var_scai_base = factor(val_scai, levels = liv_scai),
var_egfr_10 = val_egfr,
var_lact_base = val_lact,
tdc_mcs_type = factor(val_mcs, levels = liv_mcs),
tdc_hrt_pre = factor(val_hrt, levels = liv_hrt),
var_center = factor(centro_default, levels = liv_center)
)
}
profili_estesi <- bind_rows(
crea_profilo("65y, Male, AMI", "M", "AMI"),
crea_profilo("65y, Female, AMI", "F", "AMI"),
crea_profilo("65y, Male, HF-CS", "M", "HF-CS"),
crea_profilo("65y, Female, HF-CS", "F", "HF-CS")
)
# Incidence Extraction
estrai_incidenza <- function(m, dp) {
preds <- predict(m, newdata = dp, type = "probs")
cbind(dp[, c("Profilo", "day")], preds) %>% group_by(Profilo) %>%
mutate(Prob_Stay_Hosp = cumprod(`0_None`), S_prev = lag(Prob_Stay_Hosp, default = 1),
Inc_Home = S_prev * `1_Home`,
Inc_Rehab = S_prev * `2_Rehab`,
Inc_Hospital = S_prev * `3_Hospital`,
Inc_Hospice = S_prev * `4_Hospice/Palliative Care`,
Inc_Death_Hosp = S_prev * `5_Death`) %>%
ungroup() %>% dplyr::select(Profilo, day, Prob_Stay_Hosp, Inc_Home, Inc_Rehab, Inc_Hospital, Inc_Hospice, Inc_Death_Hosp)
}
lista_inc <- lapply(list_multinom, function(m) estrai_incidenza(m, profili_estesi))
inc_pooled <- bind_rows(lista_inc) %>% group_by(Profilo, day) %>% summarise(across(everything(), mean), .groups="drop")
# Post-Discharge Survival Extraction
profili_post_completi <- profili_estesi %>%
filter(day == 1) %>% dplyr::select(-any_of("setting_disch")) %>%
tidyr::crossing(setting_disch = c("1_Home", "2_Rehab", "3_Hospital", "4_Hospice/Palliative Care")) %>%
mutate(
status_hrt_pre = 0,
var_omt_score = val_omt,
var_diabete = val_diab,
var_fe_disch = val_fe,
var_egfr_10 = val_egfr,
storico_mcs = factor(val_mcs, levels = liv_mcs)
) %>%
group_by(Profilo, setting_disch) %>% slice(rep(1, each = max_giorni_totali)) %>% mutate(day = 1:max_giorni_totali) %>% ungroup()
estrai_sopravvivenza <- function(m, dp) {
dp$hazard <- predict(m, newdata = dp, type = "response")
if(any(is.na(dp$hazard))) dp$hazard[is.na(dp$hazard)] <- 0
dp %>% group_by(Profilo, setting_disch) %>% arrange(day) %>%
mutate(Survival = cumprod(1 - hazard)) %>% ungroup() %>% dplyr::select(Profilo, setting_disch, day, Survival)
}
lista_sur <- lapply(list_glm_post, function(m) estrai_sopravvivenza(m, profili_post_completi))
sur_pooled <- bind_rows(lista_sur) %>% group_by(Profilo, setting_disch, day) %>% summarise(Survival = mean(Survival), .groups="drop")
# Convolution
sop_finale <- inc_pooled %>% arrange(Profilo, day) %>% mutate(Occ_Home = 0, Occ_Rehab = 0, Occ_Hospice = 0, Occ_Hospital = 0)
for (p in unique(sop_finale$Profilo)) {
temp_p <- sop_finale %>% filter(Profilo == p) %>% arrange(day)
i_home <- temp_p$Inc_Home; i_rehab <- temp_p$Inc_Rehab; i_hospital <- temp_p$Inc_Hospital; i_hospice <- temp_p$Inc_Hospice
s_home <- (sur_pooled %>% filter(Profilo == p, setting_disch == "1_Home") %>% arrange(day))$Survival
s_rehab <- (sur_pooled %>% filter(Profilo == p, setting_disch == "2_Rehab") %>% arrange(day))$Survival
s_hospital <- (sur_pooled %>% filter(Profilo == p, setting_disch == "3_Hospital") %>% arrange(day))$Survival
s_hospice <- (sur_pooled %>% filter(Profilo == p, setting_disch == "4_Hospice/Palliative Care") %>% arrange(day))$Survival
s_ext_home <- c(1, s_home); s_ext_rehab <- c(1, s_rehab); s_ext_hospital <- c(1, s_hospital); s_ext_hospice <- c(1, s_hospice)
for (t in 1:max_giorni_totali) {
sop_finale$Occ_Home[sop_finale$day == t & sop_finale$Profilo == p] <- sum(i_home[1:t] * rev(s_ext_home[1:t]), na.rm=TRUE)
sop_finale$Occ_Rehab[sop_finale$day == t & sop_finale$Profilo == p] <- sum(i_rehab[1:t] * rev(s_ext_rehab[1:t]), na.rm=TRUE)
sop_finale$Occ_Hospital[sop_finale$day == t & sop_finale$Profilo == p] <- sum(i_hospital[1:t] * rev(s_ext_hospital[1:t]), na.rm=TRUE)
sop_finale$Occ_Hospice[sop_finale$day == t & sop_finale$Profilo == p] <- sum(i_hospice[1:t] * rev(s_ext_hospice[1:t]), na.rm=TRUE)
}
}
sop_finale <- sop_finale %>% mutate(Occ_Death = 1 - (Prob_Stay_Hosp + Occ_Home + Occ_Rehab + Occ_Hospital + Occ_Hospice))
plot_data_sop <- sop_finale %>%
dplyr::select(Profilo, day, Prob_Stay_Hosp, Occ_Home, Occ_Rehab, Occ_Hospital, Occ_Hospice, Occ_Death) %>%
pivot_longer(-c(Profilo, day), names_to = "Stato", values_to = "Probabilita") %>%
mutate(
Stato = factor(Stato, levels = c("Occ_Death", "Occ_Hospice", "Occ_Hospital", "Occ_Rehab", "Occ_Home", "Prob_Stay_Hosp"),
labels = c("Mortalità (Intra+Post)", "Hospice", "Trasf. Ospedale", "Riabilitazione", "A Casa", "In Ospedale"))
)
plot_data_sop <- plot_data_sop %>% mutate(Stato = case_when(Stato == 'In Ospedale'~'1. Still In Hospital',
Stato == 'A Casa'~'2. At Home',
Stato == 'Riabilitazione'~'3. Rehabilitation',
Stato == 'Trasf. Ospedale'~'4. Hospital Transfer',
Stato == 'Hospice'~'5. Hospice',
Stato == 'Mortalità (Intra+Post)'~'6. Mortality (Intra+Post)'))
plot_sop_lines <- ggplot(plot_data_sop, aes(x = day, y = Probabilita, col = Profilo, linetype = Profilo)) +
geom_line(alpha = 0.9, linewidth = 1) +
facet_wrap(~ Stato, ncol = 2, scales = 'free_y') +
ggsci::scale_color_jco() + # JCO palette for maximum discriminability
scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
scale_x_continuous(expand = c(0, 0)) +
labs(subtitle = "65 years old, SCAI C, Lactate levels 4.3, eGFR 60, No MCS in the acute phase",
x = "Days from Index Event", y = "State Occupation Probability", col = "", linetype = "") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom", strip.text = element_text(face = "bold", size = 10))
png('output/fig_SOPS.png', width = 30, height = 30, unit = 'cm', res = 1000)
print(plot_sop_lines)
dev.off()
# ==============================================================================
# 11.4 STATE OCCUPATION PROBABILITIES FOR 4 SPECIFIC PROFILES (OMT & LVEF)
# ==============================================================================
cat("Calculating S.O.P. for the 4 profiles (OMT and LVEF)...\n")
# Reusing definitions from the previous block, but now fixing Sex and Etiology
val_sesso <- "M"
val_ezio <- "AMI"
# The function now doesn't ask for sex and etiology (fixed for all 4 profiles)
crea_profilo2 <- function(nome) {
data.frame(
Profilo = nome,
day = 1:max_giorni_totali,
var_eta_10 = val_eta,
var_sesso = factor(val_sesso, levels = liv_sesso),
var_eziologia = factor(val_ezio, levels = liv_ezio),
var_scai_base = factor(val_scai, levels = liv_scai),
var_egfr_10 = val_egfr,
var_lact_base = val_lact,
tdc_mcs_type = factor(val_mcs, levels = liv_mcs),
tdc_hrt_pre = factor(val_hrt, levels = liv_hrt),
var_center = factor(centro_default, levels = liv_center)
)
}
# Creating the 4 intersections (OMT 4 vs 1 / FE >= 40% vs FE < 40%)
profili_estesi2 <- bind_rows(
crea_profilo2("OMT 4, FE ≥40"),
crea_profilo2("OMT 4, FE <40"),
crea_profilo2("OMT 1, FE ≥40"),
crea_profilo2("OMT 1, FE <40")
)
# Incidence Extraction (In-hospital curves will be identical for all 4)
lista_inc2 <- lapply(list_multinom, function(m) estrai_incidenza(m, profili_estesi2))
inc_pooled2 <- bind_rows(lista_inc2) %>% group_by(Profilo, day) %>% summarise(across(everything(), mean), .groups="drop")
# Post-Discharge Survival Extraction
profili_post_completi2 <- profili_estesi2 %>%
filter(day == 1) %>% dplyr::select(-any_of("setting_disch")) %>%
tidyr::crossing(setting_disch = c("1_Home", "2_Rehab", "3_Hospital", "4_Hospice/Palliative Care")) %>%
mutate(
status_hrt_pre = 0,
var_egfr_10 = val_egfr,
var_diabete = val_diab,
storico_mcs = factor(val_mcs, levels = liv_mcs),
# DYNAMIC ASSIGNMENT BASED ON PROFILE NAME
var_omt_score = ifelse(grepl("OMT 4", Profilo), 4, 1),
var_fe_disch = ifelse(grepl("FE ≥40", Profilo), "≥40", "<40")
) %>%
group_by(Profilo, setting_disch) %>% slice(rep(1, each = max_giorni_totali)) %>% mutate(day = 1:max_giorni_totali) %>% ungroup()
lista_sur2 <- lapply(list_glm_post, function(m) estrai_sopravvivenza(m, profili_post_completi2))
sur_pooled2 <- bind_rows(lista_sur2) %>% group_by(Profilo, setting_disch, day) %>% summarise(Survival = mean(Survival), .groups="drop")
# Convolution
sop_finale2 <- inc_pooled2 %>% arrange(Profilo, day) %>% mutate(Occ_Home = 0, Occ_Rehab = 0, Occ_Hospice = 0, Occ_Hospital = 0)
for (p in unique(sop_finale2$Profilo)) {
temp_p <- sop_finale2 %>% filter(Profilo == p) %>% arrange(day)
i_home <- temp_p$Inc_Home; i_rehab <- temp_p$Inc_Rehab; i_hospital <- temp_p$Inc_Hospital; i_hospice <- temp_p$Inc_Hospice
s_home <- (sur_pooled2 %>% filter(Profilo == p, setting_disch == "1_Home") %>% arrange(day))$Survival
s_rehab <- (sur_pooled2 %>% filter(Profilo == p, setting_disch == "2_Rehab") %>% arrange(day))$Survival
s_hospital <- (sur_pooled2 %>% filter(Profilo == p, setting_disch == "3_Hospital") %>% arrange(day))$Survival
s_hospice <- (sur_pooled2 %>% filter(Profilo == p, setting_disch == "4_Hospice/Palliative Care") %>% arrange(day))$Survival
s_ext_home <- c(1, s_home); s_ext_rehab <- c(1, s_rehab); s_ext_hospital <- c(1, s_hospital); s_ext_hospice <- c(1, s_hospice)
for (t in 1:max_giorni_totali) {
sop_finale2$Occ_Home[sop_finale2$day == t & sop_finale2$Profilo == p] <- sum(i_home[1:t] * rev(s_ext_home[1:t]), na.rm=TRUE)
sop_finale2$Occ_Rehab[sop_finale2$day == t & sop_finale2$Profilo == p] <- sum(i_rehab[1:t] * rev(s_ext_rehab[1:t]), na.rm=TRUE)
sop_finale2$Occ_Hospital[sop_finale2$day == t & sop_finale2$Profilo == p] <- sum(i_hospital[1:t] * rev(s_ext_hospital[1:t]), na.rm=TRUE)
sop_finale2$Occ_Hospice[sop_finale2$day == t & sop_finale2$Profilo == p] <- sum(i_hospice[1:t] * rev(s_ext_hospice[1:t]), na.rm=TRUE)
}
}
sop_finale2 <- sop_finale2 %>% mutate(Occ_Death = 1 - (Prob_Stay_Hosp + Occ_Home + Occ_Rehab + Occ_Hospital + Occ_Hospice))
plot_data_sop2 <- sop_finale2 %>%
dplyr::select(Profilo, day, Prob_Stay_Hosp, Occ_Home, Occ_Rehab, Occ_Hospital, Occ_Hospice, Occ_Death) %>%
pivot_longer(-c(Profilo, day), names_to = "Stato", values_to = "Probabilita") %>%
mutate(
Stato = factor(Stato, levels = c("Occ_Death", "Occ_Hospice", "Occ_Hospital", "Occ_Rehab", "Occ_Home", "Prob_Stay_Hosp"),
labels = c("Mortalità (Intra+Post)", "Hospice", "Trasf. Ospedale", "Riabilitazione", "A Casa", "In Ospedale"))
) %>%
mutate(Stato = case_when(Stato == 'In Ospedale'~'1. Still In Hospital',
Stato == 'A Casa'~'2. At Home',
Stato == 'Riabilitazione'~'3. Rehabilitation',
Stato == 'Trasf. Ospedale'~'4. Hospital Transfer',
Stato == 'Hospice'~'5. Hospice',
Stato == 'Mortalità (Intra+Post)'~'6. Mortality (Intra+Post)'))
plot_sop_lines2 <- ggplot(plot_data_sop2, aes(x = day, y = Probabilita, col = Profilo, linetype = Profilo)) +
geom_line(alpha = 0.9, linewidth = 1) +
facet_wrap(~ Stato, ncol = 2, scales = 'free_y') +
ggsci::scale_color_jco() +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), expand = c(0, 0)) +
scale_x_continuous(expand = c(0, 0)) +
labs(subtitle = "Male patient, 65 years old, AMI, SCAI C: Impact of OMT and FE after discharge",
x = "Days from Index Event", y = "State Occupation Probability", col = "", linetype = "") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom", strip.text = element_text(face = "bold", size = 10))+
coord_cartesian(xlim = c(0,180))
png('output/fig_SOPS2.png', width = 30, height = 30, unit = 'cm', res = 1000)
print(plot_sop_lines2)
dev.off()
# ==============================================================================
# 11.5 GRAPHICAL OUTPUTS
# ==============================================================================
# View the generated plots directly in the RStudio Viewer
print(plot_forest_intra)
print(plot_forest_post)
print(plot_sop_lines)
print(plot_sop_lines2)