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] 100Users 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.
## 5608 6581 7296 7175 7838 8394To 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.
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
## 0.5499991 0.0000000 0.6231056 1.2020253
## d_cdmards d_abtivmtx d_abtscmtx d_adamtx
## 0.0000000 -0.9007627 -0.9005787 -0.7743447
## d_ada d_adabiosbwwdmtx d_anamtx d_bct
## -0.5470945 -0.7678458 -0.4505512 -0.8477756
## d_bctmtx d_czp d_czpmtx d_etn
## -0.6992776 -0.7559150 -1.1912514 -0.7941323
## d_etnmtx d_etnbiosszzsmtx d_etnbiosykromtx d_gol
## -0.7628515 -0.5468121 -0.8449216 NA
## d_golmtx d_ifxmtx d_ifxbiosqbtxmtx d_placebo
## -0.8450899 -0.7734014 NA 0.3774837
## d_nbt d_rtx d_rtxmtx d_sar
## 0.0000000 -0.5237583 -0.7040334 -0.9321168
## d_sarmtx d_triple d_tcz d_tczmtx
## -0.8696777 -0.5961814 -1.0482577 -0.9883258
## d_tofmtx d_tof d_upamtx
## -0.7667502 -0.5347355 -0.7271539iviRA::nma.das28.naive$mean## A d_cdmards d_abtivmtx d_abtscmtx
## -0.9906912 0.0000000 -1.3419370 -1.2905265
## d_adamtx d_ada d_adabiosbwwdmtx d_anamtx
## -1.1957996 -0.3846224 -1.2555142 NA
## d_bct d_bctmtx d_czp d_czpmtx
## NA NA NA -2.0094080
## d_etn d_etnmtx d_etnbiosszzsmtx d_etnbiosykromtx
## -1.5127713 -1.5699353 -1.5696227 -1.6663967
## d_gol d_golmtx d_ifxmtx d_ifxbiosqbtxmtx
## NA -1.4525273 -0.9384479 NA
## d_placebo d_nbt d_rtx d_rtxmtx
## 0.4588825 0.0000000 -0.7546151 -0.9851809
## d_sar d_sarmtx d_triple d_tcz
## -1.2770095 NA -1.3012420 -1.7762413
## d_tczmtx d_tofmtx d_tof d_upamtx
## -1.9362256 -0.9513472 -0.7044028 -1.2288702iviRA::nma.haq.naive$mean## A d_cdmards d_abtivmtx d_abtscmtx
## -0.2328891 0.0000000 -0.2317792 -0.2229830
## d_adamtx d_ada d_adabiosbwwdmtx d_anamtx
## -0.3204838 -0.1601882 NA -0.1104685
## d_bct d_bctmtx d_czp d_czpmtx
## NA NA -0.3193571 -0.3883754
## d_etn d_etnmtx d_etnbiosszzsmtx d_etnbiosykromtx
## -0.1512577 -0.3369290 NA -0.2390994
## d_gol d_golmtx d_ifxmtx d_ifxbiosqbtxmtx
## NA -0.3436631 -0.2168854 -0.2400438
## d_placebo d_nbt d_rtx d_rtxmtx
## 0.1710748 0.0000000 NA -0.2510926
## d_sar d_sarmtx d_triple d_tcz
## -0.3404208 -0.2506528 -0.2673432 -0.2480185
## d_tczmtx d_tofmtx d_tof d_upamtx
## -0.2386515 -0.4050463 -0.3117164 -0.3516575In 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.
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.648000iviRA::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.08000iviRA::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.08000The 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.
iviRA::acr2haq## 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.07489When 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).
iviRA::acr2eular## 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 2When 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.
iviRA::eular2haq## eular mean se
## 1: None 0.000 0.000
## 2: Moderate -0.317 0.048
## 3: Good -0.672 0.112The 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: abtscmtx -0.001 0.001530640 -0.0040 0.0020
## 4: adamtx -0.003 0.007398095 -0.0175 0.0115
## 5: ada 0.012 0.006887882 -0.0015 0.0255
## 6: adabiosbwwdmtx -0.003 0.007398095 -0.0175 0.0115
## 7: anamtx -0.001 0.001530640 -0.0040 0.0020
## 8: bct -0.001 0.001530640 -0.0040 0.0020
## 9: bctmtx -0.001 0.001530640 -0.0040 0.0020
## 10: czp -0.001 0.001530640 -0.0040 0.0020
## 11: czpmtx -0.001 0.001530640 -0.0040 0.0020
## 12: etn -0.003 0.003826601 -0.0105 0.0045
## 13: etnmtx -0.005 0.003316387 -0.0115 0.0015
## 14: etnbiosszzsmtx -0.005 0.003316387 -0.0115 0.0015
## 15: etnbiosykromtx -0.005 0.003316387 -0.0115 0.0015
## 16: gol -0.001 0.001530640 -0.0040 0.0020
## 17: golmtx -0.001 0.001530640 -0.0040 0.0020
## 18: ifxmtx -0.001 0.003571494 -0.0080 0.0060
## 19: ifxbiosqbtxmtx -0.001 0.003571494 -0.0080 0.0060
## 20: placebo 0.031 0.002551067 0.0260 0.0360
## 21: nbt 0.031 0.002551067 0.0260 0.0360
## 22: rtx -0.001 0.001530640 -0.0040 0.0020
## 23: rtxmtx -0.001 0.001530640 -0.0040 0.0020
## 24: sar -0.001 0.001530640 -0.0040 0.0020
## 25: sarmtx -0.001 0.001530640 -0.0040 0.0020
## 26: triple -0.001 0.001530640 -0.0040 0.0020
## 27: tcz -0.001 0.001530640 -0.0040 0.0020
## 28: tczmtx -0.001 0.001530640 -0.0040 0.0020
## 29: tofmtx -0.001 0.001530640 -0.0040 0.0020
## 30: tof -0.001 0.001530640 -0.0040 0.0020
## 31: upamtx -0.001 0.001530640 -0.0040 0.0020
## sname est se lower upperThe 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.02127200The 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 seThe 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 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.7576710iviRA::ttd.all$weibull$vcov## shape scale
## shape 1.871894e-04 -5.220527e-05
## scale -5.220527e-05 3.644134e-04The 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] 2iviRA::ttd.all$weibull$anc1.index## [1] 1iviRA::ttd.all$weibull$anc2.index## [1] NASerious 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.
iviRA::ttsi## sname lograte lograte_se
## 1: cdmards -3.649659 0.135922
## 2: abtivmtx -3.352407 0.135922
## 3: abtscmtx -3.352407 0.135922
## 4: adamtx -3.352407 0.135922
## 5: ada -3.352407 0.135922
## 6: adabiosbwwdmtx -3.352407 0.135922
## 7: anamtx -3.352407 0.135922
## 8: bct -3.352407 0.135922
## 9: bctmtx -3.352407 0.135922
## 10: czp -3.352407 0.135922
## 11: czpmtx -3.352407 0.135922
## 12: etn -3.352407 0.135922
## 13: etnmtx -3.352407 0.135922
## 14: etnbiosszzsmtx -3.352407 0.135922
## 15: etnbiosykromtx -3.352407 0.135922
## 16: gol -3.352407 0.135922
## 17: golmtx -3.352407 0.135922
## 18: ifxmtx -3.352407 0.135922
## 19: ifxbiosqbtxmtx -3.352407 0.135922
## 20: placebo -23.025851 0.135922
## 21: nbt -23.025851 0.135922
## 22: rtx -3.352407 0.135922
## 23: rtxmtx -3.352407 0.135922
## 24: sar -3.352407 0.135922
## 25: sarmtx -3.352407 0.135922
## 26: triple -3.649659 0.135922
## 27: tcz -3.352407 0.135922
## 28: tczmtx -3.352407 0.135922
## 29: tofmtx -3.352407 0.135922
## 30: tof -3.352407 0.135922
## 31: upamtx -3.352407 0.135922
## sname lograte lograte_seBaseline 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.9head(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.6The 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).
iviRA::mort.or## 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.1095391The 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.
iviRA::mort.hr.haqdif## 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.03435188Two 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.
iviRA::utility.wailoo## 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.0103The 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.962219250Quality-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.156formals(sample_pars)$si_ul_range## [1] 0.2The 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.
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.
| 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 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| abatacept IV | abtiv | 750 mg IV at weeks 0, 2, 4 then Q4W | 250mg vial | 750 | 750 | mg | 8 | 13.000000 | 250.0 | mg | 1167.33000 | 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 | 1094.72250 | 0.2 | 0.3 | 164 | 1 | 0 |
| adalimumab | ada | 40 mg EOW | 40 mg/0.8 mL syringe or pen injector | 40 | 40 | mg | 13 | 26.000000 | 40.0 | mg | 2587.04500 | 0.2 | 0.3 | 0 | 0 | 0 |
| adalimumab-bwwd (biosimilar Samsung Bioepis) | adabiosbwwd | 40 mg EOW | 40 mg/0.8 mL syringe or pen injector | 40 | 40 | mg | 13 | 26.000000 | 40.0 | mg | 2069.63600 | 0.2 | 0.3 | 0 | 0 | 0 |
| anakinra | ana | 100 mg daily | 100 mg syringe | 100 | 100 | mg | 182 | 364.000000 | 100.0 | mg | 147.06714 | 0.2 | 0.3 | 0 | 0 | 0 |
| baricitinib | bct | 2 mg daily | 2 mg tablet | 2 | 2 | mg | 182 | 364.000000 | 2.0 | mg | 71.23000 | 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 | 4327.43000 | 0.2 | 0.3 | 0 | 0 | 0 |
| etanercept | etn | 50 mg QW | 50 mg/0.98 mL syringe or pen injector | 50 | 50 | mg | 26 | 52.000000 | 50.0 | mg | 1293.51500 | 0.2 | 0.3 | 0 | 0 | 0 |
| etanercept-szzs (biosimilar Sandoz) | etnbiosszzs | 50 mg QW | 50 mg/0.98 mL syringe or pen injector | 50 | 50 | mg | 26 | 52.000000 | 50.0 | mg | 1034.81200 | 0.2 | 0.3 | 0 | 0 | 0 |
| etanercept-ykro (biosimilar Samsung Bioepis) | etnbiosykro | 50 mg QW | 50 mg/0.98 mL syringe or pen injector | 50 | 50 | mg | 26 | 52.000000 | 50.0 | mg | 1034.81200 | 0.2 | 0.3 | 0 | 0 | 0 |
| golimumab | gol | 50 mg QM | 50 mg/0.5 mL syringe or pen injector | 50 | 50 | mg | 6 | 12.000000 | 50.0 | mg | 4809.02000 | 0.2 | 0.3 | 0 | 0 | 0 |
| hydroxychloroquine sulfate | hcl | 400 mg daily | 200 mg tablet | 400 | 400 | mg | 182 | 364.000000 | 200.0 | mg | 3.88200 | 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 | 1167.82000 | 0.2 | 0.3 | 164 | 0 | 1 |
| infliximab-qbtx (biosimilar Pfizer) | ifxbiosqbtx | 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 | 946.28000 | 0.2 | 0.3 | 164 | 0 | 1 |
| methotrexate | cdmards | 15mg QW | 15 mg injection | 15 | 15 | mg | 26 | 52.000000 | 15.0 | mg | 5.81000 | 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 | 4697.60000 | 0.2 | 0.3 | 164 | 0 | 0 |
| sarilumab | sar | 200 mg EOW | 200mg/1.14 mL syringe | 200 | 200 | mg | 13 | 26.000000 | 200.0 | mg | 1457.57456 | 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.24202 | 0.2 | 0.3 | 0 | 0 | 0 |
| tocilizumab | tcz | 162 mg SC EOW | 162 mg/0.9 mL syringe | 162 | 162 | mg | 13 | 26.000000 | 162.0 | mg | 1014.26000 | 0.2 | 0.3 | 0 | 0 | 0 |
| tofaticinib | tof | 5 mg BID | 5 mg tablet | 5 | 5 | mg | 364 | 728.000000 | 5.0 | mg | 74.67717 | 0.2 | 0.3 | 0 | 0 | 0 |
| upadacitinib | upa | 15 mg daily | 15 mg tablet | 15 | 15 | mg | 182 | 364.000000 | 15.0 | mg | 163.88900 | 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: abtscmtx abtsc cdmards
## 4: adamtx ada cdmards
## 5: ada ada <NA>
## 6: adabiosbwwdmtx adabiosbwwd cdmards
## 7: anamtx ana cdmards
## 8: bct bct <NA>
## 9: bctmtx bct cdmards
## 10: czp czp <NA>
## 11: czpmtx czp cdmards
## 12: etn etn <NA>
## 13: etnmtx etn cdmards
## 14: etnbiosszzsmtx etnbiosszzs cdmards
## 15: etnbiosykromtx etnbiosykro cdmards
## 16: gol gol <NA>
## 17: golmtx gol cdmards
## 18: ifxmtx ifx cdmards
## 19: ifxbiosqbtxmtx ifxbiosqbtx cdmards
## 20: placebo placebo <NA>
## 21: nbt nbt <NA>
## 22: rtx rtx <NA>
## 23: rtxmtx rtx cdmards
## 24: sar sar <NA>
## 25: sarmtx sar cdmards
## 26: triple hcl cdmards
## 27: tcz tcz <NA>
## 28: tczmtx tcz cdmards
## 29: tofmtx tof cdmards
## 30: tof tof <NA>
## 31: upamtx upa cdmards
## sname agent1 agent2The 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.
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.
iviRA::hosp.cost## haq days_mean days_se cost_pday_mean cost_pday_se
## 1: 0 to <0.5 0.26 0.5 1347.088 205.6705
## 2: 0.5 to <1 0.13 0.5 1347.088 205.6705
## 3: 1 to <1.5 0.51 0.5 1347.088 205.6705
## 4: 1.5 to <2 0.72 0.5 1347.088 205.6705
## 5: 2 to <2.5 1.86 0.5 1347.088 205.6705
## 6: >2.5 4.16 0.5 1347.088 205.6705Management costs are based on mean costs and cost ranges reported in Claxton (2016).
iviRA::mgmt.cost## service est lower upper se
## 1: chest_xray 117.37215 104.45045 130.29386 6.592828
## 2: xray_visit 57.07086 48.45639 65.68533 4.395218
## 3: outpatient_followup 201.36323 171.21259 231.51388 15.383265
## 4: tst 32.84267 32.84267 32.84267 0.000000We 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.
iviRA::prod.loss## est se
## 1: 6218.603 1643.376