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.

Complete 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
## 1:     1   1  1     6  3    1        1     0 53.88004  0.00000000   0
## 2:     1   1  1    12  5    2        1     0 54.38004  3.96800684   0
## 3:     1   1  1    18  5    2        2     0 54.88004  2.96800684   0
## 4:     1   1  1    24  5    2        3     0 55.38004  1.96800684   0
## 5:     1   1  1    30  5    2        4     0 55.88004  0.96800684   0
## 6:     1   1  1    36  5    2        5     0 56.38004 -0.03199316   0
##    eular das28 sdai cdai      haq     ttsi si yrlen  tx_cost  hosp_days
## 1:     0    NA   NA   NA 1.748439 36.06485  0   0.5 22115.24 0.66999367
## 2:     1    NA   NA   NA 1.412559 35.39017  0   0.5 22278.58 0.07371646
## 3:     1    NA   NA   NA 1.406302 34.39017  0   0.5 22278.58 0.07371646
## 4:     1    NA   NA   NA 1.400044 33.39017  0   0.5 22278.58 0.07371646
## 5:     1    NA   NA   NA 1.393786 32.39017  0   0.5 22278.58 0.07371646
## 6:     1    NA   NA   NA 1.668734 31.39017  0   0.5 22278.58 0.66999367
##    hosp_cost mgmt_cost si_cost prod_loss    utility      qalys      txseq
## 1: 877.44936  206.5151       0  4641.221 0.53897897 0.26948948 Sequence 1
## 2:  92.04092  206.5151       0  3749.631 0.52073151 0.26036576 Sequence 1
## 3:  92.04092  206.5151       0  3733.020 0.07084795 0.03542397 Sequence 1
## 4:  92.04092  206.5151       0  3716.409 0.66197195 0.33098597 Sequence 1
## 5:  92.04092  206.5151       0  3699.798 0.05287038 0.02643519 Sequence 1
## 6: 877.44936  206.5151       0  4429.647 0.58830246 0.29415123 Sequence 1

The 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.

Summary output

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.190 13.33928       1.3750       17.8150        0
## 2:     1   2 20.450 14.06886       1.2475       19.2025        0
## 3:     1   3 18.405 12.97018       1.2000       17.2050        0
## 4:     1   4 22.080 14.83969       1.3025       20.7775        0
## 5:     1   5 19.420 13.30105       1.1250       18.2950        0
## 6:     1   6 18.870 12.89980       0.9625       17.9075        0
##           dhaq   si     qalys   dqalys  tx_cost dtx_cost hosp_days
## 1: -0.01256774 0.18 10.053779 7.141130 426760.2 328453.3  14.94116
## 2:  0.15197556 0.23 10.198160 7.170670 363637.6 287592.3  16.23963
## 3:  0.07048098 0.14 10.081042 7.209824 369094.2 294109.8  12.40449
## 4:  0.07373880 0.17 11.806656 8.088598 434507.4 341601.5  21.67478
## 5:  0.18486994 0.17  9.374515 6.626901 399391.4 309583.6  21.17218
## 6:  0.19772908 0.14  9.736306 6.765677 394434.6 295464.6  10.29068
##    hosp_cost dhosp_cost mgmt_cost dmgmt_cost   si_cost dsi_cost prod_loss
## 1:  19674.73  12934.243  7926.049   5509.525  999.0275 641.9417  130789.6
## 2:  22286.10  14203.282  7866.493   5411.863 1158.2648 923.8942  163303.1
## 3:  15594.11  10526.067  6967.812   4910.285  955.3593 723.7438  111502.7
## 4:  25497.53  15910.507  8999.035   6048.138 1020.7943 863.5781  157884.1
## 5:  26560.56  17080.036  7658.567   5245.469  948.0979 704.0521  268486.4
## 6:  15224.50   9377.341  7883.200   5389.069  826.6805 628.3806  222620.7
##    dprod_loss dhc_cost dtot_cost yrs_since_approval dqalys_ann
## 1:   87873.24 347539.0  435412.3           21.73401  0.2846090
## 2:  108665.27 308131.3  416796.6           22.71512  0.2839645
## 3:   76804.50 310269.8  387074.4           22.00427  0.2890435
## 4:  102496.77 364423.8  466920.5           22.64487  0.3233713
## 5:  176261.58 332613.2  508874.7           21.61206  0.2654827
## 6:  148878.64 310859.3  459738.0           22.25244  0.2676979
##    dhc_cost_ann dprod_loss_ann
## 1:     13994.14       3457.607
## 2:     12377.06       4271.318
## 3:     12530.70       3076.920
## 4:     14814.88       4024.812
## 5:     13435.59       6960.148
## 6:     12340.11       5862.680
head(sim.out1.summary$time.means)
##    model sim month alive     qalys      haq  tx_cost hosp_cost mgmt_cost
## 1:     1   1     6   100 0.2881622 1.244936 22115.24  467.1567  206.5151
## 2:     1   1    12    99 0.2799452 1.160639 22177.94  402.1655  206.5151
## 3:     1   1    18    99 0.2821838 1.066589 20278.46  394.1646  206.5151
## 4:     1   1    24    98 0.2938425 1.044726 19628.86  382.7827  206.5151
## 5:     1   1    30    97 0.2912116 1.079306 18805.75  359.6472  206.5151
## 6:     1   1    36    96 0.2927717 1.076499 17305.43  382.4209  206.5151
##     si_cost prod_loss
## 1: 55.50153  3304.677
## 2: 56.06215  3080.909
## 3:  0.00000  2831.254
## 4:  0.00000  2773.220
## 5:  0.00000  2865.012
## 6:  0.00000  2857.560
head(sim.out1.summary$out0)
##    model sim id tx acr eular       ttd      ttsi
## 1:     1   1  1  3   1     2  2.668269  35.64467
## 2:     1   1  1  5   1     2  1.785475  93.92105
## 3:     1   1  1  8   2     1  8.051842  30.13126
## 4:     1   1  2  3   1     1  1.676008  10.30048
## 5:     1   1  2  5   2     2 15.705864 110.92480
## 6:     1   1  2  8   0     0  0.000000  46.17941