1.

Introduction

A general framework for the application of computer-aided diagnosis (CAD) to the field of diffuse optical tomography (DOT) is presented in this two-part paper. We apply the framework to the diagnosis of rheumatoid arthritis (RA) from frequency-domain DOT (FD-DOT) images of human proximal interphalangeal (PIP) joints. The data set is from a recent clinical study and consists of 219 DOT images of PIP joints (99 joints from 33 subjects with RA and 120 joints from 20 healthy subjects).¹ Absorption (${\mu }_a$) and scattering (${\mu }_s^{\prime }$) images are available for each joint.

In Part 1,² we presented a framework for processing DOT images and extracting heuristic features. The classification strength of each feature is evaluated with Kruskal–Wallis analysis of variance, Dunn’s test, and receiver operating characteristic (ROC) curve analysis. Three important observations are made. First, we observe that features of subjects with RA differ from features of healthy subjects ($p<0.05$). This implies that physiological differences between subjects with RA and healthy subjects can be captured by DOT images. Our second major finding pertains to subjects with RA who do not exhibit effusions, synovitis, or erosion on magnetic resonance imaging (MRI) and ultrasound (US) scans. The DOT images of these subjects are statistically the same as the images of subjects with RA who do exhibit effusions, synovitis, or erosion, suggesting that DOT can potentially detect the onset of RA in these joints before MRI and US could show macroscopic anatomical changes. Our third major finding shows that features from ${\mu }_s^{\prime }$ images allow for more accurate classification of each joint as affected or not affected by RA compared to features from ${\mu }_a$ images.

In this part (Part 2), we present a general framework for classifying DOT images of PIP joints as affected or not affected by RA using machine learning techniques. This approach allows for the use of multiple features in the analysis, thus going beyond the ROC curve analysis of Part 1, which was limited to the evaluation of one feature at a time. Classification of each PIP joint is performed with the best 30 features from Part 1 and five distinct classification algorithms—linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), $k$-nearest-neighbors (KNN), self-organizing maps (SOMs), and support vector machines (SVMs). We report the performance of each algorithm in terms of sensitivity (Se) and specificity (Sp).

The work in this paper goes beyond previously published material^1,3,4 in several ways. First, unlike in previous studies, we combine ${\mu }_a$ and ${\mu }_s^{\prime }$ features in the classification analysis. Second, we are substantially increasing the number of features considered from 4 to 594. In addition to basic features, such as smallest or largest ${\mu }_a$ and ${\mu }_s^{\prime }$ values, we now consider more advanced features, such as Fourier coefficients of two-dimensional (2-D) and 3-D spatial distribution of optical properties. Third, we compare the performance of five different classification algorithms to determine which scheme is most suitable for DOT imaging data. Fourth, we employ a feature-selection algorithm to determine the subset of image features that achieves highest Se and Sp. This step is essential given the large number of permutations possible when 594 features are considered together with five different classification algorithms. Here we employ an evolution strategy that identifies five to ten optimal features out of 594. Finally, we use intracluster correlation coefficients (ICCs) to compute the effective sample size (ESS) of each data group, which is used to adjust the Se and Sp values. This approach helps to account for bias that may arise as a result of treating each imaged finger as an independent sample and allows us to compute confidence intervals (CIs) for Se and Sp; a necessary treatment as our data consists of multiple fingers per subject and these images may not be statistically independent.

In the remainder of this paper we address the general format of multidimensional classification by presenting details on the five classification methods of interest. We review the mathematical and theoretical foundations of the feature-selection algorithm and present results from classification of DOT images of PIP joints as affected or not affected with RA. All classification results are validated through extensive cross-validation. The paper concludes with a discussion on the potential impact of CAD in the diagnosis of RA with DOT.

2.

Methods

Reducing the number of features from the original 594 features to a smaller subset is motivated by various factors. First, features that poorly differentiate between the two diagnostic groups in ROC analysis are unlikely to offer substantial contributions in multidimensional analysis and should be discarded. Second, to improve the generalizability of classification results, it is generally desirable for the ratio between the number of features ($l$) and data samples ($N$) to be small (generally $l/N=0.1$ to 0.2 is acceptable).⁵ Third, the complexity of classification algorithms can increase with $l$, sometimes exponentially. Ultimately, the final number of features should strike a balance between these motivating factors and the desire to include as many important features as possible.

In this work we select 30 features from the original 594 features, as it is a good compromise, resulting in $l/N=0.13$. These 30 features are selected based on Youden indices ($Y=\mathrm{Se}+\mathrm{Sp}-1$) from ROC curve analysis. The features with the 30 largest $Y$ values are selected; of these, 4 are from ${\mu }_a$ and 26 are from ${\mu }_s^{\prime }$ images. Throughout this paper we refer to features by the numbering in Table 1.

Table 1

List of features with top Youden indices from ROC curve analysis (presented in Part 1).

The original clinical data are divided into five groups (labels A, B, C, D, or E) of subjects with RA (segmented based on symptoms) and one group without RA (label H) (refer to Part 1 for the details of the clinical trial data). We established in Part 1 that the subgroups of affected subjects are not statistically different from each other based on the features we currently consider. However, each of the affected subgroups is statistically different from the cohort of healthy subjects. As a result, in the following analysis, we treat all subjects with RA as one group (affected with RA) and attempt to accurately classify an unseen data point as affected or not affected with RA.

The nomenclature used in Table 1 was established in Part 1 and follows the pattern Feature #:Projection Name:Optical Parameter. Shorthand notation is necessary because of the large number of features we considered in Part 1. For example, feature 19 is denoted by F03:GT:s, which translates to the mean (F03 or feature number 3) of the geometrically dominant transverse (GT) slice of the reduced scattering (s) reconstruction. Indices “a” and “s” denote ${\mu }_a$ and ${\mu }_s^{\prime }$ derived features, respectively. Feature numbers F01 to F05 are basic statistical features, F06 to F12 are Gaussian mixture model features, and labeling of FFT features starts from F13 for the first FFT coefficient. For 2-D images, the last FFT coefficient is F26, while for 3-D images it is F76.

As an example and to visualize the reduction of subgroups A to E into a single group, consider features 5 and 19 from Table 1. The allocation of data to subgroups A to E and H is presented in Fig. 1(a), where all six groups are visualized. The same data are presented in Fig. 1(b); however, in this plot, a single group (label RA) replaces subgroups A to E, resulting in only two groups of data. For the purpose of consistency, we use features 5 and 19 throughout this text when it is necessary to show 2-D plots.

Fig. 1

Distribution of sample features. In (a) the six distinct diagnosis groups (A, B, C, D, E, H) are identified, while in (b) the five cohorts diagnosed with RA (A, B, C, D, E) are grouped into one group (RA).

In this work we study the classification performance of various multidimensional combinations of these 30 features, starting with 2-D combinations. In Secs. 2.1 to 2.4 we briefly review the five classification algorithms. The cross-validation methodology is presented in Sec. 2.5. The feature-selection algorithm, which we use to find optimal feature combination, is presented in Sec. 2.6.

2.1.

Nearest-Neighbor Classification

The KNN algorithm is among the most basic classification algorithms because it does not require much, if any, prior knowledge of the distribution of the training or testing data. Each unseen feature vector, $\mathrm{x}$, is classified according to the density of affected or healthy data points within a spatial sphere of radius $r$ (covering $k$ neighboring data points and two distinct data classes $M$).^5,6 The rules governing the assignment of a label to each testing vector $\mathrm{x}$ are as follows:

• 1. From the training data, identify the $k$ nearest neighbors to each vector $\mathrm{x}$ using the Euclidean distance measure.

• 2. Count the number of training data vectors belonging to each class (healthy or RA).

• 3. Assign test vector $\mathrm{x}$ to the class with the maximum number of $k_i$ samples (healthy or RA).

The choice of $k$ affects classification: generally, larger values reduce the effect of noise, but make boundaries between classes less distinct [Fig. 2(a)]. The simplest version of the algorithm is when $k=1$, known as the nearest-neighbor (NN) rule. In other words, a feature vector $\mathrm{x}$ is assigned to the class of its nearest neighbor.

Fig. 2

KNN (a), QDA (b), LDA (c), and SVM (d) boundaries for a nonlinearly separable two-class problem. (a), Examples of $k=1$ and $k=13$ for KNN are shown for a new data point (▴). (d), Support vector data points are denoted by white squares and circles.

2.6.

Selection of Best Features

In our work to date, our strategy for finding the combination of features that yields the best results (Se and Sp) in the classification of RA has been to evaluate the ability to diagnose with all possible feature combinations.^1,3,4 The same analysis cannot be performed for large numbers of features as the number of combinations is $2^n-1-n$, where $n$ is the number of features. In this work $n=594$, and the number of possible combinations is astronomical.

We overcome this problem in two steps. First, we reduce the dimensionality $l$ of the feature space by considering only those 30 features with the largest Youden index ($\mathrm{Se}+\mathrm{Sp}-1$). Second, we employ an optimization algorithm to determine a subset of these 30 features that yield optimal or near-optimal classification results. The algorithm does not test all possible feature combinations; instead, it samples only a subset of these combinations while still achieving high classification accuracy.

In particular, we employ an evolution strategy algorithm generally referred to as (1, $\lambda$)-single-parent evolution strategy (ES) or greedy feature-selection rule. This is an optimization technique based on ideas of adaptation and evolution.^5,10 The ES algorithm determines which multidimensional set of features achieves optimal (or near-optimal) Se and Sp.

Beginning with a set of parent features ($p^k$), the algorithm has two steps: (1) mutation and (2) selection. In the mutation step, a total of $M$ mutants, denoted as ${\lambda }^k$, are generated to compete with their parents $p^k$. The mutation occurs in three steps: (1) $M_a$ new feature combinations are generated by adding a new feature to $p^k$; (2) $M_r$ new features are generated by replacing an existing feature in $p^k$ with one of the remaining features; (3) $M_d$ new feature combinations are generated by dropping an existing feature in $p^k$. Thus, the total number of mutants $M$ ($M_a+M_r+M_d$) in generation $k$ are obtained by adding, dropping, or replacing one feature from parent features $p^k$.

In the selection step, new parents $p^{k+1}$ are selected from the current set of parents $p^k$ and their ${\lambda }^k$ mutants. The selection of the new parent features $p^{k+1}$ is made by selecting the feature combination (with dimensionality $d$) that yields the largest augmented Youden index, defined as

Eq. (9)

$$Y^{*}(\mathrm{Se},\mathrm{Sp})=\mathrm{Se}+\mathrm{Sp}+\alpha L_{\mathrm{Se}}+\beta L_{\mathrm{Sp}}-\delta d-1,$$

where $L_{\mathrm{Se}}$ and $L_{\mathrm{Sp}}$ are the lower bounds of the CI for Se and Sp (i.e., ${\mathrm{CI}}_{\mathrm{Se}}^{\alpha }$ and ${\mathrm{CI}}_{\mathrm{Sp}}^{\alpha }$ from Part 1), respectively. The scaling factors $\alpha$, $\beta$, and $\delta$ control the contribution of the lower bound values ($L_{\mathrm{Se}}$ and $L_{\mathrm{Sp}}$) and dimensionality ($d$) of the selected feature combination on $Y^{*}$, and are all set to 0.001. In this way, feature combinations with higher lower bounds and low dimensionality are preferred.

The mutation and selection operators are applied in a loop and an iteration of the loop is called a generation, denoted by $k$. We begin the process by specifying the initial set of parents $p^0$. The sequence of generations continues until we are unable to improve $Y^{*}$, that is, $Y^{*k}=Y^{*k+1}$. The algorithm finds the feature combination with the highest Se and Sp values, favoring combinations with higher $L_{\mathrm{Se}}$ and $L_{\mathrm{Sp}}$ and lower dimensionality $d$. The process is summarized as follows:

• 1. A feature combination (may be a single feature) is chosen as the parents $p^k$ of the current generation $k$.

• 2. All $M$ possible ${\lambda }^k$ mutants of generation $k$ are obtained by adding ($M_a$), dropping ($M_d$), or replacing ($M_r$) one feature from the parent combination $p^k$ ($M=M_a+M_r+M_d$).

• 3. Of all ${\lambda }^k$ mutants and their $p^k$ parents, the combination with the largest $Y^{*}$ becomes the parents of the next generation ($p^{k+1}$).

• 4. Set $k\leftarrow k+1$ and repeat this process until the objective function $Y^{*}$ does not improve (i.e., $Y^{*k}=Y^{*k+1}$).

This procedure is formulated as an optimization problem, where the objective function, ${\varphi }_{\lambda }^k$, is defined as the winning $Y^{*}$ of generation $k$ [Fig. 3(a)]. In this sense, we seek to maximize $\varphi$ by selecting the feature combination that maximizes $Y^{*}$ at each generation and define it as

Eq. (10)

$${\varphi }_{\lambda }^k\to \underset{i}{\max }{(Y_{\lambda }^{*})}^k=\underset{i}{\max }({\mathrm{Se}}^{k,i}+{\mathrm{Sp}}^{k,i}+\alpha L_{\mathrm{Se}}^{k,i}+\beta L_{\mathrm{Sp}}^{k,i}-\delta d^{k,i}-1),$$

where the index $i$ refers to the $i$th feature combination in generation $k$. ${\mathrm{Se}}^{k,i}$ and ${\mathrm{Sp}}^{k,i}$ are the Se and Sp from the LOOCV procedure using the $i$th feature combination of the $k$th generation. Similarly, $L_{\mathrm{Se}}^{k,i}$ and $L_{\mathrm{Sp}}^{k,i}$ are the lower bounds of the 95.0% CI for ${\mathrm{Se}}^{k,i}$ and ${\mathrm{Sp}}^{k,i}$, while $d^{k,i}$ is the dimensionality of the corresponding feature combination. The algorithm guarantees that ${\varphi }_{\lambda }^{k+1}>{\varphi }_{\lambda }^k$ until convergence, where the solution converges to a near-optimal combination (a local maximum) that maximizes ${\varphi }_{\lambda }$ [Fig. 3(b)]. Se and Sp are not required to increase at each iteration.

Fig. 3

(a), Sample within-generation values of the augmented Youden index $Y^{*}$ for all possible feature combinations (mutants). (b), Sample evolution of objective function over multiple generations.

To identify the appropriate parents for the first generation ($p^0$), we first perform LOOCV using all possible 2-D feature combinations and compute $Y^{*}$ for each combination. Feature pairs that yield the five highest $Y^{*}$ values are selected as the first-generation parents $p^0$. Thus, five distinct optimization runs are executed for each combination of algorithm type (KNN, DA, SOM, SVM) and algorithm parameters (number of neighbors, discriminant type, SVM kernel, etc.), where the initial set of features $p^0$ is different for each run.

3.

Results

We start by showing three examples of typical decision boundaries. Figure 4 shows these boundaries for three different classification algorithms applied to a data set consisting of two features. The two features are the minimum and mean value of ${\mu }_s^{\prime }$ in images of all healthy subjects (blue dots) and subjects with RA (red dots). The classification algorithms that produce the decision boundaries are LDA, QDA, and SVM. In the case of LDA and QDA, all the data points influence placement of the boundary, while the support vectors determine the SVM boundary only (markers with white face). Here, 34 support vectors ($<16\%$ of the data) are identified, 20 from healthy subjects and 14 from subjects with RA. In these three examples, the entire data are used to train (i.e., no cross-validation).

Fig. 4

LDA (a), QDA (b), and SVM (c) decision boundaries separating affected from healthy data using image features 5 and 19. Support vector data points are identified by circles and squares with white face.

From these plots we determine TN, TP, FN, and FP values. In these example cases we find that LDA achieves 91% Se and 91% Sp ($\mathrm{FP}=11$, $\mathrm{FN}=9$), QDA results in 93% Se and 95% Sp ($\mathrm{FP}=6$, $\mathrm{FN}=7$), and SVM achieves 98% Se and 93% Sp ($\mathrm{FP}=9$, $\mathrm{FN}=2$). A similar analysis is performed for all other feature combinations identified by the feature-selection algorithm; however, in those cases, we do perform extensive cross-validation as described in Sec. 2.5.

Results from the evolution algorithm are summarized for each classification method in Secs. 3.1 to 3.4. For brevity, we report only the results obtained with three distinct seed parents used in the optimization algorithms (out of a possible five). The results of nonreported optimization runs are well within the trend presented by the reported cases.

In general, the optimization algorithm always converges to a combination of two to six features in one to six generations and always improves on the initial Se and Sp by 5% to 15% for all classification methods (Fig. 5). The optimal set of features typically includes basic statistical features derived from the raw reconstruction data, basic features from projections of the 3-D data set, as well as coefficients from the Fourier transform of the data (refer to Table 1 for definitions). The algorithm fails to improve on the original set of features only twice (KNN, 1 neighbor, runs 1 and 3).

Fig. 5

Sample Se and Sp evolution paths obtained with (a) KNN, (b) DA, (c) SOM, and (d) SVM classifiers. The lower bounds (Lb) for the 95% confidence interval are shown for Se and Sp, respectively. The upper bound in these cases is 100.0%.

We report the Se and Sp values to which the algorithm converges, the initial and final feature combinations, and the lower and upper bounds of the 95.0% CI for Se and Sp. The following convention is used: Se ($L_{\mathrm{Se}}$,$U_{\mathrm{Se}}$), where $L_{\mathrm{Se}}$ and $U_{\mathrm{Se}}$ are the lower ($L$) and upper ($U$) bounds around the computed Se value within which we have 95.0% confidence that the true value of Se lies. This interval is computed using the ESS and ICC as presented in Part 1 and corrects for correlation between fingers from the same subject. The ESS and ICC values used throughout the analysis are presented in Table 2.

Table 2

ICC and ESS values for each of the 30 features listed in Table 1.

3.5.

Best Feature Selection

The frequency with which each of the original 30 features appears as an optimal classifier is presented in Fig. 6. Features 8, 15, and 19 are chosen as optimal features 73.0%, 40.0%, and 73.0% of the time across all KNN optimization runs [Fig. 6(a)]. Features 5 (100.0%), 15 (80.0%), and 19 (80.0%) are chosen as optimal features most often across all DA iterations [Fig. 6(b)]. Similarly, features 3 (20.0%), 6 (20.0%), 8 (33.0%), and 19 (33.0%) are chosen as optimal features most often across all SOM iterations [Fig. 6(c)]. Finally, features 5 (88.0%), 6 (60.0%), 19 (48.0%), and 30 (48.0%) are chosen as optimal features most often across all SVM iterations [Fig. 6(d)].

Fig. 6

Frequency with which all features appear as optimal classifiers using (a) KNN, (b) DA, (c) SOM, and (d) SVM.

Feature 5 is chosen most often by DA (100.0%) and SVM (88.0%), while feature 19 is chosen most often by KNN (73.0%), and SOM ((33.0%). Feature 19 occurs as an optimal feature at least 20.0% of the time in all classification methods. Feature 15 appears as an optimal feature $>20.0\%$ of the time for three classification methods. Features 3, 5, 6, and 8 appear as an optimal feature $>20.0\%$ of the time for two classification methods. Finally, features 7, 11, 21, and 30 are chosen at least 20.0% of the time by one classification method.

4.

Discussion and Conclusions

In this two-part paper, we present a general framework for the application of CAD techniques to DOT. In Part 1 we focus on feature-extraction methods, while in Part 2, image classification with machine learning techniques takes center stage. As a specific example, the framework is applied to the classification of FD-DOT images of 219 PIP joints ($N=53$ subjects). The goal is to classify each PIP joint as affected or not affected with RA based on the analysis of single image features (Part 1) and combinations of multiple image features (Part 2). For image classification with multiple features, we compare the performances of five different algorithms, including KNN, LDA, QDA, SOM, and SVM. Given the large number of possible permutations of the 594 features we extract, it is necessary to implement a feature-selection algorithm to determine the subset of image features that achieves the highest Se and Sp in combination with each of the five classification algorithms. Results are validated through extensive cross-validation, where the ground truth is the clinical diagnosis made according to the American College of Rheumatology (ACR) 2010 criteria.¹¹

The procedure for training and testing the algorithms is a two-step process using a modified version of the LOOCV procedure. The first phase uses data from $N-1$ subjects to train the classifier and then retrospective classification is performed on the data from the remaining one subject (testing phase). This procedure is repeated so that data from each subject (three images for subjects with RA and six images for subjects without RA) is left out one time each.

We find that all five classification algorithms achieve clinically relevant sensitivities and specificities $>90\%$. The best combination of sensitivity and specificity is obtained with SVM using a polynomial kernel of degree 3. The optimal features corresponding to these results are {4,5,6,12,15,30}. For this case, the sensitivity is 100% with a 95% CI of (96.4%,100.0%). The specificity is 97.8% with a 95% CI of (93.8%,100.0%).

Features 4 and 5 capture the range and mean ${\mu }_a$ and ${\mu }_s^{\prime }$ values in the PIP joint, while features 6, 12, 15, and 30 all capture the variation of ${\mu }_s^{\prime }$ across the joint. Feature 4 corresponds to the ratio of maximum ${\mu }_a$ and minimum ${\mu }_a$ values in the transverse slice across the middle of the PIP joint (F05:GT:a); feature 5 is the minimum value of the ${\mu }_s^{\prime }$ unstructured reconstruction (F02:UV:s); feature 6 is the variance value of the unstructured reconstruction ${\mu }_s^{\prime }$ data (F04:UV:s); feature 12 is the mean value of the variation among all ${\mu }_s^{\prime }$ sagittal slices (F03:VS:s); feature 15 is the variance value of the variation among all ${\mu }_s^{\prime }$ coronal slices (F04:VC:s); and feature 30 is the absolute value of the sixth coefficient of the 2-D-FFT of the central sagittal slices of ${\mu }_s^{\prime }$ (F16:VC:s).

Features that most often achieve accurate classification are associated with global absolute values of the absorption and scattering data and their spatial variation near the PIP joint. We see evidence that features that quantify spatial variation across PIP joints are smaller for subjects with RA compared to healthy subjects. This is in concordance with our earlier findings and the images shown in Fig. 2 of Part 2. The synovial fluid as well as the surrounding tissue experience changes in optical properties in subjects with RA.^12,13 The inflammatory process starts in the synovium, leading to changes in the cell and tissue structure. Cell proliferation can be observed and the appearance of the synovial fluid changes from a clear, yellowish substance to a turbid, gray-yellowish substance. The number of leukocytes per mL increases from 100–200 in healthy conditions to 1000–100,000 during stages 1 and 2 of the disease. Leukocytes have a diameter of approximately 7–20 μm and therefore have an effect on the scattering coefficient. Furthermore, the protein content in the synovial fluid approximately triples from 10–$20\mathrm{g}/\mathrm{L}$ to 30–$60\mathrm{g}/\mathrm{L}$.^14,15 In addition, neovascularization in the surrounding tissue has been related to synovitis,¹⁶ which leads to an increase in the absorption coefficient.

Overall, the net effect of these changes is an increase in absorption and scattering in the finger joint affected by RA, resulting in an optical profile similar to the physiology surrounding the joint. Thus the spatial variation of optical properties decreases in subjects with RA. Although similar results have been described before, it appears that our CAD algorithms can extract these features with higher accuracy and produce sensitivities and specificities at levels not reported before. All five classification algorithms achieve clinically relevant sensitivities and specificities $>90.0\%$, some even $>97.0\%$. The computed 95.0% CIs for our results offer further validation that our results are robust. These observations warrant further multicenter prospective clinical studies.

One particular focus of these trials should further investigate an intriguing finding in our current study. We observe that joints of subjects with RA but without radiological evidence of effusion, erosion, or synovitis (as determined by MRI or US) are statistically identical to joints of subjects with RA that have detectable effusion, erosion, or synovitis; indeed, these joints are also statistically different from joints of healthy subjects. These joints are clear examples of cases where DOT imaging is sensitive to symptoms associated with RA before they become detectable by MRI or US. A longitudinal study is necessary to determine if these joints eventually evolve to exhibit evidence of effusion, erosion, or synovitis on MRI or US scans so as to meet the ACR criteria for RA. If, indeed, that turns out to be the case, one could prove that DOT can detect the presence of symptoms associated with RA at an earlier stage than other imaging methods.

In addition to diagnosis of RA, we expect that the framework established in this work (CAD with DOT) can be deployed in the diagnosis and monitoring of other diseases, including breast cancer and peripheral artery disease.