## 2 Materials and Methods
### Feature Selection  
Two different techniques were used and compared to select the best subset of radiomic features: LASSO regression [19] and genetic algorithm [20]. We chose the LASSO regression as it is a standard method for the feature selection procedure. The LASSO regression selects the features, minimizing the L1-loss function. Contrarily, the genetic algorithm allows a customization of the loss function, which can be tailored for the specific application.  
We used a Cox proportional hazards model with a LASSO penalty, also known as an \(1_{1}\) penalty, to fit the training data. The number of features to preserve was parametrized by the penalization coefficient of the LASSO regression: a higher penalty leads to fewer features being preserved. We used 100 penalization values evenly spaced in [10\({}^{-4}\), 1]. For each value, a Cox model was trained and evaluated using 5-fold cross-validation. We repeated the cross-validation procedure 5 times for each Cox model to give a more reliable estimation of the model performance. The parameter yielding the best cross-validation result was selected as the optimal one. Then, the best Cox model found was used for the estimation of \(\beta\) parameters: the LASSO model filters out all of the features associated to a corresponding zero value of \(\beta\).  
We developed a genetic algorithm (GA) for the feature selection procedure. Each genome identified a putative set of features, expressed as a binary pattern of the considered features. Features associated to no-null genes were preserved; null genes identified features that were excluded. Each genome was assessed using a Cox model. The resulting model was evaluated using a metric function, and the scores of the models were used as fitness values. The first model evaluation was performed using the Akaike information criterion (AIC) [21] as a metric (Figure 2). The AIC score penalizes models with a high number of features. We estimated the best model, and the corresponding subset of features, as the one with the lowest AIC score. Using the subset of features obtained by the first filtering, a second reduction was performed according to the maximization of the fitness function given by:  
\[Fitness=\mu(CV_{5})-\sigma(CV_{5}) \tag{1}\]  
where \(\mu(CV_{5})\) and \(\sigma(CV_{5})\) are the average and standard deviation scores obtained by the model in a 5-fold cross-validation repeated five times, respectively. We used the concordance index (_CI_) metric for the evaluation of the cross-validation. The concordance index [22] is defined as:  
\[CI=\frac{n^{\circ}concordantpairs}{n^{\circ}concordantpairs+n^{\circ}discordantpairs} \tag{2}\]  
Two patients were correctly ordered if the patient with the higher predicted risk experienced the death event (binary score) before the patient with the lower predicted risk. The concordant pairs were the pairs of patients correctly ordered by the model, while the discordant pairs were those that were incorrectly ordered. This fitness function aimed to maximize the overall score obtained in the repeated cross-validation, while achieving a similar score for all the cross-validation folds, by minimizing the standard deviation.  
We performed a 5-fold cross-validation repeated five times for the evaluation of the Cox model, i.e., the identification of the best features subset, using the concordance index as an evaluation metric.  
CT and PET features were independently selected and used to build independent models to evaluate their different predictive power. The two selected feature subsets were merged, further reduced by maximizing (1) with the GA and used to build a comprehensive Cox model to estimate the predictive capability of their combination. All the selected features were included in the final models, except for the ones with a \(p\)-value greater than 0.2, which were manually excluded.