Stable prediction with radiomics data
Carel F. W. Peeters Caroline Übelhör Steven W. Mes Roland Martens Thomas Koopman Pim de Graaf Floris H. P. van Velden Ronald Boellaard Jonas A. Castelijns Dennis E. te Beest Martijn W. Heymans Mark A. van de Wiel









                        
                        
 	                    
                             
                      
  









 










[]



















plain
thmTheorem
lemmaLemma
assumeAssumption
propositionProposition
corollaryCorollary
definitionDefinition
definition
exampleExample
*example_contdExample 1
remark
remarkRemark
BconBiological conjecture



⊥⊥


    
    
               







R
language=R,
basicstyle=,
numbers=none,
breaklines=true,
keywordstyle=,
commentstyle=,
stringstyle=,
otherkeywords=, {, }, 
    ,
,
frame=none,
tabsize=2,
backgroundcolor=,
caption=R code


.8in
.7in

plainnat

Stable Prediction with Radiomics Data]
Stable Prediction with Radiomics DataC.F.W. Peeters et al.]Carel F.W. Peeters**Corresponding author.[Carel F.W. Peeters]
Dept. of Epidemiology & Biostatistics 

Amsterdam Public Health research institute 

Amsterdam University medical centers, location VUmc 

Amsterdam

The Netherlandscf.peeters@vumc.nl]Caroline Übelhör[Caroline Übelhör]
Dept. of Epidemiology & Biostatistics 

Amsterdam Public Health research institute 

Amsterdam University medical centers, location VUmc 

Amsterdam

The Netherlands]
Steven W. Mes[Steven W. Mes]
Dept. Otolaryngology/Head & Neck Surgery 

Amsterdam University medical centers, location VUmc 

VUmc Cancer Center 

Amsterdam

The Netherlands]Roland Martens[Roland Martens]
Dept. of Radiology & Nuclear Medicine 

Amsterdam University medical centers, location VUmc 

VUmc Cancer Center 

Amsterdam

The Netherlands]Thomas Koopman[Thomas Koopman]
Dept. of Radiology & Nuclear Medicine 

Amsterdam University medical centers, location VUmc 

VUmc Cancer Center 

Amsterdam

The Netherlands]Pim de Graaf[Pim de Graaf]
Dept. of Radiology & Nuclear Medicine 

Amsterdam University medical centers, location VUmc 

VUmc Cancer Center 

Amsterdam

The Netherlands]Floris H.P. van Velden[Floris H.P. van Velden]
Dept. of Radiology 

Leiden University medical center

Leiden

The Netherlands]Ronald Boellaard[Ronald Boellaard]
Dept. of Radiology & Nuclear Medicine 

Amsterdam University medical centers, location VUmc 

VUmc Cancer Center 

Amsterdam

The Netherlands]Jonas A. Castelijns[Jonas A. Castelijns]
Dept. of Radiology & Nuclear Medicine 

Amsterdam University medical centers, location VUmc 

VUmc Cancer Center 

Amsterdam

The Netherlands]
Dennis E. te Beest[Dennis E. te Beest]
Biometris 

Wageningen University & Research 

Wageningen 

The Netherlands]Martijn W. Heymans[Martijn W. Heymans]
Dept. of Epidemiology & Biostatistics 

Amsterdam Public Health research institute 

Amsterdam University medical centers, location VUmc 

Amsterdam

The Netherlands]Mark A. van de Wiel[Mark A. van de Wiel]
Dept. of Epidemiology & Biostatistics 

Amsterdam Public Health research institute 

Amsterdam University medical centers, location VUmc 

Amsterdam

The Netherlands; 
MRC Biostatistics Unit 

Cambridge University 

Cambridge 

United Kingdom


 
Motivation: Radiomics refers to the high-throughput mining of quantitative features from radiographic images.
It is a promising field in that it may provide a non-invasive solution for screening and classification.
Standard machine learning classification and feature selection techniques, however, tend to display inferior performance in terms of (the stability of) predictive performance.
This is due to the heavy multicollinearity present in radiomic data.
We set out to provide an easy-to-use approach that deals with this problem.

Results:
We developed a four-step approach that projects the original high-dimensional feature space onto a lower-dimensional latent-feature space, while retaining most of the covariation in the data.
It consists of (i) penalized maximum likelihood estimation of a redundancy filtered correlation matrix.
The resulting matrix (ii) is the input for a maximum likelihood factor analysis procedure.
This two-stage maximum-likelihood approach can be used to (iii) produce a compact set of stable features that (iv) can be directly used in any (regression-based) classifier or predictor.
It outperforms other classification (and feature selection) techniques in both external and internal validation settings regarding survival in squamous cell cancers.

Availability: The  implementation of the pipeline () as well as scripts for the data experiments are available at <https://github.com/CFWP/FMradio>.

Contact: cf.peeters@vumc.nl

Supplementary information: Supplementary information is available at XXX online.

 Key words:
Classification;
Factor analysis;
High-dimensional data;
Prediction models;
Radiomics
[
    [
    Received ; accepted 
========================








§ INTRODUCTION


 §.§ Focus

Radiomics is a relatively recent addition to the omics-scene.
It encompasses the high-throughput mining of quantitative features from radiographic images.
Such images stem from technology such as computed tomography (CT), positron emission tomography (PET), and magnetic resonance imaging (MRI).
Introduced by <cit.> and matured by <cit.>, it is rapidly finding widespread application, especially in oncology research <cit.>.
The promise of radiomic data is that it may be a basis, through the information contained in standard-of-care images, for non-invasive medical decision support at low (additional) cost.
As such, it is often viewed as an addition to or alternative for screening and classification based on molecular omics.


Irrespective of the medical imaging technology used, radiomics has a delineated workflow, roughly consisting of four steps.
The first is the acquisition of images and the identification of the volumes or regions of interest (VOIs).
The second step consists of image segmentation, i.e., the manual, computer-assisted or automated delineation of the borders of the VOIs.
The third step is feature extraction: the elicitation of features that characterize the segmented VOIs.
While this may be done in an unsupervised manner, the dominant approach is to use tailored extraction algorithms that rely on expert input.
Extracted features can range in the thousands <cit.> and enrich standard radiological lexicon (e.g., tumor extension, location, cellularity) with descriptors of distributions, textures, and morphology <cit.>.
For an overview of (issues related to) steps one to three, see <cit.>

We will concern ourselves with the fourth step: downstream analysis.
While there are some committed informatical resources for radiomics, these deal with feature extraction mostly <cit.>.
Hence, to date, there has been no dedicated workflow for the downstream analysis of radiomic data.
Here, we will focus on providing such a workflow for classification and prediction problems.
This focus answers recent calls for more robust methods of analysis <cit.> and aids in harnessing radiomics' translational potential.



 §.§ The collinearity problem

The main problem to overcome in radiomic-based classification and prediction is multicollinearity in the face of high-dimensional data.
Radiomic data are often high-dimensional, in the sense that the number of features, say p, exceeds the number of observations n.
Moreover, radiomic data typically deal with a heavy multicollinearity burden: (subsets of) features that are highly correlated or linearly dependent on other (subsets of) features.
Hence, there are two sources of ill-conditioning or singularity of the feature space that hamper estimation (stability) and feature selection.
The radiomics literature contains several approaches to deal with this problem.

The first is to use classifiers amended with a penalty on the feature space <cit.>.
Arguably the most popular such penalty is the ℓ_1 (lasso) penalty <cit.> which performs automatic feature selection.
While attractive, ℓ_1-regularized classifiers will not have a unique minimizer in the case of a heavy multicollinearity burden.
Moreover, its selection capability is constrained by the sample size n.
This combination can lead to unstable selection: multiple (non or only partly overlapping) sets of features with similar predictive performance.
The elastic net, which linearly combines an ℓ_1 and an ℓ_2 (ridge) penalty <cit.>, was developed to deal with lasso's shortcomings.
The effect of the ridge penalty is to force “the estimated coefficients of highly correlated predictors to be close to each other" <cit.>, after which the lasso (de)selects these features alltogether.
However, this behavior introduces heavy bias and inadequate prediction performance when faced with highly correlated features with a negative sign <cit.>, a situation common in radiomics data.
Summarizing, the penalized classification approach deals with the problem of high-dimensionality rather than the problem of strong collinearity.

A second approach that is often encountered is to perform heuristic feature selection based on correlation redundancy <cit.>.
From a cluster of highly correlated features either a representative feature is chosen or the cluster is collapsed into a representative feature.
This approach often implies a loss of information as more features are removed or collapsed than necessary.
Moreover, it does not guarantee to alleviate the collinearity burden to such a degree that standard classifiers can be expected to give adequate estimates and prediction performance.

A third approach (sometimes used in conjunction with the second approach) is to select features on the basis of bivariate analysis between the individual features and the outcome of interest <cit.>.
This approach consists of producing a ranking of the individual features on the basis of some test or measure of their bivariate association with the outcome.
A subset of features is then subsequently selected using a (often arbitrary) cut-off.
This approach has several problems.

As with the second approach, it does not guarantee to alleviate the collinearity burden to an acceptable degree.
Moreover, it does not reflect the behavior of features once they become part of a multi-feature classifier, as the dependencies between features are not considered.
This often results in unwanted feature inclusions and exclusions <cit.>.
Lastly, bivariate filtering will only produce stable and reliable feature-sets when the features are uncorrelated, which is not the situation with radiomics data.



 §.§ Proposed approach

What we desire is an approach that simultaneously deals with both the high-dimensionality as well as the collinearity burden of radiomics data.
Moreover, we desire an approach that produces stable and robust predictors that transfer well in validation settings.
We provide, through a four-step pipeline, such an approach based on projecting the radiomic feature-space onto a lower-dimensional (near-)orthogonal characteristic feature-space that retains most of the information contained in the full data set.
It consists of (i) penalized maximum likelihood (ML) estimation of a redundancy filtered correlation matrix.
The resulting matrix (ii) is the input for projection by a ML factor analysis procedure.
This two-stage maximum-likelihood (2SML) approach can be used to (iii) produce a compact set of stable characteristic features that (iv) can be directly used in any classifier or survival model.
The benefits of the approach are (a) computational speed and efficiency, (b) the ability to use standard downstream modeling and performance evaluation, (c) robustness of predictors, and (d) a potentially substantial increase in predictive performance.

The approach bears some flavor of principal component regression, in which first a principal component analysis (PCA) is performed after which a selection of retained components serve as regressors <cit.>.
However, PCA differs conceptually from common factor analysis (FA) <cit.>:
PCA provides a deterministic transformation of the data by identifying weighted linear combinations of observed features that maximize the explanation of (solely the) observed variance.
FA, in contrast, seeks to explain the covariation of and between the observed features through a small number of explanatory latent features (known as common factors).
Hence, FA provides a model-based transformation that includes measurement error.
As such, we favor FA over PCA.
Our approach to FA for high-dimensional data also differs from the mainstream.
Instead of inducing sparsity constraints on the (latent) parameter space <cit.>, we regularize the sufficient statistic for usage in the (regular) maximization of the model-based likelihood.
As such, the proposed approach resembles the road taken by <cit.>.



 §.§ Overview

Section <ref> details the individual steps of the proposed approach.
Section <ref> applies this approach to several real data sets.
It outperforms other classification (and feature selection) techniques in an MRI-imaging based external validation setting regarding survival in oropharyngeal cancer.
In addition, the approach gives the most stable performance in a PET/CT-imaging based internal validation setting regarding survival in head and neck cancer.
Section <ref> concludes with a discussion.





§ MATERIALS AND METHODS


 §.§ Setting

Let 𝐙^T≡[z_1,…,z_n] define p-variate vectors of radiomic scores on i = 1, …, n subjects with z_i^T≡[z_i1,…,z_ip]∈ℝ^p as a realization of the random vector Z_i^T≡[Z_i1,…,Z_ip]∈ℝ^p.
We will assume that the data-matrix 𝐙 is column-wise centered and scaled.
The dependent variable of interest is contained in the vector y and may designate either a continuous, a categorical (binary, nominal, or ordinal), a count, or a survival-type outcome.
Interest then lies in predicting or classifying y using information contained in 𝐙.
The proposed pipeline is given in the next section.




 §.§ Pipeline overview

In general, the radiomic features will be highly collinear.
Moreover, in many clinical situations the data setting will be high-dimensional in the sense that p > n.
Hence, any standard matrix summarization of the independent data, such as the correlation matrix, will be either ill-behaved or singular.
The first step (i) is then to obtain a regularized estimator of the correlation matrix that is well-behaved.
Next, (ii) this regularized correlation matrix is the basis for factor-analytic data compression.
That is, the original feature-space is projected onto a (much) lower-dimensional latent feature-space.
Subsequently, (iii) factor-scores are obtained, expressing the score each person obtains on each latent feature.
These latent features can be understood as radiomic meta-features.
These factor-scores can then (iv) be used as predictors in standard (generalized linear) predictive modeling.
Below, each step is explained in more detail.



  §.§.§ Regularized correlation matrix estimation

The sample correlation matrix 𝐑 = 𝐙^T𝐙/(n-1) will likely contain redundant information, in the sense that some entries will approach perfect (negative) correlation.
This implies that the information contained in one feature is almost completely represented by another feature.
Hence, only one of the redundant features needs to be retained.
Algorithm 1 from Section <ref> of the Supplementary Material (SM) contains a procedure for removing the minimal number of redundant features under absolute marginal correlation threshold τ.
In practice, we recommend setting τ∈ [.9,.95].


Denote the correlation matrix after redundancy filtering by 𝐑̇∈ℝ^p^*× p^*, where p^* is the number of features that passes filtering.
This matrix will still be ill-behaved or singular.
Hence, we will employ a penalized ML representation of 𝐑̇<cit.>:

    𝐑̇(ϑ) = (1 - ϑ)𝐑̇ + ϑ𝐈_p^*

where ϑ denotes a strictly positive penalty parameter.
This estimator is always positive definite, has well-defined limiting behavior, and is asymptotically consistent w.r.t. 𝐑̇ (see Section <ref> of the SM).
For appropriate choices of ϑ, this estimator is also well-conditioned.

Choosing the optimal value of ϑ is done in a data-driven way using K-fold cross-validation (CV).
The K-CV procedure is computationally efficient when noticing that (i) the estimator in (<ref>) is rotation equivariant and (ii) one may employ a root-finding procedure.
This implies that an optimal value ϑ^ of ϑ can be found in the (worst-case) computational complexity of a single spectral decomposition: 𝒪(p^*3).
See Section <ref> of the SM for details.
Then, 𝐑̇(ϑ^) is used as the basis for the next step.



  §.§.§ Factor analytic data compression

FA assumes that observed variables can be grouped on the basis of their covariation or correlation into a lower-dimensional linear combination of latent features.
Let z_i^* denote the radiomic observation vector that remains after removing the features implicated in redundancy filtering.
Then this linear combination can be expressed as:

    z_i^* = Λξ_i + ϵ_i,

where Λ is a (p^*× m)-dimensional matrix of factor loadings in which each element λ_jk is the loading of the jth variable on the kth factor, j=1,…,p^*, k=1,…,m, and where ϵ_i denotes the error measurements for person i.
The ξ_i then represent realizations of a latent variable of dimension m, with m < p^*, whose elements are referred to as common factors.
Under the assumptions on the model (see Section <ref> of the SM) z_i^*∼𝒩_p^*(0, ΛΛ^T + Ψ).
Hence, the correlation matrix among the (remaining) radiomic features is implied to be decomposable into common components (ΛΛ^T) and unique components (𝔼(ϵ_iϵ_i^T) = Ψ≡(ψ_11,…, ψ_pp)).
The overarching assumption is thus that the population correlation matrix Σ_R equals the model-implied population correlation Σ(Θ) = ΛΛ^T + Ψ, where Θ = {Λ,Ψ}.
From a likelihood perspective, the sample correlation matrix 𝐑̇ is a sufficient statistic (see Section <ref> of the SM).
This insight has led to the covariance structure modeling approach that would use the following discrepancy function to find estimate Θ̂<cit.>:

    F[Σ(Θ);𝐑̇] = ln|Σ(Θ)| + tr[𝐑̇Σ(Θ)^-1] - ln|𝐑̇| - p^*.

Minimizing (<ref>) requires a positive definite and, ideally, well-behaved sample correlation matrix.
Hence, we propose to replace 𝐑̇ by 𝐑̇(ϑ^) in (<ref>) in order to obtain estimate Θ̂ = {Λ̂,Ψ̂}.
The algorithm for minimizing (<ref>) based on the ML normal equations ∂ F[Σ(Θ);𝐑̇(ϑ^)]/∂Λ = 0 and ∂ F[Σ(Θ);𝐑̇(ϑ^)]/∂Ψ = 0 is the Fletcher-Powell algorithm <cit.> introduced into FA by <cit.>.
This algorithm is widely used in many standard platforms performing FA.


Imperative is the selection of a value for m.

From the model-implied population correlation matrix and the assumptions on the factor model it follows that:

    Σ(Θ) - Ψ = ΛΛ^T,

meaning that the reduced population correlation matrix is Gramian and of rank m<cit.>.
The number of common factors can then be determined by assessing the rank of Σ(Θ) - Ψ when replacing its constituents with appropriate estimates.
We use 𝐑̇(ϑ^) as the sampling counterpart to the population correlation matrix and 𝐈_p as a conservative estimate of the unique variance matrix.
Let d(𝐀)_j denote the jth eigenvalue of the matrix 𝐀.
Some ready algebra on 𝐑̇(ϑ^) - 𝐈_p will then show that its eigenvalues are of the form (1 - ϑ^)[d(𝐑̇)_j - 1].
The basic decision rule is then to choose an optimal value of m, say m̃, by:

    m̃ := (A),      A≡{j:(1 - ϑ^)[d(𝐑̇)_j - 1] > 0}.

Hence, we are simply determining the cardinality of the set of positive eigenvalues of 𝐑̇(ϑ^) - 𝐈_p.
From the perspective of 𝐑̇(ϑ^) - 𝐈_p, all variance is considered unique or error variance, implying that the average eigenvalue is 0.
Thus, a positive eigenvalue (1 - ϑ^)[d(𝐑̇)_j - 1] indicates a latent factor dominating information content, as its contribution to variance-explanation is above and beyond mere unique variance.
We then retain all such factors.
This approach concurs with the Guttman-Kaiser rule <cit.>.
Section <ref> of the SM contains additional information on the approach to dimensionality selection.
Section <ref> of the SM contains an extensive comparative simulation indicating that the decision rule in (<ref>) can reliably function as an upper-bound estimate to the true generating latent dimension m.
We suggest that this procedure always be accompanied by substantive considerations (see Section <ref> of the SM).
Now let Θ̂_m̃ = {Λ̂,Ψ̂}_m̃ denote the ML estimate of the parameters under m̃ factors
In the remainder we will refer to the elements of Θ̂_m̃ as Λ̂ and Ψ̂ in order to avoid notational clutter.

The orthogonal factor model considered copes with an inherent indeterminacy.
Consider an arbitrary orthogonal matrix 𝐇∈ℝ^m × m.
Then Σ(Θ_m) = ΛΛ^T + Ψ = (Λ 𝐇)(Λ 𝐇)^T + Ψ.
Any method of estimation under m factors thus requires m(m - 1)/2 restrictions on Λ.
In the ML-procedure this is usually achieved by requiring that Λ^𝐓Ψ^-1Λ is a diagonal matrix and that its diagonal elements are ordered.
These restrictions are convenient from an estimation perspective, but do not necessarily carry substantive meaning.
Hence, after estimation, rotational mappings may be employed that satisfy certain criteria for interpretation.
We opt, for our solution under m̃ factors, for a rotation to orthogonal simple structure <cit.>.
We do so by maximizing the normalized Varimax criterion (see Section <ref> of the SM).
The Varimax-rotated solution, Λ̂_V, is conducive in the next step.



  §.§.§ Obtaining factor scores

After projection of the original variable-space onto the lower-dimensional factor-space, we desire factor scores:
the score each individual obtains on each of the latent factors.
There are several methods for obtaining such scores.
Here, we regress Ξ∈ℝ^n ×m̃ on 𝐙^* to obtain <cit.>:

    Ξ̂ = 𝐙^*Ψ̂^-1Λ̂_V( 𝐈_m̃ + Λ̂_V^TΨ̂^-1Λ̂_V)^-1.

Note that these scores are orthogonal by construction.
They will be used as predictors in a prediction rule.



  §.§.§ Prediction rules

The original variable-dimension is now projected onto a lower-dimensional space.
In most situations m̃ < n.
Hence, we can estimate and evaluate any low-dimensional prediction rule.
We focus on regression-based prediction such as Cox regression for survival responses and the general linear model for categorical and linear responses, with:

    g{𝔼(y)} = Ξ̂β,

where Ξ̂β denotes the predictor and where g denotes a link function.
Note that our approach allows the usage of standard estimation (such as ML) and evaluation (such as ROC) methods.





§ RESULTS


 §.§ Data

We have two data settings, both with a right-censored survival (time to event) outcome.
The response concerns realizations T_i of a positive random variable 𝒯 representing the time from starting point t = 0 to an event of interest.
That is, the response for person i is (T̃_i,δ_i).
In this response T̃_i represents min(T_i,C_i), the minimum of event time T_i and censoring time C_i.
Then δ_i = I{T_i≤ C_i}, a status indicator that is 1 when the event of interest has occurred during the study-time and 0 when the event of interest did not occur during the study-time or when the subject was lost to follow-up (right-censoring).
At t = 0 predictor (prognostic) features z_i are available which, in our case, represent p radiomic features.
The general aim would then be to predict survival probabilities based on the prognostic features available at baseline.

The first data setting concerns oropharyngeal squamous cell carcinoma (OPSCC), a common malignancy affecting the oropharynx.
The event of interest is death.
It is a validation setting in the sense that independent test and validation sets are available.
The test set comprises n = 89 malignancies and 47 events and originated at the Amsterdam University medical centers, location VUmc, Amsterdam, the Netherlands.
The validation set comprises n = 56 malignancies and 28 events and originated at the University Medical Center Utrecht, the Netherlands.
The feature-dimension concerns p = 89 radiomic features extracted from axial T1-weighted MRI scans.
Additional information on data acquisition and feature extraction can be found in <cit.> as well as Section <ref> of the SM.

The second data setting concerns primary (hypopharyngeal, oropharyngeal, and laryngeal) tumors in head and neck squamous cell carcinoma (HNSCC) and the event of interest is again death.
In this case only one data set is available for model building and model evaluation.
We have n = 174 primary tumors and counted 55 events.
The feature-dimension concerns p = 432 radiomic features extracted from low-dose ^18F-FDG-PET/CT based images.
Additional information on data acquisition and radiomic feature extraction for this setting can be found in <cit.> as well as Section <ref> of the SM.




 §.§ Model comparison

In both settings we take interest in the predictive prowess of our pipeline.
The (training) data is used to obtain the projected factor scores.
These can then be directly used as the predictors in the Cox proportional hazards model <cit.>.
Note that the survival times are not used to guide the construction of the factor scores.
Hence, the projection approach is fully separated from the fitting of the model.
We take interest in comparing the pipeline approach with several other state-of-the-art ensemble and regularized variable selection methods.

For ensemble methods we consider survival forests.
These are nonparametric learning methods where the predictor is an ensemble formed by combining many survival decision trees.
There are, in general, 2 different types of survival forest and we will consider both: random survival forests <cit.> and conditional inference forests for survival analysis or conditional survival forests for short <cit.>.
They differ in how the ensemble is constructed.
RSFs construct the ensemble “by aggregating tree-based Nelson-Aalen estimators" <cit.>.
CSFs, on the other hand, use weighted Kaplan-Meier (KM) estimates, which tend to put more weight on terminal nodes with large numbers of subjects at risk <cit.>.
<cit.> has previously found that, in a radiomic prediction setting, RFs perform well in comparison to other methods, including PCA and partial least squares regression (which bear some of the flavor of our pipeline).
The main tuning parameter for RFs is the number of trees to grow, which we set to 1,000 for both the RSF and CSF.
We set no restriction on the depth of the trees.
The remaining tuning parameters are used `factory fresh' (i.e., default settings) as implemented in the <cit.> and <cit.> packages in <cit.>.

For the variable selection methods we consider a componentwise likelihood-based offset boosting approach for the Cox proportional hazards model <cit.>.
This can be seen as a regularized regression technique.
It starts with the null model and repeatedly fits (boosting) a base-learner in order to arrive at an adequate model.
At each boosting step only a single regression coefficient is updated.
The optimal updates in each step are determined by an approximation to (the maximization of) the penalized partial log-likelihood for Cox models.
As the non-updated regression coefficients are left at 0 the approach effectively results in variable selection and a sparse final model.
More information can be found in <cit.> and <cit.>.
This Cox boosting approach was previously shown to outperform the Lasso in survival prediction settings by <cit.> and <cit.>.
The method relies on two tuning parameters: the shrinkage parameter in the penalized partial likelihood and the number of boosting steps.
The former can be chosen rather coarsely while the latter must not be too small or too large.
Both parameters are set according to the `factory fresh' settings as implemented in the  package <cit.> (implying that the optimal number of boosting steps is chosen by 10-fold CV).

Hence, we are considering a = 1, …, A = 4 models.
A cox model where the predictors are the projected latent metafeatures ξ_i (the pipeline approach), two types of RF and a Cox boosting model.
The latter three use the z_i as predictors.
Let p_i be a generic designation for either ξ_i or z_i.
Now, let π̂^a(t|p_i) denote the predicted survival probability for individual i at time t under model a.
The (ensemble) survival functions π^a(t|p_i) for Cox regression and both forest-type models can be found, e.g., in <cit.>.
They are central to the model evaluation approach discussed next.




 §.§ Model evaluation

For model evaluation we focus on prediction error through the time-dependent Brier score <cit.>.
Let ỹ_i(t) = I{T̃_i≥ t} denote the observed survival status of subject i at time t.
The empirical Brier score may then be seen as a mean square error of prediction when π̂^a(t|p_i) ∈ [0,1] is taken to be a prediction of the survival status ỹ_i(t) ∈{0,1}<cit.>.
Hence, it is a measure of inaccuracy.
See <cit.> for reasoning regarding the appropriateness of the Brier score over the c-index in survival settings.
Section <ref> of the SM contains details on the calculation of the (integrated) Brier score in both our external (data setting 1) and internal (data setting 2) validation settings.

The models under consideration will then be compared w.r.t. (integrated) Brier scores in the following manner.
First, they will be compared to each other noting that, as a measure of inaccuracy, lower Brier scores are better.
Second, they will be compared w.r.t. their respective apparent errors.
The apparent error is the prediction error obtained when the data used from training are re-used for validation, which will result in optimistic prediction errors.
The smaller the gap between apparent and (cross-)validated error, the more stable a model (approach) can be considered to be.
Third, they will be compared to the benchmark value of .25, which corresponds to predicting 50% risk constantly over time <cit.>.
One desires a model to do better than chance over the lion's share of t ∈ [0, τ], where τ≤max(T_i).
We take τ to be the median follow-up time <cit.>.
Fourth, they will be compared to the prediction error based on the KM prediction rule, which ignores all feature information and as such acts as a null model <cit.>.
One desires a model to do considerably better than the null model.
These considerations are all incorporated in prediction error curves, where the time-dependent Brier scores are plotted against time.
Section <ref> shows how the (integrated) Brier score can be used to calculate a measure of (overall) explained residual variation: R^2.
The higher the R^2 measure, the better, giving a fifth ground for comparison of the models under consideration.




 §.§ External validation performance

The redundancy filtering algorithm was run on the training data with τ set to .95 and it retained p^* = 51 out of the original p = 89 features.
Subsequently, using 5-fold CV of the log-likelihood function, the optimal penalty-value for the regularized correlation matrix was determined (Step 1).
For the factor-analytic data compression, m̃ was set to 7.
Subsequently, the (rotated) ML estimates of the factor loadings and error covariance matrices were determined for the training data (Step 2), say Λ̂_V and Ψ̂.
Factor scores were then determined (Step 3) as ξ̂_i = (𝐈_7 + Λ̂^T_VΨ̂^-1Λ̂_V)^-1Λ̂^T_VΨ̂^-1z_i^*.
Let D be a data (subset) indicator and let D_T indicate the training set while D_V indicates the external validation set.
Using i ∈ D_T will then give the factor scores for the training data with which apparent performance can be assessed, while using i ∈ D_V gives the factor scores for the validation data (based on the factor solution of the training data) with which validation performance is assessed.
Section <ref> of the SM gives additional information on Steps 1–3.

Table <ref> gives an overview of the apparent and validated integrated Brier scores and explained residual variations for the proposed approach as well as the models mentioned in Section <ref>.
All methods perform better than chance and the KM reference model.
The  approach has the lowest validation prediction error.
This translates into the  approach having the highest integrated R^2.
Its performance is also amongst the most stable performances as indicated by the relatively small difference between its apparent and validation prediction errors.
Figure <ref> in Section <ref> of the SM visualizes these results.
Especially the RSF behaves unstably, with great apparent performance but severely diminished validation performance.



The  approach is subject to the strictest evaluation of validation performance: not only are the validation data evaluated at the Cox parameter-estimates stemming from the training sample, also the validation data projection is informed by the training-based factor solution.
In determining the new factor scores through the training factor-solution, we are recalibrating the positioning of the validation observations w.r.t. the latent dimensions determined in the training.
In this respect one could argue that, with shifts in latent-trait positioning, one should also recalibrate the Cox parameter-estimates.
If we would do so the recalibrated integrated prediction error of the  approach would drop to .134 while its integrated R^2 would raise to .2 (Table <ref>).
For comparison purposes we also recalibrated the one other approach that allows for ready recalibration: Cox boosting.
The validation data were subsetted to retain those features selected by the boosting approach in the training phase.
Subsequently a Cox proportional hazards model was fitted on these data.
Its performance was then assessed on again the validation data.
Note that the feature-subsetting in the Cox Boosting recalibration makes use of the outcome information from the training set, whereas the factor-projection is fully independent from the outcome measurements.
Hence, this gives the Cox boosting recalibration a certain advantage.
Nevertheless, the recalibrated  approach retains its superior performance over the recalibrated boosting approach (Table <ref>).



 §.§ Internal validation performance

The redundancy filtering algorithm in the second data setting was also run with τ set to .95.
It retained p^* = 124 out of the original p = 432 features.
For the factor-analytic data compression, m̃ was set to 8.
Subsequently, the factor scores were determined based on the ML factor solution.
Section <ref> of the SM gives additional information on Steps 1–3 in this second data setting.

Table <ref> gives an overview of the apparent and cross-validated integrated Brier scores and explained residual variations.
These results were obtained with K = 5 and B = 500 (see Section <ref> of the SM).
Again, all methods perform better than chance and the KM reference model.
Also, again the  approach has the lowest validation prediction error and the most stable performance (smallest difference between its apparent and validation prediction errors) which again translates into the highest (and most stable) integrated R^2.
Figure <ref> visualizes these results.



One could argue that the performance of  vis-à-vis the other approaches is due to the redundancy filtering, which removed many features.
To assess this the analysis was repeated with the RSF, CSF, and Cox boosting approaches based on the same redundancy-filtered feature-set as the  approach.
Section <ref> of the SM contains the results which are qualitatively the same as above: The most stable approach resulting in the lowest prediction error remains the proposed  approach.





§ DISCUSSION

We present a new approach for classification and prediction with radiomic data.
A 2-stage ML factor-analytic projection is proposed that summarizes the high-dimensional collinear feature-space in a low-dimensional orthogonal meta-feature space.
This enables standard down-stream estimation and evaluation of radiomic-feature-based classification and prediction models.
Moreover, it results in stable prediction models by directly dealing with the strong collinearity problem characteristic of radiomic data.
This approach then outperforms, in terms of (the stability of) prediction error, the most popular current approaches in radiomic data analysis.
While we have focussed on survival settings in our examples, we stress that any prediction model can use the meta-features obtained by our projection.

We also stress the distinct advantages of our approach.
First, as the meta-features are projections of amalgams of individual features, they are more generalizable within the radiomic context.
This generalizability will—for a given set of latent projections—improve as the feature-dimension grows.
Second, the projection enables the prediction itself to be performed in a low-dimensional context.
Besides stability, this ensures low computation times and allows, in contrast to the other methods, for the quantification of uncertainty (such as through confidence intervals on the parameter estimates).
Third, as the projection is fully independent from the outcome measurements it is possible to recalibrate the parameter estimates without overfitting.
This may lead to considerable further improvements in performance.
Fourth, the meta-features are often interpretable.
This gives the proposed approach an advantage over ensemble-type methods, as interpretability in the clinical setting is often desired.

We see several inroads for further research and improvements.
First, we have worked under the assumption that the (standardized) radiomic features are normally distributed.
In case this assumption is untenable one could base 𝐑, and subsequently 𝐑̇, on either a rank-type correlation or on a semiparametric Gaussian copula model <cit.>.
Second, the penalized correlation matrix may be extended with substantive prior information.
In a sense, the estimator 𝐑̇(ϑ) is a Bayesian estimator, balancing the ill-behaved 𝐑̇ with the well-behaved 𝐈_p^*.
Instead of 𝐈_p^* one could also use more informative matrices (on the correlation scale) representing more elaborate prior information, possibly stemming from previous radiomic studies.
These considerations might further improve the predictive performance of the proposed approach.

We also note several immediate extensions to the approach.
While we focussed on using all extracted meta-features in the prediction process one could indeed also embark on feature-selection amongst these meta-features.
Moreover, we have limited ourselves to using only radiomic meta-features for prediction.
However, one can immediately combine these features with other clinical variables to possibly further improve predictive performance.
One would then be evaluating

    g{𝔼(y)} = Ξ̂β + 𝐂β'.

Lastly, the approach could also be of interest for other data settings in which the features are subject to a natural grouping, such as, for example, in metabolomics.













 §.§ Acknowledgements

We thank Ruud Brakenhoff and René Leemans of the department of Otolaryngology/Head & Neck Surgery of the Amsterdam University medical centers as well as Remco de Bree of the department of Head and Neck Surgical Oncology of the University Medical Center Utrecht, for sharing their data with us.
We also thank the study participants.
The data cohorts in the external validation study have been obtained in accordance with the Declaration of Helsinki and were approved by the Medical Ethics Committees of the associated medical centers.
All patients provided a signed informed consent.
The trials for the internal validation study have also been performed in accordance with the Declaration of Helsinki, approved by the Medical Ethics Committee of the VU University Medical Center.
A written informed consent was waived for the retrospective cohort (103 patients, reference: 2016.498), whereas for the prospective cohort (72 patients, reference: 2013.191) a written informed consent was obtained from all patients.




 §.§ Funding

This work was financially supported by the Netherlands Organisation for Health Research and Development (ZonMw), grant 10-10400-98-14002.



  §.§.§ Conflict of interest:

None declared.





[Radiomics_Main]












tocsectionSupplementary Material
 SUPPLEMENTARY MATERIAL
 

Stable Prediction with Radiomics Data




This supplement contains text, figures, and  code in support of the Main Text.
Section <ref> contains additional technical information on various steps of the proposed pipeline.
Section <ref> contains details on the evaluation of predictive performance.
Section <ref> contains a description of simulation efforts regarding the selection of an appropriate latent-dimensionality.
Sections <ref> and <ref> contain additional information on the data used in the Main Text.
 code with which all analyzes and simulations can be repeated are found on <https://github.com/CFWP/FMradio>.





§ PIPELINE DETAILS

Section <ref> contains the redundancy filtering algorithm.
Section <ref> gives details on the regularized estimator for the correlation matrix.
Subsequently, Section <ref> explains how the regularization parameter is chosen.
Section <ref> then delves into the assumptions underlying the factor analysis model.
Sections <ref> and <ref> show how to select the dimension of the latent vector and how to obtain factor scores.
Most of the development will be in the general notation p and m for the feature and latent dimensions, respectively.
These may be replaced, without loss of generality, by p^* and m̃.





 §.§ Redundancy filtering

Algorithm <ref> contains pseudocode for redundancy filtering of a correlation matrix.



Consider the following toy example explaining the workings of Algorithm <ref>.


Say we have the following 4 × 4 correlation matrix on variables A, B, C, and D:

    ccccc    A     B     C     D
    c[cccc]
      A     1     .95     .95     .30 
    
      B     .95     1     .30     .30 
    
      C     .95     .30     1     .95 
    
      D     .30     .30     .95     1 
    .

Assume that we consider absolute correlations ≥ .95 to indicate redundancy.
It is clear that, if we would blindly remove all columns and corresponding rows in which an (absolute) correlation appears that abides the redundancy criterion, we would be left with no matrix at all.
Also, if we would treat this matrix as representing a collinear block, then choosing one representative feature would also imply a loss of information as no feature represents all other features in terms of information.
For example, if we would choose feature A as the representative feature then we would loose much information regarding feature D: A is highly collinear with features B and C, but not with D.

Hence, we apply the RF procedure of Algorithm <ref>.
It starts by evaluating, for each row, the number of times the corresponding feature exhibits an absolute correlation ≥ .95.
This information is collected in the vector v.
For our correlation matrix this evaluation would amount to:

    ccccccc    A     B     C     D       v[j] ∑_j'1_{|𝐑[j,j']| ≥τ}
    c[cccc]cc
      A     1       .95     .95     .30          v[1] = 3
    
      B     .95     1       .30     .30          v[2] = 2
    
      C     .95     .30     1       .95          v[3] = 3
    
      D     .30     .30     .95       1          v[4] = 2
    .

The feature(s) that have the most absolute correlations exceeding the threshold can be thought of a having their information represented by the most alternative features.
Hence, such features are natural candidates for removal.
The algorithm then evaluates for which indices (corresponding to row and column numbers) v attains a maximum: (v = max(v)).
In this case (v = max(v)) = [1,3], as both feature A and feature C have 3 (absolute) correlations exceeding our threshold.
In the situation where (v = max(v)) is a vector of length greater than 1, the first index in this vector is retained, i.e.: c (v = max(v))[1].
Thus, at current, c = 1.
Then, as at this point, max(v) = 3 > 2 (i) an updated 𝐑 is formed by removing the first row and column (𝐑𝐑[-1,-1]), corresponding to feature A, (ii) v is emptied, and (iii) a new round of evaluation is started.

The new round of evaluation starts with the updated correlation matrix.
The vector v is now evaluated as:

    cccccc    B     C     D       v[j] ∑_j'1_{|𝐑[j,j']| ≥τ}
    c[ccc]cc
      B     1       .30     .30          v[1] = 1
    
      C     .30     1       .95          v[2] = 2
    
      D     .30     .95       1          v[3] = 2
    .

Now, (v = max(v)) = [2,3] and c (v = max(v))[1] = 2.
As max(v) = 2 (i) an updated 𝐑 is formed by removing the second row and column (𝐑𝐑[-2,-2]), corresponding to feature C, (ii) v is emptied, and (iii) a new round of evaluation is started.
Again, the subsequent round of evaluation starts with the updated correlation matrix.
This matrix and the corresponding vector v are now given as:

    ccccc    B     D       v[j] ∑_j'1_{|𝐑[j,j']| ≥τ}
    c[cc]cc
      B     1        .30          v[1] = 1
    
      D     .30        1          v[2] = 1
    .

As now max(v) < 2, we have that, logically, only the diagonal elements of the remaining matrix exceed the threshold.
Hence, the evaluation stops and the correlation matrix in (<ref>) is returned as the final, redundancy filtered correlation matrix.

In this final matrix retained feature B also represents most of the information contained in removed feature A, and retained feature D also represents most of the information contained in removed feature C.
Thus, in a sense, the final matrix represents all original features.
Note that, while c (v = max(v))[1] may be termed arbitrary in case (v = max(v)) is a vector of length greater than 1, alternative choices would lead to retainment of roughly the same total information.
For example, say we would have used c (v = max(v))[2] in each case where ((v = max(v))) > 1.
This would have led to the removal of feature C in the first step and feature B in the second, giving the final correlation matrix:

    ccc    A     D
    c[cc]
      A     1      .30 
    
      D     .30      1 
    .

In this alternative final matrix removed feature B is represented by retained feature A and removed feature C is represented by retained features A and D.


The algorithm can be expanded by including the sum of absolute correlations exceeding the threshold in the assignment of index c whenever ((v = max(v))) > 1.
This would imply bookkeeping of an extra vector.





 §.§ Regularized correlation matrix estimation

Here we give some details on the regularized correlation matrix estimation used in the Main Text.
Section <ref> explains how the estimator may be understood as a penalized maximum likelihood estimator.
Section <ref> details on selecting a value for the penalty parameter.



  §.§.§ Penalized maximum likelihood

We are considering sample correlation matrices.
Let x_i^T≡[x_i1,…,x_ip]∈ℝ^p define a realization of the random vector X_i^T≡[X_i1,…,X_ip]∈ℝ^p.
Then the sample covariance matrix can be obtained as

    𝐒 = (n - 1)^-1∑_i=1^n(x_i - x̅)(x_i - x̅)^T,

where x̅ = n^-1∑_i=1^nx_i.
The correlation matrix can then be obtained directly from 𝐒 or from the standardized sample realizations (𝐒∘𝐈_p)^-1/2(x_i - x̅) ≡z_i as:

    𝐑 = (𝐒∘𝐈_p)^-1/2𝐒(𝐒∘𝐈_p)^-1/2 =
        (n-1)^-1∑_i=1^nz_iz_i^T.

If x_i ∼𝒩_p(μ, Σ) then, under the stated data transformation, z_i ∼𝒩_p(0, Σ_R), where μ, Σ, and Σ_R respectively refer to the population mean, the population covariance matrix, and the population correlation matrix.
It is well-known that, if all realizations x_i ∼𝒩_p(μ, Σ), then 𝐒 has a Wishart distribution with scale matrix Σ/(n-1) and n-1 degrees of freedom, i.e., 𝐒∼𝒲_p(n-1,Σ/(n-1)).
While the exact distribution of 𝐑 is unknown, it can be well-approximated by the Wishart distribution <cit.>.

Hence, we assume 𝐑∼𝒲_p(n-1,Σ_R /(n-1)).
The negative log-likelihood of 𝐑 (up to proportionality) is then:

    ℒ(Σ_R; 𝐑) ∝ln|Σ_R| + [𝐑Σ_R^-1].

If we would minimize (<ref>) w.r.t. Σ_R we would find that 𝐑 is the maximum likelihood (ML) estimator of Σ̂_R.
This estimator will be ill-conditioned when p → n - 1 or in the presence of strong multicollinearity and it will be singular when p > n-1.
These issues can be resolved if we adapt (<ref>) by (a) rescaling 𝐑 with (1 - ϑ) where ϑ∈ [0,1], and (b) adding the penalty ϑ(Σ_R)^-1.
We then have a penalized log-likelihood

    ℒ^p(Σ_R; 𝐑, ϑ) ∝ln|Σ_R| + [(1-ϑ)𝐑Σ_R^-1] + ϑ(Σ_R^-1),

whose minimization w.r.t. Σ_R conditional on ϑ results in:

    Σ̂_R(ϑ) :=
        (1 - ϑ)𝐑 + ϑ𝐈_p ≡𝐑(ϑ),

a convex combination of 𝐑 and 𝐈_p.

The estimator (<ref>) has certain appreciable properties.
First, it is on the correlation scale (i.e., its diagonal elements are always unity).
Second, whenever the penalty parameter is strictly positive the resulting estimate is positive definite (p.d.).
It can also be seen that the right-hand and left-hand limits are the sample correlation matrix and the identity matrix, respectively.
This latter observation enables the understanding of the estimator: It balances the unbiased but highly variable matrix 𝐑 with the biased but perfectly stable identity matrix.
We note that <cit.> arrives at the same correlation matrix estimator via a somewhat different road.




  §.§.§ Choosing the penalty parameter

From the above we get that 𝐑(ϑ) is always p.d. when the penalty parameter is strictly positive, irrespective of the ratio p/(n - 1).
It is important to make a well-informed choice regarding its value.
Choosing it too small can lead to an ill-conditioned estimate when p ⪆ n - 1, meaning that the estimate is unstable and that its use in matrix operations may lead to large error propagation <cit.>.
Choosing it too large, however, may suppress relevant data signal.
Here, we choose a data-driven approach that seeks loss efficiency.
That is, we will use K-fold cross-validation (kCV) of the log-likelihood function which, asymptotically, can be explained in terms of minimizing Kullback-Leibler divergence.
The kCV score for 𝐑(ϑ) based on a fixed choice of ϑ can be stated as:

    φ^K(ϑ) := 1/K∑_k = 1^Kn_k {ln|𝐑(ϑ)_ k| + [𝐑_k(𝐑(ϑ)_ k)^-1]},

where n_k denotes the sample size of subset k, for k = 1, …, K disjoint subsets and where 𝐑_k denotes the sample correlation matrix based on subset k, while 𝐑(ϑ)_ k denotes the regularized correlation matrix based on all samples not in k.
Using (<ref>), we choose an optimal ϑ, denoted ϑ^ such that:

    ϑ^ := _ϑ∈ (0,1]φ^K(ϑ).

The conditioning of the resulting estimate 𝐑(ϑ^) can be evaluated along the lines laid out in <cit.>.
The procedure is fast for practically accepted choices for K, such as K ∈{5,10}, when realizing that the optimization problem in (<ref>) can be combined with root-finding procedures such as the Brent algorithm <cit.>.





 §.§ Factor analytic modeling

The assumptions considered for the factor model are <cit.>:
(i) z_iz_i',∀ i≠ i'; (ii) rank(Λ)=m; (iii) ϵ_i∼𝒩_p(0, Ψ), with Ψ≡[ψ_11,…, ψ_pp], and ψ_jj>0,∀ j; (iv) ξ_i∼𝒩_m(0, 𝐈_m); and (v) ξ_iϵ_i',∀ i,i'.
The likelihood for the observations conditional on the realization of Ξ can then be expressed as:

    L(Λ, Ξ, Ψ; 𝐙) =
      ∏_i=1^nf(z_i|Λ,ξ_i,Ψ)
      =
      ∏_i=1^n(2π)^-p/2|Ψ|^-1/2exp{-1/2ϵ_i^TΨ^-1ϵ_i},

where
ϵ_i=z_i-Λξ_i.
The likelihood of the observed data can be obtained by marginalizing over ξ_i:

    L(Λ, Ψ; 𝐙)   =
      ∏_i=1^n∫ f(z_i|Λ,ξ_i,Ψ)
      g(ξ_i|𝐈_m) ∂ξ_i
       =∏_i=1^n(2π)^-p/2|ΛΛ^T+Ψ|^-1/2exp{-1/2z_i^T(ΛΛ^T+Ψ)^-1z_i}
       =(2π)^-np/2|ΛΛ^T+Ψ|^-n/2exp{-n/2[𝐑(ΛΛ^T+Ψ)^-1]
      }
      ,

giving that the factor decomposition constrains the correlation
structure of the z_i to be a function of Θ = {Λ,Ψ}, i.e.:

    Σ(Θ)=ΛΛ^T+Ψ.

Then, for existence (vi), generally (p-m)^2-p-m⩾0, simply stating that the number of nonredundant elements in the sample correlation matrix 𝐑 must be greater than or equal
to the number of freely estimable parameters in Σ, which places an upper bound on m.

Hence, under the assumptions of the factor model z_i∼𝒩_p(0, Σ(Θ)).
Moreover, from (<ref>) we understand that 𝐑 is a sufficient statistic.
Thus, we may assume that under the model 𝐑∼𝒲_p(n-1,Σ(Θ)/(n-1)).
The negative log-likelihood of 𝐑 (up to proportionality) under this correlation-structure model is then:

    ℒ(Σ(Θ); 𝐑) ∝ln|Σ(Θ)| + [𝐑Σ(Θ)^-1].
<cit.> then proposed to use the following discrepancy function to find an estimate Θ̂ = {Λ̂,Ψ̂} of Θ:

    F[Σ(Θ);𝐑] = ln|Σ(Θ)| + tr[𝐑Σ(Θ)^-1] - ln|𝐑| - p.

Minimizing (<ref>) is equivalent to minimizing (<ref>), but computationally more convenient.
Minimizing (<ref>) does require a p.d. and, ideally, well-conditioned sample correlation matrix for the estimate Θ̂ to be meaningful and stable.
Hence, we propose to replace the ill-behaved 𝐑 by the well-behaved 𝐑(ϑ^) in (<ref>) in order to obtain estimate Θ̂:

    Θ̂ := _Λ,Ψ F[Σ(Θ);𝐑(ϑ^)].

The algorithm for minimizing (<ref>) based on the ML normal equations ∂ F[Σ(Θ);𝐑(ϑ^)]/∂Λ = 0 and ∂ F[Σ(Θ);𝐑(ϑ^)]/∂Ψ = 0 is the Fletcher-Powell algorithm <cit.> introduced into FA by <cit.>.
This algorithm is widely used in many standard platforms performing FA.

For model determinacy (see Section 2.2.2. of the Main Text), the minimization is performed under the requirement that Λ^𝐓Ψ^-1Λ be diagonal with ordered diagonal elements.
This gives the canonical solution to the minimization problem which is convenient and efficient from an estimation perspective.
For interpretation purposes the resulting estimate Λ̂ may be subjected to a post-hoc rotations.
In our case this will be the normalized Varimax rotation <cit.> giving a rotation to orthogonal simple structure:

    Γ := max_𝐇{∑_k = 1^m̃∑_j = 1^p^*(Λ̂𝐇)_jk^4/ĉ_j^2 -
         1/p^*∑_k = 1^m̃[ ∑_j = 1^p^*(Λ̂𝐇)_jk^2/ĉ_j]^2
         : 𝐇^T𝐇 = 𝐈_m̃},

where ĉ_j = ∑_k = 1^m̃(Λ̂)_jk^2 denotes the retrieved communality: the amount of variance of item j explained by the latent features (also see Section <ref>).
Maximizing this argument comes down to solving quartic equations and is well-implemented in most standard packages and platforms.
The Varimax-rotated solution is then Λ̂Γ≡Λ̂_V.
In the remaining text Λ̂ may refer to the loadings matrix either before or after post-hoc rotation in order to avoid notational clutter.





 §.§ Latent dimensionality selection

The estimation process detailed in the previous section hinges upon a choice for the latent dimensionality m.
Here we give some details on the latent dimensionality selection approach used in the Main Text.
Section <ref> explains the basic decisional approach.
Section <ref> details on how this basic approach may be supplemented for additional decision support.



  §.§.§ Decision approach

From (<ref>) we get that the variance of the jth feature can be expressed as:

    ∑_k = 1^mλ_jk^2 + ψ_jj,

where ψ_jj represents the unique variance and where ∑_k = 1^mλ_jk^2≡ c_j represents the common variance for feature j.
The common variance c_j can be understood as the amount of variance of observed feature j that is explained by the m latent features.
In the literature c_j is known as the communality.
For the population correlation case, we have the consequence: c_j = 1 - ψ_jj.
This forms the basis of our decisional approach.

From a matrix perspective we have that:

    Σ(Θ) - Ψ = ΛΛ^T.

That is, the correlation matrix with communalities in the diagonal (the reduced correlation matrix) is Gramian and of rank m<cit.>.
The number of common factors can then be determined by assessing the rank of Σ(Θ) - Ψ when replacing its constituents with appropriate estimates.
We use 𝐑(ϑ^) as the sampling counterpart to the population correlation matrix.
We use 𝐈_p as a conservative estimate of the unique variance matrix.Hence, we then obtain 0 as a lower-bound estimate of the communalities, such that all variance is considered unique or error variance.
This information can be used to understand the decisional approach.

Some ready algebra will show that:

    𝐑(ϑ^) - 𝐈_p = (1 - ϑ^)(𝐑 - 𝐈_p),

giving the matrix we assume Gramian and whose rank we want to determine.
Now, let 𝐕𝐃(𝐑)𝐕^T be the spectral decomposition of 𝐑 with 𝐃(𝐑) denoting a diagonal matrix with the eigenvalues of 𝐑 on the diagonal and where 𝐕 denotes the matrix that contains the corresponding eigenvectors as columns.
Note that 𝐕𝐕^T = 𝐕^T𝐕 = 𝐈_p.
Then (<ref>) can be written as

    𝐕{(1 - ϑ^)[𝐃(𝐑) - 𝐈_p]}𝐕^T,

giving that the eigenvalues of interest are of the form (1 - ϑ^)[d(𝐑)_j - 1].
The basic decisional rule is then to choose an optimal value of m, say m̃, by:

    m̃ := (A),      A≡{j:(1 - ϑ^)[d(𝐑)_j - 1] > 0}.

Our communality estimate concurs with the weakest lower-bound estimate <cit.>, which corresponds to all variance being unique or error variance.
From this perspective, a positive eigenvalue (1 - ϑ^)[d(𝐑)_j - 1] indicates a latent factor whose contribution to variance-explanation is above and beyond mere unique variance.
We then retain all such factors.

The decisional rule in (<ref>) can be understood as a regularized version of Guttman's weakest lower-bound approach to minimum rank determination <cit.>.
Hence, in the remainder we will refer to this rule as the Guttman bound (GB).
It also concurs with a regularized version of what is known as the Kaiser-rule <cit.>.
While especially the Guttman approach has historically been used as a lower-bound estimate of the latent dimensionality, we consider the decisional rule in (<ref>) to give an upper-bound.
Section <ref> contains an extensive simulation study that shows that, in challenging situations, the rule (<ref>) provides a reliable upper-bound.
This section also explains why formal testing (such as likelihood ratio testing) is not available in the situations we consider.
In such situations the approach above gives a computationally inexpensive and reliable probe regarding the number of factors that dominate information content.




  §.§.§ Additional decision support

Here we describe tools and strategies for assessing the quality of the factor solution and the choice of latent dimensionality.
These revolve around the assessment of (i) factorability, (ii) extracted communalities, (iii) proportion of explained variance, and (iv) interpretability of extracted factors.
These assessments may inform if the result of the decisional rule above should be accepted or be treated as an upper-bound.

Factorability.
Factorability refers to the assessment if factor analysis is an appropriate tool for the data at hand.
Where appropriateness is taken to mean the ability to identify coherent common latent factors.
The basic premise of common factor analysis, from a conditional independence perspective, is that the observed features are independent given the common latent features, i.e.: z_jz_j'|{ξ_k}_k = 1^m,∀ j≠ j'.
Guttman <cit.> then proved that the inverse of the marginal (population and sample) correlation matrix should also be near diagonal for factor analysis to be appropriate.
Moreover, it necessarily approaches a diagonal matrix when the number of observed features grows to infinity.
Hence, a basic inspection of the off-diagonal elements of 𝐑(ϑ^)^-1 would be a basic evaluation of factorability in our case.
One practical option for such an inspection is to assess the Kaiser-Meyer-Olkin (KMO) measure of feature-sampling adequacy <cit.>.
For our setting the KMO index amounts to:

    ∑_j∑_j'[𝐑(ϑ^) - 𝐈_p]^2_jj'/∑_j∑_j'[𝐑(ϑ^) - 𝐈_p]^2_jj' + ∑_j∑_j'[𝐏(ϑ^) - 𝐈_p]^2_jj',

where 𝐏(ϑ^) = [𝐑(ϑ^)^-1∘𝐈_p]^-1/2𝐑(ϑ^)^-1[𝐑(ϑ^)^-1∘𝐈_p]^-1/2, the standardized regularized precision matrix.
The index basically compares the sizes of the off-diagonal entries of the regularized correlation matrix to the sizes of the off-diagonal entries of its scaled inverse (proportional to the partial correlation matrix) and takes values in [0,1].
In general the KMO index will be larger when the marginal correlations get larger and the partial correlations get smaller.
It will indeed equal unity when 𝐏(ϑ^) = 𝐈_p.
In general a KMO index between .9 and 1 is considered to be indicative of marvelous factorability <cit.>.

Comparison of communalities.
The GB tends to be very reliable when the communalities of the generating data mechanism are high, irrespective of the obtained sample size (see Section <ref> below).
Hence, it is of considerable interest to compare lower-bound estimates of the communalities to the extracted communalities.
The strongest lower-bound estimate of the communality for any feature j then is <cit.>:

    SMC_j = 1 - 1/[𝐑(ϑ^)^-1∘𝐈_p]_jj,

which is essentially the squared multiple correlation for predicting feature j from the remaining p - 1 features.
It is the best possible lower-bound estimate for the communality <cit.>.
When, overall, these are high, the generating data mechanism (from a factor analytic perspective) most likely has high communalities and, hence, the GB has reliable performance.
It is also of interest to juxtapose the SMC_j with the retrieved estimated communalities ĉ_j = ∑_k = 1^m̃(Λ̂)_jk^2 = (Λ̂Λ̂^T)_jj.
When the chosen value m̃ is sufficient then one would expect, for almost all j, that ĉ_j⪆SMC_j.
If this is not the case then one might have extracted too few factors.

Proportion of variance explained.
We might also look at the proportion of explained variance.
The proportion of variance explained by any factor k can be obtained as:

    v_k = (Λ̂^TΛ̂)_kk/p.

The proportion of variance explained by all factors is then ∑_k v_k.
We want the proportion of variance explained by all factors to be appreciable (say, in excess of 70%).
Moreover, one would want the proportion of variance explained by the m̃th factor in relation to the (m̃ - 1)th factor to be appreciable and the proportion of variance of the (m̃ + 1)th factor in relation to the m̃th factor to be negligible.
These are all qualitative judgements.
It is often informative to graph ∑_k v_k against the (possible) number of factors.
The point at which the graph flattens out is indicative of a formative number of latent factors.

Substantiveness of factors.
If a factor does not have at least 3 significant loadings (and thus indicators), than this factor is considered weak or redundant <cit.>.
Such weak factors should be considered for removal.
This can be evaluated by thresholding the loadings matrix:

    λ̂_jk^† =
    {[             0,  if |λ̂_jk| ≤ω;         λ̂_jk, if |λ̂_jk| > ω;                ].
    .

In practice ω is often set to .3.
Hence, it is important to assess the interpretability (from a meta-radiomic feature perspective) and practical significance of the extracted latent features.
Again, these are qualitative judgments.

Factor score determinacy.
The appropriateness of factor analysis and the quality of the factor extraction can also be probed through the squared multiple correlation between the observed features and the latent common factors. In our setting this squared multiple correlation for latent factor k is given as <cit.>:

    SMC_k^ξ  = (Λ̂^T𝐑(ϑ^)^-1Λ̂)_kk.

It indicates how well common factor k can be predicted by the observed features or, from another perspective, to what extent the factor scores (Section <ref>) are determinate.
The closer SMC_k^ξ is to unity, the more determinate the factor scores for latent factor k can be considered to be.
In practice a SMC_k^ξ≥ .9 would be considered (very) adequate.
Observing SMC_k^ξ that are relatively low might indicate that one has retained (a) weak latent factor(s).
We will turn to estimating factor scores next.





 §.§ Factor scores

We desire an estimate of the score each individual would obtain on each of the latent factors, i.e., we desire factor scores.
These will be obtained through the expectation of the conditional distribution ξ_i|z_i.
First, we use standard covariance algebra to obtain the joint distribution of z_i and ξ_i:

    [ z_i; ξ_i ]∼𝒩_p + m(0, Σ(Θ)_ξ^z),

with

    Σ(Θ)_ξ^z
        =
        [  Σ(Θ)_zz  Σ(Θ)_zξ; Σ(Θ)_ξ z  Σ(Θ)_ξξ ]
        =
        [ ΛΛ^T + Ψ        Λ;      Λ^T      𝐈_m ].

Here we indeed recognize the marginal distribution of the observed data as obtained via an alternative route in Section <ref>.
As the joint distribution is normal, a standard result states that the conditional expectation 𝔼(ξ_i|z_i) can be obtained as Σ(Θ)_ξ zΣ(Θ)_zz^-1z_i<cit.>.
The matrix Σ(Θ)_ξ zΣ(Θ)_zz^-1 is the matrix of (OLS) regression coefficients of ξ_i on z_i and Σ(Θ)_ξ zΣ(Θ)_zz^-1z_i is the regression function.
Hence, scores obtained through this route can be thought of as regression-type scores.
Now,

    Σ(Θ)_ξ zΣ(Θ)_zz^-1   = Λ^T(ΛΛ^T + Ψ)^-1
       = [𝐈_m - Λ^TΨ^-1Λ(𝐈_m + Λ^TΨ^-1Λ)^-1]Λ^TΨ^-1
       = (𝐈_m + Λ^TΨ^-1Λ)^-1Λ^TΨ^-1,

by the Woodbury matrix identity <cit.>.
Hence, 𝔼(ξ_i|z_i) = (𝐈_m + Λ^TΨ^-1Λ)^-1Λ^TΨ^-1z_i such that an estimate of obtained score ξ_i given an observed z_i amounts to:

    ξ̂_i = (𝐈_m + Λ̂^TΨ̂^-1Λ̂)^-1Λ̂^TΨ̂^-1z_i,

where Λ̂ and Ψ̂ are (possibly rotated) estimates of Λ and Ψ.
To collect all i = 1, …, n estimated scores from (<ref>) in the (n × m)-dimensional matrix Ξ̂ we use the observed data matrix 𝐙 and the (properties of the) matrix transpose to arrive at the expression given in Section 2.2.3 of the Main Text.

The scores obtained through (<ref>) are known as Thomson scores.
<cit.> obtained factor scores by minimizing 𝔼(𝐁z_i - ξ_i)^2.
This minimization leads to 𝐁 equalling (<ref>).
Hence, this estimator minimizes the mean squared prediction error.
There exist alternative consistent estimators of the factor scores.
<cit.> devised an unbiased estimator that, in practice, is very close to Thomson scores.
Under the orthogonal model the latent factors are orthogonal in the population and, hence, the Thomson and Bartlett-type factor scores will be near orthogonal in the sample.
<cit.> constructed an alternative estimator for the factor scores that enforces their orthogonality in the sample.
This estimator considers the ξ_i to be fixed rather than random.
For an overview of factor score estimation we confine by referring to <cit.>.





§ MODEL EVALUATION

In this section we provide details on the measures of model performance used in the Main Text.
For model evaluation we focus on prediction error through the time-dependent Brier score <cit.>.
Let ỹ_i(t) = I{T̃_i≥ t} denote the observed survival status of subject i at time t.
The empirical Brier score may then be seen as a mean square error of prediction when π̂^a(t|p_i) ∈ [0,1] is taken to be a prediction of the survival status ỹ_i(t) ∈{0,1}<cit.>.
Hence, it is a measure of inaccuracy.
Section <ref> provides details on the Brier-based evaluation of prediction error for both our internal and external validation settings.
Section <ref> explain how the Brier score may be used to calculate a measure of explained variation.



 §.§ Prediction error

Let D be a data (subset) indicator.
And say we have a training set D_T and validation set D_V available (data setting 1) of respective sample sizes n_T and n_V.
In this case the estimated Brier score is <cit.>:

    ℬ_V(t, π̂^a) = 1/n_V∑_i ∈ D_VŴ_i(t)[ỹ_i(t) - π̂_T^a(t|p_i)]^2,

where π̂_T^a denotes the survival function under model a based on training data D_T and where

    Ŵ_i(t) = [1-ỹ_i(t)]δ_i/Ĝ(T̃_i) + ỹ_i(t)/Ĝ(t),

denotes an inverse probability of censoring weighting scheme to cope with information loss due to censoring.
The Ĝ in (<ref>) denotes the marginal Kaplan-Meier (KM) estimate of the censoring distribution.

In case only one data set is available (as in data setting 2) we have to harness ourselves form overoptimistic prediction errors.
In such situations we will use a cross-validated prediction error.
Let the available data be randomly split into K roughly equally sized subsets k = 1, …, K.
The prediction rule π̂^a is then trained on D\ D_k (π̂^a_k) and the data in D_k are used as an internal validation set.
The cross-validated prediction error is then <cit.>:

    ℬ_CV(t, π̂^a) = 1/K∑_k = 1^K1/n_D_k∑_i ∈ D_kŴ_i(t)[ỹ_i(t) - π̂_k^a(t|p_i)]^2.

Note that (<ref>) is dependent on the choice of folds.
We mitigate this dependency by repeating the K-fold cross-validation (K-CV) many times and averaging the corresponding Brier scores over these repeats:

    ℬ̃_CV(t, π̂^a) = 1/B∑_b=1^Bℬ_CV(t, π̂^a)_b,

where ℬ_CV(t, π̂^a)_b represents the score (<ref>) under random K-foldage b.
In our applications we set K = 5 and B = 500.
Also, we make use of the facilities provided by the  package <cit.>, ensuring that all models under consideration are evaluated using the exact same folds.

Now, (<ref>) can be integrated for t ∈ [0, τ], with τ≤max(T_i), to arrive at the overall prediction error over the chosen time-period <cit.>:

    ℬ_a|CV^I≡1/τ∫_0^τℬ̃_CV(t, π̂^a)  ∂ t.

We take τ to be the median follow-up time <cit.>.
Note that (<ref>) can be integrated similarly to arrive at ℬ_a|V^I

 §.§ Explained residual variation

Let the survival probability under the KM null-model be denoted by π̂^0(t|p_i).
Also, let ℬ̃_CV(t, π̂^0) be the averaged cross-validated prediction error under this null model.
Then the additional explained residual variation for model a at time t vis-à-vis the KM reference model is <cit.>:

    R^2_a|CV(t) = 1 - ℬ̃_CV(t, π̂^a)/ℬ̃_CV(t, π̂^0).

To arrive at the overall explained residual variation we could use the integrated Brier scores ℬ_a|CV^I and ℬ_0|CV^I.
Note that analogous considerations can be given to the (integrated) Brier scores in the validation setting.





§ SIMULATIONS DIMENSIONALITY SELECTION

In this section we describe the numerical evaluation of the latent dimensionality selection procedure outlined in Section <ref>.
Section <ref> contains an overview of the general setup of the simulation study.
Section <ref> described the other methods of dimensionality selection to which the GB will be compared.
The results can then be found in Section <ref>.




 §.§ Setup

We seek to evaluate and compare the performance of the basic decision rule XXX in settings that emulate the complexities of radiomic data.
In particular we take interest in varying (i) the p/n ratio, (ii) the true value for m, (iii) the communality-strength, and (iv) the factor-structure complexity.

We will consider p ∈{100,200}.
For each of these feature-dimensions we will consider a true generating m ∈{5,12,20}.
For each of these possible combinations we consider various values for the generating communalities.
As stated previously, the communality c_j for feature j can be understood as the amount of variance of observed feature j explained by the m latent features.
We will consider c_j ∈{.7,.8,.9}, indicating moderate to high communality values.
We will do so by letting each observed feature be an indicator feature for a single latent feature by setting its corresponding loading to .6, indicating a moderately high loading.
The remaining m - 1 loadings for that feature are then determined as

    λ_m ≡√(c_j - .6^2/m - 1),

such that c_j = .6 + (m-1)λ_m^2.
Hence, we are explicitly not considering factor structures that are factorially pure.
For each combination of p, m, and c_j we will consider both balanced and unbalanced factor structures.
A balanced factor structure is considered to be a structure where each latent feature has (approximately) the same number of indicator features.
An unbalanced factor structure is considered to be a structure where this is not the case.
In the unbalanced setting where p = 100 we use the following numbers of indicator features:

    [40  20  15  15  10], 
        [20  10  10  10  10  10  5  5  5  5  5  5], 
        [10  10  5  5  5  5  5  5  5  5  5  5  5  5  5  3  3  3  3  3],

for the true m = 5, m = 12, and m = 20 solutions, respectively.
For the unbalanced p = 200 setting these vectors are multiplied by 2.
With these settings we can generate for each possible combination a population loadings matrix, say Λ^⋆.
We can then, according to the model, determine the population uniqueness matrix as Ψ^⋆≡𝐈_p - Λ^⋆(Λ^⋆)^T∘𝐈_p.
For each possible combination of p, m, c_j, and balancedness we then generate 100 data sets according to 𝒩_p(0, Λ^⋆(Λ^⋆)^T + Ψ^⋆) for each of the sample sizes n ∈{50, 75, 100, 150, 250}.
For each individual data set we then determine, after standardization, 𝐑(ϑ^) based on 5-fold CV (Section <ref>).
For each individual 𝐑(ϑ^) we record, for each considered approach (see Section <ref> below), the value of m deemed or selected as optimal in accordance with the decision rule.
Figure <ref> contains a schematic depiction of the simulation setup.



 §.§ Comparison

The performance of the GB will be compared to testing and model selection strategies that, in the literature, are often applied to the problem of selecting m.
In particular, it will be compared to likelihood ratio testing and various information criteria.
These are described below.



  §.§.§ Likelihood ratio test

The most formal approach to factor analytic dimensionality assessment is through likelihood ratio (LR) testing.
Let Θ̂_m = {Λ̂,Ψ̂}_m denote the ML estimate of the parameters under m factors and let Θ_m denote its population counterpart.
Testing an m-factor model against the saturated model then amounts to assessing:

    H_0:    Σ = Σ(Θ_m)
    H_a:    Σ≻ 0.

The corresponding LR criterion then converges, under the standard correlation matrix and corresponding parameter estimates under m-factors, to (n - 1) times the discrepancy function evaluated at the ML-parameters, and is approximately χ^2-distributed under certain regularity conditions <cit.>.
That is,

    𝒯_m ≡ (n-1)F[Σ(Θ̂_m);𝐑] χ^2_κ_m,

with κ_m = [(p - m)^2 - (p + m)]/2, the number of parameters in the saturated model (i.e., the unstructured sample correlation) minus the number of freely estimable parameters in the m-factor model.
Note that our interest lies with H_0.
The suggested strategy would then be to sequentially test solutions of increasing dimensionality m = 1, …, m_max≤⌊[2p + 1 - (8p + 1)^1/2]/2⌋ until H_0 is not rejected at Type-I error level α, i.e., until 𝒯_m < χ^2_κ_m(α).
This sequential testing approach is intended to avoid overfactoring (retaining too many factors) due to irregularities in the parameter space <cit.>.

Note that (<ref>) represents an asymptotic result.
In our setting, however, n is usually small relative to p.
Moreover, 𝐑 might even be unstable when p > n - 1 due to multicollinearity.
As 𝐑(ϑ) is an asymptotically unbiased estimator of the population correlation matrix 𝐑(ϑ) and 𝐑 will converge in the sample limit.
Hence, it is tempting to evaluate the quantity 𝒯_m at 𝐑(ϑ) and corresponding estimate Θ̂_m.
Although κ_m will in a sense overestimate the degrees of freedom when p > n - 1.

There might be several ways to cope with the last observation.
One might be to scale the likelihood ratio test statistic under 𝐑(ϑ) such that it is more in compliance with χ^2_κ_m.
This road was explored by <cit.>.
Their method, however, implicates duplication matrices and Kronecker products which pose serious computational bottlenecks for large p.
Another simple option would be to adapt the degrees of freedom to incorporate the rank deficiency of 𝐑.
One could do this by setting

    κ_m' = p'(p' + 1)/2 - [pm + p - m(m - 1)/2],

where p' = rank(𝐑).
In general, when p ≥ n then p' = rank(𝐑) = n - 1.
Some ready algebra on κ_m' > 0 under this situation will give that:

    m_max≤⌊-1/2[√(-4n^2 + 4n + 4p^2 + 12p + 1) - 2p - 1]⌋,

which puts serious constraints on the maximum value for m and the p to n ratio.
It is then very well possible that the true factor solution is not represented in the feasible domain for m.

Hence, we will refer 𝒯_m,ϑ≡ (n-1)F[Σ(Θ̂_m);𝐑(ϑ)] to χ^2_κ_m in a sequential testing approach as outlined above.
So, while we are aware that we cannot use 𝒯_m,ϑ∼χ^2_κ_m for inference in general, we are interested to assess its performance in the combinations indicated in Section <ref>.



  §.§.§ Information criteria

Information criteria are often used in selecting a value for m.
These are model selection criteria that try to balance model fit with model complexity.
Their general setup evaluates (minus 2 times) the maximized value of the (model-dependent) likelihood function offset with a penalty function dependent on the free parameters in the model.
We are interested in evaluating the Akaike Information Criterion <cit.> and Bayesian Information Criterion <cit.>, as these belong to the most oft used criteria <cit.>.
Now, let Σ(Θ̂_m) = Λ̂_mΛ̂_m^T + Ψ̂_m denote the 𝐑(ϑ^)-based MLE solution under m factors and let η = p(m + 1) - m(m - 1)/2 indicate the number of free parameters in the model.
Then, for our problem at hand, the general IC amounts to:

    n { p ln(2π) + ln|Σ(Θ̂_m)| + ( Σ(Θ̂_m)^-1𝐑(ϑ^) ) } + φ,

where the first term equals -2ln[L(Λ̂_m, Ψ̂_m; 𝐙)] (see equation (<ref>)) and where φ denotes the penalty.
For the AIC φ = 2η.
For the BIC φ = ln(n)η.
The strategy would then be to determine the AIC and BIC for a range of consecutive values of m.
The solution with the lowest IC score is deemed optimal.

In general ICs are found to overfactor <cit.>.
This behavior, however, is observed mostly in (numerical studies with) small scale models with relatively large n.
Here, we are interested to evaluate their performance in larger scale models with more challenging p to n ratios, such as outlined in Section <ref>.




 §.§ Results

The results can be found in Tables <ref>–<ref>.
As a reading example, consider Table <ref>, indicating results for the balanced factor structure situation where p = 100 and where the true m = 5.
The upper third gives the results when c_j = .7, the middle third represents c_j = .8 and, lastly, the lower third c_j = .9.
We then see that, at c_j = .7 and for a sample size of n = 150, the GB selects the 5-factor solution 94 times, the 6-factor solution 5 times, and, finally, the 7-factor solution 1 time.
All tables can be read in a similar manner.



  §.§.§ True m = 5

The results for the true m = 5 runs can be found in Tables <ref>–<ref>.
When p = 100 and the structure is balanced (Table <ref>) the AIC performs the best overall, in the sense that it selects the 5-factor solution most often.
The BIC and GB are more on a par in terms of selecting the true factor-dimension.
The AIC and BIC tend to underfactor (select too few factors) when not selecting the true factor-dimension.
The GB tends to (slightly) overfactor when not selecting the true factor-dimension.
The AIC, BIC and GB all give perfect performance, regardless of sample size, when the communality is high.
The LRT consistently underfactors unless both the communality is medium to high and the sample size is (relatively) large.
This picture persists when the factor-structure is unbalanced (see Table <ref>).
This picture also largely persists when p grows to 200 (Tables <ref> and <ref>).



  §.§.§ True m = 12

The results for the true m = 12 runs can be found in Tables <ref>–<ref>.
When p = 100 and the structure is balanced (Table <ref>) the GB dominates performance at all combinations of communality-value and sample size.
Again, the GB tends to slightly overfactor when not selecting the true factor-dimension.
Indeed, the AIC and BIC again tend to underfactor when not selecting the true factor-dimension, especially for the lower sample sizes when c_j ∈{.7,.8}.
This behavior is more pronounced for the BIC which generally seems to need higher sample sizes for adequate performance.
The LRT consistently underfactors.
It only selects the true factor-dimensionality adequately at the highest communality-value and sample size.
This picture persists when the factor-structure is unbalanced (see Table <ref>) although all approaches except the GB seem to diminish in performance somewhat.
When p grows to 200 (Tables <ref> and <ref>) the AIC fares somewhat better than GB when c_j ∈{.7,.8}.
The GB, however, gives perfect performance, regardless of sample size, when c_j = .9.
Again, the GB tends to overfactor under low sample size and communality-value combinations, whereas the ICs tend to underfactor in this situation.
The LRT never selects the true factor-dimensionality when p = 200 and seems most affected by unbalancedness of the factor structure.



  §.§.§ True m = 20

The results for the true m = 20 runs can be found in Tables <ref>–<ref>.
When p = 100 and the structure is balanced (Table <ref>) the GB dominates performance at all combinations of communality-value and sample size.
It only slightly overfactors and underfactors at the lower sample sizes.
The AIC, BIC, and LRT all grossly underfactor when c_j ∈{.7,.8} at all sample sizes, although the AIC fares somewhat better at the highest samples sizes.
The AIC and BIC pick up their performance when c_j = .9 and when n ≥ p.
The LRT gives its best performance when c_j = .9 at the highest sample size, but it never selects the true factor-dimensionality as optimal.
This picture persists when the factor-structure is unbalanced (see Table <ref>)
This picture also largely persists when p grows to 200 (Tables <ref> and <ref>).



  §.§.§ Discussion

The LRT gives, expectedly, the worst performance.
While the regularization approach gives the opportunity to compute the LRT statistic even when p > n, at low sample-sizes the expectation of 𝒯_m,ϑ^ tends to be below that of χ^2_κ_m; an observation also made by <cit.>.
The LRT only performs well at the combination of the highest considered sample size and the lowest considered feature size when model complexity (in terms of the generating latent dimensionality m) is low.
For example, the LRT performs well when p = 100, c_j = .9, and n ∈{150,250} when the true m = 5.
When the true m = 12 (and p and c_j as above) the LRT only performs adequate when n = 250.
And when the true m = 20 (and p and c_j as above) the LRT gives inadequate performance at all sample sizes.
We then make an observation that, to our knowledge, has not been made before: As, ceteris paribus, the complexity of the generating factor structure increases, the sample size should increase for the LRT to have adequate performance.

The ICs tend to give good performance when the complexity of the generating factor structure is low.
Their performance diminishes for more complex generating factor structures.
An effect that is more pronounced at the lower communality-values and lower sample sizes.
This effect is also more pronounced for the BIC than for the AIC.
This should be no surprise as the penalty-term for the BIC is dependent on n whereas for the AIC, it is not.
ICs have been observed to overfactor when n gets large <cit.>.
We observe, for the more challenging p to n ratios considered here, that the ICs tend to underfactor.
This tendency gets more severe when the complexity of the generating factor structure increases.
This tendency is also more pronounced for the BIC in comparison to the AIC.
Hence, the BIC and especially the AIC seem to balance model fit and complexity well at low generating-model complexities.
As this complexity increases, both the generating communality values and the sample size should be (relatively) high for adequate performance.

The GB is most directly dependent on the regularization of the correlation matrix.
It can be shown that 𝐑(ϑ^) is rotation equivariant, meaning that it operates on the eigenvalues of 𝐑 only.
When p ⪆ n the empirical eigenvalues of 𝐑 will be distorted, in the sense that the largest eigenvalues are too large whereas the smallest are too small.
The estimator 𝐑(ϑ^) shrinks the eigenvalues towards the central value of unity (1).
The decision rule regarding the GB is directly based on the eigenvalues of 𝐑(ϑ^).
We notice that, for generating factor structures of lower factor-dimension, the GB tends to overfactor when the communalities are low to medium in combination with lower sample sizes.
At these combinations some eigenvalues are (still) overestimated after shrinkage (driven by large sampling fluctuation).
When the generating factor-structure gets more complex, however, the shrunken eigenvalues seem better spread over the latent dimensions that matter and the GB performs well.
Even when the communalities are low and the sample size is small relative to p.
When it doesn't select the true generating m as optimal it has the (general) tendency to slightly overfactor.

In our performance assessment the AIC seems the only contender for the GB.
Our preference lies with the GB for the following reasons:

  *  Its considerable better performance when the generating factor-structure is more complex;
  
  *  Its tendency to (slightly) overfactor in cases in which it does not select the true generating value for m, as opposed to the other methods' tendency for underfactoring. Overfactoring leads to far less error in factor loading estimates and loading-structure representation (so central to the computation of the desired factor scores) than underfactoring <cit.>;
  
  *  Its considerable lower computational complexity. Its computation requires little more than a single spectral decomposition, which is of worst case computational complexity 𝒪(p^3).

In conclusion, we will use the GB as an upper-bound to the problem of selecting a value for m.
The methods for additional decision support mentioned in Section <ref> may then be combined with the GB to inform a final decision.








§ ADDITIONAL INFORMATION ORAL SQUAMOUS CELL CARCINOMA DATA

In this section we give additional information on the oral squamous cell carcinoma data as used in the external validation setting.
Section <ref> contains an overview of data acquisition.
Section <ref> contains information on the factor solution obtained in the training projection.
Lastly, Section <ref> contains visualizations of the analysis results.




 §.§ Image processing and radiomic feature extraction

Two independent cohorts of human papillomavirus-negative oropharyngeal squamous cell carcinoma patients were available.
Data on the Amsterdam University medical centers (location VUmc) training cohort was gathered from 2008 to 2012.
Data on the University Medical Center Utrecht validation cohort was gathered from 2010 to 2013.
The oropharyngeal tumors were imaged with axial T1-weighted MRI scans without gadolinium enhancement.
Different MRI vendors and protocols were used in image acquisition.
The images were segmented using VelocityAI 3.1 (Velocity Medical Solutions, Atlanta GA, USA).
Identification of the volumes of interest was performed manually by a senior radiologist.
The delineated regions were exported as Digital Imaging and Communications in Medicine Radiotherapy Structure Sets <cit.>.

A total of 545 raw radiomic features were extracted from the MRI-images using in-house build software <cit.> relating to, e.g., morphology, texture, and intensity.




Texture features were derived from grey-level co-occurrence matrices or from grey-level run-length matrices, calculated from either the x-y plane (2D) or volumetrically (3D), either per direction (and then averaged) or after merging over all possible directions (combined).
In addition, images were discretized to bin sizes of 32, 64 and 128.
The texture features were subsequently averaged over the 12 representations (each combination of bin size, dimensionality, and direction).
This left a total of p = 89 radiomic core features.
Table <ref> gives a listing of all final features.
Additional information on the data acquisition and radiomic features can be found in <cit.>.
The processed data object is available from the Authors upon reasonable request.



 §.§ Factor solution training data

The first step in the factor modeling pipeline is to get a regularized estimate of the correlation matrix on the nonredundant radiomic features.
Figure <ref> visualizes the redundancy pattern among the full set of 89 radiomic features.
The redundancy filtering algorithm (Algorithm <ref>) was run with τ set to .95 and it retained p^* = 51 out of the original p = 89 features.
The remaining correlation matrix was then subjected to penalized ML estimation.
The optimal value for the penalty parameter was determined by 5-fold CV of the log-likelihood function and was found to be .04950291.
The resulting regularized correlation matrix was well-conditioned, as indicated by a low condition number (292.97).



The second step consists of factor analytic data compression.
The GB indicates 10 factors.
In accordance with its purported usage (see Section <ref>) we treat 10 as an upper-bound to the optimal dimension of the latent vector.
As the proportion of additional variance explained flattens after 7 common factors and as solutions in excess of 7 factors produce weak factors, we set m̃ = 7.
The Varimax-rotated maximum-likelihood estimate of the factor loadings matrix is given in Table <ref>.
The 7 common factors explain approximately 76% of the variance in the observables.



In the third step, factor scores are obtained through the approach outlined in Section <ref> as well as Section 2.2.3. of the Main Text.
Using the estimates of the factor solution in conjunction with the training data will give the factor scores for the training data.
Using the estimates of the factor solution in conjunction with the validation data will give the factor scores for the validation data.
The factor scores are highly determinate with SMC_k^ξ > .9 for all factors.




 §.§ Prediction results

Figure <ref> gives a visualization of the prediction results provided in Section 3.4. of the Main Text.






§ ADDITIONAL INFORMATION HEAD AND NECK SQUAMOUS CELL CARCINOMA DATA

In this section we give additional information on the head and neck squamous cell carcinoma data as used in the internal validation setting.
Section <ref> contains an overview of data acquisition.
Section <ref> contains information on the factor solution obtained in the training projection.
Lastly, Section <ref> gives the analysis results when redundancy-filtering is applied for all methods.




 §.§ Image processing and radiomic feature extraction

A total of 175 head and neck squamous cell carcinoma patients were available from the Amsterdam University medical centers.
On these patients low-dose ^18F-FDG-PET/CT was performed according to the European Association of Nuclear Medicine (EANM) guidelines 2.0 on a Gemini-TF or Ingenuity TF-64 PET/CT (Philips Heatlthcare, Best, The Netherlands) with EARL (EANM Research Ltd.) accreditation.
The PET images were reconstructed using iterative ordered subsets expectation maximization with the parameters optimized for the head and neck area (i.e., 4 iterations, 16 subsets, 53-dimensional Gaussian filter) with photon attenuation correction.
Reconstructed images had an image matrix size of 144 × 144, voxel size of 4 × 4 × 4.
Low-dose-CT was collected using a beam current of 50 at 120 for anatomical correlation of ^18F-FDG uptake and attenuation correction.
CT-scans were reconstructed using an image matrix size of 512 × 512 resulting in pixel sizes of 1.17 × 1.17 and a slice thickness of 5 .

Radiomics features were extracted from an in-house build Accurate tool <cit.>.
It provides a 3D implementation of feature extraction methods for shape, intensity, texture and wavelet-type features.
Wavelet-type features were obtained through a wavelet transform using a Coiflet scaling function.
For the texture and wavelet analysis, images were discretized to a fixed bin size of 64 bins and 0.25 standardized uptake value (SUV) in the CT and PET, respectively <cit.>.
This choice of bin sizes resulted in similar numbers of bins in both modalities.
On average 64 bins were analyzed.
For each patient a total of p = 436 radiomic features were extracted on the primary tumors.
Table <ref> gives a listing of all extracted features.

The ngldmFeatures3D_Dependence count percentage, ngldmFeatures2Dmrg_Dependence count percentage, gldzmFeatures2Davg_Grey level variance GLDZM, and Intensity histogram_minimum features were removed as they were behaving as de facto constants (no variation over the respective patients).
In addition, 1 person with a follow-up time of 0 was removed from the data.
Hence, the final data object concerned p = 432 radiomic features extracted from the primary tumors of n = 174 patients.

Additional information on the data acquisition and radiomic features can be found in <cit.>.
The processed data object is available from the Authors upon reasonable request.




 §.§ Factor pattern training data

Again, the first step in the factor modeling pipeline is to get a regularized estimate of the correlation matrix on the nonredundant radiomic features.
Figure <ref> visualizes the redundancy pattern among the full set of 432 radiomic features.
The redundancy filtering algorithm (Algorithm <ref>) was run with τ set to .95 and it retained p^* = 124 out of the original p = 432 features.
The remaining correlation matrix was then subjected to penalized ML estimation.
The optimal value for the penalty parameter was determined by 5-fold CV of the log-likelihood function and was found to be .0591859.
The resulting regularized correlation matrix was well-conditioned, as indicated by a relatively low condition number (805.10).



The second step consists of factor analytic data compression.
The GB indicates 13 factors as an upper-bound to the optimal dimension of the latent vector.
As the proportion of additional variance explained flattens after 8 common factors and as solutions in excess of 8 factors produce weak factors, we set m̃ = 8.
The Varimax-rotated maximum-likelihood estimate of the factor loadings matrix is given in Table <ref>.
The 8 common factors explain approximately 79% of the variance in the observables.
The factor scores obtained in the third step are highly determinate with SMC_k^ξ > .94 for all factors.



 §.§ Results under redundancy filtering for all methods

Figure <ref> and Table <ref> contain the results when the RSF, CSF, and Cox boosting approaches are based on the same redundancy-filtered feature-set as the  approach.
Table <ref> contains the apparent and averaged cross-validated prediction errors (integrated Brier scores) as well as the overall time-dependent explained variation based on integrated Brier scores with the KM estimator as the reference model.
We see that the three alternative methods all fare somewhat better when redundancy filtering is applied.
Nonetheless, the overall results remain qualitatively the same as discussed in Section 3.5. of the Main Text.
The most stable approach resulting in the lowest prediction error is the proposed  approach.







[Radiomics_SM]