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 49.34167  0.00000000   0
## 2:     1   1  1    12  5    2        1     0 49.84167  1.01616606   1
## 3:     1   1  1    18  5    2        2     0 50.34167  0.01616606   1
## 4:     1   1  1    24  5    2        3     0 50.84167 -0.98383394   1
## 5:     1   1  1    30  8    3        1     0 51.34167  0.00000000   0
## 6:     1   1  1    36 14    4        1     0 51.84167          NA   0
##    eular das28 sdai cdai      haq      ttsi si yrlen    tx_cost hosp_days
## 1:     0    NA   NA   NA 2.061876  25.58617  0   0.5 22554.8552 0.7550633
## 2:     2    NA   NA   NA 1.372245  37.78928  0   0.5 22339.1997 0.6093987
## 3:     2    NA   NA   NA 1.366207  36.78928  0   0.5 22339.1997 0.6093987
## 4:     2    NA   NA   NA 1.947925  35.78928  0   0.5 22339.1997 0.2276890
## 5:     0    NA   NA   NA 1.947925 101.26948  0   0.5 14550.2739 0.2276890
## 6:     0    NA   NA   NA 1.947925        NA  0   0.5   615.8625 0.2276890
##    hosp_cost mgmt_cost si_cost prod_loss    utility      qalys      txseq
## 1: 1018.1458   179.451       0  5896.261 0.20580269 0.10290134 Sequence 1
## 2:  836.5431   179.451       0  3924.151 0.59441250 0.29720625 Sequence 1
## 3:  836.5431   179.451       0  3906.885 0.66479710 0.33239855 Sequence 1
## 4:  302.4408   179.451       0  5570.399 0.55065326 0.27532663 Sequence 1
## 5:  302.4408   179.451       0  5570.399 0.03168531 0.01584266 Sequence 1
## 6:  302.4408   179.451       0  5570.399 0.69352011 0.34676006 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 dhaq ## 1: 1 1 21.630 14.80567 1.0175 20.6125 0 0.2277345 ## 2: 1 2 23.265 15.43161 1.0675 22.1975 0 0.1895127 ## 3: 1 3 17.400 12.19669 0.8650 16.5350 0 0.1369388 ## 4: 1 4 20.610 14.12631 1.0625 19.5475 0 0.2013628 ## 5: 1 5 20.320 13.90575 0.8050 19.5150 0 0.2217152 ## 6: 1 6 21.235 14.41850 1.0825 20.1525 0 0.2513651 ## si qalys dqalys tx_cost dtx_cost hosp_days hosp_cost ## 1: 0.24 10.781622 7.595900 404620.3 319132.0 21.251140 26756.43 ## 2: 0.26 12.268406 8.253489 503911.6 378024.5 20.856204 26019.48 ## 3: 0.13 9.480336 6.772433 315813.9 263269.5 11.214770 16389.78 ## 4: 0.11 10.599809 7.426626 427604.8 323790.1 15.030416 20373.70 ## 5: 0.18 9.791653 6.946507 350917.6 281202.6 14.792480 19951.80 ## 6: 0.15 10.967505 7.600056 381010.6 298179.4 9.644956 12744.59 ## dhosp_cost mgmt_cost dmgmt_cost si_cost dsi_cost prod_loss dprod_loss ## 1: 17279.423 7763.052 5313.785 1176.1537 913.8319 172889.1 114264.52 ## 2: 16688.794 8837.976 5862.203 1453.5467 983.4124 212493.7 137498.08 ## 3: 10738.498 7123.387 4993.203 635.7526 504.5201 182292.2 123132.51 ## 4: 13082.205 8663.688 5938.181 517.9750 385.7249 136780.7 90882.26 ## 5: 11980.865 7173.839 4909.330 982.4709 750.3847 187826.4 122587.31 ## 6: 7705.313 8794.608 5971.514 964.6748 737.1856 218022.0 143544.78 ## dhc_cost dtot_cost yrs_since_approval dqalys_ann dhc_cost_ann ## 1: 342639.0 456903.6 22.72163 0.2960171 13672.26 ## 2: 401558.9 539057.0 22.51578 0.3209304 15814.13 ## 3: 279505.7 402638.2 22.28904 0.2650214 11158.08 ## 4: 343196.2 434078.5 22.38295 0.2906845 13617.82 ## 5: 298843.2 421430.5 22.49300 0.2691092 11782.40 ## 6: 312593.4 456138.2 22.70989 0.2950108 12305.85 ## dprod_loss_ann ## 1: 4463.188 ## 2: 5377.726 ## 3: 4809.389 ## 4: 3573.674 ## 5: 4774.171 ## 6: 5567.512 head(sim.out1.summary$time.means)
##    model sim month alive     qalys      haq  tx_cost hosp_cost mgmt_cost
## 1:     1   1     6   100 0.2645210 1.309635 22554.86  571.5816   179.451
## 2:     1   1    12    98 0.2580874 1.212526 22447.03  474.8873   179.451
## 3:     1   1    18    98 0.2707105 1.161609 20685.29  448.5548   179.451
## 4:     1   1    24    97 0.2976020 1.148461 20127.85  432.2242   179.451
## 5:     1   1    30    96 0.2763700 1.161930 19131.92  445.0664   179.451
## 6:     1   1    36    96 0.2852212 1.158026 17965.50  454.0132   179.451
##      si_cost prod_loss
## 1:  98.01281  3745.108
## 2:  50.00653  3467.411
## 3:  50.00653  3321.804
## 4:   0.00000  3284.206
## 5: 102.09667  3322.722
## 6:   0.00000  3311.559
head(sim.out1.summary\$out0)
##    model sim id tx acr eular       ttd      ttsi
## 1:     1   1  1  3   1     1 17.019789  27.61814
## 2:     1   1  1  5   0     0  0.000000  14.01862
## 3:     1   1  1  8   2     2 47.005763  18.30061
## 4:     1   1  2  3   0     0  0.000000 129.20949
## 5:     1   1  2  5   1     1  2.850435 105.22817
## 6:     1   1  2  8   0     0  0.000000  10.96291