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 58.26744 15.58959   2     2
## 2:     1   1  1    12  3    1        2     0 58.76744 14.58959   2     2
## 3:     1   1  1    18  3    1        3     0 59.26744 13.58959   2     2
## 4:     1   1  1    24  3    1        4     0 59.76744 12.58959   2     2
## 5:     1   1  1    30  3    1        5     0 60.26744 11.58959   2     2
## 6:     1   1  1    36  3    1        6     0 60.76744 10.58959   2     2
##    das28 sdai cdai      haq    ttsi si yrlen tx_cost hosp_days hosp_cost
## 1:    NA   NA   NA 1.570994 20.1043  0   0.5 23059.2 0.1288273   164.905
## 2:    NA   NA   NA 1.559978 19.1043  0   0.5 23059.2 0.1288273   164.905
## 3:    NA   NA   NA 1.548961 18.1043  0   0.5 23059.2 0.1288273   164.905
## 4:    NA   NA   NA 1.537945 17.1043  0   0.5 23059.2 0.1288273   164.905
## 5:    NA   NA   NA 1.526928 16.1043  0   0.5 23059.2 0.1288273   164.905
## 6:    NA   NA   NA 1.515912 15.1043  0   0.5 23059.2 0.1288273   164.905
##    mgmt_cost si_cost prod_loss   utility      qalys      txseq
## 1:  192.6214       0  5857.471 0.2200887 0.11004436 Sequence 1
## 2:  192.6214       0  5816.395 0.1478795 0.07393977 Sequence 1
## 3:  192.6214       0  5775.320 0.7014200 0.35071001 Sequence 1
## 4:  192.6214       0  5734.245 0.7219599 0.36097995 Sequence 1
## 5:  192.6214       0  5693.170 0.7828545 0.39142723 Sequence 1
## 6:  192.6214       0  5652.095 0.6811429 0.34057145 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 17.520 12.25783       0.9400       16.5800        0
## 2:     1   2 18.595 13.24302       1.3600       17.2350        0
## 3:     1   3 19.260 13.59348       1.1875       18.0725        0
## 4:     1   4 20.615 14.10622       1.4150       19.2000        0
## 5:     1   5 21.640 14.78227       1.0825       20.5575        0
## 6:     1   6 19.810 13.61937       1.6400       18.1700        0
##          dhaq   si     qalys   dqalys  tx_cost dtx_cost hosp_days
## 1: 0.10450258 0.28  9.585609 6.746590 371082.8 293629.9  10.90902
## 2: 0.10011371 0.15 10.287727 7.376698 436209.8 340951.9  16.74565
## 3: 0.28765800 0.17  9.252847 6.658584 351492.4 278175.1  18.60475
## 4: 0.25043886 0.11 10.640150 7.493608 382385.2 302934.8  15.17841
## 5: 0.28773956 0.34 11.135726 7.861864 420931.9 335965.8  15.41210
## 6: 0.01583112 0.21 11.081349 7.660244 403132.8 316866.1  13.07036
##    hosp_cost dhosp_cost mgmt_cost dmgmt_cost   si_cost  dsi_cost prod_loss
## 1:  16430.04   11256.38  6749.454   4722.242 1765.0196 1370.4585 160670.50
## 2:  19722.96   13682.33  7113.040   5065.778  788.0281  616.0245 118662.07
## 3:  25176.70   16997.81  7752.903   5471.908 1187.0986  941.0270 113284.02
## 4:  21664.12   13382.18  8138.437   5568.886  565.8199  358.6323  96160.05
## 5:  22515.46   13657.18  8167.234   5579.033 2392.2745 1839.1954 199655.50
## 6:  15709.33   10593.67  7768.706   5340.982 1438.3266 1139.9077 229599.33
##    dprod_loss dhc_cost dtot_cost yrs_since_approval dqalys_ann
## 1:  110655.18 310979.0  421634.2           21.84676  0.2638491
## 2:   83703.59 360316.1  444019.7           22.21022  0.2933768
## 3:   78210.97 301585.8  379796.8           22.78544  0.2653046
## 4:   63369.25 322244.5  385613.7           22.76552  0.2961682
## 5:  130297.98 357041.2  487339.2           22.45231  0.3125644
## 6:  154923.62 333940.6  488864.3           22.03369  0.3005091
##    dhc_cost_ann dprod_loss_ann
## 1:     12256.39       4370.069
## 2:     14513.44       3319.824
## 3:     12317.51       3094.495
## 4:     12874.07       2496.981
## 5:     14378.64       5158.309
## 6:     13379.72       6098.195
head(sim.out1.summary$time.means)
##    model sim month alive     qalys      haq  tx_cost hosp_cost mgmt_cost
## 1:     1   1     6   100 0.2682705 1.281949 23059.20  453.2312  192.6214
## 2:     1   1    12    97 0.2805215 1.160986 23165.62  426.8784  192.6214
## 3:     1   1    18    93 0.2914691 1.106799 21552.65  422.8916  192.6214
## 4:     1   1    24    92 0.2838542 1.153716 20268.94  444.1992  192.6214
## 5:     1   1    30    91 0.2845286 1.108191 19330.88  392.3578  192.6214
## 6:     1   1    36    88 0.2995234 1.085653 17874.57  386.0777  192.6214
##      si_cost prod_loss
## 1:  63.03642  4779.763
## 2:  64.98600  4328.750
## 3: 203.34328  4126.713
## 4: 342.58922  4301.642
## 5:   0.00000  4131.901
## 6:   0.00000  4047.870
head(sim.out1.summary$out0)
##    model sim id tx acr eular      ttd         ttsi
## 1:     1   1  1  3   0     0 0.000000 1.006781e+02
## 2:     1   1  1  5   0     0 0.000000 7.882900e+01
## 3:     1   1  1  8   0     0 0.000000 1.787979e+01
## 4:     1   1  1 14   0     0 0.000000 5.588959e+08
## 5:     1   1  2  3   0     2 5.997217 5.392255e+01
## 6:     1   1  2  5   0     1 8.107790 3.231932e+01