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