Identification of prognostic and predictive biomarkers in high-dimensional data with PPLasso

Let $\mathbf{y}$ be a continuous response or endpoint and $t_{1}$, $t_{2}$ two treatments. Let also $\mathbf{X}_{1}$ (resp. $\mathbf{X}_{2}$) denote the design matrix for the $n_{1}$ (resp. $n_{2}$) patients with treatment $t_{1}$ (resp. $t_{2}$), each containing measurements on $p$ candidate biomarkers:

To take into account the potential correlation that may exist between the biomarkers in the different treatments, we shall assume that the rows of $\mathbf{X}_{1}$ (resp. $\mathbf{X}_{2}$) are independent centered Gaussian random vectors with a covariance matrice equal to $\boldsymbol{\Sigma_{1}}$ (resp. $\boldsymbol{\Sigma_{2}}$).

To model the link that exists between $\mathbf{y}$ and the different types of biomarkers we propose using the following model:

where $(y_{i1},\dots,y_{in_{i}})$ corresponds to the response of patients with treatment $t_{i}$, $i$ being equal to 1 or 2,

with $\alpha_{1}$ (resp. $\alpha_{2}$) corresponding to the effects of treatment $t_{1}$ (resp. $t_{2}$). Moreover, $\boldsymbol{\beta}_{1}=(\beta_{11},\beta_{12},\ldots,\beta_{1p})^{\prime}$ (resp. $\boldsymbol{\beta}_{2}=(\beta_{21},\beta_{22},\ldots,\beta_{2p})^{\prime}$) are the coefficients associated to each of the $p$ biomarkers in treatment $t_{1}$ (resp. $t_{2}$) group, ^′ denoting the matrix transposition and $\epsilon_{11},\dots,\epsilon_{2n_{2}}$ are standard independent Gaussian random variables independent of $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$. When $t_{1}$ stands for the standard treatment or placebo, prognostic biomarkers are defined as those having non-zero coefficients in $\boldsymbol{\beta}_{1}$. According to the definition of prognostic biomarkers, their effect should indeed be demonstrated in the absence of therapy or with a standard therapy that patients are likely to receive. On the other hand, predictive biomarkers are defined as those having non-zero coefficients in $\boldsymbol{\beta}_{2}-\boldsymbol{\beta}_{1}$ because they aim to highlight different effects between two different treatments.

Model (2) can be written as:

with $\boldsymbol{\gamma}=(\alpha_{1},\alpha_{2},\boldsymbol{\beta}_{1}^{\prime},\boldsymbol{\beta}_{2}^{\prime})^{\prime}$. The Lasso penalty is a well-known approach to estimate coefficients with a sparsity enforcing constraint allowing variable selection by estimating some coefficients by zero. It consists in minimizing the following penalized least-squares criterion (Tibshirani (1996)):

where $\left\lVert\textbf{u}\right\rVert_{2}^{2}=\sum_{i=1}^{n}u_{i}^{2}$ and $\left\lVert\textbf{u}\right\rVert_{1}=\sum_{i=1}^{n}|u_{i}|$ for $\textbf{u}=(u_{1},\dots,u_{n})$. A different sparsity constraint was applied to $\boldsymbol{\beta}_{1}$ and $\boldsymbol{\beta}_{2}-\boldsymbol{\beta}_{1}$ to allow different sparsity levels. Hence we propose to replace the penalty $\lambda\left\lVert\boldsymbol{\gamma}\right\rVert_{1}$ in (4) by

Thus, a first estimator of $\boldsymbol{\gamma}$ could be found by minimizing the following criterion with respect to $\boldsymbol{\gamma}$:

where $D_{1}=[\mathbf{\textrm{Id}}_{p},\mathbf{0}_{p,p}]$ and $D_{2}=[-\mathbf{\textrm{Id}}_{p},\mathbf{\textrm{Id}}_{p}]$, with $\mathbf{\textrm{Id}}_{p}$ denoting the identity matrix of size $p$ and $\mathbf{0}_{i,j}$ denoting a matrix having $i$ rows and $j$ columns and containing only zeros. However, since the inconsistency of Lasso biomarker selection is originated from the correlations between the biomarkers, we propose to remove the correlation by “whitening” the matrix $\mathbf{X}$. More precisely, we consider $\widetilde{\mathbf{X}}=\mathbf{X}\boldsymbol{\Sigma}^{-1/2}$, where

and define $\boldsymbol{\Sigma}^{-1/2}$ by replacing in (7) $\boldsymbol{\Sigma}_{i}$ by $\boldsymbol{\Sigma}_{i}^{-1/2}$, where $\boldsymbol{\Sigma}_{i}^{-1/2}=\mathbf{U}_{i}\mathbf{D}_{i}^{-1/2}\mathbf{U}_{i}^{T}$, $\mathbf{U}_{i}$ and $\mathbf{D}_{i}$ being the matrices involved in the spectral decomposition of $\boldsymbol{\Sigma}_{i}$ for $i=1$ or 2. With such a transformation the columns of $\widetilde{\mathbf{X}}$ are decorrelated and Model (3) can be rewritten as follows:

where $\widetilde{\boldsymbol{\gamma}}=\boldsymbol{\Sigma}^{1/2}\boldsymbol{\gamma}$. The objective function (6) thus becomes:

Let us define a first estimator of $\widetilde{\boldsymbol{\gamma}}=(\widetilde{\alpha}_{1},\widetilde{\alpha}_{2},\widetilde{\boldsymbol{\beta}}_{1}^{\prime},\widetilde{\boldsymbol{\beta}}_{2}^{\prime})$ as follows:

for each fixed $\lambda_{1}$ and $\lambda_{2}$. To better estimate $\widetilde{\boldsymbol{\beta}}_{1}$ and $\widetilde{\boldsymbol{\beta}}_{2}$, a thresholding was applied to $\widehat{\widetilde{\boldsymbol{\beta}}}_{0}(\lambda_{1},\lambda_{2})=(\widehat{\widetilde{\boldsymbol{\beta}}}_{10}(\lambda_{1},\lambda_{2})^{\prime},\widehat{\widetilde{\boldsymbol{\beta}}}_{20}(\lambda_{1},\lambda_{2})^{\prime})^{\prime}$. For $K_{1}$ (resp. $K_{2}$) in $\{1,\ldots,p\}$, let $\textrm{Top}_{K_{1}}$ (resp. $\textrm{Top}_{K_{2}}$) be the set of indices corresponding to the $K_{1}$ (resp. $K_{2}$) largest values of the components of $|\widehat{\widetilde{\boldsymbol{\beta}}}_{10}(\lambda_{1},\lambda_{2})|$ (resp. $|\widehat{\widetilde{\boldsymbol{\beta}}}_{20}(\lambda_{1},\lambda_{2})|$), then the estimator of $\widetilde{\boldsymbol{\beta}}=(\widetilde{\boldsymbol{\beta}}_{1}^{\prime},\widetilde{\boldsymbol{\beta}}_{2}^{\prime})$ after the correction is denoted by $\widehat{\widetilde{\boldsymbol{\beta}}}(\lambda_{1},\lambda_{2})=(\widehat{\widetilde{\boldsymbol{\beta}}}_{1}^{(\widehat{K}_{1})}(\lambda_{1},\lambda_{2}),\widehat{\widetilde{\boldsymbol{\beta}}}_{2}^{(\widehat{K}_{2})}(\lambda_{1},\lambda_{2}))$ where the $j$th component of $\widehat{\widetilde{\boldsymbol{\beta}}}_{i}^{(K_{i})}(\lambda_{1},\lambda_{2})$, for $i=1$ or 2, is defined by:

Note that the corrections are only performed on $\widehat{\widetilde{\boldsymbol{\beta}}}_{0}$, the estimators $\widehat{\widetilde{\alpha}}_{1}$ and $\widehat{\widetilde{\alpha}}_{2}$ were not modified. The choice of $K_{1}$ and $K_{2}$ will be explained in Section 2.4.

To illustrate the interest of using a thresholding step, we generated a dataset based on Model 3 with parameters described in Section 3.1 and $p=500$. Moreover, to simplify the graphical illustrations, we focus on the case where $\lambda_{1}=\lambda_{2}=\lambda$. Figure 1 displays the estimation error associated to the estimators of $\widetilde{\boldsymbol{\beta}}(\lambda)$ before and after the thresholding. We can see from this figure that the estimation of $\widetilde{\boldsymbol{\beta}}(\lambda)$ is less biased after the correction. Moreover, we observed that this thresholding strongly improves the final estimation of $\boldsymbol{\gamma}$ and the variable selection performance of our method.

With $\widehat{\widetilde{\boldsymbol{\beta}}}=(\widehat{\widetilde{\boldsymbol{\beta}}}_{1}^{\prime},\widehat{\widetilde{\boldsymbol{\beta}}}_{2}^{\prime})$, the estimators of $\boldsymbol{\beta}_{1}$ and $\boldsymbol{\beta}_{2}-\boldsymbol{\beta}_{1}$ can be obtained by $\widehat{\boldsymbol{\beta}}_{10}=\boldsymbol{\Sigma}_{1}^{-1/2}\widehat{\widetilde{\boldsymbol{\beta}}}_{1}$ and $(\widehat{\boldsymbol{\beta}}_{20}-\widehat{\boldsymbol{\beta}}_{10})=\boldsymbol{\Sigma}_{2}^{-1/2}\widehat{\widetilde{\boldsymbol{\beta}}}_{2}-\boldsymbol{\Sigma}_{1}^{-1/2}\widehat{\widetilde{\boldsymbol{\beta}}}_{1}$. As previously, another thresholding was applied to $\widehat{\boldsymbol{\beta}}_{10}$ and $\widehat{\boldsymbol{\beta}}_{20}$: for $i=1$ or 2,

for each fixed $\lambda_{1}$ and $\lambda_{2}$. The biomarkers with non-zero coefficients in $\widehat{\boldsymbol{\beta}}_{1}=\widehat{\boldsymbol{\beta}}_{1}^{(M_{1})}$ (resp. $\widehat{\boldsymbol{\beta}}_{2}^{(M_{2})}-\widehat{\boldsymbol{\beta}}_{1}^{(M_{1})}$) are considered as prognostic (resp. predictive) biomarkers, where the choice of $M_{1}$ and $M_{2}$ is explained in Section 2.4.

To illustrate the benefits of using an additional thresholding step, we used the dataset described in Section 2.2. Moreover, to simplify the graphical illustrations, we also focus on the case where $\lambda_{1}=\lambda_{2}=\lambda$. Figure 8 in the Supplementary material displays the number of True Positive (TP) and False Positive (FP) in prognostic and predictive biomarker identification with and without the second thresholding. We can see from this figure that the thresholding stage limits the number of false positives. Note that $\alpha_{1}$ and $\alpha_{2}$ are estimated by $\widehat{\widetilde{\alpha}}_{1}$ and $\widehat{\widetilde{\alpha}}_{2}$ defined in (10).

For each $(\lambda_{1},\lambda_{2})$ and each $K_{1}$, we computed:

where $\widehat{\widetilde{\boldsymbol{\gamma}}}^{(K1,K2)}(\lambda_{1},\lambda_{2})=(\widehat{\widetilde{\alpha}}_{1},\widehat{\widetilde{\alpha}}_{2},\widehat{\widetilde{\boldsymbol{\beta}}}_{1}^{(K_{1})^{\prime}},\widehat{\widetilde{\boldsymbol{\beta}}}_{2}^{(K_{2})^{\prime}})$ defined in (10) and in (11). It is displayed in the left part of Figure 2.

For each $\lambda_{1}$, $\lambda_{2}$ and a given $\delta\in(0,1)$, the parameter $\widehat{K_{2}}$ is then chosen as follows for each $K_{1}$:

The $\widehat{K}_{2}$ associated to each $K_{1}$ are displayed with ’*’ in the left part of Figure 2. Then $\widehat{K_{1}}$ is chosen by using a similar criterion:

The values of $\widetilde{\textrm{MSE}}_{(K_{1},\widehat{K}_{2})}(\lambda_{1},\lambda_{2})$ are displayed in the right part of Figure 2 in the particular case where $\lambda_{1}=\lambda_{2}=\lambda$, $\delta=0.95$ and with the same dataset as the one used in Section 2.2. $\widehat{K}_{1}$ is displayed with a red star. This value of $\delta$ will be used in the following sections. However, choosing $\delta$ in the range (0.9,0.99) does not have a strong impact on the variable selection performance of our approach.

The parameters $\widehat{M}_{1}$ and $\widehat{M}_{2}$ are chosen in a similar way except that $\widetilde{\textrm{MSE}}_{K_{1},K_{2}}(\lambda_{1},\lambda_{2})$ is replaced by $\widehat{\textrm{MSE}}_{M_{1},M_{2}}(\lambda_{1},\lambda_{2})$ where:

with $\widehat{\boldsymbol{\gamma}}^{(M_{1},M_{2})}(\lambda_{1},\lambda_{2})=(\widehat{\widetilde{\alpha}}_{1},\widehat{\widetilde{\alpha}}_{2},\widehat{\boldsymbol{\beta}}_{1}^{(M_{1})^{\prime}},\widehat{\boldsymbol{\beta}}_{2}^{(M_{2})^{\prime}})$ defined in (10) and (12). In the following, $\widehat{\boldsymbol{\gamma}}(\lambda_{1},\lambda_{2})=\widehat{\boldsymbol{\gamma}}^{(\widehat{M}_{1},\widehat{M}_{2})}(\lambda_{1},\lambda_{2})$.

As the empirical correlation matrix is known to be a non accurate estimator of $\boldsymbol{\Sigma}$ when $p$ is larger than $n$, a new estimator has to be used. Thus, for estimating $\boldsymbol{\Sigma}$ we adopted a cross-validation based method designed by Boileau et al. (2021) and implemented in the cvCovEst R package (Boileau et al., 2021). This method chooses the estimator having the smallest estimation error among several compared methods (sample correlation matrix, POET (Fan et al. (2013)) and Tapering (Cai et al. (2010)) as examples). Since the samples in treatments $t_{1}$ and $t_{2}$ are assumed to be collected from the same population, $\boldsymbol{\Sigma}_{1}$ and $\boldsymbol{\Sigma}_{2}$ are assumed to be equal.

For the sake of simplicity, we limit ourselves to the case where $\lambda_{1}=\lambda_{2}=\lambda$. For choosing $\lambda$ we used BIC (Bayesian Information Criterion) which is widely used in the variable selection field and which consists in minimizing the following criterion with respect to $\lambda$:

where $n$ is the total number of samples, $\textrm{MSE}(\lambda)=\|\mathbf{y}-\mathbf{X}\widehat{\boldsymbol{\gamma}}(\lambda)\|_{2}^{2}$ and $k(\lambda)$ is the number of non null coefficients in the OLS estimator $\widehat{\boldsymbol{\gamma}}$ obtained by re-estimating only the non null components of $\widehat{\boldsymbol{\beta}}_{1}$ and $\widehat{\boldsymbol{\beta}}_{2}-\widehat{\boldsymbol{\beta}}_{1}$. The values of the BIC criterion as well as those of the MSE obtained from the dataset described in Section 2.2 are displayed in Figure 3.

Table 2 in the supplementary material provides the True Positive Rate (TPR) and False Positive Rate (FPR) when $\lambda$ is chosen either by minimizing the MSE or the BIC criterion for this dataset. We can see from this table that both of them have TPR=1 (all true positives are identified). However, the FPR based on the BIC criterion is smaller than the one obtained by using the MSE.

All simulated datasets were generated from Model (3) where the $n_{1}$ ($n_{2}$) rows of $\mathbf{X}_{1}$ ($\mathbf{X}_{2}$) are assumed to be independent Gaussian random vectors with a covariance matrix $\boldsymbol{\Sigma}_{1}=\boldsymbol{\Sigma}_{2}=\boldsymbol{\Sigma}_{bm}$, and $\boldsymbol{\epsilon}$ is a standard Gaussian random vector independent of $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$. We defined $\boldsymbol{\Sigma}_{bm}$ as:

where $\boldsymbol{\Sigma}_{11}$ (resp. $\boldsymbol{\Sigma}_{22}$) are the correlation matrix of prognostic (resp. non-prognostic) biomarkers with off-diagonal entries equal to $a_{1}$ (resp. $a_{3}$). Morever, $\boldsymbol{\Sigma}_{12}$ is the correlation matrix between prognostic and non-prognostic variables with entries equal to $a_{2}$. In our simulations $(a_{1},a_{2},a_{3})=(0.3,0.5,0.7)$, which is a framework proposed by Xue and Qu (2017). We checked that the Irrepresentable Condition (IC) of Zhao and Yu (2006) is violated and thus the standard Lasso cannot recover the positions of the null and non null variables. For each dataset we assumed randomized treatment allocation between standard and experimental arm with a 1:1 ratio, i.e. $n_{1}=n_{2}=50$. We further assume a relative treatment effect of 1 ($\alpha_{1}=0$ and $\alpha_{2}=1$). The number of biomarkers $p$ varies from 200 to 2000. The number of active biomarkers was set to 10 (i.e. 5 purely prognostic biomarkers with $\boldsymbol{\beta}_{1j}=\boldsymbol{\beta}_{2j}=b_{1}=1$ $(j=1,...,5)$ and 5 biomarkers both prognostic and predictive with $\boldsymbol{\beta}_{1j}=b_{1}$ and $\boldsymbol{\beta}_{2j}=b_{2}=2$ $(j=6,...,10)$).

We considered several evaluation criteria to assess the performance of the methods in selecting the prognostic and predictive biomarkers: the $\textrm{TPR}_{\textrm{prog}}$ as the true positive rate (i.e. rate of active biomarkers selected) and $\textrm{FPR}_{\textrm{prog}}$ the false positive rate (i.e. rate of inactive biomarkers selected) of the selection of prognostic biomarkers, and similarly for predictive biomarkers with $\textrm{TPR}_{\textrm{pred}}$ and $\textrm{FPR}_{\textrm{pred}}$. We further note $\textrm{TPR}_{\textrm{all}}$ and $\textrm{FPR}_{\textrm{all}}$ the criterion of overall selection among all candidate biomarkers regardless their prognostic or predictive effect. The objective of the selection is to maximize the $\textrm{TPR}_{\textrm{all}}$ and minimize the $\textrm{FPR}_{\textrm{all}}$. All metrics were calculated by averaging the results of 100 replications for each scenario.

For the proposed method, different results are presented. $\textrm{PPLasso}_{\Sigma}$ (resp. PPLasso) corresponds to the results of the method by considering the true (resp. estimated) matrix $\boldsymbol{\Sigma}_{bm}$. For estimating $\boldsymbol{\Sigma}_{bm}$, we used the approach explained in Section 2.5. Two choices of parameters are also given: “optimal” and “min(bic)”. The former uses as a value for $\lambda$ the one that maximizes $(\textrm{TPR}_{\textrm{all}}-\textrm{FPR}_{\textrm{all}})$ and the latter uses the approach presented in Section 2.6. In order to compare the performance of our approach to the best performance that could be reached by Elastic Net, Lasso and Adaptive Lasso, we used for these methods the “optimal” parameters namely those maximizing $(\textrm{TPR}_{\textrm{all}}-\textrm{FPR}_{\textrm{all}})$. All these three methods were implemented with the glmnet R package, the best parameter $\alpha$ involved in Elastic Net was chosen in the set $\{0.1,0.2,\dots,0.9\}$. The choice of “min(bic)” is only applied to our method and corresponds to a choice of $\lambda$ that could be used in practical situations. For ease of presentation, the abbreviation EN (resp. AdLasso) refers to Elastic Net (resp. Adaptive Lasso) in the following.

Figure 4 shows the selection performance of PPLasso and other compared methods in the simulation scenario presented in Section 3.1. PPLasso achieved to select all prognostic biomarkers ($\textrm{TPR}_{\textrm{prog}}$ almost 1) even for large $p$, with limited false positive prognostic biomarkers selected. As compared to the optimal $\lambda$ maximizing $(\textrm{TPR}_{\textrm{all}}-\textrm{FPR}_{\textrm{all}})$, the one selected with the BIC tends to select some false positives (average: 33 ($\textrm{FPR}_{\textrm{prog}}=0.17$) for $p=200$ and 10 ($\textrm{FPR}_{\textrm{prog}}=0.005$) for $p=2000$). The results obtained from the oracle and estimated $\boldsymbol{\Sigma}_{bm}$ are comparable. Selection performance of predictive biomarkers is slightly lowered as compared to prognostic biomarkers. Even if the false positive selection is quite similar between prognostic and predictive biomarkers, PPLasso missed some true predictive biomarkers when $\lambda$ is selected with the BIC criterion (average $\textrm{TPR}_{\textrm{pred}}=$ 0.98 and 0.80 for oracle and estimated $\boldsymbol{\Sigma}_{bm}$, respectively, with $p=2000$). In this scenario where the IC is violated, PPLasso globally outperforms Lasso, Elastic Net and Adaptive Lasso. Although Elastic Net showed higher TPR than Lasso and Adaptive Lasso, they all failed in selecting all truly prognostic and predictive biomarkers, and the number of missed active biomarkers increased with the dimension $p$. For example, for Elastic Net, $\textrm{TPR}_{\textrm{prog}}$ = 0.85 and 0.53, $\textrm{TPR}_{\textrm{pred}}$ = 0.81 and 0.61 for $p=200$ and $2000$, respectively.