Heterogeneous Disease Trajectories Explain Variable Radiographic, Function and Quality of Life Outcomes in the Canadian Early Arthritis Cohort (CATCH)

Our objective was to identify distinct trajectories of disease activity state (DAS) and assess variation in radiographic progression, function and quality of life over the first two years of early rheumatoid arthritis (ERA). The CATCH (Canadian early ArThritis CoHort) is a prospective study recruiting ERA patients from academic and community rheumatology clinics in Canada. Sequential DAS28 scores were used to identify five mutually exclusive groups in the cohort (n = 1,586) using growth-based trajectory modeling. Distinguishing baseline sociodemographic and disease variables, treatment required, and differences in radiographic progression and quality of life measures over two years were assessed. The trajectory groups are characterized as: Group 1 (20%) initial high DAS improving rapidly to remission (REM); Group 2 (21%) initial moderate DAS improving rapidly to REM; Group 3 (30%) initial moderate DAS improving gradually to low DAS; Group 4 (19%) initial high DAS improving continuously to low DAS; and Group 5 (10%) initial high DAS improving gradually only to moderate DAS. Groups differed significantly in age, sex, race, education, employment, income and presence of comorbidities. Group 5 had persistent steroid requirements and the highest biologic therapy use. Group 2 had lower odds (OR 0.22, 95%CI 0.09 to 0.58) and Group 4 higher odds (OR 1.94, 95%CI 0.90 to 4.20) of radiographic progression compared to Group 1. Group 1 had the best improvement in physical function (Health Assessment Questionnaire 1.08 (SD 0.68) units), Physical Component Score (16.4 (SD 10.2) units), Mental Component Score (9.7 (SD 12.5) units) and fatigue (4.1 (SD 3.3) units). In conclusion, distinct disease activity state trajectories explain variable outcomes in ERA. Early prediction of disease course to tailor therapy and addressing social determinants of health could optimize outcomes.


Introduction
High variability has been noted within and between individuals with early rheumatoid arthritis (ERA), not only in their initial presentation, but also in the clinical progression of their disease [1]. This suggests the possibility that discrete disease trajectories exist at presentation, which may be explained by personal variables including sociodemographic status, disease characteristics, and/or health status, but also treatment received including initial treatment choice and dosing, escalation schemes, timing of therapy initiation, use of steroids, and tight control algorithms [2]. Understanding which baseline characteristics influence patient trajectories and important outcomes such as radiographic progression, physical function and quality of life measures in ERA could provide additional prognostic information to inform treatment strategy.
Recently, an ERA cohort study reported differing trajectories based on patient DAS28 scores over one year, showing patients could have a fast or slow response or a poor outcome, despite a treat-to-target strategy, determined with growth mixture modelling [3]. Analogous to studies defining 'responders' and 'non-responders', this analytic strategy allows further refinement of outcome ascertainment, as homogeneous groups of patients similar to each other at baseline and over time for a given outcome can be characterized longitudinally in a dataset, with the important predictors of the various treatment courses determined. Building on this important work, we applied a group-based trajectory modelling strategy to data from the CATCH (Canadian Early Arthritis Cohort) Study. Our objectives were to i) evaluate whether patients could be clustered over time based on their disease activity; ii) identify whether group-based modifiable and non-modifiable sociodemographic and disease related factors differentiate group membership at disease onset; and iii) determine if such clinical distinctions are concordant with differences in radiographic progression, physical function and quality of life measures during the first two years of disease.

Methods Subjects
Data from the CATCH study (January 2007 to May 2014) were used for this study. CATCH is a prospective observational cohort study of adults with ERA from 21 academic and community-based clinics across Canada that has been previously described [4]. Briefly, patients were included with 12 months of symptoms, 2 swollen joints or 1 swollen metacarpophalangeal (MCP) or proximal interphalangeal (PIP) joint, plus at least one RA related disease factor. The primary analysis included all CATCH subjects observed in at least 2 study visits and meeting either 2010 [5] or 1987 [6] RA Classification Criteria at the baseline visit.

Data Elements
Standardized study visits include a baseline assessment at the time of enrolment, quarterly reviews in the first year, biannually in year 2, and annually thereafter. As this is an observational study performed in clinical settings, not all subjects have anti-cyclic citrullinated peptide (anti-CCP) tested (as it is not reimbursed by all local health authorities). Physician measured joint counts, assessment of disease activity, medication received in a treat-to-target strategy, C-Reactive Protein (CRP) and Erythrocyte Sedimentation Rate (ESR), are collected at all visits. Patient self-reported function (Health Assessment Questionnaire, HAQ [7]) and visual analogue scales for global evaluation of well-being and pain in the past week (rated from 0-10) are collected every 3 months in the first year, and every 6 months in the second year. Subjects are classified for their disease activity state at each visit based on the Disease Activity Score using 28 joints (DAS28-ESR) with established cutpoints: remission (REM, DAS28<2.6), low disease activity state (LDAS, DAS28 2.6-3.2), moderate disease activity state (MDAS, DAS28 3.2-5.1) or high disease activity state (HDAS, DAS28 >5.1) [8].

Radiographic Measures
Sequential-order van der Heijde modification of Sharp scores (vdHSS) were assigned by one trained reader (intra-class correlation coefficient for intra-observer variability 0.74) in a subset of patients (n = 488, 31%) who had annual serial radiographs performed and at least 2 study visits [9]. Significant radiographic progression was defined as an annual increase above the smallest detectible change (SDC), which was 3.5 units/year in this cohort [10]. Radiographs performed 3-6 months prior to study entry were not repeated but were used as baseline films.

Quality of Life Measures
The Veteran's RAND 12 Item Health Survey (VR-12) [11] are completed annually. From twelve items in the questionnaire, a Physical Health Summary Measure (PCS-physical component score) and a Mental Health Summary Measure (MCS-mental component score) are summarized. The minimum clinically important difference (MCID) for these scores is 2.5 [12]. Fatigue was also considered as a quality of life measure in this study, captured by a 0-10 numeric rating scale, with an MCID of 2 [13].

Statistical Methods
The group-based trajectory modelling (GBTM) strategy proposed by Nagin [14] was applied to the data to test for the existence of, and optimal number of, distinct groups for analysis. In the first step, the Bayesian information criterion (BIC) statistic [15] is calculated to provide guidance on the number of distinct trajectory groups that can be defined from the data. In the second stage of the process, a cubic polynomial is used to specify the shape of each trajectory. Clinical acumen is applied to summarize the distinctive features of the data in as parsimonious a fashion as possible. Subjects are then assigned to the group where their posterior membership probability was highest (i.e. to which they most likely belong). The analysis was carried out in SAS (version 9.2), using PROC TRAJ, which uses a general quasi-Newton procedure to estimate parameters that maximize the likelihood function.
Descriptive statistics for continuous variables are presented as means (standard deviation, SD) and categorical variables as frequencies with percentages. Sociodemographic, disease and treatment variables distinguishing the groups were identified by chi-square tests for categorical variables and by Kruskal-Wallis one-way ANOVA for continuous variables. We adjusted for multiple comparisons using the Bonferroni correction.
Due to the limitations of radiographic data analysis using last-observation-carried-forward or linear imputation we used the multiple imputation method described by Baron [16].

Ethics
Individual consent was obtained, following ethics approval through each site's Research Ethics Board.

Subject Cohort
As of May 2014, a total of 2,524 subjects were enrolled in CATCH. Of these, 1,586 subjects met RA classification criteria, consented to ongoing study participation, and had a minimum of 2 study visits over 24 months (Fig 1).

Trajectory Groups
Five unique trajectory groups characterized by sequential DAS28 scores were predicted (Fig 2). Subjects were assigned to the predicted groups as follows: Group 1 (20%) initial high disease activity state (HDAS) improving rapidly to remission (REM); Group 2 (21%) initial moderate disease activity state (MDAS) improving rapidly to REM; Group 3 (30%) initial MDAS improving gradually to low disease activity state (LDAS); Group 4 (19%) initial HDAS   Table 2 summarizes the sociodemographic, health status and disease characteristics of subjects in the five identified groups. Of the sociodemographic variables, age, race, education, income, and employment status were significantly different between groups (all between-groups Table 1. Inception characteristics of subjects in CATCH and of patients with available serial radiographs. Results are reported as a Mean (SD) unless otherwise noted. RF, rheumatoid factor; Anti-CCP, anti-cyclic citrullinated peptide; ESR, erythrocyte sedimentation rate; DAS28, disease activity score using a 28-joint count; HAQ Health Assessment Questionnaire; DMARD, disease modifying anti-rheumatic drug adjusted p-values <0.001). A higher proportion of subjects in Groups 4 or 5 were >50 years of age (n = 224, 74% and n = 108, 68% respectively) relative to Groups 1, 2 or 3, and 28% (n = 44) of Group 5 were non-Caucasian. Subjects in Groups 4 or 5 had less frequently achieved higher than high school education (n = 119, 39% and n = 65, 41%, compared to >52% in Groups 1, 2 or 3). In Groups 4 or 5 <44% of subjects were employed compared to >55% in the other groups. With regards to health status, subjects in Groups 4 or 5 had a higher mean number of comorbidities (2.6 (SD 2.1) each, p<0.001). Despite Group 5 having the longest symptom duration at diagnosis (mean 203 days (SD 111)), no differences in disease characteristics at baseline for serology, baseline radiographic scores or the presence of erosions were seen.

Treatment by Group
The proportion of patients in each group receiving methotrexate, methotrexate in combination with another DMARD, biologics and corticosteroids at each study visit appears in Table 2 and Fig 3. The majority of subjects were already receiving therapy at their baseline visit, with only 14% of Group 5 patients being treatment-naïve, and with >70% of subjects on methotrexatebased regimes (Group 1 82%, Group 2 73%, Group 3 67%, Group 4 85%, Group 5 71%) in a mean dose exceeding 19 mg weekly. Differences emerge during the course of observation. In particular, Group 5 has the lowest proportion of patients remaining on methotrexate and/or combination therapy, but with the highest proportion (nearly 50%) moving on to biologic therapy, yet only achieving MDAS by two years. Group 4, which has the highest proportion of subjects on methotrexate and/or combination therapy, has approximately 30% on biologic therapy by two years, and achieves LDAS by two years. Both these groups also have the highest proportion of subjects on corticosteroids throughout the observation period. In contrast, Group 2 begins in MDAS and achieves REM, has declining methotrexate use over two years of followup and <7% prevalence of biologics throughout.

Radiographic Outcomes
There were no significant differences between groups at baseline in the proportion with erosions or vdHSS (Table 2). However, significant differences in radiographic damage progression are evident by the first year of follow-up (between-group chi-squared p<0.001; Fig 4). Radiographic progression was most frequent in subjects assigned to Groups 4 and 5 at 34% (n = 21/ 61) and 33% (n = 10/30) respectively, compared to 5% (n = 4/11) in Group 2, 19% (n = 16/85) in Group 3 and 21% (n = 14/68) in Group 1. Using Group 1 as the reference in unadjusted

Physical Function
Baseline values by trajectory group are reported in Table 2. Groups 1, 4 and 5 have the worst baseline HAQ score. However, Group 1 has the greatest improvement in HAQ score (mean change of 1.08 (SD 0.68), with Groups 2, 3 and 4 also improving significantly (mean change 0.34 (SD 0.51), 0.34 (SD 0.61) and 0.71 (SD 0.78)) respectively (Table 3). Group 5 did not improve during follow-up (mean change 0.13, SD 0.66), and only 42% improved more than the minimal clinically important difference.

Quality of Life Measures
Baseline values by trajectory group are reported in Table 2

Discussion
Our study uses group-based trajectory modelling to demonstrate distinct heterogeneity in the course of disease in ERA, which is ultimately reflected in variable degrees of radiographic damage progression and improvements in physical function and quality of life measures. We detected five discrete disease trajectories, categorized by the disease activity state at inception to the cohort and then observed over two years. These trajectories are characterized by variations in sociodemographic and health status factors at baseline and less so disease characteristics, and which emerge despite the initial aggressive treatment received. Patients in the group experiencing the largest improvement in function and quality of life measures (Group 1, HDAS to REM) and the group with the best prognosis (Group 2, MDAS to REM, with the least radiographic progression and treatment needs) are more often <50 years old, with less medical comorbidity. Patients with the worst prognosis (Group 4, beginning in HDAS and achieving LDAS with high odds for radiographic progression; Group 5 beginning in HDAS and only achieving MDAS) are more likely to have lower levels of education, are less frequently employed and have lower total household incomes. Group 5 is additionally characterized by a high proportion of non-Caucasian ethnicities. These findings speak to the importance of addressing social determinants of health to improve ERA outcomes, as even aggressive treatment could not prevent radiographic progression nor provide robust improvements in quality of life measures.
Although the relationship of treatment with these trajectories is complicated, the majority of patients were treated aggressively with methotrexate-based regimens in doses exceeding 19 mg weekly. The group with the worst disease activity trajectory (Group 5) had less frequent use of methotrexate (albeit in similar doses if used compared to the other groups) and other Baseline Visit Biologic 11 (3%) 6 (2%) 4 (1%) 8 (3%) 6 (4%) 1.0 * Bonferroni Method. Table does not display some variables included in the Bonferroni correction, including site size the patient was recruited from, the proportion meeting ACR 1987 or 2010 criteria, all methotrexate exposed (yes/no), the proportion on non-methotrexate DMARDs, and oral steroids. No significant differences were seen in these variables between groups. ** On proportion of subjects with radiographs available: Group 1 n = 102, Group 2 n = 114, Group 3 n = 131, Group 4 n = 91, Group 5 n = 49 DMARDs and more corticosteroid and biologic use during follow-up, without realizing the treatment goal of LDAS or REM. This group may represent patients who have had contraindications to usual first-line therapies, or differences in tolerating, adhering to, or accepting use of DMARD medication at higher doses, possibly influenced by cultural differences or health literacy [17]. Innate heterogeneous biologic differences are also likely to contribute to disease course, reflected in the very different responses between groups 1 and 4, in which methotrexate therapy is used similarly. To illustrate, Groups 1 and 4 both started in high disease activity and received similar treatment, yet only patients in Group 1 were rapid responders to methotrexate (DAS28 reduction of 2.96 within 3 months). Given these observed differences in treatment  response, there is an opportunity to identify unique molecular signatures in these discrete groups that could inform patient-specific initial therapeutic choices.
Our study complements the study by Siemons et al [3] in which trajectories of 'Fast Response', 'Slow Response' and 'Poor Outcome' groups were identified over 12 months despite a treat to target strategy, although notably with lower doses of methotrexate used compared to the population we describe. That study compared baseline characteristics for disease activity (DAS28, tender and swollen joint counts, inflammatory markers), disease characteristics (rheumatoid factor positive), limited sociodemographic factors (body mass index, age, sex), and selected patient-reported and well-being measures (patient global, pain, SF-36 physical, SF-36 mental) and only male sex was predictive of fast response. They concluded that future studies should attempt to identify more specific risk factors for poor outcome, to enable early identification of patients in need of alternative therapeutic approaches. Norton has examined heterogeneity in ERA by defining functional trajectories using scores from the Health Assessment Questionnaire, and examining the relationship with mortality [18] and disability progression [19]. Norton also described sociodemographic differences between trajectory groups, notably that females, higher medical comorbidity scores, low education, low social class and low employment levels predominate in the group with the highest level of HAQ disability. This group also experienced higher mortality after adjustment for important confounders. Although our study did not use a mortality endpoint, we identified that significant radiographic progression, which is tightly associated with disability [20], was observed in one-third of patients in the groups with HDAS at treatment initiation and a flat response curve, in just one year of observation. Thus, we might expect to see findings similar to that of Norton in our cohort in the long-term. Taken together, these data demonstrating disease heterogeneity underscore the need for early identification of patient-specific factors to predict a patient's destined trajectory.
The GBTM strategy had its origins in criminology theory, but has gained a foothold in psychology, mental health, epidemiology, chronic disease and medication adherence studies [21,22]. The advantages of this approach in longitudinal data analysis are the ability to statistically derive whether true heterogeneity exists in a dataset, while accounting for possible chance differences across individuals [23]. Patterns of outcomes are identified, and the selection of the number of trajectories are informed statistically, thus providing more objective analysis rather than just applying a simple dichotomous outcome selected by the researcher. As an example of application from the modelling to the clinical setting, a recent Canadian retrospective study of swollen joint count trajectories in juvenile inflammatory arthritis was performed [24]. Berard et al. identified five latent classes, which were clinically and statistically distinct from the International League of Associations for Rheumatology (ILAR) categorizations. Thus this type of analytic approach will allow us to redefine outcome trajectories that reflect heterogeneity within diseases. An additional advantage of this modelling strategy is that it allows analysis of baseline predictors for heterogeneous outcomes, thus providing a structure to explore personalized treatment selection.
Limitations of the study were addressed where possible. Longitudinal studies are susceptible to bias from patients lost to follow-up, who may differ systematically from patients remaining in the study. To assess the potential impact of incomplete follow-up, we performed the analysis on 343 patients with complete 2-year follow-up data. The model continued to identify 5 trajectories with proportionately similar numbers of patients in each group, providing further credence to the observed heterogeneity in this population. Similarly, when we performed the analysis on treatment naïve subjects or those minimally exposed to methotrexate our results did not change (S1 Fig). Moreover, this analytical method allows inclusion of all data when missing data are at random. Employment status at entry to the cohort may already reflect consequences of the disease rather than a sociodemographic factor that influences prognosis. Additionally, unmeasured confounders associated with disease trajectories such as therapeutic intolerances and adherence may not be evenly distributed across groups.

Conclusion
In the era of cost constraints, means to identify patients destined to have a severe disease course using novel statistical methods, additional variables, or biomarkers are desired, which will allow provision of the right therapy for the right patient resulting in maximized function, gainful employment and less destruction from this disease. These data confirm that heterogeneous disease trajectories exist in ERA, and that emerge despite initial aggressive therapy. This may reflect limitations in the use of available therapies due to patient comorbidities, health literacy, or individual response to standard therapies.
Supporting Information S1 Fig. Predicted Group Trajectories in Treatment Naïve or Minimally Exposed Early Rheumatoid Arthritis based on DAS28 with 95% CI (n = 343). Five predicted group trajectories (solid or dashed lines) and 95% confidence interval limits (shaded) are depicted from the group-based trajectory modelling. Percentages reflect the predicted proportion of subjects in each group, which differs marginally from the actual group characterization in the dataset. (TIFF)