sensitivity_figures

Document compiled on July 08, 2024.

Helper functions

Dynamic figure (code)

Static figure (code)

ONE-WAYS

Ranges

intervals <- 15
treatment_effect_range <- seq(0,
                              control_effect,
                              length.out = intervals) 
ac_lec_range <- seq(1,
                     10,
                     length.out = intervals) 
hr_ac_use_range <- seq(1,
                       hr_ac_use * 3,
                       length.out = intervals) 
hr_ich_lec_range <- seq(1,
                       hr_ich_lec * 3,
                       length.out = intervals)
ry_ich_base_range <- seq(0,
                       ry_ich_base * 2,
                       length.out = intervals)
parameters <- uCDR_estiamtes <- rbind(c(1,0.80,0.74,0.59,0.36), # Landeriro et al. community-reported
                                      c(1,0.86,0.85,0.86,0.82), # Landeriro et al. patient-reported
                                      c(1,0.80,0.74,0.59,0.36), # Landeriro et al. care-giver proxy
                                      c(1,0.88,0.80,0.71,0.63), # Sullivan/Hyle et al.
                                      c(1,0.73,0.69,0.53,0.38)) # Neumann et al.

ICH base rate

The ICH base rate for an individual depends on their characteristics, including sex, age, and medical history. This sensitivity analysis applies a common multiplier across all individuals.

lec x AC interaction on ICH

Sensitivity to the additional interaction effect of apixiban and lecanemab (1-10, base case = 2.21). x-axis below is the total ICH hazard ratio for a patient on both apixiban and lecanemab, i.e. the interaction effect (varies) x lecanemab HR on ICH (~2.02) x apixiban HR on ICH (~2.21) such that in the base case, the total hazard ratio is ~9.92. This is a misleading axis for the blue strategy (Apixiban alone), which does not include the interaction effect.

###For abstract

Lecanemab hazard ratio on ICH

Sensitivity to the lecanemab-specific hazard on ICH. Here again the x-axis displays the total hazard, including apixiban and the interaction. Again, this is a misleading axis for the blue strategy (Apixiban alone), which does not include the interaction effect.

# ggplot
data_hr_ich_lec <- read.csv('/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/hr_ich_lec_plotdata.csv')
plot_data$variable <- plot_data$variable 

p <- generate_ggplot(data_hr_ich_lec)

p + geom_vline(xintercept = hr_ich_lec , color = 'grey', linetype = 'dashed') +
    labs(title = 'Sensitivity to lecanmeab HR on ICH') +
    xlab('hazard ratio')

Lecanemab treatment effect (CDR)

Sensitivity to the treatment effect of lecanemab on cognition as measured by the CDR-SB score. The control group experiences a CDR-SB score increase of 1.66 points over 18-months. I vary the CDR-SB score increase over 18-months in the treatment group from 0 (perfect efficacy) to 1.66 (no efficacy).

# ggplot
data_treatment_effect <- read.csv('/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/treatment_effect_plotdata.csv')

p <- generate_ggplot(data_treatment_effect)

p + geom_vline(xintercept = treatment_effect, color = 'grey', linetype = 'dashed') +
    labs(title = 'Lecanemab treatment effect') +
    xlab('treatment effect')

Control CDR-SB rate of change

Sensitivity of model to the base rate of cognitive decline, i.e. the rate of change in CDR-SB score in the control group. The rate of cognitive decline in the treated group is then modeled using a 1) absolute treatment effect of ~0.5 CDR-SB points per year or 2) a relative treatment effect of a ~37% decline in the rate of CDR-SB progression compared to control.

# ggplot

# absolute TE
data_treatment_effect <- read.csv('/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/control_effect_abs_plotdata.csv')

p <- generate_ggplot(data_treatment_effect)

p + geom_vline(xintercept = treatment_effect, color = 'grey', linetype = 'dashed') +
    labs(title = 'Control group CDR-SB score change (absolute lecanemab TE)') +
    xlab('CDR-SB increase in control group')

# relative TE
data_treatment_effect <- read.csv('/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/control_effect_rel_plotdata.csv')

p <- generate_ggplot(data_treatment_effect)

p + geom_vline(xintercept = treatment_effect, color = 'grey', linetype = 'dashed') +
    labs(title = 'Control group CDR-SB score change (relative lecanemab TE)') +
    xlab('CDR-SB increase in control group')

Cognition-related HR on background mortality

data_cog_bmort <- read.csv('/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/cog_bmort_hr_plotdata.csv')

p <- generate_ggplot(data_cog_bmort)

p + geom_vline(xintercept = 1, color = 'grey', linetype = 'dashed') +
    labs(title = 'MCI- and dementia-related hazard ratio on background mortality') +
    xlab('multiplier')

BY AGE

helper functions

ICH base rate plots

plot_data_65_74 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ry_ich_base_plotdata_65_74.csv")
plot_data_75_84 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ry_ich_base_plotdata_75_84.csv")
plot_data_85_90 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ry_ich_base_plotdata_85_90.csv")

# process each dataset
plot_QALE_1 <- process_plot_data(plot_data_65_74,'QALE',0.0001,0.0012,11.5,13.5,basecase = ry_ich_base)
plot_QALE_2 <- process_plot_data(plot_data_75_84,'QALE',0.0001,0.0012,11.5,13.5,basecase = ry_ich_base)
plot_QALE_3 <- process_plot_data(plot_data_85_90,'QALE',0.0001,0.0012,11.5,13.5,basecase = ry_ich_base)

plot_ich_1 <- process_plot_data(plot_data_65_74,'ICH.count.per.1000.pys',0.0001,0.0012,0,250,basecase = ry_ich_base)
plot_ich_2 <- process_plot_data(plot_data_75_84,'ICH.count.per.1000.pys',0.0001,0.0012,0,250,basecase = ry_ich_base)
plot_ich_3 <- process_plot_data(plot_data_85_90,'ICH.count.per.1000.pys',0.0001,0.0012,0,250,basecase = ry_ich_base)

plot_ischemic_1 <- process_plot_data(plot_data_65_74,'Ischemic.per.1000.pys',0.0001,0.0012,15,23,basecase = ry_ich_base)
plot_ischemic_2 <- process_plot_data(plot_data_75_84,'Ischemic.per.1000.pys',0.0001,0.0012,15,23,basecase = ry_ich_base)
plot_ischemic_3 <- process_plot_data(plot_data_85_90,'Ischemic.per.1000.pys',0.0001,0.0012,15,23,basecase = ry_ich_base)

plot_cdr_1 <- process_plot_data(plot_data_65_74,'CDR.SB.increase',0.0001,0.0012,0.8,1.6,basecase = ry_ich_base)
plot_cdr_2 <- process_plot_data(plot_data_75_84,'CDR.SB.increase',0.0001,0.0012,0.8,1.6,basecase = ry_ich_base)
plot_cdr_3 <- process_plot_data(plot_data_85_90,'CDR.SB.increase',0.0001,0.0012,0.8,1.6,basecase = ry_ich_base)

plot_data_65_74$total.deaths <- plot_data_65_74$ICH.deaths.per.1000.pys + plot_data_65_74$Ischemic.deaths.per.1000.pys + plot_data_65_74$background.death.per.1000.pys
plot_data_75_84$total.deaths <- plot_data_75_84$ICH.deaths.per.1000.pys + plot_data_75_84$Ischemic.deaths.per.1000.pys + plot_data_75_84$background.death.per.1000.pys
plot_data_85_90$total.deaths <- plot_data_85_90$ICH.deaths.per.1000.pys + plot_data_85_90$Ischemic.deaths.per.1000.pys + plot_data_85_90$background.death.per.1000.pys

plot_deaths_1 <- process_plot_data(plot_data_65_74,'total.deaths',0.0001,0.0012,0,500,basecase = ry_ich_base)
plot_deaths_2 <- process_plot_data(plot_data_75_84,'total.deaths',0.0001,0.0012,0,500,basecase = ry_ich_base)
plot_deaths_3 <- process_plot_data(plot_data_85_90,'total.deaths',0.0001,0.0012,0,500,basecase = ry_ich_base)

# create plots
QALE_plot <- create_combined_plot(plot_QALE_1, plot_QALE_2, plot_QALE_3, 
                                  dvar = "QALMs",
                                  ivar = "ICH base rate",
                                  xlab = "ICH base rate", 
                                  ylab = "Quality-adjusted life months")

print(QALE_plot)

ich_plot <- create_combined_plot(plot_ich_1, plot_ich_2, plot_ich_3, 
                                 dvar = 'ICH',
                                 ivar = "ICH base rate",
                                 xlab = "ICH base rate", 
                                 ylab = "ICH count per 1000 pys")
print(ich_plot)

ischemic_plot <- create_combined_plot(plot_ischemic_1, plot_ischemic_2, plot_ischemic_3, 
                                      dvar = 'stroke',
                                      ivar = "ICH base rate",
                                      xlab = "ICH base rate", 
                                      ylab = "Stroke count per 1000 pys")
print(ischemic_plot)

cdr_plot <- create_combined_plot(plot_cdr_1, plot_cdr_2, plot_cdr_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "ICH base rate",
                                 xlab = "ICH base rate", 
                                 ylab = "Increase in CDR-SB score")
print(cdr_plot)

deaths_plot <- create_combined_plot(plot_deaths_1, plot_deaths_2, plot_deaths_3, 
                                    dvar = 'CDR-SB',
                                    ivar = "ICH base rate",
                                    xlab = "ICH base rate", 
                                    ylab = "Total deaths per 1000 pys")
print(deaths_plot)

AC x lec interaction

Sensitivity of the interaction effect on QALE estimates, broken down by age group. Figure included in abstract

plot_data_65_74 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ac_lec_plotdata_65_74.csv")
plot_data_75_84 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ac_lec_plotdata_75_84.csv")
plot_data_85_90 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ac_lec_plotdata_85_90.csv")
plot_data_65_74$variable <- plot_data_65_74$variable * hr_ich_lec * hr_ac_use
plot_data_75_84$variable <- plot_data_75_84$variable * hr_ich_lec * hr_ac_use
plot_data_85_90$variable <- plot_data_85_90$variable * hr_ich_lec * hr_ac_use

# process each dataset
plot_QALE_1 <- process_plot_data(plot_data_65_74,'QALE',4.5,15,11.5,13.5,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_QALE_2 <- process_plot_data(plot_data_75_84,'QALE',4.5,15,11.5,13.5,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_QALE_3 <- process_plot_data(plot_data_85_90,'QALE',4.5,15,11.5,13.5,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)

plot_ich_1 <- process_plot_data(plot_data_65_74,'ICH.count.per.1000.pys',4.5,15,0,250,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_ich_2 <- process_plot_data(plot_data_75_84,'ICH.count.per.1000.pys',4.5,15,0,250,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_ich_3 <- process_plot_data(plot_data_85_90,'ICH.count.per.1000.pys',4.5,15,0,250,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)

plot_ischemic_1 <- process_plot_data(plot_data_65_74,'Ischemic.per.1000.pys',4.5,15,15,23,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_ischemic_2 <- process_plot_data(plot_data_75_84,'Ischemic.per.1000.pys',4.5,15,15,23,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_ischemic_3 <- process_plot_data(plot_data_85_90,'Ischemic.per.1000.pys',4.5,15,15,23,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)

plot_cdr_1 <- process_plot_data(plot_data_65_74,'CDR.SB.increase',4.5,15,0.8,1.6,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_cdr_2 <- process_plot_data(plot_data_75_84,'CDR.SB.increase',4.5,15,0.8,1.6,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_cdr_3 <- process_plot_data(plot_data_85_90,'CDR.SB.increase',4.5,15,0.8,1.6,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)

plot_data_65_74$total.deaths <- plot_data_65_74$ICH.deaths.per.1000.pys + plot_data_65_74$Ischemic.deaths.per.1000.pys + plot_data_65_74$background.death.per.1000.pys
plot_data_75_84$total.deaths <- plot_data_75_84$ICH.deaths.per.1000.pys + plot_data_75_84$Ischemic.deaths.per.1000.pys + plot_data_75_84$background.death.per.1000.pys
plot_data_85_90$total.deaths <- plot_data_85_90$ICH.deaths.per.1000.pys + plot_data_85_90$Ischemic.deaths.per.1000.pys + plot_data_85_90$background.death.per.1000.pys

plot_deaths_1 <- process_plot_data(plot_data_65_74,'total.deaths',4.5,15,0,500,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_deaths_2 <- process_plot_data(plot_data_75_84,'total.deaths',4.5,15,0,500,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)
plot_deaths_3 <- process_plot_data(plot_data_85_90,'total.deaths',4.5,15,0,500,basecase = ac_lec_interaction * hr_ich_lec * hr_ac_use)

# create plots
QALE_plot <- create_combined_plot(plot_QALE_1, plot_QALE_2, plot_QALE_3, 
                                  dvar = "QALMs",
                                  ivar = "lec x apixaban interaction effect",
                                  xlab = "Hazard of ICH with lecanemab and apixaban together", 
                                  ylab = "Quality-adjusted life months")

print(QALE_plot)

ich_plot <- create_combined_plot(plot_ich_1, plot_ich_2, plot_ich_3, 
                                 dvar = 'ICH',
                                 ivar = "lec x apixaban interaction effect",
                                 xlab = "Hazard of ICH with lecanemab and apixaban together", 
                                 ylab = "ICH count per 1000 pys")
print(ich_plot)

ischemic_plot <- create_combined_plot(plot_ischemic_1, plot_ischemic_2, plot_ischemic_3, 
                                      dvar = 'stroke',
                                      ivar = "lec x apixaban interaction effect",
                                      xlab = "Hazard of ICH with lecanemab and apixaban together", 
                                      ylab = "Stroke count per 1000 pys")
print(ischemic_plot)

cdr_plot <- create_combined_plot(plot_cdr_1, plot_cdr_2, plot_cdr_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lec x apixaban interaction effect",
                                 xlab = "Hazard of ICH with lecanemab and apixaban together", 
                                 ylab = "Increase in CDR-SB score")
print(cdr_plot)

deaths_plot <- create_combined_plot(plot_deaths_1, plot_deaths_2, plot_deaths_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lec x apixaban interaction effect",
                                 xlab = "Hazard of ICH with lecanemab and apixaban together", 
                                 ylab = "Total deaths per 1000 pys")
print(deaths_plot)

Lec ICH hazard rate plots

plot_data_65_74 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/hr_ich_lec_plotdata_65_74.csv")
plot_data_75_84 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/hr_ich_lec_plotdata_75_84.csv")
plot_data_85_90 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/hr_ich_lec_plotdata_85_90.csv")

# process each dataset
plot_QALE_1 <- process_plot_data(plot_data_65_74,'QALE',1,10,11.5,13.5,basecase = hr_ich_lec)
plot_QALE_2 <- process_plot_data(plot_data_75_84,'QALE',1,10,11.5,13.5,basecase = hr_ich_lec)
plot_QALE_3 <- process_plot_data(plot_data_85_90,'QALE',1,10,11.5,13.5,basecase = hr_ich_lec)

plot_ich_1 <- process_plot_data(plot_data_65_74,'ICH.count.per.1000.pys',1,10,0,250,basecase = hr_ich_lec)
plot_ich_2 <- process_plot_data(plot_data_75_84,'ICH.count.per.1000.pys',1,10,0,250,basecase = hr_ich_lec)
plot_ich_3 <- process_plot_data(plot_data_85_90,'ICH.count.per.1000.pys',1,10,0,250,basecase = hr_ich_lec)

plot_ischemic_1 <- process_plot_data(plot_data_65_74,'Ischemic.per.1000.pys',1,10,15,23,basecase = hr_ich_lec)
plot_ischemic_2 <- process_plot_data(plot_data_75_84,'Ischemic.per.1000.pys',1,10,15,23,basecase = hr_ich_lec)
plot_ischemic_3 <- process_plot_data(plot_data_85_90,'Ischemic.per.1000.pys',1,10,15,23,basecase = hr_ich_lec)

plot_cdr_1 <- process_plot_data(plot_data_65_74,'CDR.SB.increase',1,10,0.8,1.6,basecase = hr_ich_lec)
plot_cdr_2 <- process_plot_data(plot_data_75_84,'CDR.SB.increase',1,10,0.8,1.6,basecase = hr_ich_lec)
plot_cdr_3 <- process_plot_data(plot_data_85_90,'CDR.SB.increase',1,10,0.8,1.6,basecase = hr_ich_lec)

plot_data_65_74$total.deaths <- plot_data_65_74$ICH.deaths.per.1000.pys + plot_data_65_74$Ischemic.deaths.per.1000.pys + plot_data_65_74$background.death.per.1000.pys
plot_data_75_84$total.deaths <- plot_data_75_84$ICH.deaths.per.1000.pys + plot_data_75_84$Ischemic.deaths.per.1000.pys + plot_data_75_84$background.death.per.1000.pys
plot_data_85_90$total.deaths <- plot_data_85_90$ICH.deaths.per.1000.pys + plot_data_85_90$Ischemic.deaths.per.1000.pys + plot_data_85_90$background.death.per.1000.pys

plot_deaths_1 <- process_plot_data(plot_data_65_74,'total.deaths',1,10,0,500,basecase = hr_ich_lec)
plot_deaths_2 <- process_plot_data(plot_data_75_84,'total.deaths',1,10,0,500,basecase = hr_ich_lec)
plot_deaths_3 <- process_plot_data(plot_data_85_90,'total.deaths',1,10,0,500,basecase = hr_ich_lec)

# create plots
QALE_plot <- create_combined_plot(plot_QALE_1, plot_QALE_2, plot_QALE_3, 
                                  dvar = "QALMs",
                                  ivar = "lecanemab ICH hazard",
                                  xlab = "Lecanemab ICH hazard", 
                                  ylab = "Quality-adjusted life months")

print(QALE_plot)

ich_plot <- create_combined_plot(plot_ich_1, plot_ich_2, plot_ich_3, 
                                 dvar = 'ICH',
                                 ivar = "lecanemab ICH hazard",
                                 xlab = "Lecanemab ICH hazard", 
                                 ylab = "ICH count per 1000 pys")
print(ich_plot)

ischemic_plot <- create_combined_plot(plot_ischemic_1, plot_ischemic_2, plot_ischemic_3, 
                                      dvar = 'stroke',
                                      ivar = "lecanemab ICH hazard",
                                      xlab = "Lecanemab ICH hazard", 
                                      ylab = "Stroke count per 1000 pys")
print(ischemic_plot)

cdr_plot <- create_combined_plot(plot_cdr_1, plot_cdr_2, plot_cdr_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lecanemab ICH hazard",
                                 xlab = "Lecanemab ICH hazard", 
                                 ylab = "Increase in CDR-SB score")
print(cdr_plot)

deaths_plot <- create_combined_plot(plot_deaths_1, plot_deaths_2, plot_deaths_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lecanemab ICH hazard",
                                 xlab = "Lecanemab ICH hazard", 
                                 ylab = "Total deaths per 1000 pys")
print(deaths_plot)

Apixaban ICH hazard rate plots

plot_data_65_74 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/hr_ac_use_plotdata_65_74.csv")
plot_data_75_84 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/hr_ac_use_plotdata_75_84.csv")
plot_data_85_90 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/hr_ac_use_plotdata_85_90.csv")

# process each dataset
plot_QALE_1 <- process_plot_data(plot_data_65_74,'QALE',1,10,11.5,13.5,basecase = hr_ac_use)
plot_QALE_2 <- process_plot_data(plot_data_75_84,'QALE',1,10,11.5,13.5,basecase = hr_ac_use)
plot_QALE_3 <- process_plot_data(plot_data_85_90,'QALE',1,10,11.5,13.5,basecase = hr_ac_use)

plot_ich_1 <- process_plot_data(plot_data_65_74,'ICH.count.per.1000.pys',1,10,0,250,basecase = hr_ac_use)
plot_ich_2 <- process_plot_data(plot_data_75_84,'ICH.count.per.1000.pys',1,10,0,250,basecase = hr_ac_use)
plot_ich_3 <- process_plot_data(plot_data_85_90,'ICH.count.per.1000.pys',1,10,0,250,basecase = hr_ac_use)

plot_ischemic_1 <- process_plot_data(plot_data_65_74,'Ischemic.per.1000.pys',1,10,15,23,basecase = hr_ac_use)
plot_ischemic_2 <- process_plot_data(plot_data_75_84,'Ischemic.per.1000.pys',1,10,15,23,basecase = hr_ac_use)
plot_ischemic_3 <- process_plot_data(plot_data_85_90,'Ischemic.per.1000.pys',1,10,15,23,basecase = hr_ac_use)

plot_cdr_1 <- process_plot_data(plot_data_65_74,'CDR.SB.increase',1,10,0.8,1.6,basecase = hr_ac_use)
plot_cdr_2 <- process_plot_data(plot_data_75_84,'CDR.SB.increase',1,10,0.8,1.6,basecase = hr_ac_use)
plot_cdr_3 <- process_plot_data(plot_data_85_90,'CDR.SB.increase',1,10,0.8,1.6,basecase = hr_ac_use)

plot_data_65_74$total.deaths <- plot_data_65_74$ICH.deaths.per.1000.pys + plot_data_65_74$Ischemic.deaths.per.1000.pys + plot_data_65_74$background.death.per.1000.pys
plot_data_75_84$total.deaths <- plot_data_75_84$ICH.deaths.per.1000.pys + plot_data_75_84$Ischemic.deaths.per.1000.pys + plot_data_75_84$background.death.per.1000.pys
plot_data_85_90$total.deaths <- plot_data_85_90$ICH.deaths.per.1000.pys + plot_data_85_90$Ischemic.deaths.per.1000.pys + plot_data_85_90$background.death.per.1000.pys

plot_deaths_1 <- process_plot_data(plot_data_65_74,'total.deaths',1,10,0,500,basecase = hr_ac_use)
plot_deaths_2 <- process_plot_data(plot_data_75_84,'total.deaths',1,10,0,500,basecase = hr_ac_use)
plot_deaths_3 <- process_plot_data(plot_data_85_90,'total.deaths',1,10,0,500,basecase = hr_ac_use)

# create plots
QALE_plot <- create_combined_plot(plot_QALE_1, plot_QALE_2, plot_QALE_3, 
                                  dvar = "QALMs",
                                  ivar = "apixaban ICH hazard",
                                  xlab = "Apixaban ICH hazard", 
                                  ylab = "Quality-adjusted life months")

print(QALE_plot)

ich_plot <- create_combined_plot(plot_ich_1, plot_ich_2, plot_ich_3, 
                                 dvar = 'ICH',
                                  ivar = "apixaban ICH hazard",
                                  xlab = "Apixaban ICH hazard",  
                                 ylab = "ICH count per 1000 pys")
print(ich_plot)

ischemic_plot <- create_combined_plot(plot_ischemic_1, plot_ischemic_2, plot_ischemic_3, 
                                      dvar = 'stroke',
                                  ivar = "apixaban ICH hazard",
                                  xlab = "Apixaban ICH hazard",  
                                      ylab = "Stroke count per 1000 pys")
print(ischemic_plot)

cdr_plot <- create_combined_plot(plot_cdr_1, plot_cdr_2, plot_cdr_3, 
                                 dvar = 'CDR-SB',
                                  ivar = "apixaban ICH hazard",
                                  xlab = "Apixaban ICH hazard", 
                                 ylab = "Increase in CDR-SB score")
print(cdr_plot)

deaths_plot <- create_combined_plot(plot_deaths_1, plot_deaths_2, plot_deaths_3, 
                                 dvar = 'CDR-SB',
                                  ivar = "apixaban ICH hazard",
                                  xlab = "Apixaban ICH hazard", 
                                 ylab = "Total deaths per 1000 pys")
print(deaths_plot)

Lec TE on CDR-SB score plots

plot_data_65_74 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/treatment_effect_plotdata_65_74.csv")
plot_data_75_84 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/treatment_effect_plotdata_75_84.csv")
plot_data_85_90 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/treatment_effect_plotdata_85_90.csv")

# process each dataset
plot_QALE_1 <- process_plot_data(plot_data_65_74,'QALE',0,control_effect,11.5,13.5,basecase = treatment_effect)
plot_QALE_2 <- process_plot_data(plot_data_75_84,'QALE',0,control_effect,11.5,13.5,basecase = treatment_effect)
plot_QALE_3 <- process_plot_data(plot_data_85_90,'QALE',0,control_effect,11.5,13.5,basecase = treatment_effect)

plot_ich_1 <- process_plot_data(plot_data_65_74,'ICH.count.per.1000.pys',0,control_effect,0,250,basecase = treatment_effect)
plot_ich_2 <- process_plot_data(plot_data_75_84,'ICH.count.per.1000.pys',0,control_effect,0,250,basecase = treatment_effect)
plot_ich_3 <- process_plot_data(plot_data_85_90,'ICH.count.per.1000.pys',0,control_effect,0,250,basecase = treatment_effect)

plot_ischemic_1 <- process_plot_data(plot_data_65_74,'Ischemic.per.1000.pys',0,control_effect,15,23,basecase = treatment_effect)
plot_ischemic_2 <- process_plot_data(plot_data_75_84,'Ischemic.per.1000.pys',0,control_effect,15,23,basecase = treatment_effect)
plot_ischemic_3 <- process_plot_data(plot_data_85_90,'Ischemic.per.1000.pys',0,control_effect,15,23,basecase = treatment_effect)

plot_cdr_1 <- process_plot_data(plot_data_65_74,'CDR.SB.increase',0,control_effect,0,control_effect,basecase = treatment_effect)
plot_cdr_2 <- process_plot_data(plot_data_75_84,'CDR.SB.increase',0,control_effect,0,control_effect,basecase = treatment_effect)
plot_cdr_3 <- process_plot_data(plot_data_85_90,'CDR.SB.increase',0,control_effect,0,control_effect,basecase = treatment_effect)

plot_data_65_74$total.deaths <- plot_data_65_74$ICH.deaths.per.1000.pys + plot_data_65_74$Ischemic.deaths.per.1000.pys + plot_data_65_74$background.death.per.1000.pys
plot_data_75_84$total.deaths <- plot_data_75_84$ICH.deaths.per.1000.pys + plot_data_75_84$Ischemic.deaths.per.1000.pys + plot_data_75_84$background.death.per.1000.pys
plot_data_85_90$total.deaths <- plot_data_85_90$ICH.deaths.per.1000.pys + plot_data_85_90$Ischemic.deaths.per.1000.pys + plot_data_85_90$background.death.per.1000.pys

plot_deaths_1 <- process_plot_data(plot_data_65_74,'total.deaths',0,control_effect,0,500,basecase = treatment_effect)
plot_deaths_2 <- process_plot_data(plot_data_75_84,'total.deaths',0,control_effect,0,500,basecase = treatment_effect)
plot_deaths_3 <- process_plot_data(plot_data_85_90,'total.deaths',0,control_effect,0,500,basecase = treatment_effect)

# create plots
QALE_plot <- create_combined_plot(plot_QALE_1, plot_QALE_2, plot_QALE_3, 
                                  dvar = "QALMs",
                                  ivar = "lecanemab's TE on CDR-SB score",
                                  xlab = "Lecanemab treatment effect (complete -> no efficacy)", 
                                  ylab = "Quality-adjusted life months")

print(QALE_plot)

ich_plot <- create_combined_plot(plot_ich_1, plot_ich_2, plot_ich_3, 
                                 dvar = 'ICH',
                                 ivar = "lecanemab's TE on CDR-SB score",
                                 xlab = "Lecanemab treatment effect (complete -> no efficacy)",
                                 ylab = "ICH count per 1000 pys")
print(ich_plot)

ischemic_plot <- create_combined_plot(plot_ischemic_1, plot_ischemic_2, plot_ischemic_3, 
                                      dvar = 'stroke',
                                      ivar = "lecanemab's TE on CDR-SB score",
                                      xlab = "Lecanemab treatment effect (complete -> no efficacy)", 
                                      ylab = "Stroke count per 1000 pys")
print(ischemic_plot)

cdr_plot <- create_combined_plot(plot_cdr_1, plot_cdr_2, plot_cdr_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lecanemab's TE on CDR-SB score",
                                 xlab = "Lecanemab treatment effect (complete -> no efficacy)",
                                 ylab = "Increase in CDR-SB score")
print(cdr_plot)

deaths_plot <- create_combined_plot(plot_deaths_1, plot_deaths_2, plot_deaths_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lecanemab's TE on CDR-SB score",
                                 xlab = "Lecanemab treatment effect (complete -> no efficacy)",
                                 ylab = "Total deaths per 1000 pys")
print(deaths_plot)

CDR-SB progression in control (ABSOLUTE treatment effect)

plot_data_65_74 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/control_effect_plotdata_65_74.csv")
plot_data_75_84 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/control_effect_plotdata_75_84.csv")
plot_data_85_90 <- read.csv("/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/control_effect_plotdata_85_90.csv")

# process each dataset
plot_QALE_1 <- process_plot_data(plot_data_65_74,'QALE',control_effect / 2,control_effect * 2,11.5,13.5,basecase = control_effect)
plot_QALE_2 <- process_plot_data(plot_data_75_84,'QALE',control_effect / 2,control_effect * 2,11.5,13.5,basecase = control_effect)
plot_QALE_3 <- process_plot_data(plot_data_85_90,'QALE',control_effect / 2,control_effect * 2,11.5,13.5,basecase = control_effect)

plot_ich_1 <- process_plot_data(plot_data_65_74,'ICH.count.per.1000.pys',control_effect / 2,control_effect * 2,0,250,basecase = control_effect)
plot_ich_2 <- process_plot_data(plot_data_75_84,'ICH.count.per.1000.pys',control_effect / 2,control_effect * 2,0,250,basecase = control_effect)
plot_ich_3 <- process_plot_data(plot_data_85_90,'ICH.count.per.1000.pys',control_effect / 2,control_effect * 2,0,250,basecase = control_effect)

plot_ischemic_1 <- process_plot_data(plot_data_65_74,'Ischemic.per.1000.pys',control_effect / 2,control_effect * 2,15,23,basecase = control_effect)
plot_ischemic_2 <- process_plot_data(plot_data_75_84,'Ischemic.per.1000.pys',control_effect / 2,control_effect * 2,15,23,basecase = control_effect)
plot_ischemic_3 <- process_plot_data(plot_data_85_90,'Ischemic.per.1000.pys',control_effect / 2,control_effect * 2,15,23,basecase = control_effect)

plot_cdr_1 <- process_plot_data(plot_data_65_74,'CDR.SB.increase',control_effect / 2,control_effect * 2,0,control_effect,basecase = control_effect)
plot_cdr_2 <- process_plot_data(plot_data_75_84,'CDR.SB.increase',control_effect / 2,control_effect * 2,0,control_effect,basecase = control_effect)
plot_cdr_3 <- process_plot_data(plot_data_85_90,'CDR.SB.increase',control_effect / 2,control_effect * 2,0,control_effect,basecase = control_effect)

plot_data_65_74$total.deaths <- plot_data_65_74$ICH.deaths.per.1000.pys + plot_data_65_74$Ischemic.deaths.per.1000.pys + plot_data_65_74$background.death.per.1000.pys
plot_data_75_84$total.deaths <- plot_data_75_84$ICH.deaths.per.1000.pys + plot_data_75_84$Ischemic.deaths.per.1000.pys + plot_data_75_84$background.death.per.1000.pys
plot_data_85_90$total.deaths <- plot_data_85_90$ICH.deaths.per.1000.pys + plot_data_85_90$Ischemic.deaths.per.1000.pys + plot_data_85_90$background.death.per.1000.pys

plot_deaths_1 <- process_plot_data(plot_data_65_74,'total.deaths',control_effect / 2,control_effect * 2,0,500,basecase = control_effect)
plot_deaths_2 <- process_plot_data(plot_data_75_84,'total.deaths',control_effect / 2,control_effect * 2,0,500,basecase = control_effect)
plot_deaths_3 <- process_plot_data(plot_data_85_90,'total.deaths',control_effect / 2,control_effect * 2,0,500,basecase = control_effect)

# create plots
QALE_plot <- create_combined_plot(plot_QALE_1, plot_QALE_2, plot_QALE_3, 
                                  dvar = "QALMs",
                                  ivar = "lecanemab's TE on CDR-SB score",
                                  xlab = "Lecanemab treatment effect (complete -> no efficacy)", 
                                  ylab = "Quality-adjusted life months")

print(QALE_plot)

ich_plot <- create_combined_plot(plot_ich_1, plot_ich_2, plot_ich_3, 
                                 dvar = 'ICH',
                                 ivar = "lecanemab's TE on CDR-SB score",
                                 xlab = "Lecanemab treatment effect (complete -> no efficacy)",
                                 ylab = "ICH count per 1000 pys")
print(ich_plot)

ischemic_plot <- create_combined_plot(plot_ischemic_1, plot_ischemic_2, plot_ischemic_3, 
                                      dvar = 'stroke',
                                      ivar = "lecanemab's TE on CDR-SB score",
                                      xlab = "Lecanemab treatment effect (complete -> no efficacy)", 
                                      ylab = "Stroke count per 1000 pys")
print(ischemic_plot)

cdr_plot <- create_combined_plot(plot_cdr_1, plot_cdr_2, plot_cdr_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lecanemab's TE on CDR-SB score",
                                 xlab = "Lecanemab treatment effect (complete -> no efficacy)",
                                 ylab = "Increase in CDR-SB score")
print(cdr_plot)

deaths_plot <- create_combined_plot(plot_deaths_1, plot_deaths_2, plot_deaths_3, 
                                 dvar = 'CDR-SB',
                                 ivar = "lecanemab's TE on CDR-SB score",
                                 xlab = "Lecanemab treatment effect (complete -> no efficacy)",
                                 ylab = "Total deaths per 1000 pys")
print(deaths_plot)

TWO-WAYS

Lecanemab efficacy and apixaban x Lec ICH hazard

ac_lec_vs_lec_TE_results <- readRDS('./Results/sensitivityanalysis/ac_lec_vs_lec_TE_results.rds')

intervals <- 15
treatment_effect_range <- seq(0,
                              control_effect,
                              length.out = intervals) 

ac_lec_range <- seq(1,
                     10,
                     length.out = intervals) 

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(ac_lec_vs_lec_TE_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(ac_lec_vs_lec_TE_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(ac_lec_vs_lec_TE_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), ac_lec_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','ac_lec_interaction','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$ac_lec_interaction <- as.factor(results_df$ac_lec_interaction)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * ((control_effect - treatment_effect) - min(treatment_effect_range)) / (max(treatment_effect_range) - min(treatment_effect_range)) + 1,
                       y = (intervals - 1) * (ac_lec_interaction - min(ac_lec_range)) / (max(ac_lec_range) - min(ac_lec_range)) + 1)

# plot

ggplot(results_df, aes(x = efficacy, y = ac_lec_interaction)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lec efficacy and AC x lec ICH hazard",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "Apixaban x Lec interaction effect",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

write.csv(results_df,'/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ac_lec_vs_lec_TE_results_plotdata.csv')

ac_lec_vs_lec_TE_results <- readRDS('./Results/sensitivityanalysis/ac_lec_vs_lec_TE_results_65_74.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(ac_lec_vs_lec_TE_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(ac_lec_vs_lec_TE_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(ac_lec_vs_lec_TE_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), ac_lec_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','ac_lec_interaction','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$ac_lec_interaction <- as.factor(results_df$ac_lec_interaction)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * ((control_effect - treatment_effect) - min(treatment_effect_range)) / (max(treatment_effect_range) - min(treatment_effect_range)) + 1,
                       y = (intervals - 1) * (ac_lec_interaction - min(ac_lec_range)) / (max(ac_lec_range) - min(ac_lec_range)) + 1)


# plot

ggplot(results_df, aes(x = efficacy, y = ac_lec_interaction)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lec efficacy and AC x lec ICH hazard (65-74)",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "Apixaban x Lec interaction effect",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

write.csv(results_df,'/Users/taidinger/Github/lecanemab_manuscript/R/Results/plot_data/ac_lec_vs_lec_TE_results_plotdata_65_74.csv')

Lecanemab efficacy and apixaban ICH hazard

lecTE_acHR_results <- readRDS('./Results/sensitivityanalysis/lecTE_acHR_results.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecTE_acHR_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecTE_acHR_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecTE_acHR_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), hr_ac_use_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','hr_ac_use','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$hr_ac_use <- as.factor(results_df$hr_ac_use)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * ((control_effect - treatment_effect) - min(treatment_effect_range)) / (max(treatment_effect_range) - min(treatment_effect_range)) + 1,
                       y = (intervals - 1) * (hr_ac_use - min(hr_ac_use_range)) / (max(hr_ac_use_range) - min(hr_ac_use_range)) + 1)

# plot

ggplot(results_df, aes(x = efficacy, y = hr_ac_use)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lecanemab efficacy and apixaban ICH hazard",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "Apixaban ICH hazard ratio",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

lecTE_acHR_results <- readRDS('./Results/sensitivityanalysis/lecTE_acHR_results_65_74.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecTE_acHR_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecTE_acHR_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecTE_acHR_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), hr_ac_use_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','hr_ac_use','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$hr_ac_use <- as.factor(results_df$hr_ac_use)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * ((control_effect - treatment_effect) - min(treatment_effect_range)) / (max(treatment_effect_range) - min(treatment_effect_range)) + 1,
                       y = (intervals - 1) * (hr_ac_use - min(hr_ac_use_range)) / (max(hr_ac_use_range) - min(hr_ac_use_range)) + 1)

# plot

ggplot(results_df, aes(x = efficacy, y = hr_ac_use)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lecanemab efficacy and apixaban ICH hazard (65-74)",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "Apixaban ICH hazard ratio",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

Lecanemab efficacy and lecanemab ICH hazard

lecTE_lecHR_results <- readRDS('./Results/sensitivityanalysis/lecTE_lecHR_results.rds')
intervals <- 15
# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecTE_lecHR_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecTE_lecHR_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecTE_lecHR_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), hr_ich_lec_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','hr_ich_lec','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$hr_ich_lec <- as.factor(results_df$hr_ich_lec)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * ((control_effect - treatment_effect) - min(treatment_effect_range)) / (max(treatment_effect_range) - min(treatment_effect_range)) + 1,
                       y = (intervals - 1) * (hr_ich_lec - min(hr_ich_lec_range)) / (max(hr_ich_lec_range) - min(hr_ich_lec_range)) + 1)

# plot

ggplot(results_df, aes(x = efficacy, y = hr_ich_lec)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lecanemab efficacy and ICH hazard",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "Lecanemab ICH hazard ratio",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

lecTE_lecHR_results <- readRDS('./Results/sensitivityanalysis/lecTE_lecHR_results_65_74.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecTE_lecHR_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecTE_lecHR_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecTE_lecHR_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), hr_ich_lec_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','hr_ich_lec','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$hr_ich_lec <- as.factor(results_df$hr_ich_lec)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * ((control_effect - treatment_effect) - min(treatment_effect_range)) / (max(treatment_effect_range) - min(treatment_effect_range)) + 1,
                       y = (intervals - 1) * (hr_ich_lec - min(hr_ich_lec_range)) / (max(hr_ich_lec_range) - min(hr_ich_lec_range)) + 1)

# plot

ggplot(results_df, aes(x = efficacy, y = hr_ich_lec)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lecanemab efficacy and ICH hazard (65-74)",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "Lecanemab ICH hazard ratio",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

Lecanemab efficacy and ICH base rate

lecTE_ichbase_results <- readRDS('./Results/sensitivityanalysis/lecTE_ichbase_results.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecTE_ichbase_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecTE_ichbase_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecTE_ichbase_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), ry_ich_base_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','ry_ich_base','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$ry_ich_base <- factor(results_df$ry_ich_base)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals * max(treatment_effect_range) - (intervals - 1) * (max(treatment_effect_range) - (control_effect - treatment_effect))) / max(treatment_effect_range),
                       y = (intervals * max(ry_ich_base_range) - (intervals - 1) * (max(ry_ich_base_range) - ry_ich_base)) / max(ry_ich_base_range))

# plot

ggplot(results_df, aes(x = efficacy, y = ry_ich_base)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lecanemab efficacy and ICH base rate",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "ICH base rate",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

lecTE_ichbase_results <- readRDS('./Results/sensitivityanalysis/lecTE_ichbase_results_65_74.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecTE_ichbase_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecTE_ichbase_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecTE_ichbase_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(round(treatment_effect_range,2), ry_ich_base_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                         results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                         results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                         results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('treatment_effect','ry_ich_base','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$efficacy <- as.factor(abs(results_df$treatment_effect - control_effect))
results_df$ry_ich_base <- factor(results_df$ry_ich_base)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals * max(treatment_effect_range) - (intervals - 1) * (max(treatment_effect_range) - (control_effect - treatment_effect))) / max(treatment_effect_range),
                       y = (intervals * max(ry_ich_base_range) - (intervals - 1) * (max(ry_ich_base_range) - ry_ich_base)) / max(ry_ich_base_range))

# plot

ggplot(results_df, aes(x = efficacy, y = ry_ich_base)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lecanemab efficacy and ICH base rate (65-74)",
    x = "Lecanemab efficacy \n(slows cognitive decline by x points over 18 mos)",
    y = "ICH base rate",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

Lec ICH hazard and apixaban x Lec ICH hazard

lecHR_ACxLec_results <- readRDS('./Results/sensitivityanalysis/lecHR_ACxLec_results.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecHR_ACxLec_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecHR_ACxLec_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecHR_ACxLec_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(hr_ich_lec_range, ac_lec_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                              results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                              results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                              results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('hr_ich_lec','ac_lec_interaction','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$hr_ich_lec_range <- as.factor(results_df$hr_ich_lec)
results_df$ac_lec_range <- as.factor(results_df$ac_lec_interaction)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * (hr_ich_lec - min(hr_ich_lec_range)) / (max(hr_ich_lec_range) - min(hr_ich_lec_range)) + 1,
                       y = (intervals - 1) * (ac_lec_interaction - min(ac_lec_range)) / (max(ac_lec_range) - min(ac_lec_range)) + 1)

# plot

ggplot(results_df, aes(x = hr_ich_lec, y = ac_lec_interaction)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lec hazard and AC x lec interaction hazard on ICH",
    x = "Lecanemab ICH hazard ratio",
    y = "Apixaban x lec interaction hazard on ICH",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

lecHR_ACxLec_results <- readRDS('./Results/sensitivityanalysis/lecHR_ACxLec_results_65_74.rds')

# Combine the optimal strategies and their corresponding QALEs into a new array
combined_results <- array(dim = c(intervals, intervals, 3))

# function to grab the second highest QALE value
second_highest <- function(x) {
  sort(x, decreasing = TRUE)[2]  # Sort in decreasing order and take the second value
} 

# Find the optimal strategy for each parameter combination
combined_results[, , 1] <- apply(lecHR_ACxLec_results[,,'qale',], c(1, 2), which.max) 
combined_results[, , 2] <- apply(lecHR_ACxLec_results[,,'qale',], c(1, 2), second_highest) 


# Reshape the array to have dimensions for each parameter combination and strategies
combined_results[, , 3] <- apply(lecHR_ACxLec_results, c(1, 2), max)

# Expand the grid of parameter values
params_grid <- expand.grid(hr_ich_lec_range, ac_lec_range)

# Flatten the combined results to match the grid
optimal_strategy_flat <- as.vector(combined_results[, , 1])
second_strategy_flat <- as.vector(combined_results[, , 2])
max_QALE_flat <- as.vector(combined_results[, , 3])

# Combine into a data frame
results_df <- cbind(params_grid, optimal_strategy = optimal_strategy_flat, max_QALE = max_QALE_flat, second_best_QALE  = second_strategy_flat)

results_df$optimal_strategy_name <- case_when(results_df$optimal_strategy == 1 ~ 'no lec no ac',
                                              results_df$optimal_strategy == 2 ~ 'no lec yes ac',
                                              results_df$optimal_strategy == 3 ~ 'yes lec no ac',
                                              results_df$optimal_strategy == 4 ~ 'yes lec yes ac')

colnames(results_df) <- c('hr_ich_lec','ac_lec_interaction','optimal_strategy','max_QALE','second_best_QALE','optimal_strategy_name')

# 2-way sensitivity analysis tile plot

results_df$hr_ich_lec_range <- as.factor(results_df$hr_ich_lec)
results_df$ac_lec_range <- as.factor(results_df$ac_lec_interaction)
results_df$max_QALE <- results_df$max_QALE * 365 / 12
results_df$second_best_QALE <- results_df$second_best_QALE * 365 / 12

# Create the heatmap

# basecase point
source('./Scripts/parameters.R')
basecase <- data.frame(x = (intervals - 1) * (hr_ich_lec - min(hr_ich_lec_range)) / (max(hr_ich_lec_range) - min(hr_ich_lec_range)) + 1,
                       y = (intervals - 1) * (ac_lec_interaction - min(ac_lec_range)) / (max(ac_lec_range) - min(ac_lec_range)) + 1)

# plot

ggplot(results_df, aes(x = hr_ich_lec, y = ac_lec_interaction)) +
    geom_tile(aes(fill = factor(optimal_strategy_name), alpha  = max_QALE - second_best_QALE)) +
    scale_fill_manual(values = c("no lec no ac" = "tomato3",
                                 "no lec yes ac" = "skyblue",
                                 "yes lec no ac" = "green",
                                 "yes lec yes ac" = "orange")) +  # Assign red and blue to strategies
    scale_alpha_continuous(range = c(0.4, 1)) +  # Adjust alpha range for intensity variation
    geom_point(data = basecase, aes(x = x, y = y), shape = 4) +
    labs(title = "Sensitivity to lec hazard and AC x lec interaction hazard on ICH (65-74)",
    x = "Lecanemab ICH hazard ratio",
    y = "Apixaban x lec interaction hazard on ICH",
    fill = "Optimal Strategy",
    alpha = "Difference compared \nto second-best strategy \n(days in perfect health)") +
    theme_minimal()

Lec ICH hazard ratio x ICH base rate

pending