Multifaceted radiomics for distant metastasis prediction in head & neck cancer

Before training and testing, a deep learning with stacked autoencoder strategy is designed to transform the features from PET, CT features as well as clinical parameters into one representation feature set in M-radiomics. We aim to generate a comprehensive Pareto-optimal model set through sorting classifier specific Pareto-optimal models in a non-dominated way in the training stage. The goal of testing stage is to fuse the non-zero weight models through analytic ER rule and get the final output probability. The problem formulation is described as follows.

The training stage of M-radiomics is shown in figure 2. Assume that the radiomic features and clinical parameters are denoted by ${T^{\,i}} = \left\{ {T_1^{\,i},\,T_2^{\,i},\, \cdots ,T_{{M_i}}^{\,i}} \right\},\,i = 1, \cdots ,I$, where M_i represents the feature number for each modality and I is the modality number. Before training, an autoencoder model is generated using all the features from the different modalities, and the hidden layer is taken as the final feature set which is fed into the model, that is:

$$\begin{equation}{T_{train}} = SSAE\left( {{T^{\,i}}} \right),\,i = 1, \cdots ,I,\end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

where $SSAE$ represents stacked autoencoder and ${T^{\,i}}$ is the feature set generated from the hidden layer in training stage. Assume that there are N base classifiers and the model parameter set for each classifier is denoted by ${\beta _n},\,n = 1, \cdots ,N$. Since the feature selection influences the model performance, we perform feature selection and model training simultaneously for each base classifier which is:

$$\begin{equation}{f^{\,\,n}} = \mathop {\max}\limits_{T,\,{\beta _n}} \left( {f_{sim\_sen}^{\,\,n},\,f_{sim\_spe}^{\,\,n}} \right),\,n = 1, \cdots ,N,\,\end{equation} \tag{ 2 }$$

where ${f^{\,\,n}}$ is the objective functions, and $f_{sim\_sen}^{\,\,n}$, and $f_{sim\_spe}^{\,\,n}$ represent similarity-based sensitivity and specificity, respectively. After optimizing by an iterative multi-objective immune algorithm (IMIA), the classifier-specific Pareto-optimal model set ${p_n},n = 1, \cdots ,N$ is obtained and the final Pareto-optimal model set $D$ is generated as:

$$\begin{equation}D = NS\left( {{p_n}} \right),\,n = 1,, \cdots ,N,\end{equation} \tag{ 3 }$$

where $NS$ represents the non-dominated sorting strategy (Deb et al 2002). The NS strategy is used for sorting all possible solutions in the population and updating the current population to generate a new population with the same size as the initial setting.

Before performing the test as shown in figure 3, the feature set denoted by ${T_{test}}$ is generated through $SSAE$. Since we hope to obtain the balanced results, models that achieve good balance between sensitivity and specificity are remained, while the other models are pruned. Assume there are M models with non-zero weight, and probability output is denoted by ${p_m},m = 1, \cdots ,M$. The final output P can be obtained by fusing the M models with weight ${w_m},m = 1, \cdots ,M$ and reliability ${r_m},m = 1, \cdots ,M$ through ERRF, that is:

$$\begin{equation}P = ERRF\left( {{p_m},{w_m},{r_m}} \right),\,m = 1, \cdots ,M.\end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

The features from different modalities are jointed represented through SSAE in M-radiomics. An illustration of two hidden layer stacked sparse autoencoder is shown in figure 4. SSAE (Suk et al 2015, Pan et al 2016, Zhou et al 2019a) consists of a single layer sparse autoencoder (AE) based on a neural network, which consists of an encoder and a decoder. The input vector is fed into a concise representation through connections between input and hidden layers in encoding step, while the decoding step is trying to reconstruct the input vector precisely through encoded feature representation.

Zoom In Zoom Out Reset image size

Assume that x is an input vector with the D dimension containing the radiomic features and clinical parameters in our study. Then the encoder is performed by mapping x to another vector denoted by z, that is:

$$\begin{equation}{z_1} = {h_1}\left( {{W_1}x + {b_1}} \right),\end{equation} \tag{ 5 }$$

where h₁ is a transfer function for encoder, W₁ is a weight matrix and b₁ is a bias vector in the first layer. The decoder is to map the encoded representation z₁ back into an estimate of x, that is:

$$\begin{equation}\hat x = {h_2}\left( {{W_2}x + {b_2}} \right),\end{equation} \tag{ 6 }$$

where ${h_2}$ is a transfer function for decoder, ${W_2}$ is a weight matrix and ${b_2}$ is a bias vector in the second layer. The average output activation measure ${\hat \rho _i}$ of $ith$ input vector in training set is calculated as:

$$\begin{equation}{\hat \rho _i} = \frac{1}{n}\mathop \sum \nolimits_{j = 1}^{\,\,n} {z_{i,1}}\left( {{x_j}} \right) = \frac{1}{n}\mathop \sum \nolimits_{j = 1}^{\,\,n} h\left( {W_{i,1}^T{x_j} + {b_{i,1}}} \right),\end{equation} \tag{ 7 }$$

where $n$ is the number of training samples and ${x_j}$ represents $jth$ training sample. $W_{i,1}^T$ is the $ith$ row of the weight matrix ${W_1}$ and ${b_{i,1}}$ is the $ith$ entry of the bias vector. Then the sparsity denoted by ${{{\Omega }}_{sparsity}}$ is calculated by measuring the dissimilarity between ${\hat \rho _i}$ and its desired value denoted by $\rho$ through Kullback–Leibler divergence (KL), that is:

$$\begin{equation}{{{\Omega }}_{sparsity}} = \mathop \sum \nolimits_{i = 1}^D KL(\rho ||{\hat \rho _i}) = \mathop \sum \nolimits_{i = 1}^D \rho {\text{log}}\left( {\frac{\rho }{{{{\hat \rho }_i}}}} \right) + \left( {1 - \rho } \right)\log \left( {\frac{{1 - \rho }}{{1 - {{\hat \rho }_i}}}} \right).\end{equation} \tag{ 8 }$$

By minimizing the cost function, ${\hat \rho _i}$ and $\rho$ will be closer to each other. Additionally, it is possible that the sparsity regularizer become small with increased ${W_1}$ and decreased ${z_1}$. To avoid this issue, ${L_2}$ regularization is added, that is:

$$\begin{equation}{{{\Omega }}_{weight}} = \frac{1}{2}\mathop \sum \nolimits_l^{\,\,l} \mathop \sum \nolimits_j^{\,\,n} \mathop \sum \nolimits_i^{\,\,k} {\left( {w_{j,i}^{\,\,l}} \right)^2},\end{equation} \tag{ 9 }$$

where $L$ is the number of hidden layer, $n$ is the number of training samples and $k$ is the feature number. Therefore, by adjusting the mean squared error, the cost function denoted by $E$ for the training sparse autoencoder is:

$$\begin{equation}E = \frac{1}{N}\mathop \sum \nolimits_{n = 1}^{\,\,k} \mathop \sum \nolimits_{k = 1}^{\,\,k} {\left( {{x_{kn}} - {{\hat x}_{kn}}} \right)^2} + \lambda^{\ast}{{{\Omega }}_{weight}} + \beta^{\ast}{{{\Omega }}_{sparsity}},\end{equation} \tag{ 10 }$$

where $\lambda$ is the coefficient for ${L_2}$ regularization and $\beta$ is the coefficient for sparse penalty. The objective for sparse AE is to minimize the reconstruction error by optimizing above function. The SSAE is a deep learning architecture of sparse AE which can be constructed by stacking another hidden layer. Then the second hidden layer is taken as the final feature set that we will use in the next step.

As the feature selection influences the model training, feature selection and model parameter training are performed simultaneously (Zhou et al 2017). Assume we $K$ training samples with the label vector $T = \left[ {{T_1},{T_2}} \right]$ ([1,0] for patients without distant metastasis and [0,1] for patients with distant metastasis). The model parameter set is denoted by $MP$ and the feature set is denoted by $FS$. In our previous model (Zhou et al 2017), sensitivity denoted by ${f_{sen}}$ and specificity denoted by ${f_{spe}}$ are used, and calculated them as follows:

$$\begin{equation}\,{f_{sen}} = \frac{{TP}}{{TP + FN\,}},\,\,{f_{spe}} = \frac{{TN}}{{TN + FP\,}},\end{equation} \tag{ 11 }$$

where $TP$ is the number of true positives, $TN$ is the number of true negatives, $FP$ is the number of false positives, and $FN$ is the number of false negatives.

However, the above objective functions are calculated based on the label, which may lead to less reliable results (Chen et al 2018) as the predicted output probability could be further away from the true label vector. For example, we consider a model with output of [0.9, 0.1] is more reliable than a model with output of [0.55, 0.45] as the output of the former model is closer to the ideal label vector [1,0]. In the label-based multi-objective model, sensitivity and specificity are calculated by setting the probability threshold at 0.5, which does not push the solution of output probability vector to the ideal label vector [0,1] or [1,0]. When we design a probability-based objective function without setting any threshold, we can obtain more reliable results by pushing the output probability to the ideal vector.

If we consider the label vector as a probability distribution, the problem is converted as measuring the similarity between predicted probability output and the label vector. Since Jaccard distance (Cha 2007) is an effective probability density distribution measurement, it is adopted here, that is:

$$\begin{equation}sim\left( {{P^{\,\,k}},{T^{\,\,k}}} \right) = \frac{{\mathop \sum \nolimits_{i = 1}^2 P_i^{\,\,k}T_i^{\,\,k}}}{{\mathop \sum \nolimits_{i = 1}^{\,\,n} {{\left( {P_i^{\,\,k}} \right)}^2} + \mathop \sum \nolimits_{i = 1}^{\,\,n} {{\left( {T_i^{\,\,k}} \right)}^2} - \mathop \sum \nolimits_{i = 1}^{\,\,n} {{\left( {P_i^{\,\,k}T_i^{\,\,k}} \right)}^2}}},\,k = 1, \cdots ,K,\end{equation} \tag{ 12 }$$

where $sim$ represents the similarity measurement. The objective function becomes:

$$\begin{equation}\mathop {f = {\text{max}}}\limits_{MP,FS\!\!\!\!\!\!\!\!\!\!\!} \mathop \sum \nolimits_{k = 1}^{\,\,k} {\text{sim}}\left( {{P^{\,\,k}},{T^{\,\,k}}} \right),\end{equation} \tag{ 13 }$$

Then we define the similarity-based sensitivity and specificity, which are denoted by ${f_{sim\_sen}}$, ${f_{sim\_spe}}$, respectively. Assume that ${P_{tp}} = \left\{ {P_{tp}^1,P_{tp}^2,\, \cdots ,P_{tp}^{TP}} \right\}$ represents the probability output of true positives and the corresponding labelled vector is ${T_{tp}} = \left\{ {T_{tp}^{\,1},T_{tp}^{\,2},\, \cdots ,T_{tp}^{TP}} \right\}$. The similarity measure of true positives $T{P_{sim}}$ is:

$$\begin{equation}T{P_{{\text{sim}}}} = \mathop {\sum\nolimits_{k = 1}^{TP} {{\text{sim}}\left( {P_{tp}^{\,\,k},T_{tp}^{\,\,k}} \right)} }\nolimits_{}^{} = \sum\nolimits_{k = 1}^{TP} {\frac{{\mathop \sum \nolimits_{i = 1}^2 P_{tp,j}^{\,\,k}T_{tp,j}^{\,\,k}}}{{\mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{tp,i}^{\,\,k}} \right)}^2} + \mathop \sum \nolimits_{i = 1}^2 {{\left( {T_{tp,i}^{\,\,k}} \right)}^2} - \mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{tp,i}^{\,\,k}T_{tp,i}^{\,\,k}} \right)}^2}}}} .\end{equation} \tag{ 14 }$$

Similarly, we can get the similarity measure of true negatives, false positives and false negatives, which are denoted by $T{N_{sim}}$, $F{P_{sim}}$ and $F{N_{sim}}$, respectively. They are:

$$\begin{equation}T{N_{sim}} = \mathop {\sum\nolimits_{k = 1}^{TN} {\frac{{\mathop \sum \nolimits_{i = 1}^2 P_{tn,j}^{\,\,k}T_{tn,j}^{\,\,k}}}{{\mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{tn,i}^{\,\,k}} \right)}^2} + \mathop \sum \nolimits_{i = 1}^2 {{\left( {T_{tn,i}^{\,\,k}} \right)}^2} - \mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{tn,i}^{\,\,k}T_{tn,i}^{\,\,k}} \right)}^2}}}} }\nolimits_{}^{} ,\end{equation} \tag{ 15 }$$

$$\begin{equation}F{P_{sim}} = \sum\nolimits_{k = 1}^{FP} {\frac{{\mathop \sum \nolimits_{i = 1}^2 P_{fp,j}^{\,\,k}T_{fp,j}^{\,\,k}}}{{\mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{fp,i}^{\,\,k}} \right)}^2} + \mathop \sum \nolimits_{i = 1}^2 {{\left( {T_{fp,i}^{\,\,k}} \right)}^2} - \mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{fp,i}^{\,\,k}T_{fp,i}^{\,\,k}} \right)}^2}}}} ,\end{equation} \tag{ 16 }$$

$$\begin{equation}F{N_{sim}} = \mathop \sum \nolimits_{k = 1}^{FN} \frac{{\mathop \sum \nolimits_{i = 1}^2 P_{fn,j}^{\,\,k}T_{fn,j}^{\,\,k}}}{{\mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{fn,i}^{\,\,k}} \right)}^2} + \mathop \sum \nolimits_{i = 1}^2 {{\left( {T_{fn,i}^{\,\,k}} \right)}^2} - \mathop \sum \nolimits_{i = 1}^2 {{\left( {P_{fn,i}^{\,\,k}T_{fn,i}^{\,\,k}} \right)}^2}}}.\end{equation} \tag{ 17 }$$

Similar to the sensitivity and specificity, we can also calculate ${f_{sim\_sen}}$, ${f_{sim\_spe}}$:

$$\begin{equation}\,{f_{sim\_sen}} = \frac{{T{P_{sim}}}}{{T{P_{sim}} + F{N_{sim}}\,}},\,\,{f_{sim\_spe}} = \frac{{T{N_{sim}}}}{{T{N_{sim}} + F{P_{sim}}\,}}.\end{equation} \tag{ 18 }$$

Finally, the probability-based objective function is defined as:

$$\begin{equation}f = {\max _{MP,FS}}\left( {{f_{sim\_sen}},{f_{sim\_spe}}} \right).\end{equation} \tag{ 19 }$$

Iterative multi-objective immune algorithm (IMIA) is a multi-objective evolutionary algorithm which has shown improved performance in obtaining solutions with more balanced sensitivity and specificity as compared those single-objective bases approaches (Zhou et al 2017). However, IMIA can only handle the problem with one base classifier, while M-radiomics contains multiple base classifiers. Therefore, an enhanced version of IMIA termed as IMIA-II is developed in this work. The procedure of IMIA-II is shown as follows and a brief workflow in summarized in table 1.

Table 1. Procedures of IMIA-II.

Input: Initial population, set $j = 0.$

while $j > G$,

Step 1: Proportional clonal operation.

Step 2: Static mutation operation.

Step 3: Deleting operation.

Step 4: AUC based fast nondominated sorting approach for updating solution set.

End while

Step 5: Classifier-specific Pareto-optimal model combination and population updating.

Output: Pareto-optimal model set.

Step 1—Initialization: Model parameters are optimized directly as the value is continuous. We use ${G_{max}}$ to denote the maximal generation and $D\left(\, j \right) = \left\{ {{d_1}, \cdots ,{d_P}} \right\}$ denote the population, with $P$ as the number of individuals in one population and $d$ is the solution in current generation. For each base classifier, Step 2- step 5 are performed until the algorithm achieves ${G_{max}}$.

Step 2—Clonal operation: Proportional cloning is used (Gong et al 2008). To diversify the population, an individual with a larger crowding-distance is reproduced more times, and the clonal time ${q_i}$ for each individual is calculated as: ${q_i} = \left\lceil{n_c} \times \frac{{\delta \left( {{d_i},D} \right)}}{{\mathop \sum \nolimits_{j = 1}^{\left| A \right|} \delta \left( {{d_j},D} \right)}}\right\rceil,$ where ${n_c}$ is the expectant value of the clonal population and $\lceil\ \rceil$ is the ceiling operator. $\delta \left( {{d_i},D} \right)$ represents the crowding distance.

Step 3—Mutation operation: The mutation operation (Hang et al 2010) will be performed on the cloned population$C\left(\, j \right)$. For each locus in the individual, a random mutation probability ($M{P_i}$) will be generated. If $M{P_i}$ is larger than a general mutation probability ($GM{P_i}$), the mutation will occur. The mutated population is denoted by $M\left(\, j \right)$. $D\left(\, j \right)$ and $M\left(\, j \right)$ are combined to form a new population$F\left(\, j \right)$.

Step 4—Deleting operation: If there are duplicated solutions in the new population $F\left(\, j \right)$, we will only keep the unique one and delete other duplicated solutions. If the size of $F\left(\, j \right)$ is less than P, step 2 should be used to generate more mutated individuals; otherwise, step 5 should be used. When the same solutions are generated in the Pareto-optimal set, the diversity of the individuals in a population will be reduced. We perform the deleting operation to ensure that all the solutions in the Pareto-optimal set are different after executing the clonal and mutation operations.

Step 5—Updating population: The selected features and model parameters for each individual are taken as the input for model (Yu et al 2005). Then, ${f_{sen\_sim}}$ and ${f_{spe\_sim}}$ are calculated through cross-validation. The individual in all feasible solutions is sorted in the descending order using a fast non-dominated sorting approach (Deb et al 2002) according to the AUC of each solution. We will obtain the Pareto-optimal set according to AUC because it is one of the most important criteria for the performance evaluation of a predictive model.

Step 6—Termination detection. When the algorithm reaches ${G_{max}}$, go to the next step; otherwise, go to step 2.

Step 7—Classifier-specific Pareto-optimal model combination. The final Pareto-optimal model set is generated by combining all the classifier-specific Pareto-optimal models through a AUC-based non-domainted sorting method.

Compared with traditional radiomics, more base classifiers are integrated into M-radiomics. Therefore, the main difference between these two algorithms is that IMIA-II can optimize a model with multiple base classifiers, while IMIA can only train a model with single classifier. Moreover, as each base classifier generates a Pareto-optimal model set, an additional step (Step 5 in table 1) for combining classifier-specific Pareto-optimal models is introduced into IMIA-II. As such, IMIA-II can be considered as a generalized optimization algorithm which can train the model with one or more classifiers.

Before fusing the probabilities of all the individual models through ERRF, reliability and weight are calculated. The fusion is performed through evidential reasoning (ER) rule (Yang and Xu 2013) which is a generic probabilistic inference tool and takes both relative weight and reliability into account. ER rule is a generalization of the ER approach (Yang and Xu 2002) that was originally developed for multiple criteria decision analysis based on the D-S theory (Shafer 1976). It has been applied to multiple group decision analysis (Fu et al 2015), safety assessment (Zhao et al 2016), and data classification (Xu et al 2017). However, the current ER rule is a recursive algorithm, which is difficult to perform fusion directly. Therefore, we theoretically inferred a new analytic ER rule with an analytic expression in this work.

Assume that there are $N$ Pareto-optimal models denoted by $C = \left\{ {{C_1},\, \cdots ,{C_N}} \right\}$ with $M$ classes, where ${P_i} = \left\{ {p_i^1,\,p_i^2, \cdots ,p_i^M} \right\},\,i = 1, \cdots N$ is the output probability for a test sample $x$. In our case study, $M$ equals 2. The reliability and weight for each individual model are denoted by ${r_i}$ and ${w_i},\,i = 1, \cdots ,N$, and the validation set is denoted by $V$. Reliability is defined as the similarity between the output probability of the test sample $x$ and output probabilities of K nearest samples in validation set $V$, that is:

$$\begin{equation}r_i=\left\{ {\begin{array}{lr} 0&{{{\rm when}}\,\,{l_i} \ne {l_j};j = 1, \cdots ,K} \\ 1&{{{\rm when}}\,\,{l_i} = {l_j} \wedge {p_{{l_j}}} = 1;\,j = 1, \cdots ,K} \\ {0 < {r_i} < 1}&\quad{{{\rm in\,\, other\,\, situations}}} \end{array}} \right.,\end{equation} \tag{ 20 }$$

where ${l_i}$ is the output labels of test sample $x$. ${l_j}$ are the output labels of K nearest neighbor samples in $V$, and ${p_{{l_j}}}$ are the corresponding output probabilities.

To better understand the reliability, we use decision making by group experts as an example. When we evaluate the reliability of an expert, a reasonable solution is that we can find several experts who have the similar background with this expert; and the reliability can be evaluated by comparing the decision result with all of other experts. In this study, as the nearest samples (e.g. K samples) of current test sample has the similar background, the reliability of each model can be evaluated through the output probabilities of these samples. Since the output probability for model ${C_i}$ is the probability distribution, calculating the reliability of ${C_i}$ for test sample $x$ is transformed into measuring the similarity of probability distribution. Dice distance is a commonly used probability distribution measurement (Cha 2007), which is adopted here. We first define the dissimilarity among different classifiers. Assume that output probability of test sample $x$ is denoted by $p_i^m$and output probability of one validation sample ${v_k},\,k = 1, \cdots ,K$ in K nearest neighbors is denoted by $p_{{i_k}}^m$. The dissimilarity measure denoted by ${D_{i,k}}\left( {x,\,{v_k}} \right)$ between test sample $x$ and validation sample ${v_k}$ is calculated as:

$$\begin{equation}{D_{i,k}}\left( {x,\,{v_k}} \right) = \frac{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m - p_{{i_k}}^m} \right)}^2}}}{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {p_{{i_k}}^m} \right)}^2}}},i = 1, \cdots ,N;k = 1, \cdots ,K,\end{equation} \tag{ 21 }$$

where $i$ represents $ith$ model and $m$ represents the $mth$ class. $k$ is the $kth$ validation sample in in K nearest validation samples. As described in equation (20), when ${l_i} = {l_j} \wedge {p_{{l_j}}} = 1;\,j = 1, \cdots ,N,\,j \ne i$, ${r_i} = 1$. Therefore, ${D_{i,k}}\left( {x,\,{v_k}} \right)$is changed as:

$$\begin{equation}{D_{i,k}}\left( {x,\,{v_k}} \right) = \frac{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m - p_{{i_k}}^m} \right)}^2}*\mathop \sum \nolimits_{m = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}}{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {p_{{i_k}}^m} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}},\end{equation} \tag{ 22 }$$

where $p_i^{\,\,l}$ is the maximal output probability for test sample $x$. Then ${D_{i,k}}\left( {x,\,{v_k}} \right)$ is simplified as:

$$\begin{equation}{D_{i,k}}\left( {x,\,{v_k}} \right) = \frac{{\left( {\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m - p_{{i_k}}^m} \right)}^2} + 1} \right)*\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}}{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {p_{{i_k}}^m} \right)}^2} + 2\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}},\end{equation} \tag{ 23 }$$

The similarity is then calculated as:

$$\begin{equation}{S_{i,k}}\left( {x,\,{v_k}} \right) = 1 - \frac{{\left( {\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m - p_{{i_k}}^m} \right)}^2} + 1} \right)*\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}}{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {p_{{i_k}}^m} \right)}^2} + 2\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}},\end{equation} \tag{ 24 }$$

where ${S_{i,k}}\left( {x,\,{v_k}} \right)$ represents the similarity measurement. As there are $K$ validation samples, the overall similarity denoted by ${S_i}\left( x \right)$ is:

$$\begin{equation}{S_i}\left( x \right) = \frac{1}{K}\mathop \sum \nolimits_{k = 1}^{\,\,k} \left( {1 - \frac{{\left( {\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m - p_{{i_k}}^m} \right)}^2} + 1} \right)*\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}}{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {p_{{i_k}}^m} \right)}^2} + 2\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}}} \right).\end{equation} \tag{ 25 }$$

Assume that $NV$ ($NV \leqslant K$) is the number of validation samples of which the output labels are the same as the test sample. Then we add $\frac{{NV}}{K}$ into ${r_i}\left( x \right)$ so that it can satisfies all the constraints defined in equation (20). The reliability is then calculated as:

$$\begin{equation}{r_i}\left( x \right) = \frac{{NV}}{{{K^2}}}\mathop \sum \nolimits_{k = 1}^{\,\,k} \left( {1 - \frac{{\left( {\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m - p_{{i_k}}^m} \right)}^2} + 1} \right)*\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}}{{\mathop \sum \nolimits_{m = 1}^M {{\left( {p_i^m} \right)}^2} + \mathop \sum \nolimits_{m = 1}^M {{\left( {p_{{i_k}}^m} \right)}^2} + 2\mathop \sum \nolimits_{j = 1}^M {{\left( {1 - p_i^{\,\,l}} \right)}^2}}}} \right),i = 1, \cdots ,N.\end{equation} \tag{ 26 }$$

As predictive models are desired to have balanced sensitivity and specificity, the weight for the model with extremely imbalanced sensitivity and specificity is set as zero, whilst the weight for the remaining Pareto-optimal models are set as non-zero as shown in figure 5. As AUC is a good measurement for evaluating the model performance, it is considered for weight calculation as well. In summary, by considering the above two factors, the weight is calculated as:

$$\begin{equation}{w_i} = \left\{ {\begin{array}{*{20}{c}} {\frac{{\frac{{f_{sim\_sen}^{\,i}}}{{f_{sim\_spe}^{\,i}}} + AU{C_i}}}{2}}&{{\text{when}}\,\,\,\,0.5 \leqslant \frac{{f_{sim\_sen}^{\,\,l}}}{{f_{sim\_spe}^{\,\,l}}} \leqslant 1} \\ {\frac{{\frac{{f_{sim\_spe}^{\,i}}}{{f_{sim\_spe}^{\,i}}} + AU{C_l}}}{2}}&{{\text{when\,}}\,\,\,0.5 \leqslant \frac{{f_{sim\_spe}^{\,\,l}}}{{f_{sim\_sen}^{\,\,l}}} \leqslant 1} \\ 0&{{\text{Other\,situation}}} \end{array}.} \right.\end{equation} \tag{ 27 }$$

Zoom In Zoom Out Reset image size

where $f_{sim\_sen}^{\,\,i}$, $f_{sim\_spe}^{\,\,i}$ and $AU{C_i}$ represent the similarity-based sensitivity and specificity, and AUC for model $l$ in training stage.

The final output probability denoted by ${p^m}\left( x \right)$ for test sample $x$ is calculated by fusing the output probabilities of all the Pareto-optimal models through analytic ER rule, that is:

$$\begin{equation}{p^m}\left( x \right)\, = AER\left( {{P_i},{w_i},{r_i}} \right),\,i = 1, \cdots ,N,\end{equation} \tag{ 28 }$$

where $AER$ represents analytic ER rule. The details of ER rule and the inference procedure of analytic ER rule are described in appendix A and B, respectively. After introducing analytic ER rule expression, ${p^m}\left( x \right)$ is calculated as:

$$\begin{equation}{p^m}\left( x \right) = \frac{{k\left[ {\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{{w_i}p_i^m}}{{1 + {w_i} - {r_i}}} + \frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right) - \mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)} \right]}}{{1 - k\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)}},\,m = 1, \cdots ,M,\end{equation} \tag{ 29 }$$

where the normalization factor $k$ is:

$$\begin{equation}k = {\left[ {\mathop \sum \nolimits_{m = 1}^M \left( {\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{{w_i}p_i^m}}{{1 + {w_i} - {r_i}}} + \frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)} \right) - \left( {M - 1} \right)\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)} \right]^{ - 1}}.\end{equation} \tag{ 30 }$$

In this study, FDG-PET and radiotherapy planning CT (figure 6) from 188 patients are used. All these patients had pre-treatment FDG-PET/CT scans between April 2006 and November 2017, which were downloaded from the cancer imaging archive (TCIA) (Vallières et al 2017). The follow up time ranges from 6 months to 112 months and the median follow up are about 43 months. Sixteen percent (16%) of these patients had distant metastasis. We also used clinical parameters extracted from clinical charts, including gender, age, T stage, N stage, M stage, TNM group stage. There are 257 features extracted from PET and CT images, respectively, including intensity, texture and geometry features.

Zoom In Zoom Out Reset image size

The hidden layer size in SSAE is set as 100 based on the recommendation from reference (Guo et al 2016). Three base classifiers including radial basis function-based support vector machine (SVM), decision tree (DT) and K-nearest neighbor (KNN) are used. When training the model using IMIA-II, the population size is set as 100 and generation time is 50. The mutation rate is set as 0.9. The validation sample number K is set as 10 for reliability calculation. Five group experiments are designed as follows. (1) We compared M-radiomics with and without SSAE to see how much benefit SSAE brings. Same as most radiomics strategies, the features from different modalities are directly combined into a single vector and are fed into the model in the comparative model (without SSAE). (2) To show the improved performance when adopting multiple base classifiers, each base classifier was used for building a model and all the other setting was same as the M-radiomics. (3) We compared IMIA-II with IMIA to test the reliability of the new training algorithm. All the experimental settings were the same for these two optimization algorithms. (4) We also evaluated the reliability of ERRF strategy. Evidential reasoning-based fusion (ERF) (Chen et al 2019) and weighted fusion (WF) are used for comparison in this experiment. (5) To evaluate the overall performance, we compared M-radiomics with multi-objective (MO) radiomics model (Zhou et al 2017) and single-objective (SO) model. In MO, all the features from three modalities are concatenated into a long vector and one optimal model is selected from Pareto-optimal model. Meanwhile, SVM is used as base classifier in MO.

Five-fold cross-validation with 60% training, 20% validation and 20% testing was performed for all the experiments. AUC, accuracy (ACC), sensitivity (SEN) and specificity are used for model evaluation. Note that SEN and SPE are calculated when the cutoff value for probability output is set as 0.5. Besides these criteria, we also use similarity-based sensitivity (SIM-SEN) and specificity (SIM-SPE) to evaluate the model reliability.

Figure 7 shows the comparison results for M-radiomics with and without SSAE. When introducing the SSAE, M-radiomics can obtain better performance in SEN, SPE, ACC, AUC and SIM-SEN. The improved performance can be attributed to that SSAE extracts more discriminative high-level features and discovers the joint information among different modalities. The difference between SEN and SPE on the model with SSAE (0.06) is lower than the model without SSAE (0.15). Compared with the improvement on SPE and SIM-SPE, SEN and SIM-SEN are greatly improved further for the model with SSAE. These results indicate that SSAE can capture more discriminative feature such that more samples are correctly predicted for both positive and negative cases. As the dataset used in this study contains much more negative cases (84%) than positive cases (16%), greater improvement is observed in SEN and SIM-SEN than SPE and SIM-SPE.

Zoom In Zoom Out Reset image size

In this experiment, we compared the single base classifier model with M-radiomics and the results are summarized in table 2. The performance of different models is compared with the paired t-test at a significance level of 0.05. The results show that there is a statistically significant difference between M-radiomics and other single-classifier-based models. Although SPE for DT is higher than M-radiomics, SEN for DT is much lower than M-radiomics. Meanwhile, SEN and SPE are highly imbalanced in DT and KNN. SVM can obtain more balanced results, but M-radiomics obtains better performance than SVM. The overall performance of M-radiomics outperforms single classifier model. The reason is that the performance of different base classifiers varies greatly for test samples with different distributions. When fusing multiple base classifiers in M-radiomics, the model diversity is greatly increased and can handle the test data with different distributions better.

This experiment is to compare the model performance after optimizing through IMIA and IMIA-II, and the results are shown in figure 8. Compared with IMIA, the model trained by IMIA-II can obtain better performance, which shows that our newly defined probability-based objective function can get a more reliable model. Specifically, IMIA-II can obtain higher AUC and SIM-SEN, which suggests that the model trained by IMIA-II is more reliable.

Zoom In Zoom Out Reset image size

We compare the model reliability when using different fusion strategies and the results are shown in figure 9. ERRF can obtain best results in six evaluation criteria among three strategies. Compared with another two strategies, the main difference of ERRF is it introduces the model reliability during the fusion process, which makes the model more reliable. Meanwhile, both ERF and ERRF outperform WF. ER is a special case of ER rule (Yang and Xu 2013), which demonstrates that ER rule is an outstanding probabilistic inference tool for information fusion.

Zoom In Zoom Out Reset image size

The overall performance of M-radiomics is evaluated in this experiment and radiomics models with multi-objective (MO) and single objective (SO) are used for comparison. The ROC curves for SO, MO and M-radiomics on one running time are shown in figure 10. M-radiomics outperforms both SO and MO models. The average AUC for MO, SO and M-radiomics are 0.76, 0.81 and 0.84, respectively. These results demonstrate that the newly developed strategies in M-radiomics can overcome the challenges in commonly used radiomics approaches and the overall performance is improved.

Zoom In Zoom Out Reset image size

As there's no local ignorance in our problem, they are not considered in ER rule. Therefore, the BD under no local ignorance for each evidence ${e_i}$ is represented as:

$$\begin{equation}{e_i} = \left\{ {\left( {{\theta _h},{p_{{h,i}}}} \right),\,h = 1, \cdots ,H;\mathop \sum \nolimits_{{\theta _h} = 1}^M {p_{{h,i}}} = 1} \right\},\,i = 1, \cdots ,N.\end{equation} \tag{ B.1 }$$

Then WBDR is:

$$\begin{equation}{m_i} = \{ \left( {{\theta _h},{{\tilde w}_i}{p_{h,i}}} \right),h = 1, \cdots ,H;{\text{(}}{P_i}\left( \Theta \right),\left( {1 - {{\tilde w}_i}} \right){\text{\} }},\,i = 1, \cdots ,N\end{equation} \tag{ B.2 }$$

where ${\theta _h}$ represents the class and ${p_{h,i}}$ is the output probability of individual classifier $i$. ${\tilde w_i}$ is the new weight.

Since normalization in the evidence combination can be performed in the end and does not change the results, it is applied it in the end. Assume that ${\hat m_{{\theta _h},l}},\,h = 1, \cdots ,H$ and ${\hat m_{P\left( {{\Theta }} \right),l}}$ denote the WBDR generated by combining the first $l$ evidence. The combination of two evidences with $l = 2$ is considered at first. And the combined WBDR which is generated through integrating the two evidences using orthogonal sum operation are calculated as:

$$\begin{align} {{{\hat m}_{{\theta _h},2}}}& = {{m_{{\theta _h},1}}{m_{P\left( \Theta \right),2}} + {m_{{\theta _h},2}}{m_{P\left( \Theta \right),1}} + {m_{{\theta _h},1}}{m_{{\theta _h},2}}}\nonumber \\ {}& = {{m_{{\theta _h},1}}\left( {{m_{{\theta _h},2}} + {m_{P\left( \Theta \right),2}}} \right) + {m_{{\theta _h},2}}{m_{P\left( \Theta \right),1}}}\nonumber \\ {}& = {{m_{{\theta _h},1}}\left( {{m_{{\theta _h},2}} + {m_{P\left( \Theta \right),2}}} \right) + {m_{P\left( \Theta \right),1}}\left( {{m_{{\theta _h},2}} + {m_{P\left( \Theta \right),2}}} \right) - {m_{P\left( \Theta \right),1}}{m_{P\left( \Theta \right),2}}} \\ {}& = {\left( {{m_{{\theta _h},1}} + {m_{P\left( \Theta \right),1}}} \right)\left( {{m_{{\theta _h},2}} + {m_{P\left( \Theta \right),2}}} \right) - {m_{P\left( \Theta \right),1}}{m_{P\left( \Theta \right),2}}} \nonumber\\ {}& = {\mathop \prod \nolimits_{i = 1}^2 \left( {{m_{{\theta _h},i}} + {m_{P\left( {{\Theta }} \right),i}}} \right) - \mathop \prod \nolimits_{i = 1}^2 {m_{P\left( \Theta \right),i}}.} \nonumber \end{align} \tag{ B.3 }$$

And

$$\begin{equation}{\hat m_{P\left( {{\Theta }} \right),2}} = \mathop \prod \nolimits_{i = 1}^2 {m_{P\left( \Theta \right),i}}.\end{equation} \tag{ B.4 }$$

Assume that the following equations are true for the (l-1) evidences. Let ${l_1} = l - 1$ and:

$$\begin{equation}{\hat m_{{\theta _h},{l_1}}} = \mathop \prod \nolimits_{i = 1}^{l - 1} \left( {{m_{{\theta _h},i}} + {m_{P\left( {{\Theta }} \right),i}}} \right) - \mathop \prod \nolimits_{i = 1}^{l - 1} {m_{P\left( \Theta \right),i}},\end{equation} \tag{ B.5 }$$

$$\begin{equation}{\hat m_{P\left( {{\Theta }} \right),{l_1}}} = \mathop \prod \nolimits_{i = 1}^{l - 1} {m_{P\left( \Theta \right),i}}.\end{equation} \tag{ B.6 }$$

The above combined probability output is further aggregated with the lth evidence. So the combined probability masses are given as:

$$\begin{align} {{{\hat m}_{\theta ,l}}}& = {{{\hat m}_{{\theta _h},{l_1}}}{m_{{\theta _h},l}} + {{\hat m}_{{\theta _h},{l_1}}}{m_{P\left( {{\Theta }} \right),l}} + {m_{{\theta _h},l}}{{\hat m}_{P\left( {{\Theta }} \right),{l_1}}}} \nonumber\\ {}& = {{{\hat m}_{{\theta _h},{l_1}}}({m_{{\theta _h},l}} + {m_{P\left( {{\Theta }} \right),l}}) + {m_{{\theta _h},l}}{{\hat m}_{P\left( {{\Theta }} \right),{l_1}}}} \nonumber\\ {}& = {{{\hat m}_{{\theta _h},{l_1}}}({m_{{\theta _h},l}} + {m_{P\left( {{\Theta }} \right),l}}) + {{\hat m}_{P\left( {{\Theta }} \right),{l_1}}}\left( {{m_{{\theta _h},l}} + {m_{P\left( {{\Theta }} \right),l}}} \right) - {{\hat m}_{P\left( {{\Theta }} \right),{l_1}}}{m_{P\left( {{\Theta }} \right),l}}} \nonumber\\ {}& = {({{\hat m}_{{\theta _h},{l_1}}} + {{\hat m}_{P\left( {{\Theta }} \right),{l_1}}})({m_{{\theta _h},l}} + {m_{P\left( {{\Theta }} \right),l}}) - {{\hat m}_{P\left( {{\Theta }} \right),{l_1}}}{m_{P\left( {{\Theta }} \right),l}}} \\ {}& = {({m_{{\theta _h},l}} + {m_{P\left( {{\Theta }} \right),l}})\left( {\left( {\mathop \prod \nolimits_{i = 1}^{l - 1} \left( {{m_{{\theta _h},i}} + {m_{P\left( {{\Theta }} \right),i}}} \right) - \mathop \prod \nolimits_{i = 1}^{l - 1} {m_{P\left( \Theta \right),i}}} \right) + \mathop \prod \nolimits_{i = 1}^{l - 1} {m_{P\left( \Theta \right),i}}} \right)}\nonumber\\ &\quad - {{m_{P\left( {{\Theta }} \right),l}}\mathop \prod \nolimits_{i = 1}^{l - 1} {m_{P\left( \Theta \right),i}}} \nonumber\\ {}& = {\mathop \prod \nolimits_{i = 1}^{\,\,l} \left( {{m_{{\theta _h},i}} + {m_{P\left( {{\Theta }} \right),i}}} \right) - \mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( \Theta \right),i}}.}\nonumber \end{align} \tag{ B.7 }$$

And,

$$\begin{equation}{\hat m_{P\left( {{\Theta }} \right),l}} = {m_{P\left( {{\Theta }} \right),l}}{\hat m_{P\left( {{\Theta }} \right),{l_1}}} = {m_{P\left( {{\Theta }} \right),l}}\mathop \prod \nolimits_{i = 1}^{l - 1} {m_{P\left( \Theta \right),i}} = \mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( \Theta \right),i}}.\end{equation} \tag{ B.8 }$$

In the following, the combined WBDR results are normalized. Assume that $k$ is the normalization factor, therefore,

$$\begin{equation}k\left( {\mathop \sum \nolimits_{h = 1}^H {{\hat m}_{{\theta _h},l}} + {{\hat m}_{P\left( {{\Theta }} \right),l}}} \right) = 1.\end{equation} \tag{ B.9 }$$

That is:

$$\begin{equation}k\left( {\mathop \sum \nolimits_{h = 1}^H \left( {\mathop \prod \nolimits_{i = 1}^{\,\,l} \left( {{m_{{\theta _h},i}} + {m_{P\left( {{\Theta }} \right),i}}} \right) - \mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( {{\Theta }} \right),i}}} \right) + \mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( {{\Theta }} \right),i}}} \right) = 1,\end{equation} \tag{ B.10 }$$

$$\begin{equation}k\left( {\mathop \sum \nolimits_{h = 1}^H \left( {\mathop \prod \nolimits_{i = 1}^{\,\,l} \left( {{m_{{\theta _h},i}} + {m_{P\left( {{\Theta }} \right),i}}} \right)} \right)} \right) - k\left( {H - 1} \right)\mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( {{\Theta }} \right),i}} = 1.\end{equation} \tag{ B.11 }$$

So,

$$\begin{equation}k = {\left( {\mathop \sum \nolimits_{h = 1}^H \left( {\mathop \prod \nolimits_{i = 1}^{\,\,l} \left( {{m_{{\theta _h},i}} + {m_{P\left( {{\Theta }} \right),i}}} \right)} \right) - \left( {H - 1} \right)\mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( {{\Theta }} \right),i}}} \right)^{ - 1}},\end{equation} \tag{ B.12 }$$

Therefore,

$$\begin{equation}{m_{\theta ,l}} = k{\hat m_{\theta ,l}},\,\,{m_{P\left( {{\Theta }} \right),l}} = k{\hat m_{P\left( {{\Theta }} \right),l}},\end{equation} \tag{ B.13 }$$

where ${m_{\theta ,l}}$ and ${m_{P\left( {{\Theta }} \right),l}}$ are the combined WBDR after normalization. So the BD ${p_m}$after combining l evidence is:

$$\begin{equation}\,{p_h} = \frac{{{m_{\theta ,l}}}}{{1 - {m_{P\left( \Theta \right),l}}}} = \frac{{k{{\hat m}_{\theta ,l}}}}{{1 - k{{\hat m}_{P\left( {{\Theta }} \right),l}}}} = \frac{{k\left( {\mathop \prod \nolimits_{i = 1}^{\,\,l} \left( {{m_{\theta ,i}} + {m_{P\left( {{\Theta }} \right),i}}} \right) - \mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( \theta \right),i}}} \right)}}{{1 - k\mathop \prod \nolimits_{i = 1}^{\,\,l} {m_{P\left( \theta \right),i}}}},\,h = 1, \cdots ,H.\end{equation} \tag{ B.14 }$$

Based on equation (A.4) and $l = N$, the final BD is:

$$\begin{equation}{p_h} = \frac{{k\left[ {\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{{w_i}{p_{h,i}}}}{{1 + {w_i} - {r_i}}} + \frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right) - \mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)} \right]}}{{1 - k\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)}},\,h = 1, \cdots ,H,\end{equation} \tag{ B.15 }$$

and $k$ is:

$$\begin{equation}k = {\left[ {\mathop \sum \nolimits_{h = 1}^H \left( {\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{{w_i}{p_{h,i}}}}{{1 + {w_i} - {r_i}}} + \frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)} \right) - \left( {H - 1} \right)\mathop \prod \nolimits_{i = 1}^{\,\,N} \left( {\frac{{1 - {r_i}}}{{1 + {w_i} - {r_i}}}} \right)} \right]^{ - 1}}.\end{equation} \tag{ B.16 }$$