The iviRA package does not allow users to run the IPS with parameters set to fixed values, but instead recognizes that parameter values are inherently uncertain. As such, all parameters are randomly sampled from their (joint) probability distribution and the IPS is run for each randomly sampled parameter set. The parameters are sampled using the function sample_pars, which generates a probability distribution for all model parameters.

The random samples depend on the underlying statistical estimates of the distribution of the parameters. We can generate a sample of size 100 using default values.

parsamp <- sample_pars(n = 100, input_data = input.dat)

The object parsamp returned from sample_pars is a list of random draws of the parameters used in the IPS.

names(parsamp)
##  [1] "n"                 "acr"               "das28"            
##  [4] "haq"               "acr2haq"           "acr2das28"        
##  [7] "acr2sdai"          "acr2cdai"          "acr2eular"        
## [10] "eular2haq"         "rebound"           "haq.lprog.tx"     
## [13] "haq.lprog.age"     "haq.lcgm"          "lt"               
## [16] "mort.logor"        "mort.loghr.haqdif" "ttd.all"          
## [19] "ttd.da"            "ttd.eular"         "ttsi"             
## [22] "tx.cost"           "hosp.cost"         "mgmt.cost"        
## [25] "si.cost"           "utility.mixture"   "utility.wailoo"   
## [28] "si.ul"             "utility.tx.attr"   "prod.loss"

The number of samples is also returned for convenience.

parsamp$n
## [1] 100

Users can generate samples with iviRA defaults or feed custom inputs into sample_pars. For example, we can generate a sample in which the cost of a serious infection is $7,000 instead of its default value of $5,873.

parsamp2 <- sample_pars(n = 100, input_data = input.dat, si_cost = 7000)
summary(parsamp2$si.cost)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5702    6408    6985    7010    7526    8397

To view the full list of arguments in sample_pars that can changed by the user view the documentation by typing ?sample_pars. All default arguments are loaded with the package and are based on IVI’s statistical estimates and synthesis of the literature. The parameter estimates can be viewed by typing the command data(package = "iviRA") into the R console.

Below we describe these parameter estimates and separate them into 8 main categories: treatment effects estimated from a network meta-analysis (NMA); mappings between various measured of treatment response, disease activity, and/or function; long-term progression of disease, treatment duration, serious infections, mortality, utility, and costs.

Treatment effects during initial 6 month period

The effect of treatment on ACR response at 6 months, change in DAS28 at 6 months, and the change in the HAQ score at 6 months are estimated using a NMA. Default NMA parameter estimates are based on results from an NMA conduced by IVI using a Bayesian random effects approach. Each object is a list containing the mean and variance-covariance matrix of the posterior distribution of the parameters.

names(iviRA::nma.acr.naive)
## [1] "mean" "vcov"
iviRA::nma.acr.naive$mean
##          A         z1         z2         z3  d_cdmards d_abtivmtx 
##  0.6329128  0.0000000  0.6414534  1.2243365  0.0000000 -0.7924898 
##   d_adamtx      d_ada   d_etnmtx      d_etn   d_golmtx   d_ifxmtx 
## -0.7712937 -0.3653499 -0.9443291 -0.5837728 -0.8819216 -0.6760986 
##  d_placebo   d_tczmtx      d_tcz   d_czpmtx d_abtscmtx      d_nbt 
##  0.5629207 -0.8111553         NA -1.2556327 -0.8405713  0.0000000 
##   d_rtxmtx   d_tofmtx      d_rtx      d_tof      d_czp      d_gol 
## -0.7595393 -0.9111716 -0.5802354 -0.2177813 -0.6178098         NA
iviRA::nma.das28.naive$mean
##          A  d_cdmards d_abtivmtx   d_adamtx      d_ada   d_etnmtx 
## -1.1211094  0.0000000 -1.1052728         NA -0.2868123 -1.4966904 
##      d_etn   d_golmtx   d_ifxmtx  d_placebo   d_tczmtx      d_tcz 
## -0.7606214 -1.4478931 -0.7988820  0.6348693 -1.8619306 -0.8324263 
##   d_czpmtx d_abtscmtx      d_nbt   d_rtxmtx   d_tofmtx      d_rtx 
## -2.0054776 -1.1220057  0.0000000 -1.1918102         NA -0.8508762 
##      d_tof      d_czp      d_gol 
##         NA -0.2687301         NA
iviRA::nma.haq.naive$mean
##           A   d_cdmards  d_abtivmtx    d_adamtx       d_ada    d_etnmtx 
## -0.22058104  0.00000000 -0.22532161 -0.32124573 -0.07129560 -0.34065638 
##       d_etn    d_golmtx    d_ifxmtx   d_placebo    d_tczmtx       d_tcz 
## -0.24156205 -0.27247774 -0.21678447  0.20456867 -0.24777323  0.02597454 
##    d_czpmtx  d_abtscmtx       d_nbt    d_rtxmtx    d_tofmtx       d_rtx 
## -0.35480646 -0.21472087  0.00000000 -0.27057962          NA          NA 
##       d_tof       d_czp       d_gol 
## -0.28168946 -0.28382294          NA

In the PSA, we use a multivariate normal distribution to sample the NMA parameters. The advantage of this approach is that we do not have to load the entire posterior distribution with the package. The justification is the Bayesian central limit theorem, which states that under large samples, the posterior distribution will approach a multivariate normal distribution.

Treatment response mappings

There are various mapping between treatment response measures, disease activity, and functional status that can be used within the IVI-RA model. The exact mapping utilized for a given model run will depend on the model structure select.

The first set of mappings relates ACR response to the three measure of disease activity. The ACR response to SDAI mapping can either be based on the Leflunomide dataset or the inception cohort (see Aletaha and Smolen) while mapping for DAS28 and CDAI can only be based on the inception cohort. The mappings correspond to model structures where itreat_switch = "acr-das28-switch", itreat_switch = "acr-sdai-switch", and itreat_switch = "acr-cdai-switch". Here, we show parameters estimated from the inception cohort.

iviRA::acr2das28$inception
##          acr      mean     lower     upper
## 1:  ACR < 20  0.000000  0.000000  0.000000
## 2: ACR 20-50 -1.550000 -1.860000 -1.240000
## 3: ACR 50-70 -1.542727 -1.851273 -1.234182
## 4:   ACR 70+ -3.310000 -3.972000 -2.648000
iviRA::acr2sdai$inception
##          acr      mean     lower     upper
## 1:  ACR < 20   0.00000   0.00000   0.00000
## 2: ACR 20-50 -13.70000 -16.44000 -10.96000
## 3: ACR 50-70 -14.88182 -17.85818 -11.90545
## 4:   ACR 70+ -30.10000 -36.12000 -24.08000
iviRA::acr2cdai$inception
##          acr      mean     lower     upper
## 1:  ACR < 20   0.00000   0.00000   0.00000
## 2: ACR 20-50 -11.30000 -13.56000  -9.04000
## 3: ACR 50-70 -12.87273 -15.44727 -10.29818
## 4:   ACR 70+ -27.60000 -33.12000 -22.08000

The relationship between ACR response and HAQ is used when a model structure uses itreat_haq = "acr-haq". Evidence on mean HAQ by ACR response category is based on evidence reported in Carlson 2015.

##          acr  mean      se
## 1:  ACR < 20 -0.11 0.06765
## 2: ACR 20-50 -0.44 0.05657
## 3: ACR 50-70 -0.76 0.09059
## 4:   ACR 70+ -1.07 0.07489

When itreat_haq = "acr-eular-haq" is selected, the model must simulate the relationship between ACR response and EULAR response. The estimated relationship is based on the number of patients in each ACR response category for a given EULAR response category in the Veterans Affairs Rheumatoid Arthritis (VARA) database as reported in the 2016 NICE health technology assessment (HTA).

##           eular_none eular_moderate eular_good
## ACR < 20         755            136         57
## ACR 20-50          4             27         26
## ACR 50-70          2              2         10
## ACR 70+            0              2          2

When the H2 pathway is used, the change in HAQ at 6 months must be simulated after EULAR response is simulated. We base this relationship on statistical analysis of the British Society for Rheumatology Biologics Register (BSRBR) data by NICE (2016). In particular, the 2016 NICE HTA reported mean changes in HAQ by EULAR response category at 6 months for all patients.

##       eular   mean    se
## 1:     None  0.000 0.000
## 2: Moderate -0.317 0.048
## 3:     Good -0.672 0.112

Long-term HAQ progression

The progression of HAQ over time is modeled in two ways: first, assuming a constant annual rate, and second, by using a LCGM. The LCGM can currently only be used to model progression in the absence of tDMARDs while a constant annual rate can be used for all treatments.

The constant annual rate of progression for patients on cDMARDs (and NBT) are based on Wolfe and Michaud’s (2010) analysis of 3,829 RA patients who switched from non-biologic treatment to biologic treatment and participated in the National Data Bank for Rheumatic Diseases (NDB). The progression rate for patients on biologics are based on treatment-specific rates from Michaud et al. (2011) adjusted for baseline HAQ score, age, sex, education, smoking, BMI, comorbidity, and RA onset. Importantly, The average HAQ rate among patients on a biologic was -0.001 in the Michaud et al. analysis, which instills confidence that the reported HAQ progression rates for different bDMARDs can be directly compared with the overall annual HAQ progression rate of 0.031 reported in the 2010 Wolfe and Michaud analysis.

iviRA::haq.lprog$tx
##        sname    est          se   lower  upper
##  1:  cdmards  0.031 0.002551067  0.0260 0.0360
##  2: abtivmtx -0.001 0.001530640 -0.0040 0.0020
##  3:   adamtx -0.003 0.007398095 -0.0175 0.0115
##  4:      ada  0.012 0.006887882 -0.0015 0.0255
##  5:   etnmtx -0.005 0.003316387 -0.0115 0.0015
##  6:      etn -0.003 0.003826601 -0.0105 0.0045
##  7:   golmtx -0.001 0.001530640 -0.0040 0.0020
##  8:   ifxmtx -0.001 0.003571494 -0.0080 0.0060
##  9:  placebo  0.031 0.002551067  0.0260 0.0360
## 10:   tczmtx -0.001 0.001530640 -0.0040 0.0020
## 11:      tcz -0.001 0.001530640 -0.0040 0.0020
## 12:   czpmtx -0.001 0.001530640 -0.0040 0.0020
## 13: abtscmtx -0.001 0.001530640 -0.0040 0.0020
## 14:      nbt  0.031 0.002551067  0.0260 0.0360
## 15:   rtxmtx -0.001 0.001530640 -0.0040 0.0020
## 16:   tofmtx -0.001 0.001530640 -0.0040 0.0020
## 17:      rtx -0.001 0.001530640 -0.0040 0.0020
## 18:      tof -0.001 0.001530640 -0.0040 0.0020
## 19:      czp -0.001 0.001530640 -0.0040 0.0020
## 20:      gol -0.001 0.001530640 -0.0040 0.0020

The constant annual rate of HAQ progression is also adjusted by age based on differences between overall and age specific rates as reported by Michaud et al. (2011).

iviRA::haq.lprog$diff.age
##       age    est          se       lower       upper
## 1:   < 40 -0.020 0.004152837 -0.02813941 -0.01186059
## 2:  40-64 -0.008 0.001082326 -0.01012132 -0.00587868
## 3:  >= 65  0.017 0.002179633  0.01272800  0.02127200

The LCGM is based on the analysis by Norton et al. (2014). The parameters from the LCGM are contained in the list iviRA::haq.lcgm. The coefficient estimates are stored in iviRA::haq.lcgm$coef.

iviRA::haq.lcgm$coef
##                    var data parameter         est    se
##  1:        intercept_1 eras     beta1   0.6377896 0.058
##  2:           linear_1 eras     beta1  -1.0086461 0.074
##  3:        quadratic_1 eras     beta1  -0.6493539 0.027
##  4:            cubic_1 eras     beta1   1.3547476 0.003
##  5:        intercept_2 eras     beta2   0.9503192 0.058
##  6:           linear_2 eras     beta2  -0.1091174 0.020
##  7:        quadratic_2 eras     beta2  -3.3681108 0.002
##  8:            cubic_2 eras     beta2   3.6994158 0.064
##  9:        intercept_3 eras     beta3   1.2652436 0.064
## 10:           linear_3 eras     beta3  -0.1316658 0.056
## 11:        quadratic_3 eras     beta3  -2.5311796 0.021
## 12:            cubic_3 eras     beta3   3.5376129 0.002
## 13:        intercept_4 eras     beta4   1.9347732 0.063
## 14:           linear_4 eras     beta4  -0.5398728 0.073
## 15:        quadratic_4 eras     beta4   1.1961146 0.027
## 16:            cubic_4 eras     beta4  -0.1087785 0.003
## 17:         constant_2 eras    delta2  -3.4960000 0.622
## 18:              age_2 eras    delta2   0.0250000 0.007
## 19:           female_2 eras    delta2   0.8410000 0.196
## 20:            das28_2 eras    delta2   0.3040000 0.080
## 21: disease_duration_2 eras    delta2   0.0320000 0.016
## 22:     rheum_factor_2 eras    delta2   0.2140000 0.237
## 23:     acr_criteria_2 eras    delta2   0.2780000 0.225
## 24:            imdq4_2 eras    delta2   0.9930000 0.366
## 25:         constant_3 eras    delta3  -6.6860000 0.660
## 26:              age_3 eras    delta3   0.0370000 0.007
## 27:           female_3 eras    delta3   1.6940000 0.214
## 28:            das28_3 eras    delta3   0.5730000 0.076
## 29: disease_duration_3 eras    delta3   0.0460000 0.017
## 30:     rheum_factor_3 eras    delta3   0.3150000 0.250
## 31:     acr_criteria_3 eras    delta3   0.4130000 0.236
## 32:            imdq4_3 eras    delta3   1.1190000 0.342
## 33:         constant_4 eras    delta4 -12.0550000 1.102
## 34:              age_4 eras    delta4   0.0820000 0.011
## 35:           female_4 eras    delta4   1.9760000 0.269
## 36:            das28_4 eras    delta4   0.8000000 0.086
## 37: disease_duration_4 eras    delta4   0.0420000 0.021
## 38:     rheum_factor_4 eras    delta4   0.2980000 0.290
## 39:     acr_criteria_4 eras    delta4   0.9390000 0.316
## 40:            imdq4_4 eras    delta4   1.4290000 0.381
##                    var data parameter         est    se

The coefficient vectors beta1, beta2, beta3, and beta4, contain coefficients from a third order polynomial regression predicting the growth of the HAQ score over time for members of class 1, 2, 3, and 4, respectively; delta2, delta3, and delta4 are coefficient vectors containing predictors of membership in classes 2-4 with class 1 being the reference group. A variance-covariance matrix, iviRA::haq.lcgm$vcov, is used to sample parameters from a multivariate normal distribution. The indices of the parameters of the variance-covariance matrix correspond to the row numbers in iviRA::haq.lcgm$coef.

Time to treatment discontinuation

Time to treatment discontinuation depends on the pathway (S1-S6) used to model treatment switching. If S1 is selected, a single treatment discontinuation curve is used for all patients. The parameters needed to simulate this curve are stored in. In S2-S5, time to treatment discontinuation is stratified by the level of disease activity, and in S6 treatment duration depends on EULAR response (st. The parameters used to simulate treatment duration based on S1, S2-S5, and S6 are stored in iviRA::ttd.all, iviRA::ttd.da, and iviRA::ttd.eular, respectively.

The treatment duration curves are estimated using parametric survival models with the flexsurv package. Parameters from the survival analyses used in the IVI-RA model are stored in two-level lists. The names in the first level correspond to the name of the parametric survival model.

names(iviRA::ttd.all)
## [1] "exponential" "weibull"     "gompertz"    "gamma"       "llogis"     
## [6] "lnorm"       "gengamma"

7 parametric distributions are currently supported including the exponential, Weibull, Gompertz, gamma, log-logistic, lognormal, and generalized gamma distributions. The parameter of the survival model are stored in the second level of the list. For example, we can examine results from a Weibull fit.

names(iviRA::ttd.all$weibull)
## [1] "est"        "vcov"       "loc.index"  "anc1.index" "anc2.index"

Point estimates are stored in the element est and the variance-covariance matrix is stored in the element vcov.

iviRA::ttd.all$weibull$est
##      shape      scale 
## -0.1028306  3.7576710
iviRA::ttd.all$weibull$vcov
##               shape         scale
## shape  1.871894e-04 -5.220527e-05
## scale -5.220527e-05  3.644134e-04

The parametric distributions are characterized according to a location parameter governing the mean of the distribution. Models can also contain ancillary parameters which describe the shape, variance, or higher moments of a distribution. The location and ancillary parameters associated with each distribution are shown in the table below.

Distribution Location Ancillary 1 Ancillary 2 iviRA name
exp Exponential rate - - exponential
weibull Weibull scale shape - weibull
gompertz Gompertz rate shape - gompertz
gamma Gamma rate shape - gamma
llogis Log-logistic scale shape - llogis
lnorm Lognormal meanlog sdlog - lnorm
gengamma Generalized gamma mu sigma Q gengamma

The iviRA package currently only allows for models in which the location parameter, but not the ancillary parameters, depends on covariates. The indices of the location parameter, first ancillary parameter, and second ancillary parameter in the est and vcov objects are contained in the loc.index, anc1.index, and anc2.index objects, respectively. For example, in a regression model with no covariates, the index of the location parameter is 2 and the index of the first ancillary parameter is 1.

iviRA::ttd.all$weibull$loc.index
## [1] 2
iviRA::ttd.all$weibull$anc1.index
## [1] 1
iviRA::ttd.all$weibull$anc2.index
## [1] NA

Serious infections

Serious infection rates are based on the NMA and Cochrane review by Singh et al. (2011). In accordance with the 2016 NICE report we assume a constant serious infection rate of 0.026 for patients using cDMARDs and a constant rate of exp(iviRA::ttsi[sname == "etnmtx", lograte]) for patient using tDMARDs. A such we constantly assume that the time to a serious infection is exponentially distributed and that it does not vary across tDMARDs.

##        sname    lograte lograte_se
##  1:  cdmards  -3.649659   0.135922
##  2: abtivmtx  -3.352407   0.135922
##  3:   adamtx  -3.352407   0.135922
##  4:      ada  -3.352407   0.135922
##  5:   etnmtx  -3.352407   0.135922
##  6:      etn  -3.352407   0.135922
##  7:   golmtx  -3.352407   0.135922
##  8:   ifxmtx  -3.352407   0.135922
##  9:  placebo -23.025851   0.135922
## 10:   tczmtx  -3.352407   0.135922
## 11:      tcz  -3.352407   0.135922
## 12:   czpmtx  -3.352407   0.135922
## 13: abtscmtx  -3.352407   0.135922
## 14:      nbt -23.025851   0.135922
## 15:   rtxmtx  -3.352407   0.135922
## 16:   tofmtx  -3.352407   0.135922
## 17:      rtx  -3.352407   0.135922
## 18:      tof  -3.352407   0.135922
## 19:      czp  -3.352407   0.135922
## 20:      gol  -3.352407   0.135922

Mortality

Baseline mortality for patients without RA are based on gender-specific US life tables. More specifically, we use the probability of mortality, qx, between ages \(x\) and \(x+1\).

head(iviRA::lifetable.male)
##    age       qx     lx  dx     L      Tx   ex
## 1:   0 0.006575 100000 658 99427 7629389 76.3
## 2:   1 0.000445  99342  44 99320 7529961 75.8
## 3:   2 0.000301  99298  30 99283 7430641 74.8
## 4:   3 0.000240  99268  24 99256 7331358 73.9
## 5:   4 0.000183  99244  18 99235 7232101 72.9
## 6:   5 0.000169  99226  17 99218 7132866 71.9
head(iviRA::lifetable.female)
##    age       qx     lx  dx     L      Tx   ex
## 1:   0 0.005516 100000 552 99514 8105137 81.1
## 2:   1 0.000382  99448  38 99429 8005623 80.5
## 3:   2 0.000225  99410  22 99399 7906193 79.5
## 4:   3 0.000175  99388  17 99379 7806794 78.5
## 5:   4 0.000151  99371  15 99363 7707415 77.6
## 6:   5 0.000132  99356  13 99349 7608052 76.6

The probability of mortality is adjusted according to a patient’s baseline HAQ score as well as changes in the HAQ score from baseline. We use the log odds ratio of the effect of baseline HAQ on mortality of 0.798 from Wolfe et al. (2003).

##    var   or or_se or_lower or_upper     logor logor_lower logor_upper
## 1: haq 2.22  0.24     1.79     2.75 0.7975072   0.5822156    1.011601
##     logor_se
## 1: 0.1095391

The effect of a change in HAQ from baseline on mortality is based on the effect of 0.25 change in the HAQ score on log hazard of mortality as reported in Michaud et al. (2012). The log hazard ratios are stratified by time-period, as reported in the original paper.

##    month   hr hr_lower hr_upper     loghr loghr_lower loghr_upper
## 1:     6 1.12     1.08     1.17 0.1133287  0.07696104   0.1570037
## 2:    12 1.16     1.11     1.21 0.1484200  0.10436002   0.1906204
## 3:    24 1.16     1.10     1.21 0.1484200  0.09531018   0.1906204
## 4:    36 1.21     1.14     1.28 0.1906204  0.13102826   0.2468601
## 5:    48 1.19     1.11     1.27 0.1739533  0.10436002   0.2390169
##      loghr_se
## 1: 0.02041943
## 2: 0.02200559
## 3: 0.02431427
## 4: 0.02954948
## 5: 0.03435188

Utility

Two models can be used to simulate utility. If the model structure utility_model = "mixure" is chosen, then the Hernandez Alava et al. (2013) model is used; otherwise, if utility_model = "wailoo", the logistic regression equation reported in Wailoo et al. (2006) is used.

Point estimates (i.e., log odds ratios) and standard errors associated with each variable in the Wailoo model are available in the iviRA::utility.wailoo object.

##            var     est     se
## 1:         int  2.0734 0.0263
## 2:         age  0.0058 0.0004
## 3:     dis_dur  0.0023 0.0004
## 4:        haq0 -0.2004 0.0101
## 5:        male -0.2914 0.0118
## 6: prev_dmards  0.0249 0.0028
## 7:         haq -0.8647 0.0103

The information needed to sample from the mixture model are stored in the list iviRA::utility.mixture. The setup is the same as with the LCGM in that the list contains a table of coefficient estimates as well as the corresponding variance covariance matrix.

names(iviRA::utility.mixture)
## [1] "coef" "vcov"
iviRA::utility.mixture$coef
##            var parameter           est
## 1          HAQ     beta1  -0.089819885
## 2         HAQ2     beta1   0.000525916
## 3     Pain/100     beta1  -0.057952582
## 4      Age/10m     beta1   0.004871839
## 5     Age/10m2     beta1   0.000277046
## 6          HAQ     beta2   0.054383137
## 7         HAQ2     beta2  -0.050915914
## 8     Pain/100     beta2  -0.384088121
## 9      Age/10m     beta2   0.029093294
## 10    Age/10m2     beta2   0.002262836
## 11         HAQ     beta3  -0.141482101
## 12        HAQ2     beta3   0.015465378
## 13    Pain/100     beta3  -0.083925893
## 14     Age/10m     beta3   0.003728168
## 15    Age/10m2     beta3   0.000718860
## 16         HAQ     beta4  -0.195798978
## 17        HAQ2     beta4   0.034715726
## 18    Pain/100     beta4  -0.012745796
## 19     Age/10m     beta4  -0.004299891
## 20    Age/10m2     beta4   0.000227100
## 21  Intercept1    alpha1   0.814127363
## 22  Intercept2    alpha2   0.426581217
## 23  Intercept3    alpha3   0.329688793
## 24  Intercept4    alpha4   1.022019574
## 25        Male     alpha  -0.026542757
## 26  Variance 1  epsilon1   0.002517113
## 27  Variance 2  epsilon2   0.023997163
## 28  Variance 3  epsilon3   0.002224878
## 29  Variance 4  epsilon4   0.004379523
## 30    Variance        mu   0.002605149
## 31 Intercept 1    delta1  -1.274627693
## 32         HAQ    delta1   0.242049855
## 33    Pain/100    delta1  23.467307110
## 34   Pain/1002    delta1 -21.551251860
## 35 Intercept 2    delta2  -6.631010968
## 36         HAQ    delta2   2.193565335
## 37    Pain/100    delta2  18.371864330
## 38   Pain/1002    delta2 -13.800080430
## 39 Intercept 3    delta3  -7.476836039
## 40         HAQ    delta3   1.051656593
## 41    Pain/100    delta3  25.339552780
## 42   Pain/1002    delta3 -16.962219250

Quality-adjusted life-years are also decreased in each period that there is a serious infection and is assumed to vary by 20 percent in either direction.

formals(sample_pars)$si_ul
## [1] 0.156
formals(sample_pars)$si_ul_range
## [1] 0.2

Costs

The model incorporates 3 main types of health care sector costs: (1) drug acquisition and administration costs, (2) hospitalization costs, and (3) general management costs. Productivity losses can also be included in the analysis if the analyst would like to take a societal perspective.

Treatment costs

Drug acquisition and administration costs are separated into costs associated with the first 6 months of treatment and costs after the first 6 months. Costs are a function of the FDA recommended dose and frequency of administration, strength and dosage form, the price, and infusion. Prices with the simulation are reduced according to assumed discounts and rebates. If dosing is weight based, then the dose is calculated by multiplying the dose value (e.g., init_dose_val and ann_dose_val) by each patients weight in kilograms; otherwise the dose is simply equal to the dose value. For example, the dose for infliximab during the first six months is equal to 3 multiplied by the weight of the patient.

Information related to treatment costs are stored in the list iviRA::tx.cost. Data related to the costs of treatment for each medication are in the data table iviRA::tx.cost$cost.

library(knitr)
kable(iviRA::tx.cost$cost)
name sname dosage strength_dosage_form init_dose_val ann_dose_val dose_unit init_num_doses ann_num_doses strength_val strength_unit price_per_unit discount_lower discount_upper infusion_cost loading_dose weight_based
etanercept etn 50 mg QW 50 mg/0.98 mL syringe or pen injector 50 50 mg 26 52.000000 50.0 mg 1110.50000 0.2 0.3 0 0 0
adalimumab ada 40 mg EOW 40 mg/0.8 mL syringe or pen injector 40 40 mg 13 26.000000 40.0 mg 2220.62000 0.2 0.3 0 0 0
infliximab ifx 3 mg/kg at 0, 2, and 6 weeks, 3mg/kg Q8W, 6 mg/kg Q6W after 6 months 100 mg vial 3 6 mg/kg 5 8.666667 100.0 mg 1113.27000 0.2 0.3 164 0 1
golimumab gol 50 mg QM 50 mg/0.5 mL syringe or pen injector 50 50 mg 6 12.000000 50.0 mg 3811.18000 0.2 0.3 0 0 0
certolizumab pegol czp 400 mg at weeks 0, 2, 4 then 200 mg Q2W 400 mg kit or syringe kit (200 mg 2) 200 200 mg 16 26.000000 200.0 mg 1839.94000 0.2 0.3 0 0 0
abatacept IV abtiv 750 mg IV at weeks 0, 2, 4 then Q4W 250mg vial 750 750 mg 8 13.000000 250.0 mg 931.16000 0.2 0.3 164 0 0
abatacept SC abtsc 125 mg SC QW with IV loading dose 125mg/ml syringe 125 125 mg 26 52.000000 125.0 mg 957.14000 0.2 0.3 164 1 0
tocilizumab tcz 162 mg SC EOW 162 mg/0.9 mL syringe 162 162 mg 13 26.000000 162.0 mg 898.31000 0.2 0.3 0 0 0
rituximab rtx 1000 mg at weeks 0, 2; then Q24 W 500 mg/50ml vial 1000 1000 mg 4 4.333333 500.0 mg 4176.10000 0.2 0.3 164 0 0
tofaticinib citrate tof 5 mg BID 5mg tablet 5 5 mg 364 728.000000 5.0 mg 63.26000 0.2 0.3 0 0 0
methotrexate monotherapy cdmards 15mg QW 15 mg injection 15 15 mg 26 52.000000 15.0 mg 32.41875 0.2 0.3 0 0 0
hydroxychlorquine sulfate hcl 400mg daily 200 mg tablet 400 400 mg 182 364.000000 200.0 mg 3.18000 0.2 0.3 0 0 0
sulfazalazine ssz 1-2 g daily 500 mg tablet 2 2 g 182 364.000000 0.5 g 0.47000 0.2 0.3 0 0 0

The variable init_num_doses is the number of doses during the first 6 months after beginning a new treatment and the variable ann_num_doses is the number of doses per year after the first 6 months. The variables discount_lower and discount_upper are lower and upper bounds for the amount by which prices (price_per_unit) are decreased after discounts and rebates.

Costs for individual therapeutic agents are used to calculate costs for the (combination) treatments in a treatment sequence by using the lookup data table iviRA::tx.cost$lookup.

iviRA::tx.cost$lookup
##        sname  agent1  agent2
##  1:  cdmards cdmards      NA
##  2: abtivmtx   abtiv cdmards
##  3:   adamtx     ada cdmards
##  4:      ada     ada      NA
##  5:   etnmtx     etn cdmards
##  6:      etn     etn      NA
##  7:   golmtx     gol cdmards
##  8:   ifxmtx     ifx cdmards
##  9:  placebo placebo      NA
## 10:   tczmtx     tcz cdmards
## 11:      tcz     tcz      NA
## 12:   czpmtx     czp cdmards
## 13: abtscmtx   abtsc cdmards
## 14:      nbt cdmards      NA
## 15:   rtxmtx     rtx cdmards
## 16:   tofmtx     tof cdmards
## 17:      rtx     rtx      NA
## 18:      tof     tof      NA
## 19:      czp     czp      NA
## 20:      gol     gol      NA

The treatment names in the first column (sname) must match the names of the treatments selected in a treatment sequence. All subsequent columns refer to the names of the therapeutic agents used for a particular treatment. For example, the treatment etnmtx contains both etanercept and methotrexate/conventional cDMARDs.

Other health care sector costs

Hospitalization costs are estimated as a function of the HAQ score. In particular, Carlson 2015, report estimates of the number of days in the hospital for a given HAQ score range as well as the cost of a single day in the hospital.

##          haq days_mean days_se cost_pday_mean cost_pday_se
## 1: 0 to <0.5      0.26     0.5       1292.358     197.3144
## 2: 0.5 to <1      0.13     0.5       1292.358     197.3144
## 3: 1 to <1.5      0.51     0.5       1292.358     197.3144
## 4: 1.5 to <2      0.72     0.5       1292.358     197.3144
## 5: 2 to <2.5      1.86     0.5       1292.358     197.3144
## 6:      >2.5      4.16     0.5       1292.358     197.3144

Management costs are based on mean costs and cost ranges reported in Claxton (2016).

##                service       est     lower     upper        se
## 1:          chest_xray 112.60353 100.20681 125.00025  6.324973
## 2:          xray_visit  54.75218  46.48770  63.01666  4.216649
## 3: outpatient_followup 193.18221 164.25653 222.10788 14.758270
## 4:                 tst  31.50833  31.50833  31.50833  0.000000

Productivity losses

We use estimates of productivity loss from Wolfe et al. (2005). The authors provide estimates of the annual loss in earnings associated with a 1-unit increase in HAQ, which we convert to 2017 dollars.

##         est       se
## 1: 5965.952 1576.608