library("iviRA")
load("output.RData")

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)

The number of samples is also returned for convenience.

parsamp$n

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)

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)
iviRA::nma.acr.naive$mean
iviRA::nma.das28.naive$mean
iviRA::nma.haq.naive$mean

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
iviRA::acr2sdai$inception
iviRA::acr2cdai$inception

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.

iviRA::acr2haq

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).

iviRA::acr2eular

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.

iviRA::eular2haq

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

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

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

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)

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)

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

iviRA::ttd.all$weibull$est
iviRA::ttd.all$weibull$vcov

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.

library(knitr)
exp <- c("Exponential", "rate", "-", "-", "exponential")
weibull <- c("Weibull", "scale", "shape", "-", "weibull")
gompertz <- c("Gompertz", "rate", "shape", "-", "gompertz")
gamma <- c("Gamma", "rate", "shape", "-", "gamma")
llogis <- c("Log-logistic", "scale", "shape", "-", "llogis")
lnorm <- c("Lognormal", "meanlog", "sdlog", "-", "lnorm")
gengamma <- c("Generalized gamma", "mu", "sigma", "Q", "gengamma")
survdist.table <- data.frame(rbind(exp, weibull, gompertz, gamma, llogis, lnorm, gengamma))
colnames(survdist.table) <- c("Distribution", "Location", "Ancillary 1", "Ancillary 2",
                              "iviRA name")
kable(survdist.table)

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 r iviRA::ttd.all$weibull$loc.index and the index of the first ancillary parameter is r iviRA::ttd.all$weibull$anc1.index.

iviRA::ttd.all$weibull$loc.index
iviRA::ttd.all$weibull$anc1.index
iviRA::ttd.all$weibull$anc2.index

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 r exp(iviRA::ttsi[sname == "cdmards", lograte]) 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

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)
head(iviRA::lifetable.female)

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 r round(iviRA::mort.or$logor, 3) from Wolfe et al. (2003).

iviRA::mort.or

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.

iviRA::mort.hr.haqdif

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.

iviRA::utility.wailoo

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)
iviRA::utility.mixture$coef

Quality-adjusted life-years are also decreased in each period that there is a serious infection and is assumed to vary by r formatC(formals(sample_pars)$si_ul_range * 100, digits = 2) percent in either direction.

formals(sample_pars)$si_ul
formals(sample_pars)$si_ul_range

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)

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

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.

iviRA::hosp.cost

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

iviRA::mgmt.cost

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.

iviRA::prod.loss
save.image("output.RData")


InnovationValueInitiative/IVI-RA documentation built on Oct. 20, 2020, 10:02 p.m.