Data sources

We performed cohort analyses on new-onset RA patients with public prescription drug coverage identified in the Quebec Health Insurance Program (RAMQ) databases from 2002–2011. In principle, as in each province in Canada, all Quebec residents (approximately 7.6 million persons) have access to comprehensive health care in terms of physician visits and hospitalization. Administrative databases record hospitalization and physician billing data for all residents of Quebec. The data include hospitalization discharge diagnoses (a primary diagnosis and up to 15 non-primary diagnoses per hospitalization, abstracted by medical records clerks), and physician-visit billing-claim diagnostic codes (a single diagnostic code is allowed per visit). Those residents who are beneficiaries of the provincial drug plan (which includes all seniors, and any non-seniors without private drug insurance, representing approximately 40 % of residents) can also be linked to the province’s dispensed prescription claims database (which includes information about drug, dose, duration, and dates).

We received ethics approval from the Quebec Commission for Access to Information. All data were de-nominalized, and identifying health card numbers were scrambled by the Quebec Health Insurance Board. As the data are anonymous, no informed consent is required from individual patients.

Study population

To establish a population-based incident RA sample, all physician visits with an RA diagnosis code between 1 January 2002 and 31 December 2011 were identified. To increase the positive predictive value and specificity of the RA case definition, cases required at least three visit billing codes using International Classification of Diseases (ICD)-9, code 714 over a 3-year period, at least one of which was by an internist or a rheumatologist. This represents an adaptation of the RA diagnosis that was recently validated by Widdifield et al. (2013) using Ontario administrative data [7]. To identify truly incident (rather than prevalent) RA cases, we removed any patient who had any billing codes for RA prior to 1 January 2002, as data were available from 1989. Patients were followed from cohort entry (time when they fulfilled the RA case definition) until their first joint replacement surgery, or were censored at death date, or the end of study period (31 December 2011), whichever came first. We restricted our analyses to those patients who had drug coverage by the public drug program at the cohort entry and during at least 80 % of their follow-up time and who had follow up longer than 1 year.

Exposure assessment

For each prescription of MTX, or other DMARDs (sulfasalazine, chloroquine, hydroxychloroquine, leflunomide, cyclosporine, minocycline, penicillamine, and cyclophosphamide), anti-TNF inhibitors, other biologic DMARDs (anakinra, rituximab, abatacept) cyclooxygenase-2 inhibitors (COXIBs), nonselective nonsteroidal anti-inflammatory drugs (NSAIDs) and systemic steroids, the start date, number of pills, dosage, and days supplied were retrieved from the prescription claims and used to construct the daily drug exposure matrix [8]. The daily exposure matrix was then used to calculate time-dependent measures of cumulative duration of use of a specific drug, or class of drugs until a given day during the follow up [9]. For overlapping prescriptions of the same drug, the individual was assumed to have had prescriptions refilled early and completed the first prescription before starting the second. The same rule was applied to overlapping prescriptions of the same drug but with different doses and overlapping prescriptions for different drugs within the same drug class (anti-TNF agents, systemic steroids, COXIBs or NSAIDs). However, given that combination use is common for DMARDs, prescriptions for each class of DMARD were treated separately. When there was a gap of 7 days or less between two prescriptions of the same drug, or of different drugs within the same drug class (excluding DMARDs), it was assumed that the drug was taken continuously and the gap was filled with the daily dose of the second prescription. Exposures to our main anti-rheumatic drugs included 1) MTX and 2) other DMARDs besides MTX.

Outcomes

The outcome of interest was the time from cohort entry to the first joint replacement surgery (for any joint), defined using the Canadian Classification of Health Intervention (CCI) and the Canadian Classification of Diagnostic, Therapeutic and Surgical Procedures (CCP) procedure codes for joint replacement (see Appendix 1).

Covariates

Variables that were available in the administrative databases and considered potential confounders for the association between the drugs of interest and joint replacement were selected a priori and adjusted for in all the multivariable models. These included sex, age at cohort entry (in years and with a squared-age term added to account for non-linear effects), place of residence (urban or rural, defined from postal codes) social assistance status at cohort entry, and ecological measures based on census data, on income, education level, and employment rate in the area. We used diagnostic codes from all outpatient physician and/or hospital visits during the 3 years before cohort entry to assess comorbidities (including osteoarthritis (OA), myocardial infarction, diabetes, osteoporosis, cerebrovascular disease, acute renal failure, chronic renal failure, coronary artery disease, chronic obstructive pulmonary disease (COPD), asthma, cancer and the Charlson index). To discriminate between high and low users of the health care system, we adjusted for a binary indicator of high users, defined as people with at least 20 physician visits in at least 1 year of the 3 years before baseline. To adjust for disease severity, we used a time-varying variable capturing the number of rheumatologist visits during the follow-up period (log-transformed). We also controlled for a time-varying indicator of the presence of extra-articular manifestations of the disease during the follow up; this included rheumatoid lung, Felty’s syndrome, rheumatoid carditis, eye involvement, dermatological complications (vasculitis, pyoderma gangrenosum), neuropathies and amyloidosis. Finally, we adjusted for time-dependent variables reflecting cumulative use of other drugs, anti-TNF inhibitors, other biologics, COXIBs, NSAIDs, and systemic steroids during the follow up, and binary indicators of anti-TNF, MTX, other DMARDs, COXIBs, NSAIDs, and systemic steroid use during the period of 1 year before cohort entry.

Statistical analyses

Descriptive statistics were used to characterize the study population. We used a Cox proportional hazards (PH) regression model with time-dependent variables measuring drug use for 1) MTX and 2) other DMARDs. As mentioned, the model also controlled for concomitant drug exposure (anti-TNF inhibitors, other biologics, COXIBs, NSAIDs, steroids). Our primary analyses considered the effects of early use of MTX and/or DMARDs, that is, in the first year of follow up only. Different time windows for cumulative drug use were considered in alternate models [10]: 1) during the second year of follow up, 2) throughout the entire follow up, and 3) throughout follow up but not including the year prior to the index time for each event. The cumulative use was obtained by summing the duration of all prescriptions for the relevant drugs, up to a given day over the relevant time period. We also tested for an interaction term between the cumulative effects of MTX and other DMARDs to account for possibly increased or decreased risks among users of both drug classes. Adjusted hazard ratios (HRs) with 95 % CIs were generated. The fit to the data of the different models was compared with the Akaike Information Criterion (AIC) [11].

In sensitivity analyses, we weighted the cumulative exposure to evaluate if weighting past exposure by recent use would improve the prediction of joint replacement surgery [12]. This method estimates from the data the relative weights of the timing of past exposure of each of MTX and other DMARDs on risk of outcome. We also performed sensitivity analysis excluding subjects with previous diagnosis of OA or excluding the cases of hip replacements, as these are less likely to be attributed to RA in the early years. Finally in additional sensitivity analysis, we used the propensity scores to adjust for potential differences associated with the exposure of primary interest, i.e., between the characteristics of patients who received different treatments during the first year after the cohort entry. To this end, we first used logistic regression to estimate separate propensity scores (PS) for the first-year treatment with 1) MTX and 2) other DMARDs. Both PS estimated the probability of receiving the respective treatment as a function of all time-fixed characteristics available at the entry into the cohort. In two separate analyses, we then included both PS in the multivariable Cox models that also included the cumulative first-year exposures to 1) MTX, 2) other DMARDs, and 3) their interaction, and all time-dependent potential confounders. The difference between the two PS-adjusted models was that the first also included all time-fixed variables, while the second excluded these time-fixed variables.