Weighted Concordance Index Loss-based Multimodal Survival Modeling for Radiation Encephalopathy Assessment in Nasopharyngeal Carcinoma Radiotherapy

We acquire eligible 4,816 patients from a cancer center. The eligibility of enrolling in the NPC-REP dataset is that diagnosis and treatment data in NPC radiotherapy and follow-up are traceable and recorded. We binarize a label for each patient according to whether the REP confirmation by MR image diagnosis during follow-up. If diagnosing out RTLI, the label of NPC patient is REP, otherwise non-REP. At the same time, we record the time intervals (in Months) between NPC radiotherapy and REP confirmation or the last follow-up of non-REP. Then we trace back the diagnosis and treatment data during NPC radiotherapy for each patient. RD images show the radiation dose distribution to the human body during radiation therapy. Clinicians sketch the tumor outline in closed curve coordinate form to yield the RS image. Then we apply the mask of RS images to generate the input ROIs of RD and CT images. As Fig 1 shows, we crop out the identical regions of interest (ROI) for CT images for diagnosis and RD images according to regions of the bilateral temporal lobe in RS images. Then we feed the three 3D ROIs into the multimodal learning network for extracting visual features. It is confirmed that the imaging dose delivered during NPC radiotherapy is positively related to REP risk [24]. Hence, in addition to considering image data, we introduce 18 features manually calculated from DVH, which relates radiation dose to tissue volume in radiation therapy planning [8]. Moreover, we also observe the affection of demographic information (age, gender), clinic stages (T/N/TNM), and treatment options (radio-chemotherapy and radiotherapy) to REP risk. The detailed statistics of demographic and clinical characteristics for the NPC-REP dataset are shown in Appendix A1.

We notice that most NPC patients with REP are diagnosed after radiotherapy for about 37 months. And it is almost impossible to occur REP when the follow-up time is less than two months and more than 73 months. Such a statistics observation about the onset time of REP illustrates that the risk probability of suffering REP is highly related to some factors in NPC radiotherapy. These factors cause the onset of REP in a relatively definite period. At the same time, we can cast REP assessment as a 3-years survival analysis task according to the average REP confirmation time intervals. Thus we build a deep multimodal survival network (MSN) to identify potential factors from the NPC-REP dataset, which consists of a binarized label, time intervals, image data (CT, RD, and RS images), and non-image data (24 features in total). As Fig 2 shows, one feature extractor learns visual features of 3D ROIs intercepted from CT, RD, and RS images, and the other explores the contributions of non-image data to REP risk. The 24 features of non-image data are directly fed into the MLP module to output the REP risk $R_{nv}$. In the branch of visual feature extractor, we first apply convolution network as a backbone, followed by non-local attention capturing the global inter-dependencies [23]. Next, we use a global average pooling (GAP) to merge feature maps and two fully-connected layers (FCs) to generate the REP risk $R_{v}$. After generating the two risk probabilities for two modality data, we further utilize our WCI loss function to train the network. The combination of the two predicted probabilities is defined as follows:

where $R_{nv},R_{v}\in[0,1]$ and $w_{nv}$ and $w_{v}$ are their corresponding weights. Both are initially set to 0.5. We can tune the two weights to observe the contribution of different modal data on assessing the REP risk.

CI is a standard performance metric for survival models that corresponds to the probability that the model correctly orders a randomly chosen pair of patients in terms of event time [4, 25]. To calculate this indicator, we first pair all samples in the dataset of sample size $n$ with each other to get $n(n-1)/{2}$ sample pairs. Then, we filter out evaluation sample pairs from all sample pairs, as follows:

where $T_{i}$ and $T_{j}$ are the time intervals of $i^{th}$ and $j^{th}$ NPC patients respectively, and $O_{i}=1$ denotes that $i^{th}$ NPC patient is observed as REP. $\mathbb{E}$ is the set of sample pairs $(i,j)$ based on the ordering of survival times, i.e., the time intervals of $i^{th}$ NPC patient is shorter than $j^{th}$ NPC patient. Assuming the REP risks predicted by models for sample pairs $(i,j)$ are $R_{i}$ and $R_{j}$, we assert the prediction outcomes are correct if $R_{i}>R_{j}$ and $T_{i}<T_{j}$. We mark those sample pairs whose predicted results coincide with the actual as a set $\mathbb{E}_{t}$. The CI value is indicated as the number ratio of $\mathbb{E}_{t}/{\mathbb{E}}$. The higher CI explains away the better competence of survival models in predicting REP risks.

Because the assumptions that proportional hazards are invariant over time and log-linear limit the Cox loss function (See Appendix A2) used for survival models, BCI loss function (See Appendix A3) directly formulates the optimizing problem of survival models as maximizing the CI. The Cox depends only on the ranks of the observed survival time rather than on their actual numerical value [21]. Given the CI as a performance measure in survival analysis, BCI can effectively exploit the changes of proportional hazards with survival time. However, BCI is sensitive to uneven sampling, thus yielding poor predictive accuracy. Uneven sampling refers to the sample pairs meeting the constraints of Equation 2 in each batch failing to equally cover REP samples ($O_{i}=1$) which are the key learning objectives. BCI balances CI loss for all REP samples in each batch by adopting a global average operation, not paying more attention to those REP samples with fewer sample pairs.

By referring to the weighted cross-entropy loss function, our WCI enhances the contributions of those REP samples with fewer sample pairs in each batch by a dual average operation, as follows:

where $\tau$ is a temperature hyper-parameter [10, 26] and initially set as 0.1. $N_{O=1}$ (outer) is the number of REP samples, and $N_{T_{j}\geq T_{i}}$ (inner) is the number of sample pairs of $i^{th}$ REP sample. If the relation of predicted risk probabilities is $R_{i}<R_{j}$ (wrong prediction), the CI loss is amplified, and if $R_{i}>R_{j}$ (correct prediction), the CI loss infinitely approximates to zero. After tuning the risk difference, we use an inner average operation $1/{N_{T_{j}\geq T_{i}}}$ and an outer average operation $1/{N_{O=1}}$ to prevent fluctuation in the optimization process caused by uneven sampling per batch. The outer average operation is balanced for all REP samples, while the inner average operation assigns different weights for REP samples according to their number of sample pairs.

Dataset. We split the NPC patients in the NPC-REP dataset as a training set with 3,253 patients (264 REP and 2,989 non-REP), a validation set with 575 patients (43 REP and 532 non-REP), and a test set with 988 patients (61 REP and 927 non-REP). Moreover, we view left and right temporal lobes as independent samples during model training and inference.

Loss functions. We experiment with four loss functions to demonstrate the gains of our WCI, including standard cross-entropy (CE) [27], censored cross-entropy (CCE) [25], BCI [21], and Cox [16]. We adopt the survival time of 3-years (36 Months) for training and evaluation in terms of the time intervals in the NPC-REP dataset. For the CE, we define the label as 1 if REP confirmation and as 0 if the last follow-up time of non-REP is beyond 36 months. And the CCE extends CE used for classification models to train survival predictions with right-censored data by discretizing event time into intervals $(0,36],(36,\infty)$. The BCI maximizes a lower bound on CI, explaining the high CI-scores of proportional hazard models observed in practice. The Cox for every sample is a function of all samples in the training data and is usually used for fitting Cox proportional hazard models [2, 7].

Evaluation metrics. CI is the most frequently used evaluation metric of survival models [1, 11]. It is defined as the ratio of correctly ordered (concordant) pairs to comparable pairs. Two samples $i$ and $j$ are comparable if the sample with lower observed time $T$ confirmed REP during follow-up, i.e., if $T_{j}>T_{i}$ and $O_{i}=1$, where $O_{i}=1$ is a binary event indicator. A comparable pair $(i,j)\in\mathbb{E}$ is concordant if the estimated risk $R$ by an assessment model is higher for NPC patients with lower REP confirmation time, i.e., $R_{i}>R_{j}\wedge T_{j}>T_{i}$, otherwise the pair is discordant. We also report the time-dependent AUC [15, 20] for all loss functions, apart from CI.

Optimizer. Our MSN model is trained from scratch using the SGD optimizer for all compared loss functions. Specifically, the learning rate is decay scheduled by cosine annealing [18] setting 2e-4 as an initial value and 1e-3 from 5th epochs. We set the batch size as a multiple of 8 for all tasks, i.e., 128 (8 GPUs with 16 input images per GPU). The parameters of networks are optimized in 60 epochs with a weight decay of 3e-5 and a momentum of 0.8.

Multimodal data Analysis. The intention of acquiring treatment and diagnosis data in NPC radiotherapy launches on the confirmation that the imaging dose delivered during NPC radiotherapy is positively related to REP risk. Hence, we employ our MSN to jointly leverage non-image data (DVH+DT6) and images (CT+RD+RS) to predict REP risk. We report the CI and 3-years AUC in Fig 3 for compared loss functions used for training our MSN. We can observe that combining DVH and DT6 (orange line) brings more gains than DVH (blue line) alone or CT+RD+RS (green line), and using all three data types (red line) achieves the best predictive accuracy. Such experimental results can demonstrate the contributions of multimodal data of NPC radiotherapy to the improvement of REP prediction and give supportive evidence for the above confirmation. Next, we analyze the contribution proportion of two modality data: images and non-image data. We conduct an ablation study for Equation 1 to find the best weight ratio for two probabilities predicted by the two branches in our MSN. Table 1 shows the best CI (bold) achieved by setting the proportion of $w_{nv}$ and $w_{v}$ as 0.7:0.3. The more contributions of non-image data than images coincide with the reports of Fig 3 in which non-image data (orange line) yields more gains than image data (green line).

Loss Function Analysis. After confirming the utility of multimodal data and analyzing the contributions of different modalities, we further discuss the benefits of our WCI according to the report in Fig 3. We interpret CI as the probability of correct ranking of REP occurring probability by casting predicting REP risk as a survival analysis task. Hence, CE used for learning classification models can not achieve sound accuracy without considering event times. Compared to CE, Cox can significantly improve the performance by exploiting the ranks of the observed survival time. On the other hand, Cox is beat by BCI due to ignoring the changes of proportional hazards with survival time. Although BCI enjoys the gain of directly maximizing CI, its performance is slightly lower than CCE due to sensitivity to uneven sampling. Identical to BCI, CCE also considers the contributions of the proportional hazard of sample pairs in learning the survival model. Based on the comparable performance of BCI and CCE, WCI without $\tau$ can further obtain better performance by encouraging our MSN to effectively learn REP samples on the condition of uneven sampling per batch. By using $\tau$ help model convergence, our WCI achieves the highest CI (82.52) and 3-years AUC (83.75) on two modality data (DVH+DT6+CT+RD+RS). Although WCI is designed initially for the CI indicator used for evaluating survival models, the best performance on 3-years AUC can also prove the excellent generalization of our MSN trained by our WCI. We also conduct McNemar tests [17] on CI between our WCI and other loss functions, including CE, Cox, BCI, CCE, and WCI w/o $\tau$. All comparisons in a p-value of less than 0.01 can support that our WCI brings significant improvement to CI.

Ablation studies of $\tau$. By casting REP assessment as a survival analysis task, our WCI exhibits manifest advantage over other loss functions by introducing the dual average operation and temperature hyper-parameter $\tau$. Apart from the dual average operation preventing our MSN from being sensitive to uneven sampling, we also utilize $\tau$ to sharpen the risk difference of sample pairs to help model convergence. We conduct an ablation study for $\tau$ to observe the heating and cooling effect of tuning the CI loss for our WCI. Table 2 shows that the best value of $\tau$ used to train our MSN is 0.1 (bold). And when $\tau=1.0$, our WCI without imposing the tuning on CI achieves the second-highest CI of 82.24 (underline) over two modality data (DVH+DT6+CT+RD+RS). When $\tau>1.0$, it occurs the opposite effect. At the same time, the adjustment force should not be too large from the performance of setting $\tau=0.05$ or $\tau=0.02$. The ablation studies in Table 2 confirm that $\tau$ can help model convergence by sharpening the risk difference. Theoretically, the significant risk difference helps improve the predicted certainty probability, i.e., reduces the entropy of probability distribution. Benefiting from the insensitive to uneven sampling and sharpening risk difference, our WCI exhibits more stability and convergence than BCI by comparing standard deviations of batch loss (See Appendix A4).