TCE analysis
Danilo Lofaro
04 giugno 2019
library(readxl)
library(tidyverse)
library(ggsci)
library(gridExtra)
library(scales)
library(tableone)
library(knitr)
library(kableExtra)
library(ggpubr)
library(cowplot)
library(plotROC)
library(extrafont)
#font_import()
loadfonts(device = "win")
windowsFonts(OpenSans = windowsFont("Open Sans"))
fct_case_when <- function(...) {
args <- as.list(match.call())
levels <- sapply(args[-1], function(f) f[[3]]) # extract RHS of formula
levels <- levels[!is.na(levels)]
factor(dplyr::case_when(...), levels=levels)
}
theme_Publication_void <- function(base_size=14, base_family="sans") {
library(grid)
library(ggthemes)
(theme_foundation(base_size=base_size, base_family=base_family)
+ theme(plot.title = element_text(face = "bold",
size = rel(1.2)),
text = element_text(),
panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour = NA),
axis.title = element_text(size = rel(1)),
axis.title.y = element_text(angle=90,vjust =2),
axis.title.x = element_text(vjust = -0.2),
axis.text = element_text(),
axis.line = element_blank(),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.key = element_rect(colour = NA),
legend.position = "bottom",
legend.direction = "horizontal",
legend.key.size= unit(0.2, "cm"),
legend.margin = unit(0, "cm"),
legend.title = element_text(face="italic"),
plot.margin=margin(10,5,5,5,"mm"),
strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
strip.text = element_text(face="bold")
))
}ds_raw <- read_excel("../data/Data base TCE 2018 TOT 16 MARZO.xlsx", range = "A5:EA106", col_names = F)ds_raw2 <- ds_raw %>%
dplyr::select(-c(34,35,90,91,128,129)) %>%
rename_all(~ c("surname","name","dob","sex", "date_gca","date_adm","primary.brain.injury", "antiaggr_anticoag", "neuroimaging", "cranial_fractures","facial_fractures", "extremities_fractures","toracic_trauma","abdominal_trauma","vertebral_trauma","myelic_trauma","vascular_trauma", "icu_isp", "diagnosis_adm", "crs_adm","lcf_adm","drs_adm","tracheo_adm","breath_adm","feed_adm", "urinary_cateth_adm","bedsore_adm","asworth_adm","poa_adm","craniolac_adm","hydrocef_adm","erbi_sub_a_adm","erbi_sub_b_adm","cirs_cardiac_adm", "cirs_hypertension_adm", "cirs_vascular_adm", "cirs_respiratory_adm", "cirs_visual_otorhinolaryngo_adm", "cirs_gastro_upper_adm", "cirs_gastro_lower_adm", "cirs_hepatic_adm", "cirs_renal_adm", "cirs_genito_urinary_adm", "cirs_muscle_skeletal_skin_adm", "cirs_neurologic_adm", "cirs_endocrine_adm", "cirs_psichiatric_behave_adm", "cirs_cancers_adm", "leukocytes", "hemoglobin", "platelets", "pct", "weight", "bmi", "respiratory_inf_adm", "blood_inf_adm", "uti_adm", "stool_inf_adm", "csf_inf_adm", "brain_inf_adm", "skin_inf_adm", "ear_eye_inf_adm", "mrsa_adm", "pseudomonas_resistant_adm", "acinetobacter_resistant_adm", "vre_adm", "clostridium_difficilis_adm", "klebsiella_esbl_adm", "proteus_mirabilis_adm", "ecoli_ebls_adm", "other_inf_adm", "cirs_cardiac_rehab", "cirs_hypertension_rehab", "cirs_vascular_rehab", "cirs_respiratory_rehab", "cirs_visual_otorhinolaryngo_rehab", "cirs_gastro_upper_rehab", "cirs_gastro_lower_rehab", "cirs_hepatic_rehab", "cirs_renal_rehab", "cirs_genito_urinary_rehab", "cirs_muscle_skeletal_skin_rehab", "cirs_neurologic_rehab", "cirs_endocrine_rehab", "cirs_psichiatric_behave_rehab", "nsh_rehab", "other_compl", "respiratory_inf_rehab", "blood_inf_rehab", "uti_rehab", "stool_inf_rehab", "csf_inf_rehab", "brain_inf_rehab", "skin_inf_rehab", "ear_eye_inf_rehab", "mrsa_rehab", "pseudomonas_resistant_rehab", "acinetobacter_resistant_rehab", "vre_rehab", "clostridium_difficilis_rehab", "klebsiella_esbl_rehab", "proteus_mirabilis_rehab", "ecoli_ebls_rehab", "other_inf_rehab","btx_interv","tracheo_interv", "date_dism", "type_dism", "diagnosis_dism", "crs_dism","lcf_dism","drs_dism","tracheo_dism","breath_dism","feed_dism", "urinary_cateth_dism","bedsore_dism","asworth_dism","poa_dism","craniolac_dism","hydrocef_dism","erbi_sub_a_dism","erbi_sub_b_dism", "therapy_dism", "gos_dism")) %>%
mutate(cirs_psichiatric_behave_adm = str_replace(string = cirs_psichiatric_behave_adm, pattern = "O", replacement = "0"))ds_preproc <- ds_raw2 %>%
mutate_at(.vars = vars(dob, date_gca, date_adm, date_dism),
.funs = list(as.Date.POSIXct)) %>%
mutate_at(.vars = vars(crs_adm, tracheo_adm, cirs_psichiatric_behave_adm, pct, weight, bmi, cirs_respiratory_rehab, tracheo_dism),
.funs = list(as.numeric)) %>%
mutate_at(.vars = vars(lcf_adm,lcf_dism),
.funs = list(~as.numeric(as.roman(.)))) %>%
mutate_at(.vars = vars(breath_adm, breath_dism),
.funs = list(~fct_case_when(
. == 1 ~ "Autonomous",
. == 2 ~ "Autonomous + O2",
. == 3 ~ "Mechanical"))) %>%
mutate_at(.vars = vars(feed_adm, feed_dism),
.funs = list(~fct_case_when(
. == 4 ~ "Oral",
. == 1 ~ "Parenteral",
. == 2 ~ "NG-tube",
. == 3 ~ "PEG"))) %>%
mutate_at(.vars = vars(hydrocef_adm, hydrocef_dism),
.funs = list(~fct_case_when(
. == 0 ~ "No",
. == 1 ~ "Yes",
. == 2 ~ "DVE",
. == 3 ~ "DVP"))) %>%
mutate(id = paste("pz",str_pad(row_number(), 4, pad = "0"),sep ="_"),
primary.brain.injury = case_when(
primary.brain.injury == 1 ~ "CAR STR AUTO GUIDA",
primary.brain.injury == 2 ~ "INC STR AUTO PASSEGGERO",
primary.brain.injury == 3 ~ "INC STR MOTO GUIDA",
primary.brain.injury == 4 ~ "INC STR MOTO PASSEGGERO",
primary.brain.injury == 5 ~ "PEDONE",
primary.brain.injury == 6 ~ "CADUTA ALTO",
primary.brain.injury == 7 ~ "CADUTA ACCIDENTALE",
primary.brain.injury == 8 ~ "FAF",
primary.brain.injury == 9 ~ "CORPO PENETRANTE",
primary.brain.injury == 10 ~ "VIOLENZA",
primary.brain.injury == 11 ~ "SPORT",
primary.brain.injury == 12 ~ "ALTRO"),
urinary_cateth_adm = case_when(
urinary_cateth_adm > 1 ~ 1,
TRUE ~ urinary_cateth_adm),
antiplatelet_adm = fct_case_when(
antiaggr_anticoag == 1 ~ 1,
TRUE ~ 0),
blood_thinner_adm = fct_case_when(
antiaggr_anticoag == 2 ~ 1,
TRUE ~ 0),
antiaggr_anticoag = as.numeric(!antiaggr_anticoag==0),
neuroimaging = str_replace_all(coalesce(str_replace(as.character(as.numeric(neuroimaging)),"\\.",","),neuroimaging)," ",""),
neuroimaging = str_replace_all(neuroimaging, ",,",","),
neuroimaging = str_replace_all(neuroimaging, "(.*)(,$)", "\\1"),
diagnosis_adm = factor(recode(diagnosis_adm, "EMERSIONE" = "Emersion", "EMERS" = "Emersion", "sv" = "VS", "SV" = "VS", "SMC" = "MCS")),
other_compl_cat = as.numeric(!other_compl %in% c(0,NA)),
pct_cat = case_when(
pct < .1 ~ "no infection",
pct >= .1 & pct < .25 ~ "unlikely infection",
pct >= .25 & pct < .5 ~ "possible infection",
TRUE ~ "infection"),
diagnosis_dism = fct_relevel(factor(recode(diagnosis_dism, "EMERSIONE" = "Emersion", "MERSIONE" = "Emersion", "SV" = "VS", "SMC" = "MCS", "DECESSO" = "Death")),"Death", after = 3),
therapy_dism = str_replace_all(coalesce(str_replace(as.character(as.numeric(therapy_dism)),"\\.",","),therapy_dism)," ","")) %>%
select(-c(other_compl))ds_neuroimaging <- ds_preproc %>%
dplyr::select(id, neuroimaging) %>%
separate(neuroimaging, into = paste("neuroimg", 1:5, sep = "_")) %>%
gather(n_img, neuroimg, neuroimg_1:neuroimg_5, na.rm = T) %>%
mutate(neuroimg = fct_case_when(
neuroimg == 1 ~ "danno_assonale_diffuso",
neuroimg == 2 ~ "emorragia_subaracnoidea",
neuroimg == 3 ~ "ematoma_epidurale",
neuroimg == 4 ~ "ematoma_subdurale",
neuroimg == 5 ~ "contusioni_cerebrali",
neuroimg == 6 ~ "emorragia_ventricolare",
neuroimg == 7 ~ "ematoma_intraparenchimale",
neuroimg == 8 ~ "corpo_estraneo",
neuroimg == 9 ~ "frattura_cranica_infossata"),
val = 1) %>%
select(-n_img) %>%
filter(!is.na(neuroimg))
ds_neuroimg_wide <- ds_neuroimaging %>%
#unite(treat, treatment, n, sep = "_") %>%
spread(key = neuroimg, value = val, fill = 0)ds_therapy <- ds_preproc %>%
select(id, therapy_dism) %>%
separate(therapy_dism, into = paste("treat", 1:7, sep = "_")) %>%
gather(n_treat, treatment, treat_1:treat_7, na.rm = T) %>%
mutate(treatment = fct_case_when(
treatment == 0 ~ "nessuno",
treatment == 1 ~ "antiepilettici",
treatment == 2 ~ "antipertensivi",
treatment == 3 ~ "eparina",
treatment == 4 ~ "protez_gastrica",
treatment == 5 ~ "antispastici",
treatment == 6 ~ "antidepressivi",
treatment == 7 ~ "dopaminergici",
treatment == 8 ~ "t_orm_sost",
treatment == 9 ~ "antibiotici",
treatment == 10 ~ "cortisone",
treatment == 11 ~ "antipsicotici",
treatment == 12 ~ "anticoag",
treatment == 13 ~ "other"),
val = 1) %>%
select(-n_treat) %>%
filter(!is.na(treatment))
ds_therapy_wide <- ds_therapy %>%
#unite(treat, treatment, n, sep = "_") %>%
spread(key = treatment, value = val, fill = 0)ds_proc <- ds_preproc %>%
select(-c(neuroimaging, therapy_dism)) %>%
left_join(ds_neuroimg_wide) %>%
left_join(ds_therapy_wide) %>%
mutate(date_adm = case_when(
id == "pz_0032" ~ as.Date("2017-12-04"),
TRUE ~ date_adm),
days_icu = as.numeric(date_adm - date_gca),
days_rehab = as.numeric(date_dism - date_adm),
age_event = round(as.numeric(date_gca - dob)/365.25, 1))ds_tce <- ds_proc %>%
select(id, age_event, sex, days_icu, diagnosis_adm, primary.brain.injury, icu_isp, antiaggr_anticoag:vascular_trauma, crs_adm:other_inf_adm, cirs_cardiac_rehab:cirs_psichiatric_behave_rehab, days_rehab, diagnosis_dism:drs_dism, erbi_sub_a_dism:gos_dism) %>%
select(-c(cirs_cancers_adm, cirs_neurologic_adm, cirs_neurologic_rehab)) %>%
mutate(tot_comorb_adm = reduce(mutate_all(select(., cirs_cardiac_adm:cirs_psichiatric_behave_adm), .funs = (~replace_na(as.numeric(.>1),0))), `+`),
severe_comorb_adm = as.numeric(tot_comorb_adm>2),
tot_comorb_rehab = reduce(mutate_all(select(., c(cirs_cardiac_rehab:cirs_psichiatric_behave_rehab)), .funs = list(~as.numeric(.>2))), `+`),
severe_comorb_rehab = as.numeric(tot_comorb_rehab>2),
death = as.numeric(diagnosis_dism == "Death"),
diagnosis_adm2 = factor(diagnosis_adm, levels = c("Emersion", "MCS", "VS", "Death")),
diagnosis_chng = as.numeric(diagnosis_adm) - as.numeric(diagnosis_dism),
improvement = case_when(
diagnosis_adm == "Emersion" ~ NA_real_,
TRUE ~ as.numeric(diagnosis_chng>0)),
age_cat = cut(age_event, breaks = c(0,21,40,60,100), c("<=21","22-39","40-59", ">=60")),
nr_comorbidity_adm = factor(severe_comorb_adm, 0:1, c("<2",">=2")),
nr_comorbidity_rehab = factor(severe_comorb_rehab, 0:1, c("<2",">=2")),
primary_injury_cat = fct_case_when(
primary.brain.injury == "CADUTA ACCIDENTALE" ~ "Accidental fall",
primary.brain.injury == "CADUTA ALTO" ~ "High-level fall",
primary.brain.injury %in% c("INC STR AUTO GUIDA", "INC STR AUTO PASSEGGERO") ~ "Car accident",
primary.brain.injury %in% c("INC STR MOTO GUIDA", "INC STR MOTO PASSEGGERO") ~ "Motorcycle accident",
primary.brain.injury == "PEDONE" ~ "Pedestrian accident",
TRUE ~ "Other"),
drs_delta = drs_dism - drs_adm,
drs_adm_cat = cut(drs_adm,breaks = c(-1,0,1,3,6,11,16,21,24,29,31), labels = c("None", "Mild", "Partial", "Moderate", "Moderately severe", "Severe", "Extremely severe", "VS", "Severe VS", "Death")),
drs_dism_cat = cut(drs_dism,breaks = c(-1,0,1,3,6,11,16,21,24,29,31), labels = c("None", "Mild", "Partial", "Moderate", "Moderately severe", "Severe", "Extremely severe", "VS", "Severe VS", "Death")),
drs_adm_nr = as.numeric(drs_adm_cat),
drs_dism_nr = as.numeric(drs_dism_cat),
funct_impr = factor(drs_dism_nr < drs_adm_nr, c("FALSE","TRUE"), c("no","yes")),
funct_impr_num = as.numeric(funct_impr)-1,
funct_worse = factor(drs_dism_nr < drs_adm_nr, c("FALSE","TRUE"), c("yes","no")),
funct_worse_num = (funct_impr_num-1)*-1) %>%
filter(days_icu<88)ds_cause <- ds_tce %>%
group_by(primary_injury_cat, .drop = FALSE) %>%
summarise(n = n()) %>%
mutate(freq = n / sum(n))ds_trauma <- ds_tce %>%
select(cranial_fractures:vascular_trauma, age_cat) %>%
group_by(age_cat) %>%
#mutate(n = n()) %>%
summarise_all(.funs = funs(tot = sum, n = n())) %>%
rename(n = cranial_fractures_n) %>%
rename_all(~ str_remove(., "_tot")) %>%
select(-c(facial_fractures_n:vascular_trauma_n)) %>%
gather(trauma, total, -c(age_cat,n)) %>%
mutate(trauma = str_to_title(str_replace(trauma, "_", " ")),
trauma = factor(trauma, levels = unique(trauma)),
freq=total/n)
# ds_fractures <- ds_trauma %>%
# filter(str_detect(trauma, "Fract")) %>%
# mutate(trauma = factor(str_remove(trauma, " Fractures")))
ds_trauma <- ds_trauma %>%
filter(str_detect(trauma, "Trauma")) %>%
mutate(trauma = factor(str_remove(trauma, " Trauma")))ds_neuroimaging_plot <- ds_neuroimaging %>%
inner_join(ds_tce %>% select(id,age_cat)) %>%
mutate(neuroimg = str_to_title(str_replace_all(neuroimg,"_"," ")),
neuroimg = factor(neuroimg, unique(neuroimg), c("Other", "SAH", "DAI", "Cerebral Contusions", "SDH", "Other", "Other", "Other", "Other"))) %>%
mutate(neuroimg = fct_relevel(neuroimg, "DAI", "Cerebral Contusions", "SAH", "SDH", "Other")) %>%
group_by(neuroimg) %>%
summarise(tot = sum(val)) %>%
#inner_join(ds_trauma %>% select(age_cat, n) %>% distinct()) %>%
mutate(freq = tot/96)ds_cirs_long <- ds_preproc %>%
select_at(vars(c(id,contains("cirs")))) %>%
gather(comorbidity, cirs, -id) %>%
separate(comorbidity, into=c("comorbidity","time"), sep = "_(?=[^_]+$)") %>%
filter(!comorbidity %in% c("cirs_neurologic", "cirs_cancers")) %>%
mutate(comorbidity = str_remove(comorbidity, "cirs_"),
comorbidity = factor(comorbidity,levels = c("cardiac", "hypertension", "vascular", "respiratory", "visual_otorhinolaryngo", "gastro_upper", "gastro_lower", "hepatic", "renal", "genito_urinary", "muscle_skeletal_skin","endocrine", "psichiatric_behave"), labels = c("Cardiac", "Hypertension", "Vascular\nHematopoietic", "Respiratory", "EENTL", "Upper GI", "Lower GI", "Hepatic", "Renal", "Genitourinary", "Musculoskeletal\nSkin","Endocrine\nBreast", "Psychiatric\nillness")),
val = as.numeric(cirs>0),
time = factor(time, labels=c("Admission", "Rehab"))) %>%
inner_join(ds_tce %>% select(id,age_cat))
ds_cirs_plot <- ds_cirs_long %>%
group_by(time, comorbidity, .drop = FALSE) %>%
summarise(count = sum(val), freq = sum(val)/n()) %>%
replace_na(list(count = 0, freq = 0))ds_erbi <- ds_tce %>%
select(id, age_cat,erbi_sub_a_adm, erbi_sub_a_dism, erbi_sub_b_adm, erbi_sub_b_dism) %>%
rename_all(~ str_remove(string = ., pattern = "erbi_sub_")) %>%
gather(time, erbi, a_adm:b_dism) %>%
separate(time, c("Subscale", "time"), sep = "_") %>%
mutate(Subscale = factor(Subscale, labels=c("A","B")),
time = factor(time, labels=c("Admission", "Dismission")))ds_drs <- ds_tce %>%
select(id,age_cat, nr_comorbidity_adm, nr_comorbidity_rehab, drs_adm, drs_dism) %>%
mutate(drs_delta = drs_dism - drs_adm) %>%
rename_all(~ str_remove(string = ., pattern = "drs_")) %>%
gather(time, drs, adm:dism) %>%
mutate(time = factor(time, c("adm","dism"),c("Admission", "Dismission")),
drs_cat = cut(drs,breaks = c(-1,0,1,3,6,11,16,21,24,29,31), labels = c("None", "Mild", "Partial", "Moderate", "Moderately severe", "Severe", "Extremely severe", "VS", "Severe VS", "Death")))Results
vars_tab1 <- c("age_event", "sex", "days_icu", "primary_injury_cat", "cranial_fractures", "facial_fractures", "extremities_fractures", "toracic_trauma", "abdominal_trauma", "vertebral_trauma", "myelic_trauma", "vascular_trauma", "days_rehab", "tracheo_adm", "breath_adm", "feed_adm", "urinary_cateth_adm", "bedsore_adm", "craniolac_adm")
ds_tab1 <- as.data.frame(print(CreateTableOne(vars = vars_tab1, data = ds_tce, factorVars = c("sex", "primary_injury_cat", "cranial_fractures", "facial_fractures", "extremities_fractures", "toracic_trauma", "abdominal_trauma", "vertebral_trauma", "myelic_trauma", "vascular_trauma","tracheo_adm", "breath_adm", "feed_adm", "urinary_cateth_adm", "bedsore_adm", "craniolac_adm")), catDigits = 1, contDigits = 1, noSpaces = T, missing = F, test = F, printToggle = F, nonnormal = c("days_icu", "days_rehab")))
#ds_tab11 <- as.data.frame(print(CreateTableOne(vars = vars_tab1, data = ds_tce, factorVars = c("sex", "diagnosis_adm", "tracheo_adm", "breath_adm", "feed_adm", "urinary_cateth_adm", "bedsore_adm", "poa_adm", "craniolac_adm"), strata = "age_cat"), catDigits = 1, contDigits = 1, noSpaces = T, missing = F, test = F, printToggle = F, nonnormal = c("days_icu", "crs_adm", "drs_adm", "days_rehab")))
ds_tab1 <- ds_tab1 %>%
rownames_to_column(var = "vars") %>%
mutate(vars = str_trim(vars))
kable(ds_tab1, caption = "Table 1. Charasteristics of patients at admission in rehab centers") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
add_indent(c(8:10, 12:14))| vars | Overall |
|---|---|
| n | 96 |
| age_event (mean (SD)) | 43.6 (20.8) |
| sex = M (%) | 69 (71.9) |
| days_icu (median [IQR]) | 27.0 [20.0, 35.0] |
| primary_injury_cat (%) | |
| Accidental fall | 17 (17.7) |
| High-level fall | 14 (14.6) |
| Car accident | 13 (13.5) |
| Motorcycle accident | 20 (20.8) |
| Pedestrian accident | 7 (7.3) |
| Other | 25 (26.0) |
| cranial_fractures (%) | |
| 0 | 58 (60.4) |
| 1 | 37 (38.5) |
| 5 | 1 (1.0) |
| facial_fractures = 1 (%) | 45 (46.9) |
| extremities_fractures = 1 (%) | 43 (44.8) |
| toracic_trauma = 1 (%) | 49 (51.0) |
| abdominal_trauma = 1 (%) | 17 (17.7) |
| vertebral_trauma = 1 (%) | 23 (24.0) |
| myelic_trauma = 1 (%) | 5 (5.2) |
| vascular_trauma = 1 (%) | 3 (3.1) |
| days_rehab (median [IQR]) | 72.0 [43.8, 128.2] |
| tracheo_adm = 1 (%) | 59 (61.5) |
| breath_adm (%) | |
| Autonomous | 75 (78.9) |
| Autonomous + O2 | 18 (18.9) |
| Mechanical | 2 (2.1) |
| feed_adm (%) | |
| Oral | 36 (37.9) |
| NG-tube | 41 (43.2) |
| PEG | 18 (18.9) |
| urinary_cateth_adm = 1 (%) | 91 (94.8) |
| bedsore_adm = 1 (%) | 31 (32.3) |
| craniolac_adm = 1 (%) | 23 (24.0) |
vars_tab2 <- c("diagnosis_adm", "crs_adm", "drs_adm", "tot_comorb_adm", "diagnosis_dism", "crs_dism", "drs_dism", "tot_comorb_rehab")
ds_tab2 <- ds_tce %>%
select(id, vars_tab2) %>%
rename(totComorb_adm = tot_comorb_adm,
totComorb_dism = tot_comorb_rehab) %>%
gather(key, value, -id) %>%
separate(col = key, into = c("var","time"), sep = "_") %>%
spread(var, value, fill = 0) %>%
mutate_at(.vars = vars(crs, drs, totComorb),
.funs = list(as.numeric))
ds_tab22 <- as.data.frame(print(CreateTableOne(vars = c("diagnosis", "crs", "drs", "totComorb"), data = ds_tab2, factorVars = "diagnosis", strata = "time"), catDigits = 1, contDigits = 1, noSpaces = T, missing = F, test = T, printToggle = F, nonnormal = c("crs", "drs", "totComorb")))
ds_tab22 <- ds_tab22 %>%
rownames_to_column(var = "vars") %>%
mutate(vars = str_trim(vars)) %>%
rename(Admission = adm,
Dismission = dism) %>%
select(-test)
kable(ds_tab22, caption = "Table 2. Diagnosis, CRS, DRS and nr. of CIRS comorbidity at admission and dismission from rehab") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
add_indent(c(3:6))| vars | Admission | Dismission | p |
|---|---|---|---|
| n | 96 | 96 | |
| diagnosis (%) | <0.001 | ||
| Death | 0 (0.0) | 11 (11.5) | |
| Emersion | 55 (57.3) | 73 (76.0) | |
| MCS | 21 (21.9) | 8 (8.3) | |
| VS | 20 (20.8) | 4 (4.2) | |
| crs (median [IQR]) | 23.0 [7.0, 23.0] | 23.0 [21.0, 23.0] | 0.014 |
| drs (median [IQR]) | 18.0 [16.0, 21.0] | 9.0 [5.0, 15.2] | <0.001 |
| totComorb (median [IQR]) | 5.0 [3.8, 6.0] | 2.0 [1.0, 3.0] | <0.001 |
g <- ggplot(ds_cause, aes(x = primary_injury_cat, y = freq))
g + geom_bar(stat = "identity", colour="black", fill = "grey50") +
labs(x="TBI cause", y="Frequency") +
#scale_fill_grey(name = "Age (years)", start = 0.5, end = 1) +
#geom_text(stat = "count", aes(label=..count.., y = (..count..)/sum(..count..)), vjust = .25, color = "black", fontface = "bold", size = 5, position = position_dodge2()) +
scale_y_continuous(breaks = seq(0,.25,.05), labels = scales::percent(seq(0,.25,.05))) +
scale_x_discrete(labels = str_replace(levels(ds_cause$primary_injury_cat), " ", "\n")) +
#facet_wrap(~ age_cat) +
theme_bw(base_family="OpenSans") +
theme(legend.position = "bottom",
text = element_text(size=16),
axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 0, b = 0, l = 0)),
axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 0)))
Figure 1. Mechanisms of Traumatic Brain Injury (TBI)
g <- ggplot(ds_neuroimaging_plot, aes(x = neuroimg, y = freq))
g + geom_bar(stat="identity", colour="black", fill = "grey50") +
#scale_fill_grey(name = "Age strata", start = 0, end = 1) +
#scale_y_continuous(labels = scales::percent) +
labs(x="Primary injury", y = "Frequency") +
scale_y_continuous(breaks = seq(0,.75,.15), labels = scales::percent(seq(0,.75,.15))) +
#scale_x_discrete(labels = str_replace(levels(ds_neuroimaging_plot$neuroimg), " ", "\n")) +
#facet_wrap(~ age_cat) +
theme_bw(base_family="OpenSans") +
theme(legend.position = "bottom",
text = element_text(size=16),
axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 0, b = 0, l = 0)),
axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 0)))
Figure 2. Lesions revealed at imaging
g <- ggplot(data=ds_cirs_plot,aes(x=comorbidity, y=freq, fill=time))
g + geom_bar(data = ds_cirs_plot %>% filter(time=="Rehab"), stat = "identity", colour = "black") +
geom_bar(data = ds_cirs_plot %>% filter(time=="Admission"), aes(y=freq*(-1)), stat = "identity", colour = "black") +
scale_y_continuous(limits = c(-1,1), breaks=seq(-1,1,.25), labels = paste0(abs(seq(-100,100,25)),"%")) +
scale_x_discrete(limits=levels(fct_rev(ds_cirs_plot$comorbidity))) +
scale_fill_grey(name = "Time") +
labs(x = "Comorbidity", y = "Frequency") +
coord_flip() +
#facet_wrap(~ age_cat) +
theme_bw(base_family="OpenSans") +
theme(legend.position = "bottom",
text = element_text(size=16),
axis.text.x = element_text(size=14),
axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 0, b = 0, l = 0)),
axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 0)))
Figure 3. CIRS measured comorbidity rate at admission and after rehabilitation
g1 <- ggplot(ds_tce, aes(x=age_event, y=drs_delta)) +
geom_point(shape = 21, colour = "black", fill = "grey50", size = 5, alpha = .8) +
geom_smooth(method = "lm", colour="black", formula = y ~ x) +
labs(x="Age at TBI (years)", y = "DRS Delta", title = "A") +
#scale_x_continuous(breaks=seq(12.5,87.5,12.5)) +
theme_bw(base_family="OpenSans") +
theme(text = element_text(size=16),
plot.title = element_text(face="bold", size=24),
axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 0, b = 0, l = 0)),
axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 0)))
g2 <- ggplot(ds_tce, aes(x=tot_comorb_rehab, y=drs_delta)) +
geom_point(shape = 21, colour = "black", fill = "grey50", size = 5, alpha = .8) +
geom_smooth(method = "lm", colour="black", formula = y ~ x) +
labs(x="nr. Comorbidity in rehab", y = "", title = "B") +
lims(x=c(0,6)) +
theme_bw(base_family="OpenSans") +
theme(text = element_text(size=16),
plot.title = element_text(face="bold", size=24, hjust = 1),
axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 0, b = 0, l = 0)))
grid.arrange(g1 , g2, nrow=1)
Figure 4. Admission - Dismission DRS delta vs Age at TBI (A) and nr. of severe comorbidities (CIRS >= 2; B)
g1 <- ggplot(ds_tce, aes(y=age_event, x=drs_delta)) +
geom_point(shape = 21, colour = "black", fill = "grey50", size = 5, alpha = .8) +
geom_smooth(method = "lm", colour="black", formula = y ~ x) +
labs(y="Age at TBI (years)", x = "", title = "A") +
#scale_x_continuous(breaks=seq(12.5,87.5,12.5)) +
theme_bw(base_family="OpenSans") +
theme(text = element_text(size=16),
plot.title = element_text(face="bold", size=24),
axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 0)))
g2 <- ggplot(ds_tce, aes(y=tot_comorb_rehab, x=drs_delta)) +
geom_point(shape = 21, colour = "black", fill = "grey50", size = 5, alpha = .8) +
geom_smooth(method = "lm", colour="black", formula = y ~ x) +
labs(y="nr. Comorbidity in rehab", x = "DRS Delta", title = "B") +
#lims(x=c(0,6)) +
theme_bw(base_family="OpenSans") +
theme(text = element_text(size=16),
plot.title = element_text(face="bold", size=24),
axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 0, b = 0, l = 0)),
axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 0)))
grid.arrange(g1 , g2, nrow=2)
Figure 4. Admission - Dismission DRS delta vs Age at TBI (A) and nr. of severe comorbidities (CIRS >= 2; B)
ds_tce_out1 <- ds_tce %>%
filter(!is.na(funct_impr)) %>%
#mutate(funct_impr_neg = (funct_impr_num*-1)+1) %>%
dplyr::select(age_event:diagnosis_adm, primary_injury_cat, cranial_fractures:crs_adm, drs_adm, breath_adm, feed_adm, cirs_cardiac_rehab:cirs_psichiatric_behave_rehab, funct_worse)ds_tce_logit <- ds_tce %>%
filter(!is.na(funct_impr)) %>%
dplyr::select(age_event, age_cat, sex, diagnosis_adm, severe_comorb_adm, severe_comorb_rehab, tot_comorb_adm, tot_comorb_rehab, funct_worse, funct_worse_num) %>%
mutate(age_biv = cut(age_event, breaks = c(0,60,100), c("0","1")),
age_biv2 = cut(age_event, breaks = c(0,45,100), c("0","1")))library(epicalc)
fit <- glm(funct_worse_num ~ age_event + sex + diagnosis_adm + tot_comorb_rehab, family = "binomial", data = ds_tce_logit)
tab2 <- logistic.display(fit, simplified = T)
ds_tab2 <- as.data.frame(round(tab2$table,3)) %>%
rownames_to_column(var = "vars") %>%
mutate(CI95 = paste(lower95ci,"-",upper95ci)) %>%
dplyr::rename(p = `Pr(>|Z|)`) %>%
dplyr::select(vars, OR, CI95, p)
kable(ds_tab2, caption = "Table 3. Logistic model for association with risk of failed functional improvement") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| vars | OR | CI95 | p |
|---|---|---|---|
| age_event | 1.026 | 0.994 - 1.059 | 0.118 |
| sexM | 0.955 | 0.265 - 3.442 | 0.944 |
| diagnosis_admMCS | 2.460 | 0.558 - 10.855 | 0.235 |
| diagnosis_admVS | 2.141 | 0.555 - 8.261 | 0.269 |
| tot_comorb_rehab | 1.663 | 1.089 - 2.538 | 0.018 |
detach("package:epicalc", unload=TRUE)
#detach("package:MASS", unload=TRUE)library(caret)
set.seed(2406)
class.control <- trainControl(method = "repeatedcv", number = 5, repeats = 5, summaryFunction = twoClassSummary, classProbs = T, verboseIter = F, savePredictions = TRUE)rf.TCE <- train(funct_worse ~ ., data=ds_tce_out1, method = "rf", trControl = class.control, tuneGrid = data.frame(mtry=c(10,20,30)), preProc = c("center", "scale"), na.action = na.omit)lasso.TCE = train(funct_worse ~ ., data=ds_tce_out1, method = "glmnet", trControl = class.control,tuneGrid = expand.grid(alpha=c(.2,.3), lambda=seq(.01,.03,.01)), preProc = c("center", "scale"), na.action = na.omit)svmP.TCE <- train(funct_worse ~ ., data=ds_tce_out1, method = "svmPoly", trControl = class.control, preProc = c("center", "scale"), tuneGrid = expand.grid(C=c(1,3), scale=c(1,10), degree=2:3), na.action = na.omit)sI.rf <- as_tibble(rf.TCE$pred) %>%
filter(mtry == rf.TCE$bestTune$mtry) %>%
dplyr::select(-mtry) %>%
mutate(model = "RF")
sI.svm <- as_tibble(svmP.TCE$pred) %>%
filter(C == svmP.TCE$bestTune$C,
degree == svmP.TCE$bestTune$degree,
scale == svmP.TCE$bestTune$scale) %>%
dplyr::select(-c(C,scale,degree)) %>%
mutate(model = "SVM Polynomial")
sI.lasso <- as_tibble(lasso.TCE$pred) %>%
filter(alpha == lasso.TCE$bestTune$alpha,
lambda == lasso.TCE$bestTune$lambda) %>%
dplyr::select(-c(alpha,lambda)) %>%
mutate(model = "Lasso")
ds_roc <- sI.rf %>%
bind_rows(sI.svm,sI.lasso) %>%
dplyr::select(pred,obs, yes, model)g <- ggplot(ds_roc, aes(m=yes, d=factor(obs, levels = c("no", "yes")), colour = model, linetype = model))
g + geom_roc(n.cuts=0, size = 1.5) +
scale_linetype_manual(name = "Model", values=c("solid", "solid", "dotted"), guide = FALSE) +
scale_colour_manual(name = "Model", values = c("black", "grey60", "black")) +
coord_equal() +
style_roc() +
guides(colour = guide_legend(override.aes = list(linetype = c("solid", "solid", "dotted")))) +
theme(legend.position = "bottom",
text = element_text(size=16, family = "OpenSans"),
axis.title.x = element_text(face = "bold", margin = margin(t = 20, r = 0, b = 0, l = 0)),
axis.title.y = element_text(face = "bold", margin = margin(t = 0, r = 20, b = 0, l = 0)))
Figure 5. ROC curves for Machine Learning models tested: Random Forest (grey), Lasso regression (black solid) and SVM with Polynomial kernel (black dotted)
library(pROC)
roc_rf <- roc(sI.rf$obs, sI.rf$yes, ci = T)
rf_coords <- coords(roc_rf, "best", best.method="youden", ret=c("sens", "spec", "ppv", "npv"))
rf_coords <- round(c(roc_rf$auc, roc_rf$ci, rf_coords)[-3],3)
names(rf_coords) <- c("auc", "lowCI", "upCI", "sensitivity", "specificity", "ppv", "npv")
roc_lasso <- roc(sI.lasso$obs, sI.lasso$yes, ci = T)
lasso_coords <- coords(roc_lasso, "best", best.method="youden", ret=c("sens", "spec", "ppv", "npv"))
lasso_coords <- round(c(roc_lasso$auc, roc_lasso$ci, lasso_coords)[-3],3)
names(lasso_coords) <- c("auc", "lowCI", "upCI", "sensitivity", "specificity", "ppv", "npv")
roc_svm <- roc(sI.svm$obs, sI.svm$yes, ci = T)
svm_coords <- coords(roc_svm, "best", best.method="youden", ret=c("sens", "spec", "ppv", "npv"))
svm_coords <- round(c(roc_svm$auc, roc_svm$ci, svm_coords)[-3],3)
names(svm_coords) <- c("auc", "lowCI", "upCI", "sensitivity", "specificity", "ppv", "npv")
ds_auc <- rf_coords %>%
bind_rows(lasso_coords, svm_coords) %>%
mutate(models = c("RF", "Lasso", "SVM Poly"),
CI = paste(lowCI,"-",upCI)) %>%
dplyr::select(models, auc, CI, sensitivity, specificity, ppv, npv)
kable(ds_auc, caption = "Table 4. Machine Learning models predictive performance") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| models | auc | CI | sensitivity | specificity | ppv | npv |
|---|---|---|---|---|---|---|
| RF | 0.861 | 0.821 - 0.902 | 0.684 | 0.867 | 0.952 | 0.417 |
| Lasso | 0.819 | 0.765 - 0.873 | 0.696 | 0.867 | 0.952 | 0.426 |
| SVM Poly | 0.761 | 0.698 - 0.824 | 0.872 | 0.567 | 0.885 | 0.537 |
rf.VI <- data.frame(varImp(rf.TCE, scale=FALSE)$importance) %>%
rownames_to_column() %>%
as_tibble() %>%
arrange(-Overall) %>%
rename(variable = rowname) %>%
mutate(variable = str_replace_all(str_remove(variable, "_rehab"),"_"," "),
variable = str_replace(variable, "primary injury cat", "TBI cause: "),
variable = str_replace(variable, "breath adm", "breathing: "),
variable = str_replace(variable, "feed adm", "feeding: "),
variable = str_remove(variable, "diagnosis adm"),
variable = str_to_title(variable),
variable = str_replace(variable, "Ng","NG"),
variable = str_replace(variable, "Tbi","TBI"),
variable = str_replace(variable, "Adm","at admission"),
variable = str_replace(variable, "Event","at event"),
variable = str_replace(variable, "Behave","illness"),
variable = str_replace(variable, "Icu","in ICU"),
variable = str_replace(variable, "Crs","CRS"),
variable = str_replace(variable, "Drs","DRS"),
variable = str_replace(variable, "Smc","MCS"),
variable = str_replace(variable, "Gastro Upper","Upper GI"),
variable = str_replace(variable, "Gastro Lower","Lower GI"),
variable = str_replace(variable, "Genito Urinary","Genitourinary"),
variable = str_replace(variable, "Muscle Skeletal Skin","Musculoskeletal and skin"),
variable = str_replace(variable, "Visual Otorhinolaryngo", "EENTL"),
type = factor(str_detect(variable, "Cirs"),c("TRUE","FALSE"),c("CIRS Comorbidity","Baseline variable")),
variable = str_replace(variable, "Cirs ","")
) %>%
top_n(wt = Overall, n = 20)rf.VI %>%
arrange(Overall) %>%
mutate(variable=factor(variable,variable),
yend = Overall - .1) %>%
rename(x = variable, y = Overall) -> rf.VI2
g <- ggplot(data = rf.VI2, aes(x=x, y=y, shape = type))
g + geom_segment( aes(x=x, xend=x, y=0, yend=y-.0775), color="grey90", size=2) +
geom_point(size=9, stroke=2, color = "grey20") +
geom_text(aes(x = x, label = round(y,1)), color = "black", fontface = "bold", size=4) +
scale_shape_manual(name = "Variable Type", values = c(0,1)) +
theme_pubr(legend = "bottom", base_size = 16, base_family = "OpenSans") +
coord_flip() +
theme(
panel.grid.major.y = element_blank(),
panel.border = element_blank(),
axis.ticks.y = element_blank(),
legend.text = element_text(size = 14)
) +
xlab("") +
ylab("Mean Decrease Gini")
# g <- ggdotchart(rf.VI, x = "variable", y = "Overall",
# color = "type", # Color by groups
# #shape = 21,
# stroke = 2,
# #dot.size = 3,
# ylab = "Mean Decrease Gini",
# xlab ="",
# palette = c("grey10", "grey50"), # Custom color palette
# sorting = "descending", # Sort value in descending order
# add = "segments",
# add.params = list(color = "grey90", size = 2),
# rotate = TRUE, # Rotate vertically
# # group = "cyl", # Order by groups
# size = 10, # Large dot size
# label = round(rf.VI$Overall,1), # Add mpg values as dot labels
# font.label = list(color = "black", size = 12,
# vjust = .4, face = "bold"), # Add segments from y = 0 to dots
# legend.title = "Variable Type",
# #font.legend = c(12, "plain", "black"),
# ggtheme = theme_pubr(legend = "bottom", base_size = 16, base_family = "OpenSans") # ggplot2 theme
# )
#
# g + theme(legend.text = element_text(size = 14))