Hi everyone,
I'm currently learning how to use the heemod package in R for health economic modeling, and I'm working through one of the official vignettes to practice building a Markov model. The example compares two treatment strategies for HIV: Monotherapy and Combined Therapy.
After running the model and performing a probabilistic sensitivity analysis (PSA), I exported the results to BCEA using run_bcea() to explore additional cost-effectiveness analyses like CEACs and EVPIs.
However, I’ve encountered a confusing issue:
While the reference strategy is correctly assigned, the costs and effects from the PSA appear to be swapped between strategies in the BCEA object. That is, the values that should belong to Monotherapy are assigned to Combined Therapy, and vice versa.
This leads to misleading results in the BCEA plots and summaries, as the strategies are effectively reversed in terms of their outcomes.
You can verify this by checking the last few lines of the code, where I inspect the PSA results and the BCEA object:
View(sensibilidad_modelo$psa) # heemod PSA results
View(exportar_bcea$c) # BCEA costs
View(exportar_bcea$e) # BCEA effects
To make it easier to reproduce the issue, I’ve attached the full code I’m using — it should run without problems.
# Loading libraries:
library(heemod)
library(BCEA)
# Defining model parameters:
param <- define_parameters(
cost_zido = 2278,
cost_lami = 2086,
drug_cost_monotherapy = cost_zido,
drug_cost_combined_therapy = cost_zido + cost_lami,
other_costs_A = 2756,
other_costs_B = 3052,
other_costs_C = 9007,
other_costs_D = 0,
rr = 0.509,
p_AA_mono = 0.721,
p_AB_mono = 0.202,
p_AC_mono = 0.067,
p_AD_mono = 0.010,
p_BC_mono = 0.407,
p_BD_mono = 0.012,
p_CD_mono = 0.250,
p_AB_comb = p_AB_mono * rr,
p_AC_comb = p_AC_mono * rr,
p_AD_comb = p_AD_mono * rr,
p_BC_comb = p_BC_mono * rr,
p_BD_comb = p_BD_mono * rr,
p_CD_comb = p_CD_mono * rr,
p_AA_comb = 1 - (p_AB_comb + p_AC_comb + p_AD_comb),
utility_alive = 1,
utility_death = 0
)
# Defining transition matrix:
Transition_Matrix_Monotherapy <- define_transition(
state_names = c("A", "B", "C", "D"),
p_AA_mono, p_AB_mono, p_AC_mono, p_AD_mono,
0, C, p_BC_mono, p_BD_mono,
0, 0, C, p_CD_mono,
0, 0, 0, 1
)
Transition_Matrix_Combined_Therapy <- define_transition(
state_names = c("A", "B", "C", "D"),
p_AA_comb, p_AB_comb, p_AC_comb, p_AD_comb,
0, C, p_BC_comb, p_BD_comb,
0, 0, C, p_CD_comb,
0, 0, 0, 1
)
# Defining states:
State_A <- define_state(
drug_costs = dispatch_strategy(
Monotherapy = discount(drug_cost_monotherapy, 0.06),
Combined_Therapy = discount(drug_cost_combined_therapy, 0.06)
),
other_costs = discount(other_costs_A, 0.06),
total_costs = drug_costs + other_costs,
utility = utility_alive
)
State_B <- define_state(
drug_costs = dispatch_strategy(
Monotherapy = discount(drug_cost_monotherapy, 0.06),
Combined_Therapy = discount(drug_cost_combined_therapy, 0.06)
),
other_costs = discount(other_costs_B, 0.06),
total_costs = drug_costs + other_costs,
utility = utility_alive
)
State_C <- define_state(
drug_costs = dispatch_strategy(
Monotherapy = discount(drug_cost_monotherapy, 0.06),
Combined_Therapy = discount(drug_cost_combined_therapy, 0.06)
),
other_costs = discount(other_costs_C, 0.06),
total_costs = drug_costs + other_costs,
utility = utility_alive
)
State_D <- define_state(
drug_costs = 0,
other_costs = discount(other_costs_D, 0.06),
total_costs = drug_costs + other_costs,
utility = utility_death
)
# Defining strategies:
Monotherapy <- define_strategy(
transition = Transition_Matrix_Monotherapy,
A = State_A,
B = State_B,
C = State_C,
D = State_D
)
Combined_Therapy <- define_strategy(
transition = Transition_Matrix_Combined_Therapy,
A = State_A,
B = State_B,
C = State_C,
D = State_D
)
# Defining model:
modelo_sida <- run_model(
parameters = param,
Monotherapy = Monotherapy,
Combined_Therapy = Combined_Therapy,
cycles = 50,
cost = total_costs,
effect = utility
)
summary(modelo_sida)
# Defining probabilistic analysis:
sensibilidad_prob <- define_psa(
rr ~ lognormal(mean = 0.509, sdlog = 0.173),
other_costs_A ~ gamma(mean = 2756, sd = sqrt(2756)),
other_costs_B ~ gamma(mean = 3052, sd = sqrt(3052)),
other_costs_C ~ gamma(mean = 9007, sd = sqrt(9007)),
p_CD_mono ~ binomial(prob = 0.25, size = 40),
p_AA_mono + p_AB_mono + p_AC_mono + p_AD_mono ~ multinomial(721, 202, 67, 10)
)
# Running probabilistic analysis:
sensibilidad_modelo <- run_psa(
model = modelo_sida,
psa = sensibilidad_prob,
N = 100
)
# Results of the probabilistic analysis::
summary(sensibilidad_modelo)
# Several plots
plot(sensibilidad_modelo, type = "ce")
plot(sensibilidad_modelo, type = "ac", max_wtp = 10000, log_scale = FALSE)
plot(sensibilidad_modelo, type = "evpi", max_wtp = 10000, log_scale = FALSE)
plot(sensibilidad_modelo, type = "cov", max_wtp = 10000, log_scale = FALSE)
# Export the sensibility analysis into BCEA format:
exportar_bcea <- run_bcea(sensibilidad_modelo,
plot = TRUE,
Kmax = 10000,
ref = "Monotherapy")
# Results in BCEA format:
summary(exportar_bcea)
View(sensibilidad_modelo$psa) # View the costs and effects of the original sensibility model.
View(exportar_bcea$c) # View the costs of the BCEA sensibility model.
View(exportar_bcea$e) # View the effects of the BCEA sensibility model.
Has anyone else experienced this?
Is there a way to ensure that the PSA results are correctly mapped to the corresponding strategies when exporting from heemod to BCEA?
Any help or suggestions would be greatly appreciated!
Thanks in advance