From January 2006 to December 2018, data of patients diagnosed with SGCs and treated at Severance Hospital were retrospectively analyzed. The inclusion criteria are as follows: 1) primary malignant tumor in major or minor salivary gland and 2) diagnosed with and underwent surgery for SGCs with sufficient clinical and pathological information was available through medical records. Exclusion criteria are as follows: 1) distant metastasis at time of diagnosis, 2) surgery or radiotherapy previously performed in the head and neck area, and 3) loss during follow-up after surgery was excluded from the study. Finally, a total of 460 patients was included in the study consisting of 240 males and 220 females. The ages ranged from 14 to 99 years with a mean of 53.7 years. Tumor stage was classified based on the 8th American Joint Committee on Cancer (AJCC) staging system. This study was approved by the Institutional Review Board (IRB) of Yonsei University (3-2020-0504). Need for individual informed consent was waived as this study had a retrospective design.

Of the data set, 80% were classified as training data and used for learning, while the remaining 20% were used as test data to validate the performance of machine learning model. Survival data of cancer patients cannot be predicted by binary classification and ‘time to event’ and censored data should be considered for predicting survival. Several metrics such as accuracy, precision, recall, and area under receiver operating curve (AUC) which are commonly used to estimate the performance of binary classification are not suitable in cancer survival analysis. Therefore, Harrell’s c-index was used to analyze the performance of a survival prediction model for SGC patients [7]. The Harrell’s c-index is the most widely used index to evaluate the predictive accuracy of survival time considering death in a survival analysis model. It is a representative method for evaluating the performance of the survival prediction model and Harrell’s c-index uses a statistical method that indicates whether the predicted death time of a patient through the model is sequentially matched with the real death time of the actual patients. When the c-index value is 0.5, it corresponds to the average value of the random model and c-index value of 1 refers to a perfect match of death time ranking. C-index value over 0.8 usually indicates the strong predictive model [7,8].

We choose the following three representative models for predicting patient’s survival. The Random Survival Forest, CPH, and deep learning-based survival model (DeepSurv) were used to construct machine learning survival prediction model for SGCs patients [8-10]. The linear model CPH, the nonlinear model RFS, and the deep learning model Deepsurv were used to predict the survival of SGC patients, and to determine which model can show the best performance in predicting the survival of these patients. As Random Survival Forest uses all variables for constructing a model and assesses nonlinear effects of them, it can reduce variance and bias. To construct the Random Survival Forest models, the 0.14.0 version of scikit survival was used. The DeepSurv model devised by Katzman, et al. [8] is an open-source python module and updates weights through a feed forward network and back propagation process through a multi-layer neural network structure and presents a negative log partial likelihood. CPH model was constructed by using R program with moonBook package.

Patient demographic information, tumor location and stage, pathologic findings, recurrence, date of recurrence, recurrence site, death, date of death, and cause of death were collected and analyzed. Chi-square or Fisher’s exact test was used to evaluate differences in categorical variables between the two independent groups. An independent two-sample t-test was used to assess differences in continuous variables between the two independent groups. The Kaplan-Meier curve was used to analyze patient survival, and outcomes were assessed using a log-rank test. A p-value <0.05 was considered to indicate statistical significance. Statistical analyses were performed using R version 4.03.

A total of 460 patients was included in this study, and all patients underwent surgery as an initial treatment. After surgical treatment, 184 patients (40%) received postoperative radiotherapy, and 70 patients (15.2%) underwent concurrent chemoradiation. The primary site of cancer was the parotid gland in 388 patients (84.3%), submandibular gland in 46 patients (10.0%), sublingual gland in 20 patients (4.3%), and minor salivary gland in 6 patients (1.3%). There were 34 cases (7.4%) with facial nerve (FN) palsy before surgery. On pathologic examination of surgical specimens, 146 patients (31.7%) showed positive surgical margins, and 314 patients (68.3%) showed negative margins. LVI findings were observed in 74 patients (16.1%), PNI findings in 117 patients (25.4%), and extranodal extension (ENE) findings in 70 patients (15.2%). In pT classification, T1 was 145 (31.5%), T2 was 195 (42.4%), T3 was 70 (15.2%), and T4 was 50 (10.9%). In pN classification, N0 was 354 (77.0%), N1 was 28 (6.1%), N2 was 72 (15.7%), and N3 was 6 (1.3%). On the TNM staging system, stage I was 129 (28.0%), stage II was 153 (33.3%), stage III was 75 (16.3%), and stage IV was 103 (22.4%). On survival analysis using Kaplan-Meier curve, significant differences were observed between stage IV and other stages, but there was no significant difference between stages II and III (Fig. 1A). Other information of all patients is summarized in Table 1. In addition, baseline differences between the training and test data sets were analyzed and summarized in Table 1, and the overall survival difference between the two groups was analyzed through Kaplan-Meier survival analysis with the log-rank test (Fig. 1B). There were no significant differences between training and test data set.

The clinical prognostic factors used were: patient’s age, sex, primary site, FN palsy, adjuvant treatment, margin status, LVI, PNI, ENE, pT, pN, number of metastatic LNs, lymph node ratio (LNR), histology subtype. In the case of a linear model, multicollinearity may be a problem, so the GridSearchCV method was used to select the type of variable to be used in the CPH model. For nonlinear models, RFS and Deepsurv models, all variables were used to build the model without special variable selection because multicollinearity was not an issue.

CPH model has been constructed using R program with moonBook package. Firstly, univariate analysis was performed to evaluate prognostic factors for survival in SGCs patients. pM, LVI, ENE, FN palsy, PNI, stage, adjuvant therapy, pN, surgical margin, pT, LNR, number of metastatic LNs, age, subsite, and sex showed significant correlations with SGCs patient’s survival. On multivariate analysis, pM, stage, LVI, LNR, and age exhibited significant correlation with patient’s survival (Fig. 2). In the CPH model, when the model was built using all the parameters mentioned above, the c-index value was 0.85 for the training set and 0.81 for the test set. After tuning parameter values through the GridSearchCV function, optimal performance was achieved when 8 parameters including age, LVI, PNI, ENE, pT, pN, stage, and number of metastatic LNs were input into the model. When the CPH model was built based on these 8 parameters, the c-index value was 0.85 for the training set and 0.80 for the test set.

In the Random Survival Forest model, the c-index value for the training set was 0.86, and that for the test set was 0.82. Prediction error is calculated using OOB data (Fig. 3). The variable importance (VIMP) was obtained by measuring decrease in prediction accuracy when randomizing a particular variable. Variable with higher VIMP tends to contribute more to predictive accuracy. Feature importance of variables was shown in the order of stage, number of metastatic LNs, age, pN, pT, FN palsy, LNR, ENE, and surgical margin. Stage and age exhibited high importance in both the Random Survival Forest and CPH models. Interaction between variables is measured based on minimal depth. Stage and the number of metastatic LNs showed lowest minimal depth and are expected to be associated with other variables.

In the DeepSurv model, the c-index value was 0.72 for the training set and was 0.72 for the test set. The learning process of the DeepSurv model is visualized in Fig. 4 and overfitting was not observed in the DeepSurv model. Among the three models mentioned above, Random Survival Forest exhibited the highest performance in predicting survival of SGC patients.