Prognostic modeling in early rheumatoid arthritis: reconsidering the predictive role of disease activity scores

Objective In this prospective cohort study, we provide several prognostic models to predict functional status as measured by the modified Health Assessment Questionnaire (mHAQ). The early adoption of the treat-to-target strategy in this cohort offered a unique opportunity to identify predictive factors using longitudinal data across 20 years. Methods A cohort of 397 patients with early RA was used to develop statistical models to predict mHAQ score measured at baseline, 12 months, and 18 months post diagnosis, as well as serially measured mHAQ. Demographic data, clinical measures, autoantibodies, medication use, comorbid conditions, and baseline mHAQ were considered as predictors. Results The discriminative performance of models was comparable to previous work, with an area under the receiver operator curve ranging from 0.64 to 0.88. The most consistent predictive variable was baseline mHAQ. Patient-reported outcomes including early morning stiffness, tender joint count (TJC), fatigue, pain, and patient global assessment were positively predictive of a higher mHAQ at baseline and longitudinally, as was the physician global assessment and C-reactive protein. When considering future function, a higher TJC predicted persistent disability while a higher swollen joint count predicted functional improvements with treatment. Conclusion In our study of mHAQ prediction in RA patients receiving treat-to-target therapy, patient-reported outcomes were most consistently predictive of function. Patients with high disease activity due predominantly to tenderness scores rather than swelling may benefit from less aggressive treatment escalation and an emphasis on non-pharmacological therapies, allowing for a more personalized approach to treatment. Key Points • Long-term use of the treat-to-target strategy in this patient cohort offers a unique opportunity to develop prognostic models for functional outcomes using extensive longitudinal data. • Patient reported outcomes were more consistent predictors of function than traditional prognostic markers. • Tender joint count and swollen joint count had discordant relationships with future function, adding weight to the possibility that disease activity may better guide treatment when the components are considered separately. Supplementary Information The online version contains supplementary material available at 10.1007/s10067-024-06946-z.


Introduction
Rheumatoid arthritis (RA) is an inflammatory disease of the synovial joints that often leads to progressive joint damage and disability [1].The treat-to-target treatment (T2T) approach, wherein disease-modifying anti-rheumatic drug (DMARD) therapy are escalated until a predefined disease activity target has been achieved, has been recommended therapy for more than two decades [2].Nonetheless, some patients experience ongoing symptoms and loss of function, often with progressive joint damage [3].Improving the prediction of treatment response and outcomes is vital to guide and optimize treatment decisions [4].Prior to the T2T era, factors that predicted disability according to a 2003 systematic review included age, female sex, rheumatoid factor positivity, high pain scores, low SES, joint tenderness, and depression [5].There is currently no consensus regarding prognostic models in patients receiving the T2T approach appropriate for clinical use due to the historically small, demographically uniform cohorts, changing treatment regimens, and lack of external validation [6].
Functional outcomes in RA are frequently measured by the Health Assessment Questionnaire (HAQ), a validated and widely used instrument that was first developed in 1978.It is intended to integrate function as it relates to structural damage, disease activity, pain, and psychosocial factors [7], while the modified HAQ (mHAQ) is an abridged version to improve feasibility of use in practice [8].We conducted a case series study in a cohort of patients with early RA recruited over the period 1998 to 2021.The clinic was an early adopter of the treat-to-target approach, which poses a unique opportunity to investigate predictors of response to relatively contemporary treatment strategies.Identifying which baseline factors are amenable to treatment lends itself to more personalized treatment approaches.Four multivariable prediction models were developed to investigate the predictive factors of the mHAQ.Each model was internally validated using cross-validation [9].

Study population
The study population comprised patients enrolled in the early RA cohort at the Royal Adelaide Hospital from 1998 to 2021.Consecutive patients diagnosed with RA which met the 1987 [10] or 2010 ACR-EULAR [11] criteria for classification with RA (depending on enrolment year) were eligible.While anti-cyclic citrullinated peptide (ACPA) was only included in the 2010 classification criteria, it has been measured in this cohort in sera stored since its inception.
Patients provided informed consent for the use of their data for research (CALHN HREC approval 120,618), and subsequent approval was obtained for extraction and use of the data in this study (CALHN HREC approval 15,056).Patients were included in this cohort if the onset of symptoms of RA occurred within the preceding 12 months, they were DMARD-naive, and they were over the age of 18 [12].Patients were managed according to predefined treatment algorithms whereby they commenced triple therapy with methotrexate, sulfasalazine, and hydroxychloroquine.DMARD up-titrations were made according to previously published algorithms [13].While the treatment protocol changed in 2005 to include the use of biologic DMARDs (bDMARDs) and subsequently targeted synthetic DMARDs (tsDMARDs), the proportion of patients who commenced on DMARDs remained low due to the strict criteria for subsidized therapy in Australia [14].Data for up to 5 years post-initial diagnosis were included, given the attrition rate in a real-world population.Patients were excluded if they did not have an initial mHAQ score recorded or if data were recorded more than a month after treatment commencement (n = 3).

Outcome
The primary outcome was the mHAQ score, an extensively validated measure of self-reported function [7].Patients were assessed by a rheumatology nurse trained in metrology, and the treating rheumatologist was blinded to the mHAQ score.The resulting summed raw score ranges from 0 to 24, with a higher score indicating more dysfunction.

Candidate predictors
We considered several candidate predictors measured at baseline (defined as treatment initiation) based on subject matter knowledge and prior findings in the literature.There was a high degree of collinearity between C-reactive protein (CRP), erythrocyte sedimentation rate (ESR), and the disease activity composite measures such as the simplified disease activity index (SDAI), clinical disease activity index (CDAI), DAS28-CRP, and DAS28-ESR.We elected to include only the components of the DAS28-CRP, which combines the patient-derived or influenced measure of patient global assessment measured on a 100 mm visual analogue scale (VAS) and tender joint count (TJC) and the physician-derived measure of swollen joint count (SJC) and CRP.
DMARD use was included as a potential explanatory variable, categorized as (a) mono-or dual therapy, (b) triple therapy, (c) added leflunomide, or (d) added any other drug (i.e., cyclosporine, gold, b/tsDMARDs).Medications at 1 year were used as this allowed for the greatest separation between patients given all were managed with the same treat-to-target algorithm.As sample sizes were small and the patient cohort was predominantly of European ancestry, ethnicity was coded as a binary variable ("European ancestry" and "non-European ancestry").Socioeconomic status (SES) was based on postcode data and the Socio-economic Indexes for Areas, which divides postcodes into deciles resulting in a score from 1 to 10 with 1 indicating the most disadvantaged areas and 10 indicating the most advantaged areas [15].

Statistical analysis
Continuous variables were reported with means and standard deviations (SD) or medians and interquartile ranges (IQR) while categorical variables were reported as percentages and frequencies.Bivariate analyses were based on t-tests or Mann-Whitney U tests to investigate associations between continuous variables and binary variables.Bivariate associations between continuous variables were investigated using Pearson's correlation coefficients or Spearman's ranked correlation coefficient.Associations between two binary variables were investigated using chi-squared tests.Variables that were associated with baseline mHAQ (p < 0.10) were chosen to be included in multivariate models to reduce model complexity by limiting the degrees of freedom.The distribution of residuals was analyzed using the Kolmogorov-Smirnov (KS) test to ensure the assumption of normality was not violated for use of a standard GLM.
Four multivariate predictive models were developed that varied in the way the variables and outcomes were used: 1) A generalized linear model (GLM) [16] to predict baseline mHAQ using the variables collected at baseline 2) A GLM (Poisson) predicting mHAQ at 1 year based on variables collected at baseline 3) A GLM (Poisson) predicting mHAQ at 18 months based on variables collected at 6 months post diagnosis 4) A linear mixed effects (LME) [17] longitudinal model to predict serially measured mHAQ based on variables collected contemporaneously Variables were chosen by backward selection.Model 4 used an α of 0.05 given it had a greater a priori power as it took advantage of repeated measures.The Akaike information criteria (AIC) was used for variable selection in the smaller datasets, which corresponded to an α of 0.157 [18].The performance of each model was estimated using tenfold cross-validation, which has been shown to achieve minimal bias [9].Multiple imputation using chained equations (MICE) was performed on the training segment of the data to impute missing values (50 imputations), and resultant models were pooled using Rubin's rules [19].The missing values in the validation data were imputed using the imputation models initially developed on the training data, to avoid train-test contamination.Variables with greater than 10% of missing values were removed to avoid biasing our analysis [20].The final pooled model was tested on the validation data.The mHAQ was dichotomized at a threshold of 0.25 to construct receiver operator curves (ROC) and investigate discriminative performance.We selected this value as it approximates the minimum clinically important difference identified in previous work [21] and thus could be considered to constitute a threshold for clinically meaningful disability.
The root mean squared error (RMSE) and coefficient of determination (R 2 ) were calculated for each model.The variables that were significant across the majority of cross-validation folds (greater than 50%) were included in a final model using the entirety of the data, in order to report the variable coefficients and p values.Backward selection was conducted again at this point to remove variables with a p value above the threshold.Statistical analysis was conducted using Python 3.8.8[22], statsmodels 0.1 [23], and pyMIDAS [24].

Descriptive statistics
Patients had an initial median mHAQ of 0.6 (IQR 0.3-1.1),which improved to 0.0 (IQR = 0.0-0.375) at 1-year post-diagnosis (Table 1).At baseline, 27.9% of patients had a mHAQ of zero, 95.7% at 1 year and 93.7% at 18 months, indicating significant functional improvement with treatment.The median mHAQ score had a high initial value, with a rapid drop in mHAQ before plateauing for 5 years post-diagnosis (Fig. 1).Additionally, the median reported 5-year post-diagnosis mHAQ per calendar year did not change significantly, despite the introduction of biologics in 2005 (Supplementary Fig. 1).This suggests that treatment efficacy did not change significantly since inception of the cohort (p = 0.43).

Bivariate analysis
Bivariate analyses among variables collected at baseline were conducted to identify the variables to include in multivariate models and to investigate for collinearity that was not previously identified (Fig. 2).A total of 25 variables were included as potential predictors in our analyses (Table 1).Variables found to be associated with mHAQ score at baseline were patient global assessment, stiffness, fatigue, physician global assessment, TJC, SJC, BMI, patient pain, CRP, and socioeconomic status (SES) based on postcode.

Linear mixed effects model
A total of 397 patients with 8264 clinic appointments were used in LME longitudinal analyses.The median length of follow-up was 236 weeks (IQR 133.0-252.0weeks).The median number of clinical assessments in the cohort was 21 (IQR 14-28).

Discussion
We conducted an analysis of an early RA cohort receiving DMARDs according to a T2T approach to predict response to treatment in terms of function as measured by the mHAQ score.We developed three GLMs and one LME model to predict current and future mHAQ scores in 397 patients.The performance of the statistical models fitted ranged from an AUC of 0.64 to 0.88, R 2 ranging from 0.2 to 0.53, and RMSE ranging from 0.29 to 0.40 (Fig. 3), with the LME model being the highest performing.These results are comparable to previous work which have reported a range of AUCs from 0.78 to 0.82 with internal validation [6].
The predictive variables identified across all models were also largely consistent with previous work.Initial mHAQ has been the most consistent and strongest predictor across many studies [6], as it was in our analyses.In particular, when predicting function at 18 months, 6-month mHAQ was the only variable selected, demonstrating the value of 6-month mHAQ as a prognostic marker.Early response to DMARDs has frequently been noted as an important prognostic sign [25], which our study confirms.The components of the DAS28-CRP, TJC, patient global assessment, and CRP, were predictive of baseline mHAQ, with higher values predicting greater disability.SJC was not predictive.However, when predicting future mHAQ at 1 year, higher TJC predicted higher future mHAQ, while higher SJC predicted lower future mHAQ.This lends weight to the idea that overt signs of synovitis such as joint swelling could suggest dysfunction that is more readily modifiable with DMARDs, while tender joints alone may indicate symptoms that may not be entirely due to synovitis and hence less responsive to DMARDs.Conversely, CRP and patient global assessment were not associated with future function.
In regard to the longitudinal model, patient-derived factors made up the majority of predictors including pain, early morning stiffness, patient global assessment, TJC, and fatigue, in addition to indicators of inflammation, CRP and SJC, as well as physician global assessment.The duration of disease was also a useful predictor, with longer time since diagnosis suggesting slightly improved function.While previous work identified a "J-shaped" pattern to progression [5], with initial drop followed by slow increase in the HAQ score, at least within a 5-year time frame, this does not appear to be the case in our cohort.The identified variables were almost the same as those found in the model predicting mHAQ at baseline, the only difference being the inclusion of fatigue.Regardless of the duration of illness, the variables that predict mHAQ at the same visit remain consistent.
Female sex and age at baseline have often been identified as significant predictors of the mHAQ score in previous publications, although a recent systematic review found that 5 out of 18 (28%) prior studies did not identify an association between age and HAQ and 6 out of 21 (29%) did not identify an association between female sex and HAQ [26].Our findings were consistent with prior studies that did not find an association of age and sex with the mHAQ score.Both age and sex were excluded from further multivariate analysis in our study on the basis of initial bivariate analysis.It is also possible that our study lacked statistical power to detect modest associations of age and sex with mHAQ score.Previous studies have identified associations between marginalized ethnic groups and function [27].European vs non-European ancestry was excluded as a predictive variable from the results of initial bivariate analysis, with no association to the mHAQ identified.Our findings may reflect the small size and heterogeneity of the "non-European ancestry" group in our study cohort.Common comorbidities such as depression, osteoarthritis, chronic back pain, and fibromyalgia were also not useful predictors, despite prior evidence that comorbidities impact measurement of the mHAQ [28].However, most studies on comorbidities use a composite measure such as the Charlson Comorbidity Index [29] or the Rheumatic Disease Comorbidity Index [30], which span many more diseases.Additionally, we suspect our cohort suffered from underdiagnosis of comorbidities, particularly for fibromyalgia, in that only 6.6% of patients received a diagnosis, despite the prevalence estimate in RA patients ranging from 18 to 24% [31].We suspect because we did not have access to a more thorough measurement of comorbidities that the associations were not detected in our baseline bivariate analyses.
RF and ACPA were also not found to be predictors of mHAQ score in the current study, consistent with other studies involving contemporary patient cohorts.A 2018 systematic review found that 8 out of 11 (72.7%)papers did not identify any association between RF and HAQ [26].The 3 studies that did used data collected between 1979 and 1998 [26].Similarly, 5 out of 6 (83.3.%)papers found no association between ACPA and HAQ, although the study that found an association used a relatively large dataset (n = 1995 patients) with patients recruited from 1990 to 2009 [32].The changing predictive role of RF and ACPA may be because patients who previously would have progressed quickly now receive prompt diagnosis and treatment and no longer progress to significant dysfunction as was the norm in the past.A recent large observational study of 3251 patients suggested that RF/ACPA positivity might indicate better response to therapy [33].
Our work suggests that persistent dysfunction despite contemporary treatment may indicate a need to consider approaches other than continual treatment escalation.Recent work using the same patient cohort as analyzed in the current study has highlighted the role of non-nociceptive pain in disease activity [34].The authors note that the presence of joint pain in the absence of synovitis (i.e., a high TJC and low SJC) implies other potential mechanisms may explain high disease activity scores such as central sensitization to pain.Identifying patients who are likely to have central sensitization and thus a lesser response to treatment escalation is an important factor in guiding treatment decisions.For example, these patients could be directed down a treatment path that less aggressively escalates therapies and emphasizes alternate approaches such as physical rehabilitation, patient education, and psychosocial interventions [35].A threshold of disease activity that returns the patient to the conventional treatment path would be important, particularly where there is evidence of synovitis.While there have been discussions about the role of treatment de-escalation in certain patients, this approach would allow that to be done proactively rather than reactively, enabling personalized medicine for RA patients.
Another possibility to further improve the reliability of models to allow for eventual clinical implementation is the use of imaging of disease-related joint damage as a candidate predictor.This has not been possible until recently due to the infeasibility of systematically quantifying structural damage in the clinic.The advent of deep learning has offered methods of automatically identifying patterns in imaging data that may be more predictive than clinical and demographic data.While plain radiographs may not be sufficiently sensitive to subtle changes to improve model performance, MRI and ultrasound are being used with increasing frequency and may offer a way to bolster the ability to identify features of RA.This idea has been studied in osteoarthritis by predicting patient pain from knee X-rays, with promising results [36], but to the best of our knowledge remains unexplored in RA.
Our study has been conducted in line with best practices for prognostic modeling as specified in the TRIPOD guidelines [18].We handled missing data with multiple imputation, estimated performance with cross-validation, and used backward elimination for feature selection.While we had extensive longitudinal data, the patient group was demographically uniform and of moderate size, limiting the generalizability of our models.Our dataset sizes were sufficient based on the commonly accepted 1 to 10 rule of variables compared to degrees of freedom [37].Unfortunately, we were only able to conduct internal validation of our models due to the difficulty of accessing replication data.Independent datasets from a diverse range of institutions would allow for more thorough validation of our models.Despite this limitation, our study achieved the goal of identifying relationships between patient factors and function and confirming whether previous findings remained consistent in our cohort.

Conclusion
Accurate prediction of response to therapy in patients with RA is critical for guiding treatment decisions and offering avenues for precision therapy.Our modeling in a cohort receiving long-term T2T DMARD therapy was mostly consistent with the previously identified variables predictive of function as measured by the mHAQ.However, our study also identified that the components of the DAS28-CRP, TJC and SJC, have conflicting relationships with future function, suggesting that these factors might be better considered separately when guiding treatment decisions.

Fig. 3
Fig. 3 Receiver operator characteristic curve for each model including the AUROC, R 2 , and RMSE and their 95% confidence interval.a) GLM to predict baseline mHAQ, b) GLM predicting mHAQ at one year from baseline variables, c) GLM predicting 18-months mHAQ from variables collected at 6-months, d) a LME longitudinal model to predict serially measured mHAQ from variables collected contemporaneously

Table 1
Characteristics of study population (n = 397)

Table 2
Predictive variable for each model, coefficients, and p values