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 eular
## 1:     1   1  1     6  3    1        1     0 49.43622 12.982301   0     1
## 2:     1   1  1    12  3    1        2     0 49.93622 11.982301   0     1
## 3:     1   1  1    18  3    1        3     0 50.43622 10.982301   0     1
## 4:     1   1  1    24  3    1        4     0 50.93622  9.982301   0     1
## 5:     1   1  1    30  3    1        5     0 51.43622  8.982301   0     1
## 6:     1   1  1    36  3    1        6     0 51.93622  7.982301   0     1
##    das28 sdai cdai      haq     ttsi si yrlen  tx_cost hosp_days hosp_cost
## 1:    NA   NA   NA 1.431665 139.9466  0   0.5 20903.76 0.1401969  172.8017
## 2:    NA   NA   NA 1.427322 138.9466  0   0.5 20903.76 0.1401969  172.8017
## 3:    NA   NA   NA 1.422980 137.9466  0   0.5 20903.76 0.1401969  172.8017
## 4:    NA   NA   NA 1.418638 136.9466  0   0.5 20903.76 0.1401969  172.8017
## 5:    NA   NA   NA 1.414295 135.9466  0   0.5 20903.76 0.1401969  172.8017
## 6:    NA   NA   NA 1.409953 134.9466  0   0.5 20903.76 0.1401969  172.8017
##    mgmt_cost si_cost prod_loss   utility      qalys      txseq
## 1:  180.0968       0  4719.724 0.5971431 0.29857156 Sequence 1
## 2:  180.0968       0  4705.408 0.6357343 0.31786715 Sequence 1
## 3:  180.0968       0  4691.092 0.6341664 0.31708321 Sequence 1
## 4:  180.0968       0  4676.777 0.6894384 0.34471921 Sequence 1
## 5:  180.0968       0  4662.461 0.5915446 0.29577228 Sequence 1
## 6:  180.0968       0  4648.145 0.1065019 0.05325096 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 20.895 14.09810       1.0400       19.8550        0
## 2:     1   2 20.885 14.17093       1.1525       19.7325        0
## 3:     1   3 20.355 13.91551       0.9250       19.4300        0
## 4:     1   4 19.150 13.48714       0.6825       18.4675        0
## 5:     1   5 21.085 14.30082       1.2275       19.8575        0
## 6:     1   6 21.350 14.65675       1.6775       19.6725        0
##          dhaq   si    qalys   dqalys  tx_cost dtx_cost hosp_days hosp_cost
## 1: 0.10587768 0.26 11.35006 7.831600 437681.3 333365.8 10.738665  12461.68
## 2: 0.29679059 0.22 10.54195 7.521116 379215.8 303640.8 13.605340  17981.94
## 3: 0.27222085 0.23 10.83045 7.687801 369215.8 302111.8 10.846623  14160.47
## 4: 0.07801418 0.17 10.82399 7.723904 434832.4 334769.1  8.837646  11778.25
## 5: 0.19359357 0.15 11.40913 7.984962 466395.1 356547.5 10.998807  13810.31
## 6: 0.14586573 0.28 11.72070 8.164887 495979.5 383364.2 19.867598  28330.36
##    dhosp_cost mgmt_cost dmgmt_cost   si_cost  dsi_cost prod_loss
## 1:   7650.772  7526.247   5078.045 1652.3335 1206.1992  166015.5
## 2:  10601.706  7769.914   5272.057 1427.1379 1083.7476  200557.6
## 3:   8375.257  7838.644   5358.817 1418.4949 1089.8306  210545.0
## 4:   7888.793  7422.638   5227.683 1065.7569  817.1346  157897.1
## 5:   8124.581  7921.401   5372.659  877.7412  649.8134  199585.9
## 6:  18561.784  8469.510   5814.308 1326.7902  996.6037  100278.5
##    dprod_loss dhc_cost dtot_cost yrs_since_approval dqalys_ann
## 1:  107718.86 347300.8  455019.7           21.85011  0.3105838
## 2:  126995.33 320598.3  447593.6           22.29601  0.3012485
## 3:  136062.00 316935.8  452997.8           22.49813  0.3081736
## 4:  108624.55 348702.7  457327.2           21.60991  0.3079992
## 5:  129374.66 370694.6  500069.3           21.66821  0.3181279
## 6:   67101.71 408736.9  475838.6           21.88356  0.3268039
##    dhc_cost_ann dprod_loss_ann
## 1:     14051.53       4251.623
## 2:     13215.93       5024.046
## 3:     13045.73       5419.896
## 4:     14132.04       4311.842
## 5:     15000.58       5123.857
## 6:     16609.79       2699.484
head(sim.out1.summary$time.means)
##    model sim month alive     qalys      haq  tx_cost hosp_cost mgmt_cost
## 1:     1   1     6   100 0.2829719 1.146948 20903.76  253.4783  180.0968
## 2:     1   1    12    97 0.2862925 1.058111 21002.36  217.8435  180.0968
## 3:     1   1    18    95 0.2815945 1.068854 20202.25  227.3757  180.0968
## 4:     1   1    24    95 0.2982769 1.074813 18782.95  229.5153  180.0968
## 5:     1   1    30    95 0.2745592 1.121349 18355.88  228.3262  180.0968
## 6:     1   1    36    94 0.2920960 1.059242 17112.73  221.2133  180.0968
##      si_cost prod_loss
## 1:  63.55129  3781.105
## 2:  65.51679  3488.239
## 3:  66.89609  3523.658
## 4: 133.79219  3543.302
## 5:   0.00000  3696.714
## 6:   0.00000  3491.971
head(sim.out1.summary$out0)
##    model sim id tx acr eular        ttd         ttsi
## 1:     1   1  1  3   1     2  3.0977332 5.130062e+01
## 2:     1   1  1  5   2     2 47.4728565 1.898929e+01
## 3:     1   1  1  8   1     2  0.5384596 2.110232e+02
## 4:     1   1  1 14   0     0  0.0000000 5.457277e+09
## 5:     1   1  2  3   3     2  3.1130602 2.722444e+01
## 6:     1   1  2  5   0     0  0.0000000 3.018634e+01