We run the IPS using the function sim_iviRA. The main arguments are the selected treatment arms, the input data, the parameter samples, and the model structure. The simulation can either return all simulated output (i.e., a dataset of simulated outcomes for each model structure, sampled parameter set, individual, and time period) or summaries of the simulated output.
Here, we return complete simulation output for Arm 1 and Arm 2.
sim.out1 <- sim_iviRA(tx_seqs = txseq1, input_data = input.dat, pars = parsamp,
model_structures = mod.structs, output = "data")
sim.out1[, txseq := "Sequence 1"]
sim.out2 <- sim_iviRA(tx_seqs = txseq2, input_data = input.dat, pars = parsamp,
model_structures = mod.structs, output = "data")
sim.out2[, txseq := "Sequence 2"]
sim.out <- rbind(sim.out1, sim.out2)
head(sim.out)## model sim id month tx line tx_cycle death age ttd acr eular das28
## 1: 1 1 1 6 4 1 1 0 61.02684 0 0 0 NA
## 2: 1 1 1 12 13 2 1 0 61.52684 0 0 0 NA
## 3: 1 1 1 18 18 3 1 0 62.02684 0 0 0 NA
## 4: 1 1 1 24 21 4 1 0 62.52684 NA 0 0 NA
## 5: 1 1 1 30 21 4 2 0 63.02684 NA 0 0 NA
## 6: 1 1 1 36 21 4 3 0 63.52684 NA 0 0 NA
## sdai cdai haq ttsi si yrlen tx_cost hosp_days hosp_cost
## 1: NA NA 1.137859 53.63632 0 0.5 23983.03 0.1734779 266.7863
## 2: NA NA 1.137859 188.33893 0 0.5 26456.54 0.1734779 266.7863
## 3: NA NA 1.137859 58.40078 0 0.5 13197.05 0.1734779 266.7863
## 4: NA NA 1.137859 NA 0 0.5 0.00 0.1734779 266.7863
## 5: NA NA 1.150307 NA 0 0.5 0.00 0.1734779 266.7863
## 6: NA NA 1.162755 NA 0 0.5 0.00 0.1734779 266.7863
## mgmt_cost si_cost prod_loss utility qalys txseq
## 1: 213.4018 0 4903.448 0.65786292 0.32893146 Sequence 1
## 2: 213.4018 0 4903.448 0.15101180 0.07550590 Sequence 1
## 3: 213.4018 0 4903.448 0.06146943 0.03073472 Sequence 1
## 4: 213.4018 0 4903.448 0.73159720 0.36579860 Sequence 1
## 5: 213.4018 0 4957.090 0.78599755 0.39299877 Sequence 1
## 6: 213.4018 0 5010.733 0.69145196 0.34572598 Sequence 1The output is data.table which is an enhanced data.frame from the data.table package designed for fast data manipulation. For a given model structure (model) and sampled parameter set (sim), a sampled patient (id) remains on a given treatment until time to treatment discontinuation (ttd) becomes less than zero in a given month. Treatment discontinuation is caused by a serious infection (si = 1) if the sampled time to serious infection at treatment initiation is less than the sampled time to discontinuation. After discontinuation, a patient switches to the next line of treatment and line increments by one. Patient age (age) increases in 6-month increments.
The HAQ score (haq) during the initial 6-month period (therapy_cycle = 1), disease activity measures (das28, sdai, and cdai), and ttd depend on the model structure chosen. For example, in the first model structure, HAQ and time to treatment discontinuation depend on ACR response (acr) and in turn, EULAR response (eular). The simulation returned NA values for the three disease activity measures because treatment switching is is not a function of disease activity (i.e., S6 was chosen in the first model structure). The simulation ends when a patient dies (death = 1).
Health care sector costs consist of treatment costs (tx_cost); hospital costs (hosp_cost), which increase with the HAQ score; general management costs (mgmt_cost); and the costs of caused by serious infections (si_cost). Non-health care sector costs are those due to lost wages. Analyses performed from a societal perspective should include these productivity losses (prod_loss).
Utility (utility) depends on whether the Hernandez-Alava mixture model or Wailoo logistic regression equation is used to simulate utility. In both cases, QALYs (qalys) are calculated as a function of utility and the estimated utility loss from a serious infection.
One drawback of returning all simulation output is that it can place significant strains on computer memory. For example, if 200 model structures, 1,000 parameters sets, and 5,000 patients are simulated for 15 years (i.e., 30 model cycles), then the number of observations is 30 billion. The sim_iviRA function consequently has the option to only return summaries of simulation output.
In particular, by invoking the output = summary option, the data.table objects means, time.means, and out0 are returned. means reports simulated outcomes averaged over individuals and time and tmeans reports simulated outcomes averaged across individuals but conditional on time. out0 provides outcomes for the first model cycle including time to treatment discontinuation and time to a serious infection. Output is always conditional on the model structure and sampled parameter set.
sim.out1.summary <- sim_iviRA(tx_seqs = txseq1, input_data = input.dat, pars = parsamp,
model_structures = mod.structs, output = "summary")
head(sim.out1.summary$means)## model sim lys dlys lys_infusion lys_injection lys_oral
## 1: 1 1 19.400 13.28677 1.1350 18.2650 0
## 2: 1 2 21.310 14.38221 1.4650 19.8450 0
## 3: 1 3 19.365 13.50511 0.9825 18.3825 0
## 4: 1 4 20.525 14.06020 1.0200 19.5050 0
## 5: 1 5 17.960 12.61827 1.1025 16.8575 0
## 6: 1 6 19.560 13.58596 1.1250 18.4350 0
## dhaq si qalys dqalys tx_cost dtx_cost hosp_days
## 1: 0.16426930 0.22 10.047157 7.119037 429359.2 339349.7 17.08968
## 2: 0.20901099 0.18 10.655118 7.512766 415889.3 336534.0 54.09047
## 3: 0.13076876 0.21 10.753591 7.543713 372586.8 296686.6 16.35749
## 4: 0.16370347 0.16 10.620958 7.434538 457830.9 357063.8 12.20283
## 5: 0.14092185 0.20 9.500008 6.868745 363178.9 292775.9 12.78655
## 6: 0.05569675 0.24 10.302875 7.280763 485507.9 378836.2 14.37676
## hosp_cost dhosp_cost mgmt_cost dmgmt_cost si_cost dsi_cost prod_loss
## 1: 24382.44 15060.92 8279.990 5670.843 1254.421 895.2503 225172.44
## 2: 67490.46 45418.48 8997.288 6072.309 1195.176 898.3993 188962.97
## 3: 26057.67 17608.61 7916.669 5521.068 1162.422 966.3870 171583.54
## 4: 16179.72 10194.25 8843.203 6057.841 1086.851 779.7663 210316.54
## 5: 16511.60 10533.81 6993.642 4913.567 1162.305 929.9336 134732.97
## 6: 20256.78 12940.87 7773.689 5399.438 1139.475 848.0836 82343.05
## dprod_loss dhc_cost dtot_cost yrs_since_approval dqalys_ann
## 1: 147779.78 360976.7 508756.5 24.90399 0.2898159
## 2: 120163.44 388923.2 509086.6 25.88033 0.3080026
## 3: 117649.76 320782.7 438432.5 26.04269 0.3067010
## 4: 139459.69 374095.7 513555.4 25.65233 0.3040379
## 5: 90849.88 309153.2 400003.1 25.35447 0.2789843
## 6: 55750.27 398024.6 453774.9 24.98272 0.2970480
## dhc_cost_ann dprod_loss_ann
## 1: 14907.66 5868.013
## 2: 15984.18 4822.521
## 3: 13310.02 4705.078
## 4: 15460.91 5644.147
## 5: 12842.43 3610.262
## 6: 16338.13 2229.207head(sim.out1.summary$time.means)## model sim month alive qalys haq tx_cost hosp_cost mgmt_cost
## 1: 1 1 6 100 0.2821861 1.239386 23983.03 479.5016 213.4018
## 2: 1 1 12 100 0.2647001 1.195155 24774.55 449.2836 213.4018
## 3: 1 1 18 97 0.2786685 1.162019 23388.18 474.1959 213.4018
## 4: 1 1 24 97 0.2958194 1.144112 22568.32 464.5637 213.4018
## 5: 1 1 30 95 0.2958881 1.157702 21910.98 489.4309 213.4018
## 6: 1 1 36 95 0.2836036 1.159949 20945.67 454.0778 213.4018
## si_cost prod_loss
## 1: 57.01912 5340.962
## 2: 114.03823 5150.356
## 3: 58.78259 5007.560
## 4: 58.78259 4930.395
## 5: 60.02012 4988.957
## 6: 60.02012 4998.641head(sim.out1.summary$out0)## model sim id tx acr eular ttd ttsi
## 1: 1 1 1 4 2 2 22.6477 3.473221e+01
## 2: 1 1 1 13 0 0 0.0000 1.048569e+02
## 3: 1 1 1 18 0 0 0.0000 2.025807e+01
## 4: 1 1 1 21 0 0 0.0000 2.501323e+10
## 5: 1 1 2 4 0 0 0.0000 1.892935e+01
## 6: 1 1 2 13 0 0 0.0000 2.165614e+02