Predicting Probability of Response to Tumor Necrosis Factor Inhibitors for Individual Patients With Ankylosing Spondylitis

Key Points

Question  Is it possible to predict short-term treatment response to tumor necrosis factor inhibitors (TNFis) accurately in patients with active ankylosing spondylitis (AS)?

Findings  In a cohort study of individual participant data from 1899 patients in 10 randomized clinical trials of TNFis, models based on baseline clinical variables were developed and validated to predict the probability of having major response or having no response after 12 weeks of receiving TNFis, with moderate to high accuracy.

Meaning  These findings suggest that the probability of initial response to TNFi treatment can be accurately predicted from baseline variables, which may facilitate personalized treatment decision-making.

Importance  Tumor necrosis factor inhibitors (TNFis) have revolutionized the management of ankylosing spondylitis (AS); however, the lack of notable clinical responses in approximately one-half of patients suggests important heterogeneity in treatment response. Identifying patients likely to respond or not respond to TNFis could provide opportunities to personalize treatment strategies.

Objective  To develop models of the probability of short-term response to TNFi treatment in individual patients with active AS.

Design, Setting, and Participants  This is a retrospective cohort study using data of the TNFi group (ie, treatment group) from 10 randomized clinical trials (RCTs) of TNFi treatment among patients with active AS, conducted from 2002 to 2016. Participants were adult patients with active AS who failed nonsteroidal anti-inflammatory drugs. Included RCTs were phase 3 and 4 studies that assessed the efficacy of an originator TNFi at week 12 and/or week 24, either compared with placebo or an antirheumatic drug. The cohort was divided into a training and a testing set. Data analysis was conducted from July 1, 2019, to November 30, 2020.

Exposures  All included patients received an originator TNFi for at least 12 weeks.

Main Outcomes and Measures  Outcomes included major response and no response based on the change of AS Disease Activity Score at 12 weeks. Machine learning algorithms were applied to estimate the probability of having major response and no response for individual patients.

Results  The study included 1899 participants from 10 trials. The training set included 1207 individuals (mean [SD] age, 39 [12] years; 908 [75.2%] men), of whom 407 (33.7%) had major response and 414 (34.3%) had no response. In the reduced logistic regression models, accuracy was 0.74 for major response and 0.75 for no response. The probability of major response increased with higher C-reactive protein (CRP) level, patient global assessment (PGA), and Bath AS Disease Activity Index (BASDAI) question 2 score and decreased with higher body mass index (BMI) and Bath AS Functional Index (BASFI) score. The probability of no response increased with age and BASFI score, and decreased with higher CRP level, BASDAI question 2 score, and PGA. In the testing set (692 participants; mean [SD] age, 38 [11] years; 533 [77.0%] men), models demonstrated moderate to high accuracy.

Conclusions and Relevance  In this cohort study, the probability of initial response to TNFi was predicted from baseline variables, which may facilitate personalized treatment decision-making.

Ankylosing spondylitis (AS) also known as radiographic axial spondyloarthritis (axial SpA), is a chronic inflammatory condition characterized by inflammation in the spine, peripheral joints, and entheses and extra-articular manifestations such as uveitis, psoriasis, inflammatory bowel diseases, and aortitis.¹ The treatment options for symptom control in patients with axial SpA have expanded dramatically in the past decades with the availability of tumor necrosis factor inhibitors (TNFis), interleukin-17 (IL-17) inhibitors, and Janus kinase (JAK) inhibitors, and TNFis are recommended when patients remain with active symptoms despite maximal tolerated nonsteroidal anti-inflammatory drug therapy.^2,3 Although randomized clinical trials (RCTs) have shown that TNFis are efficacious in patients with AS, approximately one-half of patients do not achieve a notable improvement, suggesting important heterogeneity in treatment response.⁴ As clinical trials and most observational studies measure and report responses at the group level, limited guidance is available to predict treatment responses in individual patients.

In clinical practice, both patients and clinicians are interested in knowing how likely a patient is to achieve a significant response to TNFis so that treatment plans may be personalized. Several patient characteristics have been associated with a favorable response to TNFis in AS, including young age,^5-9 short disease duration,^5,7,10 male sex,^7-10 human leukocyte antigen B27 (HLA-B27) positivity,^6,7,11 and elevated inflammatory markers.^5,7,9-12 Other features, such as body mass index (BMI; calculated as weight in kilograms divided by height in meters squared), have not been investigated, even though patients with obesity and other inflammatory arthritides have been found to have poorer responses to TNFis.^13-15 More importantly, previous studies did not consider how combining these risk factors could affect the response to treatment. In addition, although previous studies have suggested factors associated with better response at the group level, none have estimated the likelihood of response using a probability score for an individual patient.

Identifying patients who are unlikely to respond to TNFis is as important as identifying responders. Such patients may be considered for only short trials of TNFi if the expectation of benefit is low or be considered for treatment with other classes of biologics, such as IL-17 inhibitors or JAK inhibitors.

The objective of this study was to develop and validate predictive models for short-term treatment response to TNFis in patients with active AS. Our goal was to provide a method to calculate probability scores of achieving major response or having no response after initiation of TNFi treatment for individual patients that was easy to apply and interpret in clinical practice.

This is a retrospective cohort study using prospectively collected data from RCTs of TNFis in patients with active AS. We requested individual participant data from trial sponsors via 2 clinical data sharing platforms (Vivli and Yale University Open Data Access [YODA]) and aggregated the data in a pooled analysis. This study followed the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD) statement. The study was exempt from institutional review board approval and the requirement for informed consent, according to 45 CFR §46.

For individual patients, the main inclusion criteria were (1) adults with AS by the modified New York criteria; (2) individuals enrolled in a randomized, double-blind trial that assessed the efficacy of an originator TNFi at week 12 and/or week 24, either compared with placebo or an antirheumatic drug, such as sulfasalazine; (3) individuals randomly assigned to the TNFi group in the trial. Based on the inclusion criteria for participants, we included RCTs based on the following criteria: (1) RCTs registered at ClinicalTrials.gov by December 31, 2018, with a study status of complete and result status of has result; (2) RCTs evaluating the efficacy of a TNFi in adult patients with active AS; (3) AS defined as fulfilling modified New York criteria; and (4) RCTs using TNFis including originator adalimumab, certolizumab, etanercept, golimumab, or infliximab. Pooling data from different medications was justified because these medications are in the same class and target the same cytokine, have similar mechanisms of action, and are used interchangeably clinically. A prior network meta-analysis reported that these medications have similar short-term efficacy in treating AS.¹⁶

To ensure the homogeneity of the study population, we excluded studies of patients with non–radiographic axial spondyloarthritis (SpA) and other SpA. This exclusion was used to remove possible concerns that poor model performance might be attributable to the inclusion of patients with different predictors of response or likelihood of response. We excluded open-label studies and open-label extensions of blinded trials. Because data from the 2 platforms could not be downloaded for potential merging due to security provisions, we developed the models using data from 1 platform (Vivli; the training set) and performed external validation using data from the other platform (YODA; the testing set).

We examined 2 outcomes at week 12: (1) major response, defined as decrease in AS Disease Activity Score (ASDAS) of 2.0 or greater and (2) no response, defined as decrease in ASDAS of less than 1.1. ASDAS is a composite score measuring disease activity. It includes assessment of patient-reported total back pain, patient global assessment (PGA) of disease activity, peripheral joint pain and swelling, duration of morning stiffness, and C-reactive protein (CRP) level or erythrocyte sedimentation rate.¹⁶ Major improvement is defined as a decrease of 2.0 or greater, and clinically important improvement is defined as a decrease of 1.1 or greater.¹⁷ Major response in this study therefore equals achieving major improvement, with ASDAS decrease of 2.0 or greater at week 12, and no response as not achieving clinically important improvement at week 12. These definitions conform with current Assessment of Spondyloarthritis International Society–European Alliance of Associations for Rheumatology (ASAS-EULAR) recommendations for continuation of TNFi or switching to a different biologic after 12 weeks of TNFi.³ It is important to note that the converse of no response is achievement of ASDAS clinically important improvement. We used binary outcomes for the study because the treatment decision is binary (ie, continue or stop treatment).

We initially considered 27 clinical characteristics at trial baseline as potential predictors (eTable 1 in the Supplement), including demographic characteristics (age, sex); BMI; laboratory tests (HLA-B27 status, CRP level), disease characteristics (disease duration; history of uveitis, inflammatory bowel disease, and/or psoriasis; tender and swollen joint counts); concurrent use of sulfasalazine, methotrexate, or corticosteroids; and prior use of a TNFi. Patient-reported measures included individual Bath ankylosing spondylitis disease activity index (BASDAI) item scores, Bath ankylosing spondylitis functional index (BASFI) score, total back pain, nocturnal back pain, and PGA. We also considered the 36-Item Short Form Health Survey (SF-36) subscales of fatigue and mental health to tap the construct of fibromyalgia. We did not include comorbidities other than BMI or SpA-related conditions because there is no evidence that common comorbidities, such as hypertension and diabetes, affect the response to TNFis in patients with AS.

After aggregating participant data from all TNFi treatment groups by data platform, we categorized participants as achieving major response (yes or no), and no response (yes or no) at week 12 and used descriptive analysis to summarize the baseline characteristics of participants in each group. For the few patients who had missing outcome data at week 12, we used their value at the most recent prior evaluation (eg, week 8, week 6).

In the training set, HLA-B27 data was missing for 1.5% of participants, which we therefore converted to 2 complimentary variables, HLA-B27 positive (yes or no) and HLA-B27 negative (yes or no). Other variables had less than 1% missing values, and we used K-nearest neighbor to impute missing data. We excluded psoriasis, inflammatory bowel disease, and prior TNFi use as predictors because their prevalence in the training set was less than 5%.

Because 4 variables (tender joint count, swollen joint count, SF-36 fatigue subscale, and SF-36 mental health subscale) were reported in only some trials, we examined their variable importance in the subset of trials in which they were included, using both logistic regression (LR) and random forest (RF) models. In these preliminary analyses, none of these variables were among the top 10 predictors of major response, while only the tender joint count was among the top 10 predictors of no response (ranked eighth). We therefore omitted these as potential predictors, which allowed for testing of 21 predictors in all available trials. We optimized prediction models using 5 different machine learning algorithms, including LR, linear discriminant analysis, support vector machine, gradient boosting tree, and RF.

We used data on the Vivli platform as a training set for model development. For each algorithm, we used 5:1 cross-validation for tuning the hyperparameters based on average accuracy (eMethods in the Supplement). We then generated the model and examined the accuracy, receiver operating characteristic area under curve (ROC AUC), sensitivity, and specificity of the optimized model (eMethods in the Supplement). We followed the same procedure for each of the 5 machine learning algorithms. We considered AUCs greater than 0.714 to represent a large effect size.¹⁸ We considered values of sensitivity and specificity of at least 0.80 as high and those between 0.40 and 0.79 as moderate.¹⁹

Because collecting data on 21 predictors has associated information costs and may not always be feasible in clinical practice, we generated reduced models for clinic use. We used permutation feature importance to rank the importance of each variable.²⁰ Permutation feature importance is defined as the decrease of model performance when the values of a given variable are randomly shuffled. Based on the importance scores, we selected the 5 most important variables for LR models and the 3 most important variables for RF models, and repeated the process of model development.

To assess how data from individual trials affected the models’ results, we recomputed the LR and RF models iteratively leaving out one of the 6 trials in the test set. We qualitatively compared model performance and variable importance rankings across these iterations.

We then externally evaluated the full and reduced models in testing set (YODA platform), and examined the accuracy, ROC AUC, sensitivity, and specificity of each model. We used bootstrapping to compute confidence intervals for the measurement statistics on 1000 random samples that included 75% of patients (with replacement) in the testing set.²¹ We also constructed calibration curves that compared the predicted probability of each response to the observed responses in 10 strata (or bins), ordered by the predicted probability. Data were extracted and analyzed using python version 3.5 and scikit-learn version 0.23.2.²²

We included data on 1899 participants from 10 eligible RCTs of TNFis (eTable 2 in the Supplement) in patients with active AS in this analysis. These participant-level data were accessible via 2 data sharing platforms.^23-32 Among these trials, 1207 participants who were assigned to the TNFi groups from 6 studies were eligible and included for model development (training set), and 692 participants from 4 studies were eligible and included for model validation (testing set).

Baseline characteristics of these participants are summarized in Table 1 and Table 2. In the training set (Table 1), the mean (SD) age was 39 (12) years, 908 (75.2%) were men, 1013 (83.9%) were HLA-B27 positive, and the median (IQR) disease duration was 5.1 (1.1-12.6) years. Mean (SD) BMI was 25.6 (5.0), and the median CRP level was 1.2 (0.4-2.6) mg/dL (to convert to milligrams per liter, multiply by 10). Major response was achieved by 407 participants (33.7%), and 414 (34.3%) had no response. Similarly, in the testing set (Table 2), the mean (SD) age was 38 (11) years, 533 (77.0%) were men, 600 (86.7%) were HLA-B27 positive, and the median (IQR) disease duration was 2.5 (0.7-7.6) years. The mean (SD) BMI was 25.6 (5.1), and the median (IQR) CRP level was 1.3 (0.5-2.8) mg/dL. Major response was achieved by 284 participants (40.1%), and 206 (29.7%) had no response.

We optimized prediction models using each of the 5 different machine learning algorithms in the training set, based on accuracy for predicting responses (that is, the proportion of patients correctly classified as improved or not by the model) using all 21 predictors. The overall performance of all 5 algorithms to predict major response and no response were similar (eTable 3 in the Supplement). In predicting Major Response, the accuracy ranged from 0.72 (gradient boosting tree) to 0.74 (LR). ROC AUC ranged from 0.79 (gradient boosting tree) to 0.81 (LR), with high specificity and moderate sensitivity. In predicting no response, the accuracy ranged from 0.73 (gradient boosting tree) to 0.75 (LR). ROC AUC ranged from 0.76 (gradient boosting tree) to 0.79 (linear discriminant analysis), with high specificity and moderate sensitivity.

Based on the overall similar performance among the methods, we continued the analysis using models based on LR and RF methods, which are more familiar to clinicians. We also used both methods to generate simplified (reduced) models for predicting major response and no response, using permutation feature importance to select a subset of variables with the greatest predictive ability. The Figure shows the ranking of the 10 most important features for each full model. For major response, the 5 most important features in the LR model were CRP level, PGA, BMI, BASFI score, and BASDAI question 2 score (eg, severity of spine pain). The probability of achieving a major response increased with higher CRP level, PGA, and BASDAI question 2 score, and decreased with higher BMI and BASFI score. Consistent with the findings in LR model, the 3 most important features in the RF model were CRP, BASDAI question 2 score, and BMI.

For no response, the 5 most important features in the LR model were CRP level, age, BASDAI question 2 score, BASFI score, and PGA (Figure). The probability of having no response increased with older age and higher BASFI and decreased with higher CRP level, BASDAI question 2 score, and PGA. Again, consistent with the findings in LR model, the 3 most important features in the RF model were CRP level, age, and BASDAI question 2 score.

We generated reduced models using these smaller subsets of variables with 5:1 cross validation. The performance of the full models and reduced models was comparable (Table 3).

In addition, we examined the performance of the full models after iteratively leaving one trial out of the training set (eTable 4 in the Supplement). Model performance was very similar across iterations (eTable 5 in the Supplement), suggesting that the results were not overly dependent on a single trial. The most important variables identified in these iterations were also quite similar to the analysis that included all trials (eFigures 1-4 and eTable 6 in the Supplement).

We then externally validated the full models and reduced models in the testing set (Table 4). Each model was tested in participants with complete data for all needed variables (ie, without imputation). Therefore, for the full models that use 21 variables for prediction, the sample size was 177, while for reduced models that use 3 or 5 variables for prediction, the sample sizes (range 625-692, depending on the model) were close to the entire testing set. The full models achieved moderate to high accuracy of 0.71 (RF model for major response) to 0.76 (RF model for no response), moderate ROC AUC of 0.65 (RF model for major response) to 0.68 (LR for no response), moderate sensitivity, and high specificity. Results for the reduced models were similar. The calibration curves demonstrated that the models predicted well at both low and high probabilities of major response and no response (eFigures 5 and 6 in the Supplement).

At a major response prevalence of 0.25, the positive predictive values (PPV) for LR-based and RF-based models ranged from 0.49 to 0.60, and negative predictive values (NPV) ranged from 0.82 to 0.84 (eTable 7 in the Supplement). At a no response prevalence of 0.25, PPVs ranged from 0.61 to 0.77, and NPVs ranged from 0.81 to 0.83.

We created online calculators based on reduced LR models (only for review purpose; need further validation, not yet for use in clinical practice), which can be used to calculate the probability scores for major response or no response. For example, a patient with CRP level of 3.0 mg/dL, BMI of 20, PGA of 8 (of 10), BASDAI question 2 (back, neck, hip pain) score of 8 (of 10), and BASFI score of 3 (of 10) before starting TNFi, the probability of having major response after 12 weeks of TNFi treatment was predicted to be 73%, meaning the patient will most likely have a favorable response to TNFi. For another patient with CRP level of 1.0 mg/dL, age of 55 years, BASDAI question 2 score of 6 (of 10), PGA of 5 (of 10), and BASFI score of 7 (of 10) before starting TNFi, the probability of having no response is predicted to be 61%, meaning the patient will most likely have no response to TNFi after 12 weeks treatment.

Response to TNFi in patients with active AS is heterogenous, highlighting a need to better tailor treatment to patients based on their likelihood of response.⁴ We developed and validated predictive models that provide probability scores for major response and no response to TNFis after 12 weeks of treatment for individual patients with active AS.

Overall, our models demonstrated moderate to high accuracy and high specificity, using only information available at the start of treatment. Both the full and reduced models provided probability scores of having major response or no response at week 12, which could help clinicians and patients make personalized treatment decisions. The reduced models only included 3 or 5 variables that can be easily collected in clinical practice (CRP level, age, BMI, BASDAI question 2 score, BASFI score, and PGA), without specialized testing or documentation, which will facilitate clinical use. The 3-variable RF models do not require the PGA, which, in contrast to the BASDAI, is often not collected in routine clinical care of patients with AS.³³

In the external testing cohort, the accuracy of different reduced models ranged from 0.70 to 0.78. Based on the sensitivity and specificity, at a prevalence of 25% for major response, the PPVs ranged from 0.49 to 0.60 and the NPVs from 0.82 to 0.84. These results indicated that, for a given patient, if the predicted probability of having major response was more than 50%, the patient may or may not have a major response to TNFis; whereas if the predicted probability was less than 50%, it is likely that the patient will not have a major response. Similarly, at a prevalence of 25% for no response, the PPVs ranged from 0.63 to 0.77 and the NPVs from 0.81 to 0.83. So, if the predicted probability of no response was more than 50%, the patient will possibly have no response, while if the predicted probability was less than 50%, the patient most likely will not have no response (ie, they would respond to treatment to some extent). Consequently, the PPVs are somewhat low when applied to a group of patients in which the prevalence of a true major response or no response is low, while the NPVs would be quite high in this scenario. However, it is important to note that these models and their applications need further investigation using practicing data in a clinical setting.

Consistent with previous findings,^5,7,9-12 in our study, CRP level, BASDAI score, BASFI score, and age were among the most important predictors. In addition, we found higher BMI associated with lower probability of major response. Importantly, and in contrast to prior studies, our models integrated information from different risk factors to estimate the probability of treatment response for individual patients.

It is possible that prediction of the ASDAS-based responses may have been aided by the major contribution of the CRP level, which is heavily weighted in the ASDAS. It is important to note that even the reduced models for the ASDAS-based responses included variables not used in the calculation of the ASDAS, such as BMI, BASFI score, and age.

This study has some limitations. Smoking has been associated with poorer response to TNFis and shorter treatment adherence.³⁴ Data on smoking were not available in some of the trials we analyzed. Second, all the participants in the training set were TNFi naive, and the models may not be similarly predictive in patients who had prior exposure to TNFis. Similarly, heterogeneity among patients included in the trials could have influenced the predictors and final models, although the iterative leave-out analysis suggests our results were robust to variations among trials. Third, we did not include nonsteroidal anti-inflammatory drug intake in our model, because data on their use were not consistent across studies. Fourth, our models have not been tested in daily clinical practice or in patients with a diagnosis of axial SpA. We focused on responses at 12 weeks because this is a common and recommended treatment decision point,³ but some patients may have delayed responses or low disease activity state. We did not include results of magnetic resonance imaging, because although these may enhance prediction, it would not be practical to obtain this imaging prior to starting TNFi treatment in all patients in clinical practice. Additionally, we do not know if pharmacogenomic data would enhance the prediction.

The models developed and validated in this study provide probability scores for achieving major response or having no response to TNFi treatment among patients with AS; they can be used to facilitate personalized decision-making in clinical practice. Confidence in choosing TNFi treatment may be enhanced with a high probability score for major response. Absence of a response in a patient predicted to have a high probability may raise a question about adherence to treatment. Conversely, a course of TNFi treatment may be terminated quickly if nonresponse occurs in a patient predicted to have a high probability of no response. It will be important to develop similar models for response to other biologic treatments so treatment options can be prioritized based on a patient’s most likely response.

Accepted for Publication: January 25, 2022.

Published: March 15, 2022. doi:10.1001/jamanetworkopen.2022.2312

Open Access: This is an open access article distributed under the terms of the CC-BY License. © 2022 Wang R et al. JAMA Network Open.

Corresponding Author: Runsheng Wang, MD, MHS, Garden State Rheumatology Consultants, 2333 Morris Ave, Ste B-111, Union, NJ 07083 (runsheng.wang2046@gmail.com).

Author Contributions: Dr Wang had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Concept and design: Wang, Ward.

Acquisition, analysis, or interpretation of data: All authors.

Drafting of the manuscript: Wang, Dasgupta.

Critical revision of the manuscript for important intellectual content: Wang, Ward.

Statistical analysis: Wang, Dasgupta.

Obtained funding: Wang, Ward.

Administrative, technical, or material support: Wang, Ward.

Supervision: Ward.

Conflict of Interest Disclosures: None reported.

Funding/Support: Dr Wang is supported by the Rheumatology Research Foundation K Bridge Award. Drs Dasgupta and Ward are supported by the Intramural Research Program, National Institute of Arthritis and Musculoskeletal and Skin Diseases, National Institutes of Health. This publication is based on research using data from data contributors (AbbVie Inc and Pfizer Inc) that has been made available through Vivli, Inc.

Role of the Funder/Sponsor: The funder had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Disclaimer: Vivli has not contributed to or approved and is not in any way responsible for the contents of this publication. This study, carried out under YODA Project 2019-3893, used data obtained from the Yale University Open Data Access Project, which has an agreement with Janssen Research & Development, LLC. The interpretation and reporting of research using this data are solely the responsibility of the authors and do not necessarily represent the official views of the Yale University Open Data Access Project or Janssen Research & Development, LLC.

Additional Contributions: The authors thank the Spondylitis Association of America for their support for Dr Wang’s research and thank Yang Liu, PhD (Virtu Financial), for discussion on statistical methods and Shuoxin Liu (City College of New York) for building the online calculator. Shuoxin Liu was compensated for his time.