Latent Class Trajectory Modeling of 2‐Component Disease Activity Score in 28 Joints Identifies Multiple Rheumatoid Arthritis Phenotypes of Response to Biologic Disease‐Modifying Antirheumatic Drugs

To determine whether using a reweighted disease activity score that better reflects joint synovitis, i.e., the 2‐component Disease Activity Score in 28 joints (DAS28) (based on swollen joint count and C‐reactive protein level), produces more clinically relevant treatment outcome trajectories compared to the standard 4‐component DAS28.


INTRODUCTION
Complexity and variability of patient trajectories in response to treatment pose significant challenges in many medical fields (1)(2)(3)(4). A current focus of medical research is facilitating a shift away from a "one-size-fits-all" approach to precision medicine approaches aimed at identifying narrower patient subgroups that would benefit from tailored interventions. However, the high variability in response to therapeutic agents over time makes prediction of response challenging (5)(6)(7).
It has previously been demonstrated that the 4-component DAS28 correlates poorly with synovitis (22). Given that the drugs used to treat RA aim to reduce synovial inflammation and prevent joint damage, it is unclear how useful the 4-component DAS28 trajectories will be as outcome measures for precision medicine studies, which aim to identify biomarkers that are predictive of therapeutic response. To address this, a revised composite disease activity score, including only the swollen joint count (SJC) and C-reactive protein (CRP) levels (2-component DAS28), has been developed and shown to better correlate with ultrasounddetected synovial inflammation compared to the original 4-component activity score (22).
We undertook this study to apply the LCMM approach in a prospective, longitudinal cohort of patients beginning biologic DMARD treatment. We used the original 4-component DAS28 and 2-component DAS28 in order to assess whether the identified trajectory groups have clinical relevance, by exploring their association with comorbidities, treatment adherence, and drug immunogenicity, all of which have previously been reported to be associated with biologic DMARD treatment response (23,24).

PATIENTS AND METHODS
Study patients. We used longitudinal data (November 2008-January 2018) from patients recruited into the Biologics in Rheumatoid Arthritis Genetics and Genomics Study Syndicate (BRAGGSS), an observational multicenter study in the UK established to study predictors of treatment response to biologic DMARDs (25). Consenting participants were eligible if they had a consultant diagnosis of RA, were >18 years old, and were about to initiate treatment with biologic DMARDs, namely tumor necrosis factor (TNF), B cell inhibitor (anti-CD20), interleukin-6 receptor (IL-6R), or T cell inhibitor (CTLA-4). Information regarding the combination with conventional DMARDs, e.g., methotrexate (MTX), was recorded. During the observational study, biologic DMARD therapy was suspended, changed, or restarted in a minority of patients. These patients were included in the analyses, and information about previous treatment with biologic DMARDs was recorded. Following the pretreatment baseline sample and data collection, patients were observed for 12 months, with follow-up visits at 3, 6, and 12 months. Participating patients provided written informed consent, and the BRAGGSS study was approved by the local ethics committee (COREC 04/Q1403/37).
The primary outcome measure to model response trajectories was the disease activity score based on the 2-component DAS28, computed as sqrt (SJC) + (0.6*ln[CRP+1]). Four-component DAS28 score and response criteria were calculated according to standard methodology (26). Variables recorded at baseline included age, sex, time since diagnosis, body mass index (BMI), smoking status (if the patient had ever smoked), functional severity measured by the Health Assessment Questionnaire (HAQ) (27), and consultant-reported comorbidities. DAS28 components (tender joint count [TJC], SJC, patient global assessment of well-being measured on a 100-mm visual analog scale [VAS], serum CRP level) were measured at baseline (i.e., before commencing any biologic DMARD treatment) and at each follow-up visit. For patients treated with TNF inhibitors, antidrug antibodies and serum drug levels were measured in serum samples collected at follow-up visits.
We included patients with ≥1 score from the 2-component DAS28 after baseline, including all measurements until the last registered follow-up. Unlike previous studies (28,29), we did not exclude patients on the basis of the number of scores, in order to avoid possible selection biases, and we also included patients who had received biologic treatment in the past. Data on patients who switched biologic DMARDs were right-censored when the change occurred. Missing data arose from noncompletion of responses to questions on demographic characteristics and missing HAQ scores. Due to the nature of the variables, we assumed they were missing at random. As in previous studies (30,31), we compared different imputation strategies (see Supplementary Methods, on the Arthritis & Rheumatology website at http://onlin elibr ary.wiley.com/doi/10.1002/art.41379/ abstract) and implemented the method with the best performance. Data were imputed with a random forest approach, with 100 trees and a maximum of 100 iterations.

LCMM.
We applied LCMM to identify patient subgroups with distinct responses to biologic DMARDs over time. Previous stratification studies in RA using LCMM focused on the 4-component DAS28 as the measure of disease activity and response to biologic treatment (18)(19)(20)(21); those analyses identified 3 groups of responders: rapid, gradual, and inadequate. In the present study, we specified 2 linear mixed-effects models: the 4-component DAS28 used in previous studies and the 2-component DAS28 used here as the dependent variable.
Mixed effects were used to account for the likely correlation of repeated measurements and included a random intercept for each individual. We fitted the model through the "lcmm" function of the R package lcmm (32). We followed the framework described by Lennon et al (33) that includes scoping model definition, refinement of the number of classes and model structure, model adequacy assessment and clinical plausibility, graphic presentations, and sensitivity analysis (Supplementary Methods, http://onlin e libr ary.wiley.com/doi/10.1002/art.41379/ abstract).
Models were adjusted for age, time since diagnosis, sex, BMI, HAQ score, biologic DMARD therapy, concomitant treatment with MTX, and previous treatment with biologic drugs. We developed scoping models for 1-10 classes to determine the optimal number of classes based on the Bayesian information criterion (BIC) and further refined the models comparing linear, quadratic, and cubic specifications of time (Supplementary Methods, http://onlin e libr ary.wiley.com/doi/10.1002/art.41379/ abstract).
To evaluate whether results were driven by differential relative effectiveness of therapies, in sensitivity analyses we repeated the analysis in patients treated with TNF inhibitors only and in those treated with any biologic DMARD in combination with MTX, comparing the discovered latent classes in terms of baseline characteristics and response variables.

Statistical analysis and comparison of latent classes.
Discovered latent classes were described and compared in terms of clinical characteristics, comorbidities, adherence to treatment, antidrug antibody positivity, and non-trough blood drug levels.
Treatment response variation in individual patients. To describe variation in treatment responses in each patient over time, we computed slopes as the difference between consecutive measurements divided by the time between the measurements (Supplementary Methods, http://onlin elibr ary.wiley.com/ doi/10.1002/art.41379/ abstract). Higher slope values indicate deterioration of a patient's condition, while negative values indicate improvements. Time to nonresponse was defined as a switch from a negative to a positive value of the slope (i.e., when the slopes start to increase). Intercept values and slopes were compared between classes using analysis of variance (ANOVA). Time to nonresponse was compared among nonresponders using Kaplan-Meier curves.
Demographic and clinical characteristics. Clinical characteristics were compared among latent classes. Categorical characteristics (sex, pretreatment smoking status, MTX cotherapy, and switch of biologic DMARDs) were compared using chi-square tests followed by Bonferroni-Holm correction, while continuous variables (age, disease duration, HAQ score, BMI, and DAS28 score) were compared using ANOVA with Tukey's post hoc test.
Association of latent classes with comorbidities and therapy (TNF inhibitor, B cell inhibitor, IL-6R inhibitor, and T cell inhibitor) was explored using logistic regression, where the latent class assignment was the regression predictor, and least squares means was used for multiple comparisons with Tukey's adjustment (from the R lsmeans package [34]). The effect of respiratory comorbidities on latent classes was assessed by stratified analyses in smokers (current or past) and nonsmokers (never).
Treatment adherence, drug levels, and antidrug antibody positivity. We examined the associations of latent class subgroup with treatment adherence, drug levels, and antidrug antibody positivity. Drug level and antidrug antibody values have previously been shown to be associated with response to biologic DMARDs (24). Self-reported adherence was recorded at 3 and 6 months following the start of therapy (35); adherent patients were defined as those who took the dose on the day agreed upon by their health care team or ≤1 day before or after, while nonadherent patients were defined as those who self-reported either missing the dose completely or taking it >1 day before or after the agreed day.
The correlation between drug levels and latent classes was investigated in terms of patient-level relative change at the next follow-up visit. We created matched samples with propensity scoring, including BMI (36), age, and sex as potential confounders without interactions and drug level changes (increasing or decreasing levels) as a binary outcome measure. Chi-square tests were then applied to assess differences in drug level changes among latent classes.
To evaluate antidrug antibody development, we used a threshold of 12 arbitrary/ml as a positive event (24). Heterogeneities of latent classes across the total cohort and the subset with antidrug antibody information were determined by Cochran's Q statistic. Antidrug antibody development events were compared in latent classes using Kaplan-Meier curves and Cox proportional hazards regression models, with adjustment for age, sex, BMI, and concomitant treatment with MTX.
Analyses were computed using R, version 3.2.3. Data are presented as the main effect with 95% confidence intervals (95% CIs), and the 5% significance level was used for main inferences.

RESULTS
Patient characteristics. Our data set included 2,991 patients with 7,567 DAS28 measurements after baseline, excluding the baseline measure. In the 12-month observation period from baseline, patients had 1-3 follow-up visits: 1,647 patients (55.07%) had 1 follow-up visit, 661 (22.10%) had 2 follow-up visits, and 683 (22.84%) had 3 follow-up visits. The first follow-up visit took place a mean ± SD of 124 ± 51 days after baseline, the second 217 ± 51 days after baseline, and the third 365 ± 50 days after baseline. Table 1   Although inclusion criteria of previous studies were different (19,21) and the mean 4-component DAS28 score at baseline was higher (20), response patterns in the present study were qualitatively comparable. We identified the following groups: a rapid responder group (n = 2,004 [67.00%]) with quick improvement in the first observation period followed by stabilization or slight increases in disease activity, a gradual responder group (n = 919 [30.73%]) with slower but consistent decrease in disease activity, and a poor responder group (n = 68 [2.27%]) (Supplementary Figure 19 and Supplementary We found no significant correlation between drug level variation and response trajectories using the 4-component DAS28 (P = 0.40). Data on antidrug antibodies were available for 475 patients receiving adalimumab or certolizumab, whose disease activity classes distributed in a similar manner to that of the whole cohort (Supplementary Table 12  Good responders demonstrated fast improvement in disease activity in the first 3 months, followed by maintenance of good response (Figure 1). The first group of good responders (n = 395 [13.2%]) and the second group (n = 1,840 [61.5%]) differed  Gradual responders demonstrated a continuous but gradual reduction in disease activity ( Figure 1). These 2 groups also had significantly different baseline scores for the 2-component DAS28,  and the low activity disease group (B) with 95% confidence intervals are shown. Score curves were determined using the loess method. In early secondary nonresponders (red), both swollen joint count (SJC) and C-reactive protein (CRP) levels followed the same pattern detected with the 2-component DAS28. In late secondary nonresponders (orange), the nonresponse pattern was mainly driven by the SJC component, while the CRP pattern was mainly flat. The trajectories of the low activity disease group (yellow) showed high tender joint counts and visual analog scale scores, and the disease activity was mainly driven by these 2 components.    Table 1 for other definitions). † Adjusted for multiple comparisons between groups.
Patients in the low disease activity group (n = 88 [2.9%]) demonstrated the lowest DAS28 scores at baseline (mean ± SD 2.24 ± 1.09), which increased modestly during follow-up (Figure 1). Initial levels of disease activity measured using the 4-component DAS28 were driven by TJC and VAS score, while SJC and CRP level showed a consistent increase in the follow-up period ( Figure 2). Latent classes were described in terms of slopes (Supplementary Figure 10 and Supplementary Table 6, http://onlin elibr ary. wiley.com/doi/10.1002/art.41379/ abstract). When we compared the 2 subgroups of secondary responders in terms of 2-component DAS28 score slopes and time to nonresponse, we found that early secondary nonresponders lost response significantly earlier (mean ± SD 212 ± 28 days after the baseline visit) than late secondary nonresponders (mean ± SD 265 ± 73 days) (P = 0.005) (Supplementary Figure 11, http://onlin elibr ary.wiley. com/doi/10.1002/art.41379/ abstract).
Patients in the low disease activity group showed modest improvements in TJC and VAS score during the observation period but demonstrated deterioration in SJC and CRP level ( Figure 2). The majority of patients in the low disease activity group were classified as nonresponders (79% by 6 months) using the 4-component DAS28-derived EULAR response criter ia (Supplementary Tables 8 and 9, http://onlin elibr ary.wiley.com/ doi/10.1002/art.41379/ abstract). Importantly, the low disease activity group included the highest proportion of patients who had received other biologic DMARD treatments in the past (36.36%).
Sensitivity analyses showed that in TNF inhibitor-treated patients (n = 2,151), the best model identified 6 classes, with an absence of the gradual responder class with low initial 2-component DAS28 score (Supplementary Figure 7, Table 2. Further details can be found in the Baseline Characteristics section of the Supplementary Methods (http://onlin elibr ary.wiley.com/doi/10.1002/art.41379/ abstract). Figure 3 illustrates comorbidity prevalence in LCMM classes. Early secondary nonresponders had a higher prevalence of chronic bronchitis/emphysema compared to good responders (good responders group 1 OR 5.38, P = 0.005; group 2 OR 5.85; P = 0.005) and gradual responders (gradual responders group 2 OR 5.54, P = 0.003). The stratified analysis of correlation between being an early secondary nonresponder and having chronic bronchitis/emphysema showed significant effects in ever-smokers (OR 2.49 [95%CI 1.14-5.42]) and never-smokers (OR 4.87 [95% CI 1.06-22.32]). Low disease activity patients had a higher prevalence of asthma compared to good responders (good responder group 1, OR 2.53 P = 0.029).
Two-component DAS28 response trajectories are associated with adherence, drug levels, and antidrug antibodies. Self-reported adherence measures, collected at follow-up visits, were available for 1,528 patients (51%). Twocomponent DAS28 slopes decreased more quickly in adherent  patients (mean ± SD −0.00128 ± 0.0581) than in nonadherent patients (mean ± SD −0.0120 ± 0.103; P = 0.002). Adherence patterns indicated that the low disease activity group, which were identified using only the 2-component DAS28, were the least adherent group (Supplementary Figure 16 The highest proportion of antidrug antibody-positive patients was observed in the early secondary nonresponders, but high antidrug antibody titers were also observed in gradual responders ( Figure 4A). Cox regression models, which included the effect of MTX co-reatment, revealed different risks for antidrug antibodies (β = 6.06 [95% CI 2.57-14.27], P = 3.71 × 10 −5 for late secondary nonresponders compared to good responders group 1; β = 1.79 [95% CI 1.00-3.22], P = 0.018 for gradual responders group 2 compared to good responders group 1) ( Figure 4B).

DISCUSSION
In the present study, we determined that using the 2-component DAS28 score as an outcome measure revealed more subgroups of patients compared to the 4-component DAS28, and the trajectories identified were clinically meaningful. While the 4-component DAS28 remains an essential tool to assess RA progression in clinical practice, the 3 trajectories identified, which can be also interpreted in terms of EULAR response criteria, may not provide sufficient granularity to help elucidate important biologic mechanisms involved in inflammation.
Strengths of this study include the fact that it is the first and largest study, based on real-world patients, to explore patientcentered factors for trajectories of drug response and that the response trajectories based on the 2-component DAS28 were identified using a structured methodologic framework (33). In addition, all classes of licensed biologic DMARD treatment regimens were considered and included in longitudinal data models, and the clinical relevance of the identified subgroups was compared using outcomes relevant to patients.
We found that the 7 subgroups demonstrated differing degrees of treatment response, mean patterns of change in disease activity, and clinical characteristics. For example, drug levels were better tracked using the 7 trajectories from the 2-component DAS28 than the 3 trajectories from the 4-component DAS28. Furthermore, antidrug antibody positivity was enriched in 2 of the nonresponder groups using the 2-component disease activity score, indicating that these trajectories are biologically relevant and more informative than those identified using the 4-component score.
Another potentially clinically useful subgroup identified using the 2-component DAS28 was the low disease activity group, in which conventional DAS28 scores were driven by the TJC and VAS. This group was less likely to respond to treatment (79% classified as EULAR nonresponders by 6 months) and less likely to adhere to treatment. It is also interesting to note that previous nonadherence, together with an increased resistance to therapy, might be related to a higher proportion of patients receiving previous biologic DMARD treatment (nearly 37% in the low disease activity group). Use of biologic drugs in this subgroup was unlikely to be clinically efficacious or cost-effective, and biomarkers to identify such patients would be useful.
While the 4-component DAS28 is an important holistic assessment that better captures depression, the 2-component DAS28 appears better able to capture biologic comorbidities, especially respiratory phenotypes. As previously demonstrated (23), antidrug antibody concentrations and treatment adherence remained the most important predictors of drug levels over time. The 2-component trajectories capture these characteristics and confirm previous findings that drug levels, treatment adherence, and antidrug antibody positivity correlate with the speed and sustainability of response.
A number of limitations should be considered. First, while the 2-component DAS28 score was developed to better correlate with synovitis and was validated by demonstrating increased association with radiographic progression in early RA, its weighting remains to be optimized in patients with established RA and those with RA receiving biologic DMARDs, and the implications for treatment decisions have yet to be established. Second, drug levels and antidrug antibody measures were available in only a subset of participants, even if similarly distributed across the 7 trajectory subgroups. Furthermore, adherence to treatments was self-reported. Third, although we assessed baseline comorbidities, adverse events leading to treatment cessation are potential confounders, especially in patients showing high rates of nonadherence, but we did not have sufficient information on those events to include them in the current analysis. A further limitation relates to the selection of participants in the BRAGGSS study, following guidelines from the UK National Institute for Care and Excellence (38). Our findings are therefore generalizable only to populations whose treatments and clinical characteristics are comparable to those outlined in these guidelines.
A final limitation pertains to the exploitation of unsupervised approaches and the exploratory hypothesis-generating nature of these approaches. However, we mitigated this effect by using clinicians' experience and the previous literature (i.e., suggesting 3 classes of responders as measured via the 4-component DAS28, which we replicate in this study), following a robust analytical framework (i.e., testing and retesting multiple models), and performing sensitivity analyses on 2 subsets to verify findings.
In pursuing the characterization of improvement in biologic changes in synovitis, relying solely on the 4-component DAS28 trajectories may result in misclassification of synovial inflammatory responses. The 2-component DAS28 trajectories provide finer granularity and may be a better measurement to reveal underlying biology. While there is currently no method implemented to predict class assignment for new patients, research to identify molecular biomarkers that can discriminate between trajectory classes at baseline will enable this. Thus, we suggest the use of the 2-component DAS28 score for biomarker discovery studies, which should encompass the understanding of adherence to treatment, immunogenicity to drugs, and general health status, including HAQ scores and comorbidities.
Better defined and physiologically contrasted phenotypes imply a better understanding of the underlying response mechanisms and are the pillar of precision medicine and well-powered biomarker discovery. Previous studies have shown that patients themselves want precision medicine approaches to receive medications that are likely to work more quickly (39). Therefore, more robust prediction models for biologic DMARD responses should be built on the basis of outcome measures reflecting the disease biology, provide evidence of larger variability in responses to biologic treatments, and allow for the identification of essential covariates to build these models. The present study shows how 2-component DAS28 trajectories capture these aspects better than traditional outcome measures and suggests that a more holistic view of patient care, considering all factors likely to influence response, should be applied in order to move toward precision medicine approaches.