Study design and setting
We conducted a population-based inception RA cohort study in Ontario. Ontario residents receive all medically necessary health services free at the point of care under a single payer healthcare system. The Ontario Drug Benefit Plan covers the cost of prescription medications for Ontarians aged 65 years and older, subject to a small copayment. Health services are recorded in administrative databases which enable comprehensive evaluations of care.
The Ontario Health Insurance Plan (OHIP) Claims History Database was used to identify physician services and associated diagnoses. Encoded physician identifiers in the claims were linked to the ICES Physician Database to identify physician specialty and practice location. Patient demographic information, place of residence and vital status were ascertained from the OHIP Registered Persons Database. Pharmacy claims from the Ontario Drug Benefit Program were available for individuals \(\ge\)66 years (allowing up to one year for coverage registration). Additional information on patient health care utilization and comorbidities were obtained from hospital discharge records in the Discharge Abstract Database and emergency department visits recorded in the National Ambulatory Care Recording System database. A complete list of databases can be found in Additional file 1: Table S1. These datasets were linked using unique encoded identifiers and analyzed at ICES (www.ices.on.ca). ICES is a prescribed entity under Sect. 45 of Ontario’s Personal Health Information Protection Act. The use of data in this study was approved by a privacy impact assessment at ICES and authorized under Sect. 45 of Ontario’s Personal Health Information Protection Act, which does not require review by a Research Ethics Board.
All incident RA patients between January 1, 2002 and December 31, 2019 were identified from the Ontario Rheumatoid Arthritis Database . Patients are included if they have at least one hospitalization, or at least two physician claims for RA over 2 years (ICD9 714.0 or ICD10 M05-M06) with at least one from a rheumatologist, internist, or orthopedic surgeon. The case definition has been validated to have a sensitivity of 78%, specificity of 100% and positive predictive value of 78% .
Individuals were excluded if they had missing demographic information, a diagnosis date before 18 years of age, were non-Ontario residents at the date of their first RA code (i.e. individuals from a another province and received care in Ontario given reciprocal agreements for care coverage between provinces), or did not have health insurance eligibility in the 5-year period prior to their first RA code (to ensure only incident RA cases were included). Individuals who did not have at least 1 year of follow-up after their first RA code were also excluded. (Flow diagram of cohort selection shown in Additional file 1: Fig. S1).
The following characteristics were determined at cohort entry: age, sex, neighborhood income quintiles derived from census data, and urban versus rural location of residence (the latter defined based on postal code and a community size of < 10,000 residents)  Fourteen Local Health Integration Networks (LHINs), as defined by the Ontario Ministry of Health, were used as geographic health service regions. Linear distance in kilometers between patients and the nearest rheumatologist was determined (calculated from the center of the patient’s postal code), and individuals residing 100 km from the nearest rheumatologist were described as living at a “remote” distance. Comorbidities in the three years prior to RA diagnosis were assessed using diagnosis codes from physician claims and hospital discharge records and applying validated health administrative algorithms [16,17,18,19,20,21] when available. In addition, the Johns Hopkins ACG® System version 11 was used to assign each patient to up to 32 Aggregated Diagnosis Groups® (ADGs) using diagnosis codes found in OHIP physician claims and hospital discharge records, using a three-year look-back period. A patient frailty indicator was also obtained using the ACG System.
System performance measure (PM) adherence
Two PMs  were operationalized as previously described [4, 5], and assessed between 2002 and 2020 for patients entering the cohort in each year. The first PM reports on access to rheumatologist care and is calculated as the percentage of incident RA patients seen by a rheumatologist within 365 days of their first RA diagnosis code by any physician.
The second PM is reported on individuals aged 66 and older who saw a rheumatologist within the first year of diagnosis (measured from the first RA code) and lived at least an additional 30 days. This PM reports on the percentage of individuals dispensed a DMARD following diagnosis confirmation at the rheumatologist visit [3, 22]. Given potential clinical challenges in starting treatment within the recommended benchmark of 14-days [3, 22] (e.g., due to delays in obtaining baseline lab investigations and/or receiving appropriate medication counselling, time for appropriate patient decision making, or patient delays in filling a prescription), we applied a 30-day benchmark. DMARDs in these analyses included conventional synthetic DMARDs, targeted synthetic DMARDs, biologic agents as well as other immunosuppressive therapies used to treat complications of RA (see Additional file 1: Table S2 for included DMARDs).
Regional rheumatology supply
Based on recommendations from the Canadian Rheumatology Association , regional rheumatologist supply was classified as optimal (at least one rheumatologist per 75,000 residents in the region of the patients’ residence) or suboptimal (less than one rheumatologist per 75,000) based on the local health region (LHIN).
We used descriptive statistics (mean or median depending on normality of the data, and frequencies) to characterize patients at cohort entry. Outcomes were assessed annually as the proportion of incident patients diagnosed each year who met each of the PMs, stratified by rheumatologist supply. Logistic regression was used to assess whether there were improvements in trends over time, and whether the improvements were associated with rheumatologist supply. Multivariable models additionally adjusted for factors which may affect access to care, including age at diagnosis, sex, income quintile, rural residence, log(distance to nearest rheumatologist), and comorbidities. Time from diagnosis to first rheumatologist visit was included in the model predicting DMARD use, and a generalized estimating equation was used to account for the clustering of patients within rheumatologist practices. Because the trends over time were not linear, the shape of the trend was characterized using fractional polynomials . In the adjusted analyses, fractional polynomials were also used to characterize the relationship between patient age and both PMs. All analyses were performed using SAS version 9.3 (SAS Institute Inc., Cary, NC). Two-tailed p-values < 0.05 were considered to be statistically significant.