library(tidyverse)── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.4.0 ✔ purrr 1.0.1
✔ tibble 3.2.1 ✔ dplyr 1.1.2
✔ tidyr 1.3.0 ✔ stringr 1.5.0
✔ readr 2.1.3 ✔ forcats 0.5.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()library(survival)
library(patchwork)
library(viridis)Loading required package: viridisLiteLiving in an MSA seeks to be consistently associated with lower negative effects of having a diagnosis for one of the ALS risk factors. One possible cause of this may be people living in MSAs being diagnosed with ALS earlier on in their disease course. One way to look at this would be to check the age at first diagnosis of ALS.
if (Sys.info()["sysname"] == "Darwin") {
model_data <- read_rds("/Volumes/lss_jsimmeri/als/model_data.rds")
} else {
model_data <- read_rds("/Shared/lss_jsimmeri/als/model_data.rds")
}
model_data <- model_data |>
mutate(
age = first_year - dobyr,
female = sex == 1,
stratum = as.numeric(as.factor(glue::glue("{first_year}-{last_year}-{dobyr}-{sex}")))
) |>
mutate(any_sx = motor | bulbar,
n_groups = motor + bulbar,
n_bulbar = speech + swallowing,
n_motor = strength + gait + involuntary_movement + pain + other + falls,
n_total = speech + swallowing + strength + gait + involuntary_movement + pain + other + falls)Living in an MSA is associated with an younger age at diagnosis:
panel_a <- survfit(
Surv(age, als) ~ lives_in_msa,
data = model_data
) |>
broom::tidy() |>
mutate(
strata = ifelse(strata == "lives_in_msa=TRUE", "Urban Area", "Non-Urban Area")
) |>
ggplot(
aes(x = time, y = estimate,
ymin = conf.low, ymax = conf.high,
color = forcats::fct_rev(strata), fill = forcats::fct_rev(strata))
) +
geom_step() +
theme_minimal() +
labs(x = "Age in Years", y = "ALS-Free Survival", color = "", fill = "") +
theme(legend.position = c(0.25, 0.25))
panel_a
Likewise, having symptoms lowers age:
panel_b <- survfit(
Surv(age, als) ~ any_sx,
data = model_data
) |>
broom::tidy() |>
mutate(
strata = ifelse(strata == "any_sx=TRUE", "Has Symptoms", "Does Not Have Symptoms")
) |>
ggplot(
aes(x = time, y = estimate,
ymin = conf.low, ymax = conf.high,
color = forcats::fct_rev(strata), fill = forcats::fct_rev(strata))
) +
geom_step() +
theme_minimal() +
labs(x = "Age in Years", y = "ALS-Free Survival", color = "", fill = "") +
theme(legend.position = c(0.25, 0.25))
panel_b
Which is interactive between MSA and symptoms:
panel_d <- survfit(
Surv(age, als) ~ any_sx + lives_in_msa,
data = model_data
) |>
broom::tidy() |>
mutate(
msa = ifelse(stringr::str_detect(strata, "msa=TRUE"), "Urban Area", "Non-Urban Area"),
sx = ifelse(stringr::str_detect(strata, "sx=TRUE"), "Has Symptoms", "Does Not Have Symptoms")
) |>
ggplot(
aes(x = time, y = estimate,
ymin = conf.low, ymax = conf.high,
color = sx, fill = sx, linetype = msa, group = paste0(sx, msa))
) +
geom_step() +
theme_minimal() +
labs(x = "Age in Years", y = "ALS-Free Survival", color = "", fill = "") +
theme(legend.position = c(0.25, 0.25))
panel_d
And there is a dose response pattern:
panel_c <- survfit(
Surv(age, als) ~ n_total_topcoded,
data = model_data |>
mutate(
n_total_topcoded = case_when(
n_total <= 4 ~ as.character(n_total),
n_total > 4 ~ "5+"
)
)
) |>
broom::tidy() |>
mutate(
strata = stringr::str_extract(strata, "[0-8]"),
strata = ifelse(strata == "5", glue::glue("{strata}+"), strata)
) |>
ggplot(
aes(x = time, y = estimate,
ymin = conf.low, ymax = conf.high,
color = strata, fill = strata)
) +
geom_step() +
theme_minimal() +
labs(x = "Age in Years", y = "ALS-Free Survival", color = "", fill = "") +
theme(legend.position = c(0.25, 0.25))
panel_c
As a single graph:
(panel_a | panel_b) / (panel_c | panel_d) +
plot_annotation(tag_levels = c("A", "B", "C", "D"))
ggsave("~/projects/als-sx/surv_fig.svg", width = 8, height = 8)Testing whether these are multiplicative/addictive or have overlap:
coxph(
Surv(age, als) ~ rate_outpatient_visits + intensity_outpatient_visits_dx +
rate_inpatient_stays + intensity_inpatient_stays_los +
elix_CHF + elix_Valvular + elix_PHTN + elix_PVD + elix_HTN +
elix_HTNcx + elix_Paralysis + elix_NeuroOther +
elix_Pulmonary + elix_DM + elix_DMcx + elix_Hypothyroid +
elix_Renal + elix_Liver + elix_PUD + elix_HIV +
elix_Lymphoma + elix_Mets + elix_Tumor + elix_Rheumatic +
elix_Coagulopathy + elix_Obesity + elix_WeightLoss +
elix_FluidsLytes + elix_BloodLoss + elix_Anemia +
elix_Alcohol + elix_Drugs + elix_Psychoses +
elix_Depression +
as.factor(any_sx) * lives_in_msa +
strata(stratum),
data = model_data,
robust = TRUE,
cluster = stratum
) |>
broom::tidy(exponentiate = TRUE, conf.int = TRUE) |>
select(term, estimate, conf.low, conf.high) |>
knitr::kable(digits = 2)| term | estimate | conf.low | conf.high |
|---|---|---|---|
| rate_outpatient_visits | 1.02 | 1.01 | 1.02 |
| intensity_outpatient_visits_dx | 1.12 | 1.09 | 1.16 |
| rate_inpatient_stays | 0.82 | 0.78 | 0.86 |
| intensity_inpatient_stays_los | 0.99 | 0.99 | 1.00 |
| elix_CHFTRUE | 0.69 | 0.65 | 0.75 |
| elix_ValvularTRUE | 0.98 | 0.93 | 1.03 |
| elix_PHTNTRUE | 0.99 | 0.91 | 1.09 |
| elix_PVDTRUE | 0.94 | 0.89 | 0.99 |
| elix_HTNTRUE | 0.99 | 0.95 | 1.03 |
| elix_HTNcxTRUE | 0.92 | 0.87 | 0.98 |
| elix_ParalysisTRUE | 2.17 | 2.05 | 2.30 |
| elix_NeuroOtherTRUE | 4.08 | 3.93 | 4.24 |
| elix_PulmonaryTRUE | 0.98 | 0.94 | 1.01 |
| elix_DMTRUE | 0.88 | 0.84 | 0.93 |
| elix_DMcxTRUE | 0.84 | 0.78 | 0.90 |
| elix_HypothyroidTRUE | 1.10 | 1.05 | 1.15 |
| elix_RenalTRUE | 0.65 | 0.59 | 0.70 |
| elix_LiverTRUE | 0.96 | 0.90 | 1.03 |
| elix_PUDTRUE | 1.04 | 0.84 | 1.28 |
| elix_HIVTRUE | 0.99 | 0.71 | 1.38 |
| elix_LymphomaTRUE | 1.18 | 1.05 | 1.33 |
| elix_MetsTRUE | 0.39 | 0.34 | 0.44 |
| elix_TumorTRUE | 0.90 | 0.85 | 0.95 |
| elix_RheumaticTRUE | 1.22 | 1.16 | 1.29 |
| elix_CoagulopathyTRUE | 0.88 | 0.81 | 0.96 |
| elix_ObesityTRUE | 0.85 | 0.81 | 0.90 |
| elix_WeightLossTRUE | 1.54 | 1.45 | 1.62 |
| elix_FluidsLytesTRUE | 0.84 | 0.80 | 0.89 |
| elix_BloodLossTRUE | 0.75 | 0.66 | 0.84 |
| elix_AnemiaTRUE | 1.04 | 1.00 | 1.09 |
| elix_AlcoholTRUE | 0.85 | 0.77 | 0.95 |
| elix_DrugsTRUE | 0.85 | 0.76 | 0.96 |
| elix_PsychosesTRUE | 0.86 | 0.82 | 0.91 |
| elix_DepressionTRUE | 0.98 | 0.94 | 1.03 |
| as.factor(any_sx)TRUE | 4.67 | 4.17 | 5.24 |
| lives_in_msaTRUE | 1.42 | 1.27 | 1.58 |
| as.factor(any_sx)TRUE:lives_in_msaTRUE | 0.82 | 0.73 | 0.92 |
And the dose-response:
coxph(
Surv(age, als) ~ rate_outpatient_visits + intensity_outpatient_visits_dx +
rate_inpatient_stays + intensity_inpatient_stays_los +
elix_CHF + elix_Valvular + elix_PHTN + elix_PVD + elix_HTN +
elix_HTNcx + elix_Paralysis + elix_NeuroOther +
elix_Pulmonary + elix_DM + elix_DMcx + elix_Hypothyroid +
elix_Renal + elix_Liver + elix_PUD + elix_HIV +
elix_Lymphoma + elix_Mets + elix_Tumor + elix_Rheumatic +
elix_Coagulopathy + elix_Obesity + elix_WeightLoss +
elix_FluidsLytes + elix_BloodLoss + elix_Anemia +
elix_Alcohol + elix_Drugs + elix_Psychoses +
elix_Depression +
as.factor(n_total_topcoded) * lives_in_msa +
strata(stratum),
data = model_data |>
mutate(
n_total_topcoded = case_when(
n_total <= 4 ~ as.character(n_total),
n_total > 4 ~ "5+"
)
),
robust = TRUE,
cluster = stratum
) |>
broom::tidy(exponentiate = TRUE, conf.int = TRUE) |>
select(term, estimate, conf.low, conf.high) |>
knitr::kable(digits = 2)| term | estimate | conf.low | conf.high |
|---|---|---|---|
| rate_outpatient_visits | 1.01 | 1.01 | 1.01 |
| intensity_outpatient_visits_dx | 1.05 | 1.01 | 1.08 |
| rate_inpatient_stays | 0.83 | 0.79 | 0.88 |
| intensity_inpatient_stays_los | 0.99 | 0.99 | 1.00 |
| elix_CHFTRUE | 0.71 | 0.66 | 0.77 |
| elix_ValvularTRUE | 0.96 | 0.91 | 1.01 |
| elix_PHTNTRUE | 1.00 | 0.92 | 1.10 |
| elix_PVDTRUE | 0.88 | 0.84 | 0.93 |
| elix_HTNTRUE | 0.99 | 0.95 | 1.03 |
| elix_HTNcxTRUE | 0.92 | 0.87 | 0.98 |
| elix_ParalysisTRUE | 1.61 | 1.52 | 1.71 |
| elix_NeuroOtherTRUE | 3.03 | 2.91 | 3.16 |
| elix_PulmonaryTRUE | 0.96 | 0.93 | 1.00 |
| elix_DMTRUE | 0.92 | 0.88 | 0.96 |
| elix_DMcxTRUE | 0.82 | 0.77 | 0.88 |
| elix_HypothyroidTRUE | 1.08 | 1.03 | 1.13 |
| elix_RenalTRUE | 0.67 | 0.62 | 0.73 |
| elix_LiverTRUE | 0.95 | 0.89 | 1.02 |
| elix_PUDTRUE | 1.00 | 0.81 | 1.24 |
| elix_HIVTRUE | 1.09 | 0.77 | 1.54 |
| elix_LymphomaTRUE | 1.19 | 1.06 | 1.35 |
| elix_MetsTRUE | 0.42 | 0.38 | 0.48 |
| elix_TumorTRUE | 0.91 | 0.86 | 0.96 |
| elix_RheumaticTRUE | 1.10 | 1.04 | 1.16 |
| elix_CoagulopathyTRUE | 0.87 | 0.80 | 0.94 |
| elix_ObesityTRUE | 0.82 | 0.78 | 0.87 |
| elix_WeightLossTRUE | 1.40 | 1.32 | 1.48 |
| elix_FluidsLytesTRUE | 0.81 | 0.77 | 0.85 |
| elix_BloodLossTRUE | 0.73 | 0.65 | 0.82 |
| elix_AnemiaTRUE | 1.02 | 0.97 | 1.06 |
| elix_AlcoholTRUE | 0.88 | 0.79 | 0.98 |
| elix_DrugsTRUE | 0.85 | 0.76 | 0.96 |
| elix_PsychosesTRUE | 0.85 | 0.81 | 0.90 |
| elix_DepressionTRUE | 0.92 | 0.88 | 0.96 |
| as.factor(n_total_topcoded)1 | 2.61 | 2.27 | 2.99 |
| as.factor(n_total_topcoded)2 | 6.82 | 5.92 | 7.85 |
| as.factor(n_total_topcoded)3 | 11.78 | 10.20 | 13.60 |
| as.factor(n_total_topcoded)4 | 17.20 | 14.68 | 20.15 |
| as.factor(n_total_topcoded)5+ | 17.64 | 14.72 | 21.13 |
| lives_in_msaTRUE | 1.44 | 1.29 | 1.60 |
| as.factor(n_total_topcoded)1:lives_in_msaTRUE | 0.96 | 0.83 | 1.11 |
| as.factor(n_total_topcoded)2:lives_in_msaTRUE | 0.77 | 0.67 | 0.90 |
| as.factor(n_total_topcoded)3:lives_in_msaTRUE | 0.76 | 0.65 | 0.88 |
| as.factor(n_total_topcoded)4:lives_in_msaTRUE | 0.74 | 0.63 | 0.88 |
| as.factor(n_total_topcoded)5+:lives_in_msaTRUE | 0.79 | 0.66 | 0.95 |