In sheejamk/MarkovModel: Decision Analysis Markov Model
In this example we will model the cost effectiveness of lamivudine/zidovudine
combination therapy in HIV infection ([Chancellor, 1997] (https://www.ncbi.nlm.nih.gov/pubmed/10169387))) and also described in Decision Modelling for Health Economic Evaluation, page 32. We have implemented the scenario in which the effect of combination therapy of lamivudine/zidovudine (as the probbaility of moving to worse heatlh states and having to accrue more cost) is only seen in first 2 years. This means that the transition probability and costs are dependent on first two cycles (As the cycle length is years).
This model aims to compare costs and utilities of two treatment
strategies, monotherapy and combined therapy.
Four states are described, from best to worst health-wise:
A: CD4 cells > 200 and < 500 cells/mm3;
B: CD4 < 200 cells/mm3, non-AIDS;
C: AIDS;
D: Death.
knitr::opts_chunk$set(
collapse = TRUE,
comment = ">"
)
library(MarkovModel)
Define health states for monotherapy - now the costs are just defined, but not defined any value.
A <- health_state("A", cost="cost_health_A+ cost_drug ",utility=1)
B <- health_state("B", cost="cost_health_B + cost_drug",utility=1)
C <- health_state("C", cost="cost_health_C + cost_drug",utility=1)
D <- health_state("D", cost=0,utility=0)
Define allowed transition probabilities and number them. The below matrix is numbered so that the maximum entry in the matrix gives the total number of allowed transitions. Column names and row names are just the names of the health states.
tmat <- rbind(c(1, 2,3,4), c(NA, 5,6,7),c(NA, NA, 8,9), c(NA,NA,NA,10))
colnames(tmat) <- rownames(tmat) <- c("A","B" ,"C","D")
Define transition matrix now, but parameters are fine now
tm <- transition_matrix(4, tmat, c("tpAtoA","tpAtoB","tpAtoC","tpAtoD",
"tpBtoB", "tpBtoC", "tpBtoD",
"tpCtoC","tpCtoD","tpDtoD" ), colnames(tmat) )
Combine the health states and define the strategy. The current strategy ie. control or intervention - here it is monotherapy
health_states <- combine_state(A,B,C,D)
mono_strategy <- strategy(tm, health_states, "mono")
Before we run the model, we need to give values to parameters, thus we define the parameter list. We need to make sure that the parameters with numerical values are given first and derived parameters later in the list. The parameters are assinged sequentially as given in the parameter list. So if there is any calcultions needed (or using functions), they are defined earlier.
mono_params=define_parameters(cost_zido = 2278,
cost_direct_med_A = 1701,
cost_comm_care_A = 1055,
cost_direct_med_B = 1774,
cost_comm_care_B = 1278,
cost_direct_med_C = 6948,
cost_comm_care_C = 2059,
tpAtoA = 1251/(1251+483),
tpAtoB = 350/(350+1384),
tpAtoC = 116/(116+1618),
tpAtoD = 17/(17+1717),
tpBtoB = 731/(731+527),
tpBtoC = 512/(512+746),
tpBtoD = 15/(15+1243),
tpCtoC = 1312/(1312+437),
tpCtoD = 437/(437+1312),
tpDtoD = 1,
cost_health_A = "cost_direct_med_A+ cost_comm_care_A",
cost_health_B = "cost_direct_med_B+ cost_comm_care_B",
cost_health_C = "cost_direct_med_C+ cost_comm_care_C",
cost_drug = "cost_zido")
We run the model with strategy mono_strategy for 20 cycles with initial state values 1,0,0,0 corresponding to A,B, C and D. Costs and qualys are discounted follwoing the values given in discount. Parameter values will have any parameter defined as health state values or tranisiton probabilties, but not their values were assgined
mono_markov <-markov_model(mono_strategy, 20,c(1,0,0,0),c(0,0,0,0),discount=c(0.06,0),mono_params)
Now for combination therapy, the transition probabilites and cost are dependent on the years or the cylce number. Thus we need to define functions that given the correct output to calculate the probabilites and costs as a function of cycle
#Define function to set the cost to be differnt for first two cycles
define_comb_cost=function(cycle,cost_lami){
if(cycle==2 || cycle ==3)
return(cost_lami)
else
return(0)
}
#Define function to set the risk ratio to be differnt for first two cycles
define_rr=function(cycle,rr){
if(cycle==2 || cycle ==3)
return(rr)
else
return(1)
}
Define health states - now the costs are just defined, but not defined any value.
A <- health_state("A", cost="cost_health_A + cost_drug",utility=1)
B <- health_state("B", cost="cost_health_B + cost_drug",utility=1)
C <- health_state("C", cost="cost_health_C + cost_drug",utility=1)
D <- health_state("D", cost=0,utility=0)
Define allowed transition probabilities and number them. The below matrix is numbered so that the maximum entry in the matrix gives the total number of allowed transitions. Column names and row names are just the names of the health states.
tmat <- rbind(c(1, 2,3,4), c(NA, 5,6,7),c(NA, NA, 8,9), c(NA,NA,NA,10))
colnames(tmat) <- rownames(tmat) <- c("A","B" ,"C","D")
Define transition matrix now, but parameters are fine now
tm <- transition_matrix(4, tmat, c("tpAtoA_rr","tpAtoB_rr","tpAtoC_rr","tpAtoD_rr",
"tpBtoB_rr", "tpBtoC_rr", "tpBtoD_rr",
"tpCtoC_rr","tpCtoD_rr","tpDtoD_rr" ), colnames(tmat) )
Before we run the model, we need to give values to parameters, thus we define the parameter list.
comb_params<-define_parameters(cost_zido = 2278,
cost_direct_med_A = 1701,
cost_comm_care_A = 1055,
cost_direct_med_B = 1774,
cost_comm_care_B = 1278,
cost_direct_med_C = 6948,
cost_comm_care_C = 2059,
tpAtoA = 1251/(1251+483),
tpAtoB = 350/(350+1384),
tpAtoC = 116/(116+1618),
tpAtoD = 17/(17+1717),
tpBtoB = 731/(731+527),
tpBtoC = 512/(512+746),
tpBtoD = 15/(15+1243),
tpCtoC = 1312/(1312+437),
tpCtoD = 437/(437+1312),
tpDtoD = 1,
rr=0.509,
cost_lami = 2086.50,
rr_cycle="define_rr(markov_cycle,rr)",
tpAtoA_rr="1-tpAtoB*rr_cycle-tpAtoC*rr_cycle-tpAtoD*rr_cycle",
tpAtoB_rr="tpAtoB*rr_cycle",
tpAtoC_rr="tpAtoC*rr_cycle",
tpAtoD_rr="tpAtoD*rr_cycle",
tpBtoB_rr="1-tpBtoC*rr_cycle-tpBtoD*rr_cycle",
tpBtoC_rr="tpBtoC*rr_cycle",
tpBtoD_rr="tpBtoD*rr_cycle",
tpCtoC_rr="1-tpCtoD*rr_cycle",
tpCtoD_rr="tpCtoD*rr_cycle",
tpDtoD_rr=1,
cost_health_A = "cost_direct_med_A + cost_comm_care_A",
cost_health_B = "cost_direct_med_B + cost_comm_care_B",
cost_health_C = "cost_direct_med_C + cost_comm_care_C",
cost_lami_cycle = "define_comb_cost(markov_cycle,cost_lami)",
cost_drug = "cost_zido + cost_lami_cycle")
# Combine the health states
health_states <- combine_state(A,B,C,D)
#The current strategy ie. control or intervention - here it is combination therapy
comb_strategy <- strategy(tm, health_states, "comb")
We run the model with strategy mono_strategy for 20 cycles with initial state values 1,0,0,0 corresponding to A,B, C and D. Costs and qualys are discounted follwoing the values given in discount. Parameter values will have any parameter defined as health state values or tranisiton probabilties, but not their values were assgined
comb_markov <-markov_model(comb_strategy, 20, c(1, 0,0,0),c(0,0,0,0),discount=c(0.06,0.0),comb_params)
Now we will define the parameters for the sensitivity analysis. The parameters are already defined as mono_params and comb_params. Their lower value and upper value are defined using definer_parameters() and create parameter table using define_parameters_sens_anal(). Then deteministic sensitvity analysis is done for the markov model of the choice here the mono_markov using the function do_sensitivity_analysis. Using the result obtained, a full report or plot (for specfic variable) can be generated using report_sensitivity_analysis() and plot_dsa() respectively. To get a report or plot of senstivity analysis on ICER on NMB, two markov models need to be passed on.
min_values<-define_parameters(cost_direct_med_B = 177.4,cost_comm_care_C = 205.9)
max_values<-define_parameters(cost_direct_med_B = 17740,cost_comm_care_C = 20590)
param_table<-define_parameters_sens_anal(mono_params, min_values, max_values)
result_dsa_control<-do_sensitivity_analysis(mono_markov,param_table)
report_sensitivity_analysis(result_dsa_control)
plot_dsa(result_dsa_control,"cost")
plot_dsa(result_dsa_control,"cost_direct_med_B",type="range")
plot_dsa(result_dsa_control,"utility",type="range")
If we generate the parameter table and deterministic sensitivityanalysis for markov model for combination treatment, we can use to generate the report and plot for both ICER and NMB. If we need
min_values<-define_parameters(cost_direct_med_B = 177.4,cost_comm_care_C = 205.9)
max_values<-define_parameters(cost_direct_med_B = 17740,cost_comm_care_C = 20590)
param_table<-define_parameters_sens_anal(comb_params, min_values, max_values)
result_dsa_treat<-do_sensitivity_analysis(comb_markov,param_table)
report_sensitivity_analysis(result_dsa_control,result_dsa_treat,1000,"mono")
plot_dsa(result_dsa_control,"cost",type="difference",result_dsa_treat,threshold=1000, comparator="mono")
plot_dsa(result_dsa_control,"NMB",type="range",result_dsa_treat,threshold=1000, comparator="mono")
```
sheejamk/MarkovModel documentation built on Jan. 23, 2020, 2:44 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
In this example we will model the cost effectiveness of lamivudine/zidovudine combination therapy in HIV infection ([Chancellor, 1997] (https://www.ncbi.nlm.nih.gov/pubmed/10169387))) and also described in Decision Modelling for Health Economic Evaluation, page 32. We have implemented the scenario in which the effect of combination therapy of lamivudine/zidovudine (as the probbaility of moving to worse heatlh states and having to accrue more cost) is only seen in first 2 years. This means that the transition probability and costs are dependent on first two cycles (As the cycle length is years).
This model aims to compare costs and utilities of two treatment strategies, monotherapy and combined therapy.
Four states are described, from best to worst health-wise:
A: CD4 cells > 200 and < 500 cells/mm3; B: CD4 < 200 cells/mm3, non-AIDS; C: AIDS; D: Death.
knitr::opts_chunk$set( collapse = TRUE, comment = ">" )
library(MarkovModel)
Define health states for monotherapy - now the costs are just defined, but not defined any value.
A <- health_state("A", cost="cost_health_A+ cost_drug ",utility=1) B <- health_state("B", cost="cost_health_B + cost_drug",utility=1) C <- health_state("C", cost="cost_health_C + cost_drug",utility=1) D <- health_state("D", cost=0,utility=0)
Define allowed transition probabilities and number them. The below matrix is numbered so that the maximum entry in the matrix gives the total number of allowed transitions. Column names and row names are just the names of the health states.
tmat <- rbind(c(1, 2,3,4), c(NA, 5,6,7),c(NA, NA, 8,9), c(NA,NA,NA,10)) colnames(tmat) <- rownames(tmat) <- c("A","B" ,"C","D")
Define transition matrix now, but parameters are fine now
tm <- transition_matrix(4, tmat, c("tpAtoA","tpAtoB","tpAtoC","tpAtoD", "tpBtoB", "tpBtoC", "tpBtoD", "tpCtoC","tpCtoD","tpDtoD" ), colnames(tmat) )
Combine the health states and define the strategy. The current strategy ie. control or intervention - here it is monotherapy
health_states <- combine_state(A,B,C,D) mono_strategy <- strategy(tm, health_states, "mono")
Before we run the model, we need to give values to parameters, thus we define the parameter list. We need to make sure that the parameters with numerical values are given first and derived parameters later in the list. The parameters are assinged sequentially as given in the parameter list. So if there is any calcultions needed (or using functions), they are defined earlier.
mono_params=define_parameters(cost_zido = 2278, cost_direct_med_A = 1701, cost_comm_care_A = 1055, cost_direct_med_B = 1774, cost_comm_care_B = 1278, cost_direct_med_C = 6948, cost_comm_care_C = 2059, tpAtoA = 1251/(1251+483), tpAtoB = 350/(350+1384), tpAtoC = 116/(116+1618), tpAtoD = 17/(17+1717), tpBtoB = 731/(731+527), tpBtoC = 512/(512+746), tpBtoD = 15/(15+1243), tpCtoC = 1312/(1312+437), tpCtoD = 437/(437+1312), tpDtoD = 1, cost_health_A = "cost_direct_med_A+ cost_comm_care_A", cost_health_B = "cost_direct_med_B+ cost_comm_care_B", cost_health_C = "cost_direct_med_C+ cost_comm_care_C", cost_drug = "cost_zido")
We run the model with strategy mono_strategy for 20 cycles with initial state values 1,0,0,0 corresponding to A,B, C and D. Costs and qualys are discounted follwoing the values given in discount. Parameter values will have any parameter defined as health state values or tranisiton probabilties, but not their values were assgined
mono_markov <-markov_model(mono_strategy, 20,c(1,0,0,0),c(0,0,0,0),discount=c(0.06,0),mono_params)
Now for combination therapy, the transition probabilites and cost are dependent on the years or the cylce number. Thus we need to define functions that given the correct output to calculate the probabilites and costs as a function of cycle
#Define function to set the cost to be differnt for first two cycles define_comb_cost=function(cycle,cost_lami){ if(cycle==2 || cycle ==3) return(cost_lami) else return(0) } #Define function to set the risk ratio to be differnt for first two cycles define_rr=function(cycle,rr){ if(cycle==2 || cycle ==3) return(rr) else return(1) }
Define health states - now the costs are just defined, but not defined any value.
A <- health_state("A", cost="cost_health_A + cost_drug",utility=1) B <- health_state("B", cost="cost_health_B + cost_drug",utility=1) C <- health_state("C", cost="cost_health_C + cost_drug",utility=1) D <- health_state("D", cost=0,utility=0)
Define allowed transition probabilities and number them. The below matrix is numbered so that the maximum entry in the matrix gives the total number of allowed transitions. Column names and row names are just the names of the health states.
tmat <- rbind(c(1, 2,3,4), c(NA, 5,6,7),c(NA, NA, 8,9), c(NA,NA,NA,10)) colnames(tmat) <- rownames(tmat) <- c("A","B" ,"C","D")
Define transition matrix now, but parameters are fine now
tm <- transition_matrix(4, tmat, c("tpAtoA_rr","tpAtoB_rr","tpAtoC_rr","tpAtoD_rr", "tpBtoB_rr", "tpBtoC_rr", "tpBtoD_rr", "tpCtoC_rr","tpCtoD_rr","tpDtoD_rr" ), colnames(tmat) )
Before we run the model, we need to give values to parameters, thus we define the parameter list.
comb_params<-define_parameters(cost_zido = 2278, cost_direct_med_A = 1701, cost_comm_care_A = 1055, cost_direct_med_B = 1774, cost_comm_care_B = 1278, cost_direct_med_C = 6948, cost_comm_care_C = 2059, tpAtoA = 1251/(1251+483), tpAtoB = 350/(350+1384), tpAtoC = 116/(116+1618), tpAtoD = 17/(17+1717), tpBtoB = 731/(731+527), tpBtoC = 512/(512+746), tpBtoD = 15/(15+1243), tpCtoC = 1312/(1312+437), tpCtoD = 437/(437+1312), tpDtoD = 1, rr=0.509, cost_lami = 2086.50, rr_cycle="define_rr(markov_cycle,rr)", tpAtoA_rr="1-tpAtoB*rr_cycle-tpAtoC*rr_cycle-tpAtoD*rr_cycle", tpAtoB_rr="tpAtoB*rr_cycle", tpAtoC_rr="tpAtoC*rr_cycle", tpAtoD_rr="tpAtoD*rr_cycle", tpBtoB_rr="1-tpBtoC*rr_cycle-tpBtoD*rr_cycle", tpBtoC_rr="tpBtoC*rr_cycle", tpBtoD_rr="tpBtoD*rr_cycle", tpCtoC_rr="1-tpCtoD*rr_cycle", tpCtoD_rr="tpCtoD*rr_cycle", tpDtoD_rr=1, cost_health_A = "cost_direct_med_A + cost_comm_care_A", cost_health_B = "cost_direct_med_B + cost_comm_care_B", cost_health_C = "cost_direct_med_C + cost_comm_care_C", cost_lami_cycle = "define_comb_cost(markov_cycle,cost_lami)", cost_drug = "cost_zido + cost_lami_cycle")
# Combine the health states health_states <- combine_state(A,B,C,D) #The current strategy ie. control or intervention - here it is combination therapy comb_strategy <- strategy(tm, health_states, "comb")
We run the model with strategy mono_strategy for 20 cycles with initial state values 1,0,0,0 corresponding to A,B, C and D. Costs and qualys are discounted follwoing the values given in discount. Parameter values will have any parameter defined as health state values or tranisiton probabilties, but not their values were assgined
comb_markov <-markov_model(comb_strategy, 20, c(1, 0,0,0),c(0,0,0,0),discount=c(0.06,0.0),comb_params)
Now we will define the parameters for the sensitivity analysis. The parameters are already defined as mono_params and comb_params. Their lower value and upper value are defined using definer_parameters() and create parameter table using define_parameters_sens_anal(). Then deteministic sensitvity analysis is done for the markov model of the choice here the mono_markov using the function do_sensitivity_analysis. Using the result obtained, a full report or plot (for specfic variable) can be generated using report_sensitivity_analysis() and plot_dsa() respectively. To get a report or plot of senstivity analysis on ICER on NMB, two markov models need to be passed on.
min_values<-define_parameters(cost_direct_med_B = 177.4,cost_comm_care_C = 205.9) max_values<-define_parameters(cost_direct_med_B = 17740,cost_comm_care_C = 20590) param_table<-define_parameters_sens_anal(mono_params, min_values, max_values) result_dsa_control<-do_sensitivity_analysis(mono_markov,param_table) report_sensitivity_analysis(result_dsa_control) plot_dsa(result_dsa_control,"cost") plot_dsa(result_dsa_control,"cost_direct_med_B",type="range") plot_dsa(result_dsa_control,"utility",type="range")
If we generate the parameter table and deterministic sensitivityanalysis for markov model for combination treatment, we can use to generate the report and plot for both ICER and NMB. If we need
min_values<-define_parameters(cost_direct_med_B = 177.4,cost_comm_care_C = 205.9) max_values<-define_parameters(cost_direct_med_B = 17740,cost_comm_care_C = 20590) param_table<-define_parameters_sens_anal(comb_params, min_values, max_values) result_dsa_treat<-do_sensitivity_analysis(comb_markov,param_table) report_sensitivity_analysis(result_dsa_control,result_dsa_treat,1000,"mono") plot_dsa(result_dsa_control,"cost",type="difference",result_dsa_treat,threshold=1000, comparator="mono") plot_dsa(result_dsa_control,"NMB",type="range",result_dsa_treat,threshold=1000, comparator="mono")
```
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.