随机森林模型在膝关节炎患者结构特征与症状定量分析中的应用

萧伊，肖峰，徐海波

Introduction

Osteoarthritis (OA) is one of the leading causes pain and disability in the aging population[1]. The main symptoms pain, stiffness, and limited function often localize to joints commonly affected by the structural changes of OA, and can be relieved by surgical procedure that targeting the damaged structure[2-3].

This fact has fostered the belief that the structural changes cause the symptoms. However, most studies of OA showed the structural abnormality on radiography was at best a weak risk factor for pain or disability[2，4].Although MRI is a more advanced imaging modality than radiography and can directly visualize the articular structure defects, MRI findings bear inconsistent and moderate relationship to the severity of the symptoms[5]. Lack of knowledge about structural determinants to OA symptoms limits specifically targeted administration strategies in clinical practice to slow disease progression and to ease suffering. Thus, there is a need to develop a method to accurately associate imaging structural findings with knee symptoms. Quantitative analysis and modeling of the structural characteristics and symptoms of the knee joint can help us better understand the pathology of the disease and help patients get effective evaluation and timely intervention.

The purpose of this study was to investigate whether these imaging based quantitative features can predict the symptoms of the knees. Further, we aimed to evaluate which of these features show the highest predictive value in each distinct knee symptom.

1 Material and Methods

1.1 Patients

This study was based on the osteoarthritis initiative[6](OAI, http://www.oai.ucsf.edu/datarelease/). OAI is a multi-center, population-based longitudinal cohort study, targeted to identify the risk factors associated with the onset and/or progression of knee osteoarthritis, and to explore potential biomarkers of the disease.

All participants were selected from the 600 subjects in a single center project in OAI, Osteoarthritis Biomarkers Consortium FNIH Project (OA FNIH)[7]. All of the included subjects were between 45 and 79 years old and consisted of different ethnic backgrounds. The exclusion criterions of this study were: (1) missing knee structural imaging evaluations or symptom scoring indicators (described in next subsections); (2) with rheumatoid or other inflammatory arthritis; (3) bilateral end stage KOA or inability to walk without aids. Finally, 551 eligible participants were retained in this study. They were randomly divided based on 7:3 ratios, such that 385 patients got assigned to the training set and the remaining 166 to the test set. Their demographic information of the participants was shown in Table 1 informed consent was obtained from all participants and the study was approved by local ethics committees.

1.2 Structure image quantitative evaluation

Five programs were used to quantitative evaluate the structure feature of the knees in this study. They were: (1) kXR FNIH bone trabecular integrity (Duke)[8-9]. Six features were included to access the vertical (compressive) and the horizontal (tensile) trabecular integrity of the knee joint. (2) kMRI FNIH QCart (Biomediq)[10-11]. Seven features were included to measure the cartilage volume and meniscal volume in several locations of the knee joint based on an automatic computer-based method. (3) kMRI FNIH Boneshape (iMorphics)[12-14]. Twelve features were included to evaluate total area of subchondral bone for 9 anatomical locations and 3D shape of the bones in the knee joint.(4) kMRI FNIH SBP (Qmetrics)[15]. Thirty-six features were included to access 9 statistical characteristics for 4 anatomical locations in the knee joint. (5) kMRI FNIH quantCartilage (Chondrometrics)[16-18]. Ninety two features were used to measure cartilage volume, thickness and other associated characteristics for main cartilage plates and their sub-regions in the knee joint.

Combining the above five sets of feature sets, for each patient, a total of 153 quantitative features were included for comprehensive evaluation of the structure of the knee joint, including bone, cartilage, meniscus and their sub-regions. Because the features came from different evaluation programs, there might be a partial correlation or redundancy between features, which should be noted in later modeling.

1.3 Symptom scoring system

The Western Ontario and McMaster Universities (WOMAC) osteoarthritis index[19]was used to evaluate the knees’ symptom of participants included in this study. In this scoring system, WOMACADL, WOMACPK and WOMACSTF were used to access the knees’ physical function, pain and stiffness extent of the patients, respectively. Besides, WOMACTS was a total score of the three scores above. A score equal to 0 in one of the above items indicates that the corresponding symptom of the knee joint is normal. In addition, a higher score means that the corresponding symptom is worse.

1.4 Features selection

Two feature selection methods were implemented to screen the features as follow steps: (1) Correlation analysis. Pearson correlation was employed to eliminate high dimensional feature redundancy between the features. If two features were highly correlated (Pearson correlation coefficient＞0.9), the one less Pearson correlated to WOMAC score was excluded. (2) Minimum-redundancy maximum-relevance (mRMR) method[20]. The remaining features were ranked using mRMR algorithm by calculating the mutual information (MI) between features and WOMAC scores. The mRMR ranks the input features by maximizing the MI with respect to the WOMAC scores and minimizing the average MI of higher ranked features. This method allows an efficient selection of relevant and non-redundant features.

1.5 Random forest modeling

Random forest (RF) was used as the predictor model of WOMAC scores based on the selected features. RF[21-22]is a nonlinear model composed of multiple decision trees. Because its advantages of noise insensitive and adjustable inner parameter, RF is adaptive in different application domains[23-25]and can effectively avoid overfitting. While RF modeling, the correlation between features and the importance of each feature can be detected, which make it suitable in the WOMAC predictors modeling.

In this study, Grid search of the RF model parameters: max_tree_num, min_samples_leaf and max_depth were employed to determine the optimal model parameters. RF modeling was conducted 10 times to obtain the averaged performance of the predictors. The metrics: R-square, mean absolute errors (MAEs) and mean squared errors (MSEs) were used to evaluate the constructed RF predictor performance in training and test sets respectively.

1.6 Statistical analysis

Statistical analysis in patients’ demographics was performed using SPSS 21.0 (IBM Corp., Armonk, New York, USA). The independent Student’s t-test was utilized to compare continuous variables, which conform to a normal distribution, while the Chi-square test was used to compare categorical variables between the two groups. A two-sided P＜0.05 was considered statistically significant.

The predictors modeling analysis were performed using R software (version 3.6.1; http://www.r-project.org) and python software (version 3.7.0; http://www.python.org). The R package “mRMRe” was used to implement the mRMR feature selection, while the python package “sklearn” was used to implement and evaluate RF predictors modeling.

2 Results

2.1 Participants characteristics

The demographic information and WOMAC scores of all 4796 participants in OAI and the 551 cases included in this study were shown in Table 1 P-value in Table 1 referred to the significance of the difference between the training set and the test set. The statistical results showed that there are no significant differences between the training set and the test set in all indicators and WOMAC scores, implying that there is no significant difference in the data distribution of the training set and the test set. In addition, the average value and composition ratio of the subjects in this study are very close to that of the entire OAI data set, indicating that the analysis results obtained based on the included data are reliable.

2.2 Features analysis

Two feature selection methods are used for feature selection. First, the original 153 features were screening through correlation analysis between the features, and highly linearly correlated features were eliminated. Then, the mRMR method was used for further screening. The criterion for this screening procedure was mutual information.

Figure 1 shows the feature selection procedure in theprediction and modeling of WOMACADL. There are a mass of highly linear correlations between the original 153 features (Pearson correlation coefficient＞0.8, Figure 1A). Highly redundant features were removed and 77 features were retained after correlation analysis. In figure 1B, it can be found that the Pearson correlation coefficients between features have decreased to below 0.9. After further screening, 20 features were obtained, and the correlation between them was further reduced (Figure 1C). Similarly, for WOMACPK, WOMACSTF and WOMACTS, 20 features were also obtained respectively, as shown in Table 2. Due to the maximizing correlation criterion of mutual information, the obtained features ranked top 20 of the mutual information scored with WOMAC.

Tab. 1 Demographics of the included patients

Tab. 2 The 20 features obtained in the WOMAC predictors after the features selection procedures

Tab. 3 Inner parameters in RF predictor models determined using grid search

Combining the features obtained based on the four WOMAC scores (Table 2), we found that some of the structural image features that affect the different symptoms (physical function, pain, stiffness) of the knee joint are highly consistent: (1) The cartilage volume in the middle femur in Chondrometrics: BMFVCL; (2) Three vectors in Imorphics: nFemurOAVector, nTibiaOAVector, andnPatellaOAVector, representing the shape of the femur, tibia, and patella, respectively. The total area of the subchondral bone of the medial femoral condyle and the lateral tibial plateau: MF_tAB and LT_tAB; (3) The cartilage volume of the medial tibia in Biomediq: MedialTibialCartilage; (4) Statistical descriptors of the curvature of the subchondral bone in medial femur in QMetics: CurAverage_MedFem, Cur_5pc_MedFem, Cur_95pc_MedFem and SCRSD_MedFem.

Tab. 4 Predictor performance in training and test sets

In addition, the structural features that affect knee pain and stiffness at the same time include statistical descriptors of the curvature of the subchondral bone in the lateral femur and medial tibia. However, it should be noted that the statistical descriptors of the curvature of the subchondral bone in lateral tibial has a smaller effect on knee symptom than that in other locations of knee joint.

2.3 RF modeling

First, grid search was performed for the RF predictor modeling to determine the internal parameters of the RF model, as shown in Table 3 respectively. Then, the RF predictor models were trained and tested with the determined internal parameters. Ten-fold cross-validation based on R-squared, MAEs and MSE were shown in Table 4 The comparison between the predicted value and the actual value was shown in Figure 2. Finally, Figure 3 showed the ranking of the importance scores of all 20 features used for modeling.

As shown in Table 4, the R-square of the four predictor test sets are all above 0.65, indicating that more than 65% of the data can be interpreted and effectively predicted by the constructed predictor model. The prediction errors were all limited to a relatively low error level: 4.5 for physical function, 1.3 for pain, and 0.6 for stiffness (MAEs). However, from the prediction results of each data point in Figure 2, it can be found that most of the data points with better predictions are concentrated in the center of the data distribution, and the prediction results for both sides of the data distribution are significantly worse. Data points with a normal symptom (score equals 0) were usually predicted to be high, while data points with a poor symptom (large score) were predicted to be low.

Figure 3 showed the contribution of input features in RF predictor model to the prediction results. It can be found that the medial cartilage volume of the medial femur (BMFVCL) contributed a lot to WOMAC physical function, pain, and total score prediction (ranking top 2), and moderately to stiffness prediction (ranking around 10). Therefore, BMFVCL was a key structural feature that deserves attention. The statistical descriptors of the curvature of the subchondral bone in medial femur also contributed greatly to physical function, pain, and total (ranking top 5), and moderately to stiffness (ranking 5-10). This is highly consistent with the contribution of BMFVCL, which might indicate that the structural characteristics of cartilage and/or subchondral bone on the medial femur may greatly affect the symptom of the knee joint.

3 Discussion

A total of 153 quantitative structural features from 5 feature sets were included in this study. After two steps of feature selection, 20 features were retained to construct and validate knee joint symptom predictor. Finally, the correlation between selected structural features and the symptom of knee joint was analyzed. This research explored the relationship between knee joint structure and symptom, and the constructed model can be effectively used to predict and evaluate knee joint symptom based on knee joint structural features.

3.1 Modeling method

A two-step feature selection method is adopted in this study. The first step of correlation analysis mainly focused on the linear correlation between features. The mutual information criterion in mRMR method in the second step further takes into account the nonlinear correlation. The optimization criterion of the greedy search was employed in mRMR[26]to search the optimal feature subset. When the number of features is large and/or their redundancy is large, this optimization method is easy to fall into the local optimum, which leads to a feature set with poor combined predictive performance. The reduction in the number and redundancy of features after the initial screening in first step can not only speed up the convergence rate of the algorithm, but also avoid local optimum features set with poor combined predictive performance.

There are two different fashions in the modeling process[27]: A simple model with a complex feature screening process or a complex model with a simple feature processing procedures. The former concentrated on feature engineering, and the obtained variables can be well explained by a constructed simple model; the latter's work focuses on the modeling process, and the resulted model performance is often better. However the results are difficult to explain. Random forest model used in this article is a good compromise. Multiple decision tree integrated learning can be used for complex nonlinear modeling[28], and applied to the case with weak linear correlation between feature variables and target dependent variables. At the same time, the importance or contribution of each feature is evaluated by calculating its corresponding Gini index. This can be used for the future explanation of the features.

3.2 Relationship between structural features and knee symptoms

The determinants of pain and mechanical dysfunction in OA keep unclear, but they are posited to involve a complex interaction between biologic, psychological, and social factors[29-31]. In our study, multiple structural biomarkers demonstrated predictive value of different OA symptoms. Of which, cartilage volume in the middle femur and abnormal shape of medial femur (represented by Cur_5pc_MedFem, Cur_95pc_MedFem) were found the greatest impact on pain, disability and total score. As the major source of disability and reason for hospital visits in patients with knee osteoarthritis (OA), pain has been proven slight correlation between the loss of cartilage volume in previous studies[4,32-33]. But some compartment-specific studies suggested that medial tibiofemoral compartment of the knee had relatively stronger relationship to pain[34-37].

It is biologically and mechanically plausible that medial tibiofemoral compartment is of more predictive value of general than lateral compartment because it transmits more load, and osteochondral changes in medial compartment contribute to increased progression of knee symptoms. Furthermore, compared with features of central medial tibia, our study showed osteochondral features of central medial femur had a much higher association with key knee symptoms such as pain and disability.

Among the medial central femur features, subchondral abnormal morphology played a more significant role in pain than cartilage volume according to the present study, it corresponds to previous studies that subchondral bone features have independent associations with structural progression,symptoms[38-39].

Several reasons could account for this: (1)Subchondral bone contains rich nerve endings and is the likely source of nociception in OA. Whereas cartilage can only generate pain through secondary mechanisms related to subchondral bone exposure, subchondral congestion and synovitis, because cartilage itself is aneural and cannot be the direct factor for the pathogenesis of joint pain[29,31]. (2)Subchondral bone, which comprises calcified cartilage and a thin cortical bone layer, has much more biomechanical function to repair, adapt, and change the shape of the joint than cartilage dose, thus contributes more to the development of osteoarthritis symptoms[31,40].

Pre-surgical WOMAC stiffness score has predictive potential for outcome of knee arthroplasty and patient satisfaction, besides, self-reported stiffness is a complex perception that cannot well represented by clinically observed range of motion[41-43].

In this study, stiffness was found highly affected by subchondral bone shape biomarkers in central lateral femur, differed a lot with other self-reported OA complaints. It is a surprising outcome that in line with no previous study, and there is a need to better understand this phenomenon.

However, several limitations should be noted. First, there are some subjective differences in the symptom scores based on the WOMAC scale. For example, the perception of pain often varies from person to person. The symptom evaluation system should raise more attention; second, the structural features we adopt are all obtained directly from the FNIH project. Most of them are onedimensional and/or two-dimensional grayscale/statistical features with limited image information. Deeper radiomics features need to be further explored in the future; third, the model's prediction results for the knee symptom score of 0 (completely normal) and large values (very poor) are poor. This suggests that there may be a difference in the correlations between these structural characteristics and the symptom. These differences should be issued in future. Preprocess of the features, different nonlinearity of the model and discretization of the WOMAC scores are potential alternatives.

In conclusion, the relationship between the structural features and the symptoms of the knee joint was explored in this study. Constructed RF model can be effectively used to predict and evaluate knee joint symptom, and the selected structural features can be future used as potential biomarkers in the knee joint symptom evaluation and treatment guiding.

Conflict of Interest: No conflict of interest to disclose.