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

from heemod.

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")

from heemod.

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

from heemod.

Github's magic ;)

from heemod.