msdes

Introduction

This shows how to run the MS discrete event simulation model in its simplest case. Refer to articles for explanation of the internal workings.

Main options

Load the required packages

library(msdes)
library(magrittr)
library(WARDEN)
library(flexsurv)
#> Loading required package: survival
library(VGAM)
#> Loading required package: stats4
#> Loading required package: splines
#> 
#> Attaching package: 'VGAM'
#> The following objects are masked from 'package:flexsurv':
#> 
#>     AICc, dgompertz, pgompertz, qgompertz, rgompertz
library(MASS)
library(utils)
library(hesim)
library(purrr)
#> 
#> Attaching package: 'purrr'
#> The following object is masked from 'package:magrittr':
#> 
#>     set_names
library(rlang)
#> 
#> Attaching package: 'rlang'
#> The following objects are masked from 'package:purrr':
#> 
#>     flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
#>     flatten_raw, invoke, splice
#> The following object is masked from 'package:magrittr':
#> 
#>     set_names
library(progressr)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:rlang':
#> 
#>     :=
#> The following object is masked from 'package:purrr':
#> 
#>     transpose
#> The following object is masked from 'package:base':
#> 
#>     %notin%
library(BCEA)
#> Registered S3 method overwritten by 'BCEA':
#>   method     from
#>   plot.evppi voi
#> The BCEA version loaded is: 2.4.83 (stable version)
#> A development version is available at https://giabaio.r-universe.dev/BCEA
#> 
#> Attaching package: 'BCEA'
#> The following object is masked from 'package:graphics':
#> 
#>     contour
library(voi)
#> 
#> Attaching package: 'voi'
#> The following object is masked from 'package:BCEA':
#> 
#>     evppi

model setup

Data is loaded with the function [rrmsdessetup()] below, all the data is self contained within this function. The user has the option to change a few of the inputs.

rrmsdessetup(  # change model settings or assumptions used in scenario analyses
  drc= 0.035, #3.5% discount rate for costs
  drq= 0.035, #3.5% discount rate for qalys
  time_horizon = 70, # number of years model is running for starting at intial age
  p_eid=0.6, # proportion of patients on external interval dosing
  jcvtest_cost= 0 #assumed free due to home service, test cost could be £247 privately but not offered in the NHS
  )

Run Model

Model is run with the function [rrmsdesmodel()] below. The user has the specify the number of patients and number of inputs.

rrmsdesmodel(n_patients=10,n_samples=10)
#> Analysis number: 1
#> Simulation number: 1
#> Time to run simulation 1: 0.59s
#> Simulation number: 2
#> Time to run simulation 2: 0.69s
#> Simulation number: 3
#> Time to run simulation 3: 0.58s
#> Simulation number: 4
#> Time to run simulation 4: 0.58s
#> Simulation number: 5
#> Time to run simulation 5: 0.68s
#> Simulation number: 6
#> Time to run simulation 6: 0.87s
#> Simulation number: 7
#> Time to run simulation 7: 0.57s
#> Simulation number: 8
#> Time to run simulation 8: 0.55s
#> Simulation number: 9
#> Time to run simulation 9: 0.58s
#> Simulation number: 10
#> Time to run simulation 10: 0.64s
#> Time to run analysis 1: 6.34s
#> Total time to run: 6.34s
#> Simulation finalized;

Plot Results

# Summarise the results

dim_results <- c(n_samples, n_treatments, n_patients)
trt_idx <- setNames(seq_along(v_treatment_names), v_treatment_names)
total_costs_ipd <- total_qalys_ipd <-  array(
  NA_real_,
  dim = dim_results,
  dimnames = list(NULL, v_treatment_names, NULL)
)

for (i_sample in seq_len(n_samples)) {
  
  df <- results[[1]][[i_sample]]$merged_df
  
  # Keep only final ("death") rows
  df <- df[df$evtname == "death", ]
  
  idx <- cbind(
    i_sample,
    trt_idx[df$arm],
    df$pat_id
  )
  
  total_costs_ipd[idx]         <- df$total_costs
  total_qalys_ipd[idx]         <- df$total_qalys
}


options(scipen = 999)
total_costs <- apply(total_costs_ipd, c(1, 2), mean)
total_qalys <- apply(total_qalys_ipd, c(1, 2), mean)


#drop natalizumab-biosimilar-SC and cladrebiene
total_costs <- total_costs[,-c(4,5)]
total_qalys <- total_qalys[,-c(4,5)]

results_table <- summarise_ce_results(wtp1 = 25000, wtp2 = 35000,
                                      total_qalys = total_qalys, 
                                      total_costs = total_costs,
                                      reference = "Natalizumab-IV")

knitr::kable(results_table)

	Total costs (95% CrI)	Total QALYs (95% CrI)	Net benefit at £25000/QALY (95% CrI)	Net benefit at £35000/QALY (95% CrI)	Incremental costs (95% CrI)	Incremental QALYs (95% CrI)	INB at £25000/QALY (95% CrI)	INB at £35000/QALY (95% CrI)	CEAC at £25000/QALY (95% CrI)	CEAC at £35000/QALY (95% CrI)
Natalizumab-IV	351,983.65 (307,662.42, 406,091.90)	9.53 (7.81, 11.37)	-113,842.77 (-174,017.68, -26,843.39)	-18,586.42 (-92,828.03, 85,484.22)	0.00 (0.00, 0.00)	0.00 (0.00, 0.00)	0.00 (0.00, 0.00)	0.00 (0.00, 0.00)	0.2	0.3
Natalizumab-SC	307,119.99 (254,825.07, 349,972.93)	8.79 (6.42, 10.81)	-87,475.73 (-130,088.51, -16,054.09)	381.97 (-57,166.74, 88,724.47)	-44,863.65 (-117,556.93, 23,377.02)	-0.74 (-3.16, 1.56)	26,367.04 (-19,743.08, 65,111.96)	18,968.39 (-49,966.19, 74,651.46)	0.3	0.2
Natalizumab biosimilar-IV	335,639.43 (300,172.95, 366,383.61)	9.37 (5.17, 13.02)	-101,346.79 (-203,893.87, -24,260.76)	-7,629.74 (-152,180.34, 104,096.81)	-16,344.21 (-63,438.11, 30,093.22)	-0.15 (-3.37, 3.64)	12,495.98 (-38,068.58, 85,145.35)	10,956.69 (-68,749.69, 121,517.63)	0.2	0.2
Ublituximab	358,531.30 (320,909.58, 382,255.80)	9.52 (6.39, 12.86)	-120,565.32 (-200,000.73, -56,702.17)	-25,378.93 (-136,137.79, 68,986.73)	6,547.66 (-39,880.26, 56,276.70)	-0.007 (-2.99, 3.24)	-6,722.55 (-77,678.40, 94,597.22)	-6,792.51 (-104,922.20, 124,893.82)	0.3	0.3
Ofatumumab	376,459.79 (333,842.01, 424,864.00)	8.88 (7.22, 10.29)	-154,457.35 (-231,609.74, -80,398.90)	-65,656.38 (-156,666.10, 21,522.09)	24,476.14 (-37,033.99, 65,102.60)	-0.65 (-1.87, 0.46)	-40,614.58 (-88,075.19, 21,261.62)	-47,069.95 (-104,895.79, 16,015.75)	0	0
Ocrelizumab	398,560.74 (342,704.38, 444,491.77)	9.44 (6.78, 12.89)	-162,498.18 (-219,503.94, -102,139.76)	-68,073.15 (-147,570.83, 21,565.49)	46,577.10 (-2,960.08, 108,924.42)	-0.083 (-1.78, 3.16)	-48,655.40 (-116,194.26, 9,154.72)	-49,486.73 (-123,248.85, 40,781.92)	0	0

#write.csv(results_table, paste0(c("results/results warden determenistic", n_samples, "x", n_patients, ".csv"), collapse = " "))

#plot results
v_results_treatment_names<- v_treatment_names[-c(4,5)]

# Create a BCEA object
des_bcea <- bcea(eff = total_qalys,
                 cost = total_costs,
                 ref = 1,
                 interventions = v_results_treatment_names)
ceplane.plot(des_bcea, wtp = 35000,
             graph = "ggplot2")


#plot a CEAC
des_multi_ce <- multi.ce(des_bcea)
ceac.plot(des_multi_ce, graph = "ggplot2",
          line = list(color = c("darkorange", "palegreen3", "blue", "grey",
                                "cyan3", "black"),
                      size=2 ),
          #ylim=c(0, 0.5),
          pos = c(1, 1)
)




#compare key events by treatment groups
modelresults_summarised <- summary_results_sim(results[[1]])
proportion_spms  <- modelresults_summarised["spms_flag",]
average_numberof_relapses <- modelresults_summarised["n_relapse",] 
average_numberof_edssincreases <- modelresults_summarised["n_edssinc_evts",] 
average_numberof_edssdecreases<- modelresults_summarised["n_edssdec_evts",] 
average_numberof_saes<- modelresults_summarised["n_sae_evts",] 
final_edss_value <- modelresults_summarised["edss_value",]
average_numberof_treatmentswitches <- modelresults_summarised["n_switch",] 
average_numberof_yearsontreatment <- modelresults_summarised["years_on_treatment",] 
average_numberof_yearson2ndline<- modelresults_summarised["t_2ndline",] 
average_numberof_yearson3rdline<- modelresults_summarised["t_3rdline",] 
average_numberof_yearson4thline<- modelresults_summarised["t_4thline",]

key_events_trt_grps <- as.data.frame(rbind(
  proportion_spms,   average_numberof_relapses, 
  average_numberof_edssincreases,   average_numberof_edssdecreases, 
  final_edss_value,   average_numberof_saes,
  average_numberof_treatmentswitches,  average_numberof_yearsontreatment,
  average_numberof_yearson2ndline,  average_numberof_yearson3rdline,
  average_numberof_yearson4thline),
  row.names=c("proportion with spms", "average no. of relapses", 
              "average no. of edss increases", "average no. of edss decreases", 
              "final edss value", "average no. of SAEs",
              "average no. of treatment switches", "average no. of years on treatment",
              "average no. of years on 2ndline", "average no. of years on 3rdline",
              "average no. of years on 4thline")
)
knitr::kable(key_events_trt_grps)

	Natalizumab-IV	Natalizumab-SC	Natalizumab biosimilar-IV	Natalizumab biosimilar-SC	Cladribine	Ublituximab	Ofatumumab	Ocrelizumab
proportion with spms	0.29 (0.1; 0.5)	0.32 (0.1; 0.6)	0.35 (0.1; 0.8)	0.26 (0; 0.6)	0.35 (0.2; 0.7)	0.34 (0.1; 0.5)	0.35 (0.2; 0.5)	0.27 (0; 0.5)
average no. of relapses	2.48 (1; 4.3)	2.13 (1.2; 3.7)	2.51 (0.8; 4.1)	2.23 (0.9; 3.2)	2.51 (1.3; 3.7)	2.25 (0.7; 3.4)	2.5 (1; 4.3)	2.57 (1.3; 5.2)
average no. of edss increases	2.66 (1.5; 4)	2.31 (1; 4.1)	2.27 (1.1; 3.8)	2.58 (1.2; 4.6)	3 (1.2; 6.3)	2.71 (1.1; 3.8)	2.39 (0.9; 4.2)	2.34 (1.4; 3.8)
average no. of edss decreases	1.74 (1.4; 2.3)	1.59 (1.2; 2)	1.5 (1; 1.9)	1.66 (1.4; 2)	1.63 (1; 2.8)	1.6 (1.1; 2.2)	1.55 (1; 2.2)	1.65 (1; 2.3)
final edss value	4.65 (3; 5.6)	4.61 (3.2; 6.2)	4.7 (3.7; 6.4)	4.69 (2.7; 6)	5.18 (4; 6.7)	4.77 (2.9; 5.9)	4.76 (3.9; 5.9)	4.6 (3.1; 6)
average no. of SAEs	4.09 (2.3; 7.2)	3.58 (2.8; 7)	4.65 (2.6; 6.6)	4.1 (3.2; 5.9)	4.12 (2.7; 6.8)	4.02 (2.5; 7.4)	7.42 (3.4; 20.8)	4.39 (2.4; 6.4)
average no. of treatment switches	1.11 (0.4; 1.9)	1.19 (0.6; 1.8)	1.12 (0.6; 1.9)	1.17 (0.6; 1.8)	1.14 (0.6; 1.7)	1.12 (0.4; 1.7)	1.25 (0.5; 2.2)	1.11 (0.6; 1.9)
average no. of years on treatment	20.75 (15.1; 30.7)	18.36 (8.7; 24.5)	20.85 (11.9; 35.7)	19.54 (12.6; 31.1)	15.02 (3.3; 23.8)	18.33 (9.9; 25.1)	19.33 (10.3; 23.1)	18.19 (11.4; 24.7)
average no. of years on 2ndline	19.43 (11.3; 32.4)	16.49 (10; 24.7)	17.65 (6.4; 32.8)	18.79 (13.1; 28.5)	14.58 (7.7; 19.7)	17.37 (7.5; 26.9)	15.58 (10.1; 21.2)	18.06 (9; 28.1)
average no. of years on 3rdline	6.97 (1.4; 12.2)	5.53 (2.2; 9.9)	5.85 (2; 14.4)	6.26 (0.6; 14.3)	5.36 (1.6; 13.9)	6.5 (2.8; 14.8)	6.17 (2.5; 15.8)	6.75 (2.9; 11.4)
average no. of years on 4thline	4.78 (0.9; 12.2)	5.27 (0; 15.9)	4 (0; 9.8)	4.83 (1.9; 8)	3.6 (0.1; 7.4)	4.26 (0; 7)	5.21 (0; 12.2)	4.34 (0.2; 7)



#  VOI analysis
#total EVPI

total_evpi <- voi::evpi(outputs = list(c = total_costs,
                                       e = total_qalys,
                                       k = c(0, 25000, 50000, 75000, 100000)))

knitr::kable(total_evpi)

k	evpi
0	6620.079
25000	16847.855
50000	74268.245
75000	118801.232
100000	163334.219

plot(x=total_evpi$k, y=total_evpi$evpi, type="l")


# Total EVPPI
# All  parameters excluding hazard ratios that are constant as it causes errors
# df_inputs <- as.data.frame(results[["input_matrix"]][,-c(73:76,81:84,89:92,79:101)])
# parameters <- as.list(colnames(df_inputs))
# 
# total_evppi <- voi::evppi(outputs = list(c = total_costs,
#                                          e = total_qalys,
#                                          k = c(25000, 35000)),
#                           inputs = df_inputs,
#                           pars = parameters,
#                           method = "bart", 
#                           verbose = TRUE)
# 
# knitr::kable(total_evppi)
# 
# #plot in base R split ober A and B panels as there are too many parameters to plot
# total_evppi_25k_a <- total_evppi$evppi[total_evppi$k ==25000 & total_evppi$evppi > mean(total_evppi$evppi)]
# names(total_evppi_25k_a) <- total_evppi$pars[total_evppi$k ==25000 & total_evppi$evppi > mean(total_evppi$evppi)]
# par(las=1, bg="white", mai=c(0.85, 4, 0.5, 0.5))
# barplot(
#   sort(total_evppi_25k_a, decreasing=FALSE), horiz = TRUE, cex.names=0.9, main="EVPPI WTP 25,000/QALY (A)", xlab= "£"
# )
# 
# total_evppi_25k_b <- total_evppi$evppi[total_evppi$k ==25000 & total_evppi$evppi < mean(total_evppi$evppi)]
# names(total_evppi_25k_b) <- total_evppi$pars[total_evppi$k ==25000 & total_evppi$evppi < mean(total_evppi$evppi)]
# barplot(
#   sort(total_evppi_25k_b, decreasing=FALSE), horiz = TRUE, cex.names=0.6, main="EVPPI WTP 25,000/QALY (B)", xlab= "£"
# )
# 
# total_evppi_35k_a <- total_evppi$evppi[total_evppi$k ==35000 & total_evppi$evppi > mean(total_evppi$evppi)]
# names(total_evppi_35k_a) <- total_evppi$pars[total_evppi$k ==35000 & total_evppi$evppi > mean(total_evppi$evppi)]
# par(las=1, bg="white", mai=c(0.85, 4, 1, 0.5))
# barplot(
#   sort(total_evppi_35k_a, decreasing=FALSE), horiz = TRUE, cex.names=0.5, main="EVPPI WTP 35,000/QALY (A)", xlab= "£"
# )
# 
# total_evppi_35k_b <- total_evppi$evppi[total_evppi$k ==35000 & total_evppi$evppi < mean(total_evppi$evppi)]
# names(total_evppi_35k_b) <- total_evppi$pars[total_evppi$k ==35000 & total_evppi$evppi < mean(total_evppi$evppi)]
# par(las=1, bg="white", mai=c(0.85, 4, 1, 0.5))
# barplot(
#   sort(total_evppi_35k_b, decreasing=FALSE), horiz = TRUE, cex.names=0.6, main="EVPPI WTP 35,000/QALY (B)",xlab= "£"
# )

Ayman Sadek

October 24, 2026

Introduction

Main options

model setup

Run Model

Plot Results