A Crowdsourcing Approach to Develop Machine Learning Models to Quantify Radiographic Joint Damage in Rheumatoid Arthritis

RA2-DREAM Challenge Data

The RA2-DREAM Challenge used radiographic images and standard SvH scores¹⁴ from 2 NIH-funded clinical studies, Bridges et al¹⁰ and Ptacek et al,¹¹ led by investigators at the UAB. A total of 674 sets of images with scores from 562 patients were used with some patients’ data at 2 time points (Figure 1). Patient demographic data, including self-reported race, were collected during enrollment. Race was categorized as Black, White, or other (including American Indian or Alaska Native, Asian, and Native Hawaiian or other Pacific Islander) and was included because it was available in each parent research study. Each radiographic set included 4 images: 1 from each foot and 1 from each hand and wrist, all taken in a posteroanterior direction. The SvH scoring system assesses JSN on 15 areas from the hand and wrist and 6 areas from the foot (range, 0-4 for each joint) and erosions from 16 joint areas of the hand and wrist (range, 0-5 for each joint) and 2 sides of 6 joints of the foot (range, 0-10 for each joint), resulting in a total score of 448. Importantly, some participants had no RA-related joint damage of the hands and wrists or feet.^10,15 Individual joints with no damage attributable to RA were used as controls (SvH score, 0) for algorithm training.

All sets of images and scores were allocated to training (367 radiographic images and SvH scores), leaderboard (119 images only), and final evaluation (188 images only) data sets. The weight-bin scheme was chosen to balance the number of patients of certain ranges of SvH scores in each database and also to decrease the amount of influence of very high or very low SvH scores would have, given that SvH scores range from 0 to 448. First, we set up 8 bins of SvH scores (or ranges: 0, 1, 2-3, 4-7, 8-20, 21-55, 56-148, and >148). The threshold of each bin (0, 1, 1.1, 1.95, 3, 4, 5, and 6) was defined by the log value of the upper limit of each bin, with the exceptions to the first 2 bins (score 0, indicating no damage, and score 1). We then use the weight (w_i) of 2^{bin_threshold} to generate equal-weighted date sets.

Challenge Procedures

The RA2-DREAM Challenge contained 3 subchallenges. Participants were tasked to develop automated methods to quickly and accurately quantify overall RA damage (subchallenge 1), JSN (subchallenge 2), and erosion (subchallenge 3) from radiographic images of hands and feet. To reach the goals, the challenge had 3 rounds: training, leaderboard, and final evaluation (Figure 1B). In the training round, teams were provided with a training data set to develop computational methods to quantify RA-related damage in typically affected joints. In the leaderboard round, teams refined their algorithms with the leaderboard data set through sequential submissions and feedback on performance scored against expert-curated SvH scores. Each team was allowed up to 3 submissions per week over 9 weeks, totaling 27 potential submissions (eFigure 1 in Supplement 1). Public leaderboards with weighted-RMSE were updated immediately after submission and evaluation. This helped teams to gain immediate feedback on their algorithms and modify them accordingly. In the final evaluation round, 2 submissions per team were allowed, and the final submission, including a written description, was evaluated to determine the team’s ranking in each subchallenge.

All teams submitted their models via the Synapse collaborative science platform (Sage Bionetworks) (eMethods in Supplement 1). Models submitted to the platform were automatically transferred to the UAB Cheaha supercomputer and executed on the challenge data to produce estimation files (Figure 1B). The challenge was finished on June 30, 2020.

Evaluation of the Performance of Algorithms

We assessed the performance of each team’s model by comparing teams’ scores of each joint (y_i) to the ground truth SvH scores (s_i) using a patient-weighted RMSE approach. To avoid minus infinity (log zero) of a healthy joint we added 1 to y_i and s_i before log transformation using the equation:

Baseline Model

A baseline model¹⁶ was developed by the RA2-DREAM Challenge organizers for comparison to the models submitted. We created the baseline model by training a 10-layer DL model without segmentation of the images to classify the damage of each joint. The subchallenge 2 and subchallenge 3 scores were summed to generate subchallenge 1 scores for the baseline model. The layered architecture followed the example described in Deep Learning with Python.¹⁷

Postchallenge Independent Validation

To validate the performance of all algorithms, we selected 50 sets of images and scores reflecting various degrees of damage from the TEAR trial.¹² In the TEAR trial, there were 2 sets of SvH scores for each set of images, generated by 2 independent readers. We ran all submitted models on the 50 sets of images and evaluated their performance against the mean SvH scores of the 2 readers. We also computed the Spearman correlation between the final evaluation and postchallenge independent validation data set to assess the reproducibility of the top-performing algorithms. The concordance index between the final evaluation and postchallenge independent validation was calculated using the concordance_index function in the Python lifelines library.¹⁸ A false discovery rate–corrected pairwise t test was also performed as a secondary comparison of performance for each team relative to the top-performing team in each subchallenge.

Reproducibility of Methods Evaluated by Bootstrapping

The reproducibility of submitted algorithms and the baseline model was determined using Bayes factor calculated for 1000 bootstrapped iterations¹⁹ of the final evaluation and postchallenge independent data sets. This analysis identified the methods that are better than, tied with, or worse than a reference method considering the performance on a random sampling of the data. In each subset, we reran all the algorithms to obtain scores for subchallenge 1, subchallenge 2, and subchallenge 3. We then used these scores to calculate Bayes factors for each subchallenge using the computeBayesFactor function of R package challengescoring²⁰ in comparison with the top-performing model using R statistical software version 4.2.1 (R Project for Statistical Computing).^{21,22,23,24,25,26}

Ensemble Modeling

As with previous DREAM Challenges, we explored the so-called wisdom of the crowds phenomenon by performing ensembling experiments using the models submitted in the final evaluation round.^{21,22,23,24,25,26} For each subchallenge, we aggregated the top performer alone, the top 2, the top 3, and so on, then evaluated these ensemble estimations by calculating the mean estimation for each joint across multiple algorithms using the previously described Bayes factor robustness analysis.

Statistical Analysis

For comparison of demographic characteristics, we reported mean (SD) or median (IQR) for continuous variables. Spearman correlation, 2-sided pairwise t test, Bayes factor calculation, and Kruskal-Wallis test were performed using Python or R packages. P values were 2-sided, and P < .05 was considered as significant. Data were analyzed from June 24, 2020, to September 7, 2021.

Selection of Top-Performing Teams for Subchallenges

All 3 subchallenges were evaluated separately based on the assumption that algorithms to score JSN may be different from those used to score erosions (Figure 1A) and that there may be an algorithm that could give an overall quantification of damage for both parameters. Patient sociodemographic and clinical characteristics used for the challenge and postchallenge independent validation are shown in Table 1. The mean (SD) age was 54.9 (13.2) years in the training data set, 51.4 (13.1) years in the leaderboard data set, and 53.5 (13.5) years in the final evaluation data set. Data sets included mostly women, with 307 (83.7%) women in the training data set, 108 (90.8%) women in the leaderboard data set, and 158 (84.0%) women in the final evaluation data set. These demographic characteristics are very consistent with typical populations of RA, which predominantly affects women.

The accuracy of each algorithm was assessed using the weighted RMSE metric, comparing the estimated scores with the expert curated SvH scores. The reproducibility of each submitted algorithm was evaluated using Bayes factor. Briefly, the final evaluation data set was subsampled, algorithms were scored on the subsampled data set, and Bayes factor was calculated after 1000 bootstrapped iterations. A Bayes factor greater than 3 was used to determine if a team’s algorithms outperformed another team’s algorithm; no ties were observed (Figure 2). Top performing teams across all subchallenges had better performance than the baseline model. Teams Shirin and HYL-YFG (Hongyang Li and Yuanfang Guan) were the top 2 ranked teams for subchallenge 1; teams HYL-YFG, Gold Therapy, and csabaibio were the top 3 for subchallenge 2; and teams Gold Therapy, csabaibio, and HYL-YFG were the top 3 for subchallenge 3 (Figure 2; eFigure 2 and eFigure 3 in Supplement 1). They were awarded a total of $50 000 in prize money (eMethods in Supplement 1). The RMSE measures how closely the estimated values from the submitted algorithm were to the known SvH scores. The RMSE of the baseline model was 4.65 for subchallenge 1, 0.71 for subchallenge 2, and 0.65 for subchallenge 3, while the lowest RMSEs were 0.44 for subchallenge 1, 0.38 for subchallenge 2, and 0.43 for subchallenge 3. Detailed descriptions of algorithm development from the teams are provided in the eMethods in Supplement 1.

Rigor and Reproducibility of Top-Performing Algorithms

After completion of the challenge, we validated the consistency of algorithms using the postchallenge independent validation data set from the TEAR study,¹² which contained lower-quality radiographs than those from the studies by Bridges et al¹⁰ and Ptacek et al,¹¹ suggesting that the postchallenge independent validation task was likely more difficult. We calculated the weighted-RMSE and compared algorithms’ performance with those in the final evaluation round (Figure 2). We also computed a concordance index ranging from 1 (same ordering) to 0.5 (random ordering) to 0 (inverse ordering) between the scores in the final evaluation round and those in postchallenge independent validation. The concordance indices were 0.71 for subchallenge 1, 0.78 for subchallenge 2, and 0.82 for subchallenge 3, suggesting that the methods were reasonably generalizable and these algorithms can be readily applied to new independent radiographic images.

There were some inconsistencies in top-performing algorithms between the final evaluation and postchallenge independent validation. For example, Team NAD in subchallenge 1 and Zbigniew Wojna in subchallenge 2 showed better performance than others in the postchallenge independent validation, suggesting these models may perform better than others in images with different resolution or quality or in data sets with differences in which joint areas had more damage, such as wrists vs proximal interphalangeal joints. These results suggest that combinations of algorithms may be more robust than individual models alone, which led us to develop and evaluate an ensemble approach.

Ensemble Models vs Individual Models

We calculated the means of estimations from the top 2 performing algorithms then added the following ranked teams (Figure 3). We performed the same bootstrap analysis as was performed for individual teams to assess the reproducibility of the ensemble models compared with the top-performing team measured by Bayes factor.

For subchallenge 1 (overall damage), the ensembled models were not significantly better than the top-performing method (Bayes factor >3) (Figure 3A). Interestingly, including estimations from team RYM substantially improved performance in subchallenge 1, although this algorithm was not one of the top-ranked models in subchallenge 1. We further explored this by evaluating the pairwise Spearman correlation for all individual final evaluation round estimations in subchallenge 1, subchallenge 2, and subchallenge 3 (eFigure 3 in Supplement 1). Notably, in subchallenge 1, team RYM had a substantially lower Spearman correlation (ρ = 0.65-0.73) with the better-performing algorithms than the pairwise Spearman correlations within the better-performing algorithms (ρ = 0.75-0.95) (eFigure 3 in Supplement 1). This suggests that the algorithm from team RYM may improve the ensemble by providing information that is orthogonal to the other top estimations, a phenomenon that has been previously observed.²⁷ After evaluating the distribution of the estimations, we found that the ensemble model that included the top 9 teams systematically performed better than the top-performing team (eFigure 2 in Supplement 1). In subchallenge 2 and subchallenge 3, estimations were less concordant between the top teams than they were in subchallenge 1 (eFigure 3 in Supplement 1). The mean of the estimations from several of the top-ranked models (up to 5 models for subchallenge 2 and up to 10 models for subchallenge 3) increased the estimation power of the algorithms substantially (Figure 3B and C; eFigure 2 and eFigure 4 in Supplement 1). The estimation tasks for JSN and erosions improved (subchallenge 2 and subchallenge 3) more with ensembling than the estimation task of overall damage (subchallenge 1). We also scored the ensemble estimations shown in Figure 3 using the Spearman correlation as an outcome metric, which showed similar results to the weighted RMSE (eFigure 4 in Supplement 1). Taken together, we found that ensembling top models yielded similar improvements in performance in all 3 subchallenges.

Impact of Individual Technical Approaches on Team Performance

To gain insights into the submitted models and aid refinements of these and other imaging algorithms, we systematically summarized the top-performing teams’ methods by conducting a postchallenge survey and evaluating the methodological approaches (Table 2).^{28,29,30,31,32,33,34} Of 13 teams included in final evaluation, 10 applied image segmentation before attempting to quantify damage. In general, teams that applied image segmentation achieved better performance for all 3 subchallenges (eMethods and eFigure 5 in Supplement 1), although the difference between using and not using segmentation was not statistically significant.

Of 13 teams in the final evaluation, 10 applied DL-based approaches. The mean scores among algorithms that used DL were not statistically significantly different (eFigure 5 in Supplement 1). Nor did we did find that teams that used ensemble approaches had more favorable scores (eFigure 5 in Supplement 1).

After the challenge, we performed visual examination of images in which there were significantly discordant scores between expert SvH and algorithm-generated scores (eMethods and eFigure 6 in Supplement 1). From individual joints in the final evaluation data set, we noted 201 outliers of 7896 JSN scores (2.5%) and 462 outliers of 8272 erosion scores (5.6%). One of the coauthors, a board-certified musculoskeletal radiologist (M.B.F.) rescored the outliers by visual inspection according to SvH methods and made corrections to 97 SvH JSN scores (1.2%) and 192 erosion scores (2.3%). The low number of revised scores suggests a high degree of accuracy in the SvH scores generated by the trained scorers. Two expert readers generated the original SvH scores in the TEAR trial.¹² The coefficients of variation between the 2 expert readers were 0.74 for overall scores (P = .15), 0.62 for JSN scores (P = .18), and 0.73 for erosion scores (P = .03). Together these results suggested that scoring erosion is more challenging than scoring JSN.

Limitations

This study has some limitations. Our results provide optimism for the future automated scoring, but overall performance is limited by the number of images used in the training. For algorithms to become truly robust, thousands of digitally acquired, annotated images from large-scale observational studies, clinical trials, and EHRs should be used for training, and models should be continuously benchmarked on new data.