Regularized estimation in sparse high-dimensional multivariate regression, with application to a DNA methylation study

1 Introduction

With the development of modern technology for data collection, high-dimensional data have become increasingly common in many scientific research fields, e.g. genome-wide studies (Lin et al., 2015), biomedical sciences (Mukherjee et al., 2015), economics and finance (Basu and Michailidis, 2015). Under these situations, the number of parameters is larger than the sample size, rendering traditional statistical procedures inappropriate. More recently, correlated data with high dimensional multivariate responses are often encountered in omics studies. Our motivating example is the Normative Aging Study, where methylation markers are taken as the multivariate outcomes. The methylation of DNA, where methyl groups are added to DNA at binding sites typically referred to as cytosine-phosphate-guanine (CpG) islands, could affect the DNA expression. DNA methylation levels measured from probes close to one another are correlated (Moen et al., 2013), resulting in high dimensional multivariate outcomes. Correlation may also exist for DNA methylation markers having the similar function, e.g. related to exposure such as smoking. It is thus necessary to account for the within-subject correlation in the estimation procedure.

Our objective is to conduct selection of regression coefficients in high dimensional multivariate DNA methylation markers. There are two challenging issues for such a study: (1) how to conduct variable selection in the high dimensional setting; and (2) how to tackle the correlation among multivariate outcomes. For the first challenge, many studies have focused on penalized methods, such as the least absolute shrinkage and selection operator (LASSO, Tibshirani, 1996), the smoothly clipped absolute deviation (SCAD, Fan and Li, 2001), the elastic net (Zou and Hastie, 2005), the adaptive LASSO (Zou, 2006), and the minimax concave penalty (MCP, Zhang, 2010). The penalized approach has been applied to many research topics, e.g. linear models (Wang and Leng, 2007; Huang et al., 2011; Fan and Lv, 2014), generalized linear models (van de Geer, 2008; Jiang et al., 2016), survival models (Fan and Li, 2002; Bradic et al., 2011; Lin and Lv, 2013). Recently, Belloni et al. (2011) proposed a pivotal square-root LASSO method, which does not rely on the knowledge of the standard deviation for the error term. Later, Belloni et al. (2014) developed a self-tuning square-root LASSO method in high-dimensional nonparametric regression analysis. Liu and Wang (2017) proposed a new procedure for optimally estimating high dimensional Gaussian graphical models using the square-root LASSO. For more topics on variable selection, please refer to Bühlmann and van de Geer (2011).

There is limited research to tackle the second challenge, especially when the responses are high-dimensional. Rothman et al. (2010) proposed an iterative algorithm for variable selection and estimation in high-dimensional multivariate regression using the LASSO penalty. Sofer et al. (2014) considered variable selection for high-dimensional multivariate regression using the penalized likelihood method, but the dimensionality of the multiple responses is still much smaller than the sample size. Liu et al. (2015) proposed a calibrated multivariate regression method for high-dimensional multivariate regression models, but they only considered the uncorrelated error structure. Li et al. (2015b) and Wilms and Croux (2017) studied the group LASSO for high-dimensional multivariate linear regression model. Most of these existing methods use cross-validation to choose tuning parameters over a full regularization path, which are computationally expensive and may potentially waste valuable training data. To deal with this problem, we will extend Belloni et al. (2011)’s (unweighted) square-root LASSO on a single outcome to multivariate outcomes with weighted square-root LASSO. The main advantage of our procedure over existing methods comes from the tuning-insensitive property, which is significantly faster than cross-validation. Another advantage is that we can use the entire dataset for variable selection, which may potentially learn a better model (Bishop et al., 2003).

The rest of the paper is organized as follows. In Section 2, we introduce the model and the weighted square-root LASSO procedure for multivariate linear regression with high-dimensional responses. In Section 3, we establish an error bound for the proposed estimator. In Section 4, we develop an efficient algorithm and conduct Monte Carlo simulations to assess the performance of our method. An empirical analysis of DNA methylation in the Normative Aging Study is presented in Section 5. Some concluding remarks are given in Section 6. All technical proofs are relegated to the Appendix.

2 Model and estimation

Consider the sparse, high-dimensional multivariate linear regression model

where Yi=(Yi1,⋯,Yip)′ is a p-dimensional response, e.g. DNA methylation markers; Xi=(Xi1,⋯,Xiq)′ is a q-dimensional covariates vector; B=(B1,⋯,Bp)′∈ℝp×q is a sparse regression coefficient matrix with Bk=(βk1,⋯,βkq)′∈ℝq, k = 1, ⋯, p; ϵi=(ϵi1,⋯,ϵip)′∈ℝp is a random error term with mean 0 and covariance matrix Σ_p. Throughout this article, we assume that the predictor’s dimension q is fixed but the dimension of response p could be larger than n. The sparse regression coefficient matrix B suggests that only a small number of coefficients are non-zero. Our interest is to estimate the coefficient matrix B and establish oracle inequalities for corresponding estimators, while account for the correlated outcomes.

Assume that (Y_i, X_i) are independently and identically distributed (i.i.d.) observations, i = 1, ⋯, n. If ϵ_i ∼ N(0, Σ_p), the negative log-likelihood function is given by

where Ωp=Σp−1 is the precision matrix (Cai et al., 2011). Denote (B1′,⋯,Bp′)′ as β=(β1,⋯,βd)′, where d = pq. Let 𝒮={j;βj≠0} be the true model with size s=|𝒮|. We use Ω_p to account for the within-subject correlation (Sofer et al., 2014), and propose the following criterion function:

In practice, the precision matrix Ω_p can be estimated using the constrained ℓ₁ minimization method (Cai et al., 2011), which has been implemented in the R package flare (Li et al., 2015a). Basically, the estimation of the sparse inverse covariance matrix (precision matrix) Ω_p can be obtained by the following optimization problem:

subject to:

where γ > 0 is a tuning parameter and Σn=1n∑i=1n(Yi−B~Xi)′(Yi−B~Xi) with B~ being a consistent estimator (e.g. ridge estimator).

Belloni et al. (2011) proposed an unweighted square-root LASSO for β with a scalar outcome. They showed that the penalty (tuning parameter) is pivotal, i.e., it does not rely on the knowledge of the error variance, nor do we need to pre-estimate it. Consequently, the estimation of β is insensitive to the tuning parameter, greatly simplifying the computation which often resorts to cross validation for tuning parameter selection. Furthermore, the square-root LASSO method achieves near-oracle performance for the estimation of β. In comparison, the ordinary LASSO needs to estimate the error variance (Belloni et al., 2011) to achieve near-oracle performance, which is a very challenging issue in high dimensional data.

Motivated by Belloni et al. (2011), the corresponding weighted square-root LASSO version for correlated outcomes in Model (1) is defined as

where λ > 0 is the tuning parameter, w_j is a known weight, j = 1, ⋯, d. We can set wj=1/|β~j| along the lines of Zou (2006) with β~_j being the ridge estimator, j = 1, ⋯, d. Of note, we choose the ridge estimator rather than the ordinary least square estimator in the adaptive LASSO since p could be larger than n in Model (1). To obtain β^ in (4), we consider the optimization problem,

For ρ ≥ 0, we have Q(β;Ωp)2ρ+ρ2≥Q(β;Ωp), so the objective function in (5) is an upper bound of that in (4). The equality is attained if and only if ρ=Q(β;Ωp). Similar to Proposition 3.1 of Liu and Wang (2017), the optimizations in (4) and (5) yield the same solution β^. This relationship between (4) and (5) provides an efficient algorithm as described below.

For given λ, we have the following procedure:

Step 0. Compute the precision matrix Ω^p (using R package flare; Li et al., 2015a) and β^ridge (using R package glmnet; Friedman et al., 2010). Set β(0)=β^ridge and ρ(0)=Q(β(0);Ω^p).

Step 1. Solve the optimization problem via coordinate descent algorithm (Friedman et al., 2010)

Step 2. Update

Step 3. Repeat Steps 1 and 2 until convergence.

Of note, Ω^_p is consistent (Cai et al., 2011) and kept unchanged in the iteratively updated procedure, so the objective function in Step 1 is convex for β, which ensures fast convergence of the algorithm. The following Lemma 1 will present an explicit expression for λ, while the optimal value for the tuning parameter in (4) will be evaluated in Section 4 via empirical studies. The estimation procedure has been implemented in R (available upon request).

3 Theoretical results

In this section, we will establish the oracle inequality for β^ defined in (4). The following lemma lays the foundation for the tuning-insensitive property of the weighted square-root LASSO procedure, which is motivated by Bickel et al. (2009)’s choice on the penalty level for LASSO. Denote Wmin=min{w1,⋯,wd}, Wmax=max{w1,⋯,wd} and N = np. We first have the following lemma.

Similar to Belloni et al. (2011), the prediction norm is given as ||e||2,N2=1Ne′𝕏′𝕏e, where 𝕏 is defined in the Appendix and e ∈ ℝ^d. To derive the oracle inequalities, we define the compatibility factor (Huang et al., 2013) and restricted eigenvalue (RE, Bickel et al., 2009) as

respectively, where 𝒞(ξ,𝒮)={e∈ℝd:||e𝒮c||1≤ξ||e𝒮||1} with ξ=cWmax+WminWmin(c−1). Hereafter, 𝒜^c denotes the complement of set 𝒜; v𝒜=(vj:j∈𝒜) for a vector v. Then, we have the following conclusion.

The following theorem gives the upper bound for the estimation error.

4 Simulation studies

In this section, we will conduct simulation studies to validate the proposed methodology. We assume that X is from N_q(0, Σ_X), where q = 2 and Σ_X = (σ_ij) is given by σ_ij = 0.7^|i−j|. The random error term is generated from N_p(0, Σ_E), we consider the following two settings for the error covariance:

Case (a). AR(1) error covariance: Σ_{E, st} = 0.8^|s−t|.

Case (b). Fractional Gaussian Noise (FGN) error covariance:

with Hurst parameter H = 0.9, so the correlation is (0.74, 0.63, 0.58, 0.55, 0.52, 0.50, 0.49, 0.48, 0.46, 0.45) for distance |s−t|=1:10. It is noted that the inverse error covariance for Case (a) is a tri-diagonal sparse matrix, while Case (b) has a dense inverse error covariance. Set β = (1, 0.6, 0.3, 1.2, 0.8, 0.5, 0, ⋯, 0)^′, i.e., the first six elements are non-zero, while the rest are all 0. We take p = 200 and 300, respectively. All simulation results are based on 100 replications with n = 100 and 200.

According to Theorem 1, the tuning parameter λ = ηWminlogdN with some η > 2. To choose the suitable λ, we consider various values of η, and the corresponding performance is plotted in Figures 1–4. We can see that the tuning parameter’s effect is limited and the procedure has better properties when η ∈ [6, 10]. Thus, we suggest η = 8 for the weighted square-root LASSO (WSR-LASSO) method in both simulation and real application. For comparison, we also consider the LASSO and square-root LASSO (SR-LASSO) with weight w_j = 1 in (4).

Figure 1

Correct model rate for WSR-LASSO with AR(1) error covariance.

Figure 2

Model error for WSR-LASSO with AR(1) error covariance.

Figure 3

Correct model rate for WSR-LASSO with FGN error covariance.

Figure 4

Model error for WSR-LASSO with FGN error covariance.

Tables 1 and 2 report the results, which include the rate that the correct model (CMR) is selected I{𝒮^=𝒮}, the false positive rate (FPR) |𝒮^\𝒮|/|𝒮^|, the false negative rate (FNR) |𝒮\𝒮^|/(d−|𝒮^|), and the model error (ME) tr[(B^−B)ΣX(B^−B)′] (Yuan and Lin, 2007). It can be seen that LASSO and SR-LASSO have similar performance, which is in line with the conclusion of Belloni et al. (2011). Moreover, the WSR-LASSO method has a higher rate of selecting the correct model and a smaller model error than LASSO and SR-LASSO. The results on the false positive and negative rates also suggest that the proposed method is preferred over LASSO and SR-LASSO in practice.

5 Application

We apply our proposed methodology to the DNA methylation (DNAm) data from the US Department of Veterans Affairs’ Normative Aging Study (NAS). We exclude participants who (i) were non-white or had missing information on race to minimize potential confounding effects of genetic ancestry, or (ii) had any cancer diagnosed and history of stroke or coronary artery disease as their blood methylation profiles could have been affected. A total of 169 individuals with samples collected at their first blood draw remain for analysis.

We are interested in the effect of smoking on DNA methylation. Gao et al. (2015) conducted a literature review on the DNA methylation change in response to active smoking exposure in adults. From it, we consider methylation markers at a total of 151 cytosine-phosphate-guanine (CpG) dinucleotides which had been reported multiple (≥2) times in the literature. In our studies, the correlations among the 151 CpG sites have a range of [−0.6440, 0.9369], manifesting high correlations among these CpGs.

We are interesting in the smoking pack year (packyr)’s effect on these DNAm markers. We also include age and BMI in the model. In total, we need to estimate 453 (151 × 3) regression coefficients. We use the proposed method in Section 2 with tuning parameter λ = 0.0026.

In Table 3, we compare our results to the original 151 CpGs listed by Gao et al. (2015). In Table 4 we report the selected CpGs and coefficient estimates. Thirty three CpGs among a total of 151 are selected by our method in the NAS data. We can see that CpGs reported more frequently in the literature are also more likely to be chosen by our method in the NAS data. For example, among 8 CpGs reported at least 7 times, 5 (62.5%) are selected by our method from the NAS data, while only 15 out of 89 (16.9%) CpGs reported two times in literature are selected by our method. The decreasing trend in Table 3 shows the consistency of our method with the literature. Of note, our method correctly identifies the top two CpGs – cg03636183 and cg05575921 (located in F2RL3 and AHRR genes), which have been reported 12 and 11 times in literature, respectively.

6 Concluding remarks

We proposed a weighted square-root LASSO method for high-dimensional multivariate regression models. We estimated the precision matrix by the CLIME method to account for the correlations between responses and obtained oracle inequalities for the estimator. Simulation studies were provided to illustrate the proposed procedure. We applied the method to study the relation between smoking and high dimensional DNA methylation markers.

There exist several topics to research in the future. First, the estimation of high-dimensional Gaussian graphical models is an active area of research (Cai et al., 2011; Cai and Yuan, 2012; Fan et al., 2013). It is of great interest to consider the joint estimation of regression coefficients and the precision matrix in (1). A possible solution is to add a penalty term in (4) for the elements of the precision matrix. Second, although it is assumed that the dimension of the covariates q is fixed, it is straightforward to extend the proposed procedure to the high-dimensional covariates setting, the main difficulty lies in the computational burden due to the ultra-high dimensional parameters. Third, statistical inference on the weighted square-root LASSO is an important and interesting topic. Fourth, since our method obviates the burden of cross-validation, our method is computationally efficient: for the DNA methylation data in Section 5, it took R software about 10.3 s to converge in a personal computer. As a reviewer suggested, it would be of interest to make the proposed algorithm scalable to genome-wide response and predictor markers, which can be implemented in high performance computing facilities. Fifth, we are interested in high-dimensional mediation analysis (Zhang et al., 2016) to determine whether high-dimensional DNA methylation markers mediate the path from intervention (e.g. diet, physical exercise) to health outcomes.