Spatial Process Decomposition for Quantitative Imaging Biomarkers Using Multiple Images of Varying Shapes
ShengLi Tzeng Jun Zhu Amy Weisman Tyler Bradshaw Robert Jeraj
stat.ME








Spatial process decomposition for quantitative imaging biomarkers
using multiple images of varying shapes
    ShengLi Tzeng^1, Jun Zhu^2,3[Corresponding author. Email address: jzhu@stat.wisc.edu],  
Amy Weisman^4, 
Tyler Bradshaw^4 and Robert Jeraj^4,5  
 
^1Department of Applied Mathematics, National Sun Yat-sen University, Taiwan
 
^2Department of Statistics,  University of Wisconsin-Madison, USA 
 
^3Department of Entomology,  University of Wisconsin-Madison, USA 
 
^4Department of Medical Physics, University of Wisconsin Madison, USA 

^5Department of Human Oncology, University of Wisconsin Madison, USA  
    Received: date / Accepted: date
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================







Quantitative imaging biomarkers (QIB) are extracted from medical images in radiomics for a variety of purposes including noninvasive disease detection, cancer monitoring, and precision medicine. The existing methods for QIB extraction tend to be ad-hoc and not reproducible. In this paper, a general and flexible statistical approach is proposed for handling up to three-dimensional medical images in an objective and principled way. In particular, a model-based spatial process decomposition is developed where the random weights are unique to individual patients for component functions common across patients. Model fitting and selection are based on maximum likelihood, while feature extractions are via optimal prediction of the underlying true image. A simulation study evaluates the properties of the proposed methodology and for illustration, a cancer image data set is analyzed and QIBs are extracted in association with a clinical endpoint. 


Keywords: 
biomarker, medical image, multi-resolution model, radiomics, spatial statistics, thin-plate splines. 








§ INTRODUCTION


With the advance of imaging techniques, radiomics has emerged as a field in medicine that concerns extraction of quantitative imaging biomarkers (QIBs) from medical images, finding a variety of applications such as noninvasive disease diagnosis, cancer monitoring, treatment response assessment, and precision medicine <cit.>.

While QIBs overcome the limitations of traditional tools including subjective visual interpretation with inter-observer variability that often occur in medical image assessment, the existing methods in radiomics could benefit more from statistical approaches in connection to the subsequent data analysis involving QIBs. 
In this paper, we propose a novel statistical method for extracting potential QIBs, providing a general framework for dealing with up to three-dimensional (3D) medical images in an objective and statistically more rigorous way.

The current procedures for developing QIBs are rather complicated, involving image acquisition, region of interest (ROI) identification and segmentation, potential biomarker quantification and extraction, and association with clinical endpoints and validation.
For example, <cit.> followed a QIB procedure that started with acquisition by digitizing dental X-ray panoramic images into 2621×1244 pixels, and then used contrast enhancement and morphological operations for segmenting jaws and teeth.
The contours of ROIs were delineated based on the segmentations and three different features were extracted (i.e., dental work performed on the individual, number of dental screws, and detection of missing teeth), followed by testing the performance of the extracted features. 

Texture analysis (TA) is a commonly used technique for potential
biomarker quantification and extraction, which systematically constructs high-throughput information from ROIs <cit.>. 
However, this technique has multiple challenges.
First, the number of extracted TA features is often much larger than the number of patients in a study, in which case popular classical statistical tools (e.g., linear regression, Cox proportional hazard regression) can not be directly applied. 
Appropriate dimension reduction methods are therefore important for not only mitigating over-fitting but also alleviating high correlations among a large number of high-throughput features. 
In practice, however, dimension reduction is not always performed. For example, <cit.>
reviewed several univariate filtering methods which associate a clinical endpoint of interest with one feature at a time, not considering the dependency between features. 
If several features have high correlations with one important feature, these filtering methods are likely to
select most of these features and the resulting strong collinearity would make many statistical tools challenging to use. 
Thus, despite their popularity, filtering methods are not satisfactory in such context. 
A simple alternative to filtering methods is to consider commonly-used features in the literature, leading
to more interpretable analysis. 
Unfortunately, commonly-used features are not necessarily the more useful ones and worse yet, focusing
on well-known features impedes the discovery of other viable biomarkers. 

Second, the existing TA techniques were developed largely for rectangular images but have been widely adopted in cancer studies where the tumor ROIs are non-rectangular. 
<cit.> identified such a shape issue arising from non-rectangular two-dimensional (2D) images and reviewed ad hoc remedies such as embedding a non-rectangular ROI in a rectangular image with blank pixels surrounding the ROI; however, this approach becomes even more problematic for 3D cases with a much larger portion of blank voxels in the cube where the ROI is embedded. 
In addition, the images to be analyzed are not conformed to the same shape or size across patients, as individual patients' ROIs are unique. 
A naive way would be to use the pixel-wise (or voxel-wise) average of TA features in a local rectangle (or cube), but would potentially lose much information in the spatial structures. 
How to extract important spatial characteristics from ROIs with different shapes and sizes is a formidable problem that seems largely unaddressed. 

A third challenge is that many TA techniques require binning the intensity values observed on a medical image into distinct categories, known as quantization, for constructing local-pattern matrices (e.g., neighboring gray-level dependence and co-occurrence) in the extraction of QIBs. 
<cit.> demonstrated by examples that quantization has a large impact on many TA features. 
Yet there is little consensus on the number of bins and the binning method to use, resulting in arbitrary quantization and hindering the reproducibility of studies. 

In light of the aforementioned challenges with TA, we propose a new statistical method applicable for extracting nearly uncorrelated features for the intensities in ROIs. 
A key strategy is model-based spatial process decomposition, which formulates the image data as linear combinations of common spatial processes, building upon previous work for similar objectives. 
For example, the wavelets and principal component analysis <cit.> developed for hyperspectral images produced multi-resolution features capturing important spatial patterns, but unfortunately was not easily applicable to non-rectangular ROIs. 
<cit.> proposed a regularized principal component analysis method for irregularly shaped 2D data, although each patient was required to share the same sampling locations. 
<cit.> used functional mixed-effects models for 2D functions by employing bivariate splines on triangulations, which handles varying shapes. 
However, the triangulation requires that each triangulated subregion contains at least a few data locations, which makes its extension to 3D images a challenge.

The novelty of our proposed methodology compared with the previous work is its ability to deal with 1D time point, 2D pixel, and 3D voxel data within a same general framework of spatial random-effects models.
Although the spatial random-effects model is a popular geostatistical tool, the purpose in geostatistics is largely spatial regression or kriging at unsampled locations <cit.>, whereas our main goal is to find the important building blocks for the observed image intensities. 
Further, for the most part, geostatistical data are 1D or 2D and there are no independent replicates, while in radiomics, the images are 2D or 3D from multiple patients.
Our innovation here is to let the random effects in our model carry weights, unique to individual patients, for the common building blocks, and an orthogonalization of these random effects produces nearly uncorrelated biomarkers. 
Finally, our proposed methodology has the advantage that the medical images require no quantization and the ROIs could have different shapes.

The remainder of this paper proceeds as follows. In Section <ref>, we propose a data model and a latent model. The methodology for parameter estimation, model selection, and feature extraction is given in Section <ref>. Section <ref> presents a simulation study evaluating the properties of
our proposed method. For illustration, a data example is presented in Section <ref>, followed by conclusions and a discussion in Section <ref>.






§ MODEL




 §.§ Data Model


Let N≥1 denote the number of patients and each patient has a 2D (or 3D) image of ROI with observed intensities across pixels (or voxels) for an index under study (e.g., fluorodeoxyglucose level in Figure <ref>). 
Let n_j denote the number of observation locations, s_j,1,…,s_j,n_j, generally the centroids of pixels or voxels, for the jth patient with j=1,…,N.
Let the union of the observation locations among all the patients be denoted as D={s_1,…,s_n}=∪_j=1^N{s_j,1,…,s_j,n_j}, where s_1,…,s_n are n distinct observation locations and D⊂ℝ^d is the collection of locations in a d-dimensional common spatial domain shared across the N patients for d=2 or 3. 
Thus, each patient has n potential observation locations, of which n_j≤ n locations have observations and the other n-n_j locations have no observations for the jth patient. 
When N=1, there is just one patient and n_j=n; that is, we only have a single realization, which is common in geostatistics but not in radiomics. 

Next, for the jth patient where j=1,…,N, we define an n_j× n incidence matrix 𝐎_j, whose (k,i)th element is 1 if s_j,k=s_i and is 0 otherwise. 
Let z_j(s_i) denote the potentially observable image intensity at location s_i for i=1,…,n and let z_j=(z_j(s_1),…,z_j(s_n))' denote the vector of image intensities at all the locations with and without observations. 
Then, 𝐎_jz_j is the vector of observed image intensity data from the jth patient. 
In the context of QIB, we define a latent variable y_j(s) to represent the true image intensity at location s∈ D and let y_j=(y_j(s_1),…,y_j(s_n))' denote the vector of true image intensities at all the locations for the jth patient.

In what follows, we assume that the image intensity data for the jth patient are observed according to z_j(s_i)=y_j(s_i)+ε_j(s_i), where ε_j(s_i) is an error corresponding to white noise, for i=1,…,n and j=1,…,N. 
Let ε_j=(ε_j(s_1),…,ε_j(s_n))^' denote a vector of errors. 
For j=1,…,N, we have 

    𝐎_jz_j=𝐎_j(y_j+ε_j),

where ε_j∼𝒩(0,σ^2I_n) is assumed to be uncorrelated with y_j, and ε_j is assumed to be uncorrelated with ε_j^* for j≠ j^*. 



 §.§ Latent Model


For j=1,…,N, we model the jth latent variable y_j(s) by

    y_j(s)=μ_j(s)+η_j(s),

where  μ_j(s) is a deterministic mean function and η_j(·) is a mean-zero spatial process with a spatial covariance function C(s,s^*)=cov(y_j(s),  y_j(s^*)) for s,s^*∈ D.  
By a truncated version of the Karhunen-Loéve expansion <cit.>, 
we model η_j(s) as a linear combination of H component functions g_h(s) such that

    η_j(s)=∑_h=1^Hθ_jhg_h(s),

where orthogonality holds for g_h(·) such that ∑_i=1^n g_h(s_i)g_h^*(s_i)=0 for h≠ h^*.
Further, let θ_j=(θ_j1,…,θ_jH)' denote a vector of random weights such that
θ_j∼𝒩(0,Λ), where Λ=diag{λ_1,…,λ_H} is an H× H diagonal matrix with variance components λ_h satisfying λ_1≥⋯≥λ_H>0, and θ_j and θ_j^* are uncorrelated for j≠ j^*. 
By (<ref>)–(<ref>), we have, for s∈ D,

    y_j(s)=μ_j(s)+∑_h=1^Hθ_jhg_h(s).

The covariance function, for s,s^*∈ D, is

    C(s,s^*)=cov(y_j(s),y_j(s^*))=∑_h=1^Hλ_hg_h(s)g_h(s^*).


For flexible forms of the component functions, we represent g_h(s) by a linear combination of K pre-specified basis functions, denoted as f_1(s),…,f_K(s), for s∈ D. That is,

    g_h(s)     =∑_k=1^Ku_hkf_k(s)=u_h^'f_K(s),

where f_K(s)=(f_1(s),…,f_K(s))' denotes a vector of K basis functions at location s and u_h=(u_h1,…,u_hK)^' denotes a vector of coefficients for the hth component function. 
With H≤ K≤ n, the representation of the component function g_h(s) with basis functions bears similarity to regression splines that approximate complex covariance functions by a set of piecewise polynomial bases <cit.>.
In particular, let  U=[u_1,…,u_H] denote a K× H matrix of coefficients and let w_j=(w_j1,…,w_jK)'=Uθ_j denote a vector of random effects.
Since θ_j∼𝒩(0,Λ) in (<ref>), we have w_j∼𝒩(0,M) where M=UΛU^'.
It is straightforward to show that the variance matrix M=UΛU^' is nonnegative-definite and w_j's are mutually uncorrelated and also uncorrelated with ε_j's
in (<ref>). 

By (<ref>) and (<ref>), we can rewrite y_j(·), for j=1,…,N, as

    y_j(s)=     μ_j(s)+∑_k=1^Kw_jkf_k(s)= μ_j(s)+w_j'f_K(s).

By (<ref>), the covariance function (<ref>), for s,s^*∈ D, now becomes 

    C(s,s^*)=cov(y_j(s),y_j(s^*))=f_K(s)'Mf_K(s^*).

The covariance function given in (<ref>) is quite flexible, as there is no need to assume stationarity and no specific parametric structure is assumed of M. 

Thus, the data model (<ref>) can be rewritten as

    z_j(s_i)= μ_j(s_i)+∑_k=1^Kw_jkf_k(s_i)+ε_j(s_i)= μ_j(s_i)+w_j'f_K(s_i)+ε_j(s_i),

such that the image intensity for a given pixel or voxel s_i∈ D within a patient's ROI has three additive elements: (1) a patient-specific mean function μ_j(s_i); (2) a linear combination of basis functions w_j'f_K(s_i) that captures spatial correlations between pixels or voxels with random effects and a class of pre-specified basis functions such that w_j∼𝒩(0,M) where M=UΛU^', and (3) a white noise ε_j(s_i) that captures measurement error and local heterogeneity such that ε_j∼𝒩(0,σ^2I_n).
Since g_h(s) and Λ are common across patients, U, f_K(s), and M are also the same for all patients. 
In other words, the image data across patients are modeled as independent realizations from (<ref>) including those components shared among patients.






§ METHODOLOGY



 §.§ Multi-resolution Basis Functions


With the observation locations {s_1,…,s_n} as the control points, we construct basis functions based on thin-plate splines (TPS) kernels denoted as, 
ψ_i(s)=(1/12)s-s_i^3 for d=1, 
ψ_i(s)=1/(8π)s-s_i^2log(s-s_i) for d=2, and
ψ_i(s)=-(1/8)s-s_i for d=3, where i=1,…,n and · denotes the L^2 norm <cit.>.
Let Ψ=(ψ_i(s_i^*))_n× n denote an n× n matrix of TPS kernels comprising ψ_i(s_i^*) for i,i^*=1,…,n.
We also define an n×(d+1) design matrix to be X=((1,s_1^')^',…,(1,s_n^')^')^' where the ith row is  (1,s_i^')^' for i=1,…,n and we let Ω=I_n-X(X'X)^-1X' be an n× n projection matrix.

With the notation above, we define the basis functions f_k(s) as follows.
Let f_1(s)=1 and f_k(s)=s_k-1 for k=2,…,d+1, where s=(s_1,…,s_d)^'.
For k=d+2,…,n, we form f_k(s) by linear combinations of the TPS kernels where the coefficients are obtained from an eigen-decomposition, which essentially extracts the orthogonal subspaces in TPS but not in the space of linear functions. 
Specifically, we pre- and post-multiply Ψ by the projection matrix Ω and then solve for the eigen-vectors V and scalars α_k that satisfy, Vdiag{α_1,…,α_n}V'=ΩΨΩ, with α_1≥⋯≥α_n.
Now, let f_k(s)=α_k-d-1^-1{ψ(s)-ΨX(X'X)^-1x}'v_k-d-1 for k=d+2,…,n,
where x=(1,s^')^'=(1,s_1,…,s_d)^' , ψ(s)=(ψ_1(s),…,ψ_n(s))', and v_k is the kth column of the eigen-vector V serving as the coefficients for the linear combination of the n TPS kernels in ψ(s). 
The rationale of the basis construction  procedure is about the smoothness for a function. Let 

    J(f)=∫_R^d∑_γ_1+⋯+γ_d=22!/γ_1!⋯γ_d!(∂^2f(s)/∂ x_1^γ_1⋯∂ x_d^γ_d)^2ds
 
measure the roughness of f. Consider two function f̃ and  f̌ satisfying ∑_i=1^nf̃^2(s_i)=∑_i=1^nf̌^2(s_i)=1. 
When f̃ represents a broader scale spatial pattern and f̌  represents a more local scale spatial pattern, it is expected that J(f̃)< J(f̌). It has been shown that 

    f_k=min_f{ J(f):∑_i=1^nf(s_i)f_1(s_i)=⋯=∑_i=1^nf(s_i)f_k-1(s_i)=0,∑_i=1^nf^2(s_i)=1};

that is, f_k(·) is the least fluctuating function under the orthogonality constraints and the unit-length identifiability condition. Hence, spatial patterns are extracted sequentially from the smoothest to the roughest in the orthogonal subspaces.  Broader-scale spatial patterns can also be easily captured by f_k, in contrast to most TAs that rely only on a few neighbors via neighboring gray-level dependence or gray-level co-occurrence.        

An example of the basis functions on an 8× 8 image of 64 pixels is given in Figure <ref>.
The multi-resolution basis functions f_k(·) are ordered from the smoothest to the roughest, corresponding to increasing k. 
Unlike most multi-resolution basis functions based on wavelets, f_k(·) can be directly applied to irregular grids of observation locations and non-rectangular spatial domains.
It has been shown that using the first K≤ n basis functions, {f_k(s);  k=1,…,K}, can approximate quite well complex, stationary or nonstationary, covariance functions (). 



 §.§ Model Estimation and Selection


We estimate the patient-specific mean function by a grand mean 

    μ̂_j(s)≡μ̂_j=n_j^-11_n_j^'𝐎_jz_j,

where 1_n_j is a vector of n_j 1's.
We let z̃_j=𝐎_jz_j-μ̂_j denote a vector of detrended observations for the jth patient; j=1,…,N. 

Given the set of basis functions {f_k(s);  k=1,…,K} in Section <ref>, we develop an expectation-maximization (EM) algorithm for obtaining the maximized likelihood (ML) estimator of σ^2 and M and denote them by σ̂^2 and M̂, respectively <cit.>.
The random vector  w_j is treated as unobserved data, and the distribution of w_j conditional on z̃_j is derived in the E-steps, while the estimates of σ^2 and M are updated in the M-steps. 

Specifically, for iteration t=0,1,2,…, given the estimated parameters σ^2(t) and M^(t) in the tth iteration, the E-step computes the conditional mean and conditional variance of w_j as

    ŵ_j^(t) =     E(w_j|z̃_1,…,z̃_n)=M^(t)F_j^'{F_jM^(t)F_j^'+σ^2(t)I_n_j}^-1z̃_j,
    Q_j^(t) =     Var(w_j|z̃_1,…,z̃_n)=M^(t)-M^(t)F_j^'{F_jM^(t)F_j^'+σ^2(t)I_n_j}^-1F_jM^(t),

where F_j=O_jF and the (i,k)th element of the n× K matrix of basis functions, F, is f_k(s_i) for i=1,…,n and k=1,…,K.
Then the M-step updates the parameters as:

    M^(t+1)    =N^-1∑_j=1^N{ŵ_j^(t)(ŵ_j^(t))^'+Q_j^(t)},
    σ^2(t+1)    =(∑_j=1^Nn_j)^-1∑_j=1^N{z̃_j^'z̃_j-2z̃_j^'F_jŵ_j^(t)+tr(F_jM^(t+1)F_j^')} .

For t=0, the estimates are initialized to be ŵ_j^(0)=0 and Q_j^(0)=I_K.
We alternate between the E- and M-steps iteratively until convergence. 

Let σ̂_K^2 and M̂_K denote the ML estimators of σ^2 and M, respectively, with K basis functions.
Also let ℓ(σ̂_K^2,M̂_K) denote the log-likelihood function for the data model (<ref>) evaluated at σ̂_K^2 and M̂_K.
We select K∈{d+1,…,K̃} for a sufficiently large K̃ based on the Akaike's information
criterion (AIC), 

    AIC(K)=-2ℓ(σ̂_K^2,M̂_K)+2df(K).

The degrees of freedom (df) in AIC(K) is df(K)=K(K+1)/2+1 if K≤ N and df(K)=KN+1-N(N-1)/2 if K≥ N; that is, the number of free parameters in M depends on the number of basis functions K and the number of patients N. 
The number of basis functions selected by AIC is denoted as K^*=min_d+1≤ K≤K̃AIC(K).
An advantage of the selection via AIC is that the computation is more efficient, in contrast to cross validation based methods that are computationally intensive and involve resampling variability <cit.>. 



 §.§ Feature Extraction


Although the spatial process η_j(·) is composed of component functions as in (<ref>), there is no explicit sample covariance matrix to decompose due to the unequal number of observation locations across patients and the large number of missing values within patients. 
Here, we let the ML estimators of M and σ^2 for the final model (<ref>) with K=K^* be denoted as M̂_K^* and σ̂_K^*^2, respectively. 
By w_j=[u_1,…,u_H]θ_j, we obtain the estimated component functions ĝ_h(s), the estimated number of component functions H^*, and the estimated eigen-values λ̂_h, for h=1,…,H^*, as follows.  

Let Ûdiag{λ̂_1,…,λ̂_K^*}Û'
be the eigen-decomposition of M̂_K^*, where Û=[û_1,…,û_K^*]
and λ̂_1≥⋯≥λ̂_H^*>0=λ̂_H^*+1=⋯=λ̂_K^*.
We only need the estimates of the first H^* components corresponding to the positive eigen-values.
Then by (<ref>), we estimate g_h(s), for h=1,…,H, by 

    ĝ_h(s)=û_h^'f_K^*(s).

In addition, an estimator of the covariance between two locations for the jth patient is 

    cov(y_j(s),y_j(s^*))=∑_h=1^H^*λ̂_hĝ_h(s)ĝ_h(s^*),
 
where ĝ_h(s) is the hth estimated component function (<ref>).

Further, let Λ̂=diag{λ̂_1,…,λ̂_H^*} denote the estimated variance matrix for θ_j.
We define an n_j× H^* matrix Ĝ_̂ĵ for the H^* component functions taken values at observation locations for the jth patient, with the (i,h)th element being ĝ_h(s_j,i) for i=1,…,n_j, h=1,…,H^*.
Thus, an empirical predictor with the smallest mean squared prediction error (MSPE) of y_j(s) is given by 

    ŷ_j(s)     =μ̂_j+cov(y_j(s),𝐎_jz_j)cov^-1(𝐎_jz_j,𝐎_jz_j)(𝐎_jz_j-μ̂_j)
        =μ̂_j+ĝ(s)^'Λ̂Ĝ_̂ĵ^'(Ĝ_jΛ̂Ĝ_̂ĵ^'+σ̂_K^*^2I_n_j)^-1z̃_j,

where  ĝ(s)=(ĝ_1(s),…,ĝ_H^*(s))^'.
It follows that the conditional expectation of θ_j given data can be obtained by equating ŷ_j(s)=μ̂_j+∑_h=1^Hθ̂_jhĝ_h(s)=μ̂_j+ĝ(s)^'θ̂_j with (<ref>) and thus,

    θ̂_j=Λ̂Ĝ_̂ĵ^'(Ĝ_jΛ̂Ĝ_̂ĵ^'+σ̂_K^*^2I_n_j)^-1z̃_j.

We then consider the vector (μ̂_j,θ̂_j^')^' to be the data-driven features for the jth patient. 

For a new patient with image intensity data denoted as z_0 observed at n_0 locations r_1,…,r_n_0, the data-driven features can be computed as follows. 
First, the patient-specific mean function is estimated to be μ̂_0=n_0^-11_n_0^'z_0 and is subtracted from z_0 to obtain z̃_0. 
We can then construct Ĝ_0, whose (i,h)th element is ĝ_h(r_i) in (<ref>) for i=1,…,n_0, h=1,…,H^*. 
Together with the estimated Λ̂ and σ̂_K^*^2, the estimated random weights θ̂_0 are obtained from (<ref>) with j=0. 






§ SIMULATION STUDY


We conducted a simulation study to investigate the properties of the proposed methodology in Section <ref>, focusing on the determination of component functions.
Four sample sizes (i.e., number of patients) N=30, 40, 50, and 60 were considered and the image for any given patient was in the unit square U=[0,1]×[0,1].
The number of pixels (or, image resolution) in the unit square was set to n=2^L× 2^L, where L=3, 4, …, 7, corresponding to the total number of pixels n=2^3× 2^3, 2^4× 2^4, …, 2^7× 2^7, respectively. 
Given a sample size N and an image resolution, different shapes of ROIs were generated within the square for N different patients, where the proportion of pixels outside the ROI (i.e., the missing rate) was selected to be between 0.25 and 0.6 mimicking the data example in Section <ref>. 
Figure <ref> illustrates the simulated ROIs for N=20 patients and n=64 × 64 pixels per image. 

The true image intensity, or signal, for the jth patient was generated according to a linear combination of cosine functions:

    y_j(s)=∑_k=1^17δ_kw_jkcos(‖s-c_k‖ ),

where the cosine functions cos(‖s-c_k‖ ) for s∈ U relied on the Euclidean distance ‖s-c_k‖ between s and 17 pre-selected center points defined to be c_k∈{(i/5,  i^*/5)^';  i,i^*=1,2,3,4}∪{ (0.5,0.5)^'}.
The random weights for the cosine functions were independently distributed as w_jk∼𝒩(0,(τ+ς_k)), where τ controlled the signal-to-noise ratio (SNR) and  ς_k introduced heterogeneity for the variances of different components.
We selected τ to be 1, 2, 3, or 4 and ς_k was drawn independently and uniformly from the interval [-0.5,0.5].
Further, δ_k=1 if the kth cosine function was included in the model and δ_k=0 otherwise. 
That is, only a portion of the cosine functions were used in generating y_j(s), determined by an image complexity level p=∑_k=1^17δ_k.
We considered seven image complexity levels, p=0,1,…,6. 
An observation of intensity z_j(s) in the pixel centered at s was then generated according to z_j(s)=y_j(s)+ε_j(s), where the errors followed the standard normal distribution ε_j(s)iid∼𝒩(0,1). 

In total, there were 560 combinations of the number of patients N, the image resolution n, the SNR factor τ, and the image complexity level p. 
For each combination, we ran 100 simulations and computed the median of the estimated numbers of extracted component functions, H^*, given in Section <ref>. 
Since the cosine functions are not orthogonal to each other, the number of component functions was not necessarily the same as the image complexity level p. 
Moreover, by (<ref>), it could happen that an underlying cosine function be approximated by more than one estimated component functions. As shown in Figure <ref>, the patterns of H^* for different sample sizes N and SNR factors τ are quite similar. 
The determining factors seem to be n and p such that as the image resolution increases and/or the image complexity increases, a larger number of component functions would be extracted.
These results are desirable, as in practice, the true eigen-functions are unknown and we expect to be able to extract more information when the underlying image structure is more complex or when the image resolution is higher with more and finer details.  






§ DATA EXAMPLE




 §.§ Cancer Image Data


In a cancer research study, <cit.> collected 3D image data from canine patients with spontaneous sinonasal tumors, investigating the association between clinical end points of interest and potential imaging biomarkers.

<cit.> adopted positron emission tomography (PET) as the imaging technique, which injects a radioactive molecule into an individual patient and records where specific molecular functions are occurring. 
PET images were also paired with computed tomography (CT) images, providing both functional and anatomic information about the patient. 
Although pre-, mid-, and post-treatment data were available, here we focus on the pre-treatment data to extract QIBs for prognostic purposes.

Three PET radiotracers were considered for each patient in this data example, namely, fluorodeoxyglucose (FDG) which measures glucose metabolism, flurothymidine (FLT) which measures proliferating cells, and Cu(II)-diacetyl-bis(N4-methylthiosemicarbazone) (Cu-ATSM) which measures hypoxia. 
All measurements were normalized by the dose injected into the patient and the body weight of the patient to achieve standardized uptake values (SUVs) for subsequent analysis. For each patient, the ROI is the contoured gross tumor volume.
The SUVs outside the ROI were regarded as missing, as it was reasonable to assume that the uptakes inside and outside the ROI were very different and only those uptakes inside the ROI were the target of interest. 
Figure <ref>(a) illustrates an example of the cancer ROI for one of the canine patients.



 §.§ Feature Extraction and Identification


We applied the methodology developed in Sections <ref>-<ref> to extract QIBs based on the 3D cancer image data described in Section <ref>.

In particular, we had three images for each of the N=22 canine patients for the three PET radiotracers, FDG, FLT, and Cu-ATSM. 
The jth patient had an ROI with n_j≤ n voxels across the three radiotracers and the images had missing values outside the ROIs. 
For each of the N=22 patients, the centroid of each ROI was transformed linearly to be the origin (0,0,0)^' in ℝ^3. 
The union of the transformed voxel centroids for all the ROIs formed the set {s_1,…,s_n} where the total number of voxels in D was n=19697. 

Model fitting was performed for each radiotracer separately based on the model given in Section <ref>, where the image data were the SUVs at transformed voxel locations within each patient's ROI. 
The model parameters were estimated by the EM algorithm and the number of basis functions K was selected by the AIC as described in Section <ref>. 
Then, the number of component functions H and the component functions g_h(·) were estimated following Section <ref>. 
We had K^*_ FDG=855, H^*_ FDG=18 for FDG, K^*_ FLT=540, H^*_ FLT=13 for FLT, and K^*_ Cu-ATSM=595, H^*_ Cu-ATSM=16 for Cu-ATSM. 
Figure <ref>(b) illustrates one of the estimated component functions ĝ_h(·).

The pairwise (Pearson) correlations between the estimated component functions ĝ_h(s) are shown in Figure <ref>. 
Within each radiotracer, the correlation between the estimated component functions was generally small, as was expected due to orthogonality.
Between radiotracers, however, there was substantial correlation with an average of 72.3% of the component functions for a given radiotracer to have more than 0.5 absolute correlation with the other two radiotracers. 
Although the metabolic mechanisms of the three radiotracers were very different, their underlying image characteristics were quite related as these correlations between different component functions reflected. 
In addition, for the jth patient, the estimated data-driven features {θ̂_jh^ tracer: h=1,…,H^*_ tracer} are plotted in the jth row, for j=1,…,22, in Figure <ref>, where the tracer is FDG, FLT, or Cu-ATSM. 
The patterns of the estimated data-driven features are similar across the three radiotracers.



 §.§ Comparison of Methods


We compared three sets of QIB features extracted from the 3D cancer image data: (1) features extracted by our proposed method based only on pre-treatment images; (2) the pre-treatment and mid-treatment features used in <cit.>; and (3) features extracted from a commonly used TA method as described in <cit.> based only on pre-treatment images. 
<cit.> found two features, the mid-treatment maximum FLT SUV and the relative change between pre- and mid-treatment mean FLT SUV, to be useful out of 26 candidate features considered in the association with a clinical end point, progression-free survival time, defined as the time from treatment initiation to any one of the events of metastasis, tumor growth, or recurrence. 
The third set by TA had 50 features, including 8 first-order features containing no information about the relative position of the voxels, 28 second-order features (with 23 features using a co-occurrence matrix and 5 features using the neighboring gray level), and 14 higher-order features (with 11 features employing a gray level run length matrix and 3 features using the neighborhood gray tone difference matrix). 
Table <ref> compares the distributions of absolute correlations between features extracted by our method and those by TA, showing clearly that our method had smaller correlations within and between features.

We also considered the association of three sets of features with progression-free survival time.
The first set based on our method had the data-driven features shown in Figure <ref>.
The second set extracted by <cit.> had the 26 features therein, while the third set by TA had 50×3=150 features from three radiotracers.
<cit.> only considered linear combinations of covariates in a Cox proportional hazards model.
To include nonlinearity, we applied a tree-based algorithm for survival data using the R package 𝚛𝚙𝚊𝚛𝚝 (). 
By the rule of thumb of needing 10 observations per covariate, we selected the best model based on a pseudo-R^2 among the Cox proportional hazards models with up to two features. 

The model fitting results are as follows. 
<cit.>'s method had a pseudo-R^2 value of 0.53 and the AIC value was 72.21. 
For the purpose of developing prognostic QIBs, we did not use the mid-treatment information for our method and the TA method. 
The TA method achieved a pseudo-R^2 value of 0.52 and the AIC value was 73.17, which were approximately the same as <cit.>, although with only pre-treatment information. 
Our proposed method had the highest pseudo-R^2 value of 0.63 and the smallest AIC of 67.16 among the three methods, indicating an improvement over the existing feature extraction methods.

As mentioned in Section <ref>, quantization could have a large impact on the features extracted by TA. 
We conducted a sensitivity analysis by recalculating the 150 features by the TA method with numbers of bins different from the default of 256 bins used above.
We examined the correlations between the features based on the default 256 bins and the recalculated features based on 128, 64, 32, and 16 bins, shown in Table <ref>.
Although the pseudo-R^2 and AIC values were similar across the numbers of bins, the recalculated features became less correlated with the features based on 256 bins, as the number of bins decreased.
These results demonstrated the challenges associated with quantization in the conventional TA method, as the results were sensitive to the choice of the bin number and it was unclear how to select an appropriate bin number. 






§ CONCLUSIONS AND DISCUSSION


In this paper, we have developed a general and flexible statistical approach for extracting QIBs from medical images in an objective and principled way. 
New statistical methodology has been developed based on spatial process decomposition where the random weights are unique to individual patients for component functions common across patients. Model fitting, model selection, and feature extractions have been demonstrated using both simulated image data and an actual cancer image data set. 
Medical images have inherent uncertainties that play an important role in the ability to develop useful QIBs. Recent works have shown that the conventional TA techniques are sensitive to image reconstruction techniques, patient motion during scanning, ROI definitions, among others <cit.>.
The method here is not a solution for all the uncertainties in PET QIB development, but does have advantages over the existing methods including approximate orthogonality among the extracted features, no need for ad-hoc quantization of features, and automatic tuning parameter selection. In addition, our method can handle non-rectangular ROIs with varying shapes among patients and applies to 1D curve, 2D pixel, and 3D voxel data in a unified way.





§ TABLES AND FIGURES



















