Sensitivity analysis - plotting dsa
vennela88 opened this issue · 4 comments
This is my model
I keep getting this error and I am not sure where the problem is?
Can someone please help me?
Hi!
Could you please create a small version of your code leading to this issue? That would allow me to fix the bug.
Best
Kevin
Hi,
Thanks for responding. This is my first model for HE. I don't know what a small version is but I have pasted the full version in case it makes it easy to understand.
😄 Vennela
1. define transitions ----
mat_aad <- define_transition(
state_names = c("af", "sr","stroke","hf", "death"),
C, p_af_sr_aad, p_af_stroke, p_af_hf, p_death_age,
p_sr_af_aad, C, p_sr_stroke, p_sr_hf, p_death_age,
0, 0, C, 0, p_stroke_dead,
0, 0, 0, C, p_hf_dead,
0, 0, 0, C, 1)
mat_ca <- define_transition(
state_names = c("af", "sr","stroke","hf", "death"),
C, p_af_sr_ca, p_af_stroke, p_af_hf, p_death_age,
p_sr_af, C, p_sr_stroke, p_sr_hf, p_death_age,
0, 0, C, 0, p_stroke_dead,
0, 0, 0, C, p_hf_dead,
0, 0, 0, C, 1)
mat_sa <- define_transition(
state_names = c("af", "sr","stroke","hf", "death"),
C, p_af_sr_sa, p_af_stroke, p_af_hf, p_death_age,
p_sr_af, C, p_sr_stroke, p_sr_hf, p_death_age,
0, 0, C, 0, p_stroke_dead,
0, 0, 0, C, p_hf_dead,
0, 0, 0, C, 1)
plot(mat_aad)
aad_sr_af_km <- read_excel("aad_sr_af_km.xlsx")
p_aad_sr_af_surv <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "weibull",
data = aad_sr_af_km)
plot(p_aad_sr_af_surv)
ca_sr_af_km <- read_excel("ca_sr_af_km.xlsx")
p_ca_sr_af_surv <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "weibull",
data = ca_sr_af_km)
plot(p_ca_sr_af_surv)
par_mod <- define_parameters(
age_base = 62,
age_cycle = model_time + age_base,
base_case = 0.30,
rrCA = 1.84,
rrSA = 2.26,
p_af_sr_aad = ifelse(markov_cycle ==1, base_case, 0),
p_af_sr_ca = ifelse(markov_cycle==1, base_caserrCA, 0),
p_af_sr_sa = ifelse(markov_cycle==1, base_caserrSA, 0),
p_sr_af_aad = ifelse(markov_cycle >= 2, (compute_surv(p_aad_sr_af_surv,
time = state_time,
km_limit = NULL)), 0),
p_sr_af = ifelse(markov_cycle >= 2, (compute_surv(p_ca_sr_af_surv,
time = state_time,
km_limit = NULL)), 0),
p_af_hf = 0.027,
p_sr_hf = 0.022,
p_af_stroke = 0.013,
p_sr_stroke = 0.009,
p_af_death = 0.025,
p_sr_death = 0.021)
hf_dead <- read_excel("hf_death_km.xlsx")
p_hf_dead_surv <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "weibull",
data = hf_dead)
plot(p_hf_dead_surv)
par_mod <- modify(par_mod,
p_aad_sr_af = compute_surv(p_aad_sr_af_surv,
time = state_time,
km_limit = NULL),
p_ca_sr_af = compute_surv(p_ca_sr_af_surv,
time = state_time,
km_limit = NULL),
p_death_age = get_who_mr(
age= age_cycle,
sex = NULL,
country = "GBR",
local = TRUE), # death based on age
stroke_mortality = 0.122,
p_stroke_dead = ifelse(state_time ==1, stroke_mortality, p_death_age),
p_hf_dead =
compute_surv(p_hf_dead_surv,
time = state_time,
km_limit = NULL))
par_mod <- modify(par_mod,
niv = 496.07,
o2_ass = 378.38,
sleep_study = 183,
lung_ft = 438.95,
hemi_para = (3niv+o2_ass+2sleep_study+lung_ft),
lrti = 1916.56,
pleu_eff = 2925.01,
ppm = 4950.09,
dccv = 2253.04,
lrti_ca = lrti*(4/60),
dccv_ca = dccv*(19/60),
comp_ca = lrti_ca+dccv_ca,
lrti_sa = lrti*(2/55),
pleu_eff_sa = pleu_eff*(1/55),
ppm_sa = ppm*(1/55),
hemi_para_sa = (hemi_para*(2/55)),
dccv_sa = dccv*(15/55),
comp_sa = lrti_sa+pleu_eff_sa+ppm_sa+hemi_para_sa+dccv_sa,
cost_sa = 19776.85,
cost_sa_cycle = ifelse(markov_cycle == 1, cost_sa+comp_sa, 0),
cost_ca= 14469.15,
cost_ca_cycle = ifelse(markov_cycle == 1, cost_ca+comp_ca, 0),
amio = 8,
amio_1yr = (amio*365),
sota = 5,
sota_1yr = (sota*2*365),
flec = 5,
flec_1yr = (flec*2*365),
drone = 113,
drone_1yr = (drone*2*365),
cost_aad = (0.6*amio_1yr)+(0.2*sota_1yr)+(0.15*flec_1yr)+(0.05*drone_1yr),
cost_aad_cycle= ifelse(markov_cycle == 1, cost_aad, 0),
cost_ae = 342,
cost_nonelective_adm = 2253,
cost_hosp_fastaf = cost_ae+cost_nonelective_adm,
cost_opd = 125,
cost_hospit_cycle_af = ifelse(state_time <=3, cost_hosp_fastaf+(cost_opd*2), cost_opd),
cost_hospit_start_stroke = 17063,
cost_hospit_end_stroke = 7225,
n_years_stroke = 2,
cost_hospit_cycle_stroke = ifelse(
state_time < n_years_stroke,
cost_hospit_start_stroke,
cost_hospit_end_stroke),
clopidogrel = 5,
clopi_1yr=(clopidogrel*365),
cost_med_stroke = clopi_1yr,
i2021 = 0.021,
i2020 = 0.011,
i2019 = 0.02,
i2018 = 0.023,
i2017 = 0.026,
i2016 = 0.009,
costhf = 13055,
cost_hospit_start_hf = (costhf+costhf*i2021+costhf*i2020+
costhf*i2019+costhf*i2018+costhf*i2017+costhf*i2016),
cost_hospit_end_hf = 3901,
n_years_hf = 2,
cost_hospit_cycle_hf = ifelse(
state_time < n_years_hf,
cost_hospit_start_hf,
(cost_hospit_end_hf*2)),
ace_arb = 2337,
beta = 1359,
aldo_a= 6603,
diuretics = 1395,
lcz = 23871,
sglti = 4770,
cost_med_hf= (ace_arb+beta+aldo_a+diuretics+lcz+sglti),
casa_qaly_bl = 0.73,
#UpersAF = 0.84,
#UpermAF = 0.80,
#qalyAF =ifelse(state_time ==1, UpersAF, UpermAF),
dr = 0.03)
state_af <- define_state(
cost_treat = dispatch_strategy(
aad = cost_aad_cycle,
ca = cost_ca_cycle,
sa = cost_sa_cycle),
cost_hospit = cost_hospit_cycle_af,
cost_med = 0,
cost_total = discount(cost_treat + cost_hospit + cost_med, r = dr),
qaly = 0.67)
state_sr <- define_state(
cost_treat = 0,
cost_hospit = 0,
cost_med = 0,
cost_total = discount(cost_treat + cost_hospit + cost_med, r = dr),
qaly = 0.75)
state_stroke <- define_state(
cost_treat = 0,
cost_hospit = cost_hospit_cycle_stroke,
cost_med = cost_med_stroke, # need to add cost of OPD for stroke patients
cost_total = discount(cost_treat + cost_hospit + cost_med, r = dr),
qaly = 0.498)
state_hf <- define_state(
cost_treat = 0,
cost_hospit = cost_hospit_cycle_hf,
cost_med = cost_med_hf,
cost_total = discount(cost_treat + cost_hospit + cost_med, r = dr),
qaly = 0.6)
state_death <- define_state(
cost_treat = 0,
cost_hospit = 0,
cost_med = 0,
cost_total = 0,
qaly = 0)
strat_aad <- define_strategy(
transition = mat_aad,
af=state_af,
sr=state_sr,
stroke=state_stroke,
hf=state_hf,
death=state_death)
strat_ca <- define_strategy(
transition = mat_ca,
af=state_af,
sr=state_sr,
stroke=state_stroke,
hf=state_hf,
death=state_death)
strat_sa <- define_strategy(
transition = mat_sa,
af=state_af,
sr=state_sr,
stroke=state_stroke,
hf=state_hf,
death=state_death)
res_mod <- run_model(
parameters = par_mod,
aad = strat_aad,
ca = strat_ca,
sa = strat_sa,
cycles = 3,
cost = cost_total,
effect = qaly,
method="beginning")
res_mod
summary(res_mod, threshold = c(10000, 20000, 30000, 40000))
plot(res_mod)
plot(res_mod, model="all")
plot(res_mod, model="all", panels="by_state")
def_dsa <- define_dsa(
age_base, 52, 72,
base_case, 0.25, 0.35,
p_af_hf, 0.012, 0.042,
p_sr_hf, 0.005, 0.039,
p_af_stroke, 0.012, 0.014,
p_sr_stroke, 0.007, 0.011,
p_af_death, 0.020, 0.030,
p_sr_death, 0.016, 0.026,
#shape, 1.4, 1.6,#shape,
#scale, 4, 6, # scale
cost_ca, 14469.15, 16639.52, #15% of patients had redo in 1st year, cost_ca + (0.15cost_ca)
cost_sa, 19776.95, 22941.15, #16% of patients had redo in 1st year, cost_sa + (0.16cost_sa)
dr, 0.025, 0.035,
n_years_stroke, 1, 3,
n_years_hf, 1, 3)
res_dsa <- run_dsa(model =res_mod, dsa = def_dsa)
res_dsa
plot(res_dsa, strategy = "ca",
result = "icer",
type = "difference")
It worked after not working all night and all morning! I am not sure how it was fixed, as not changed code.
Thank you for looking at my query though.
Vennela
Github's magic ;)