It is common to find multi-response data in many real life problems. Componentwise analysis is clearly not a good choice for multi-response data analysis, because it does not fully make use of the information available. For example, the observations of other components may contain the information for the component of interest, and such information would be completely overlooked by component-wise analysis, therefore, the resulting estimators would not be as efficient as we can expect. It is necessary to take multivariate analysis approach for multi-response data analysis. The most commonly used multivariate regression models are the multivariate linear models. The research in the multivariate linear models can be at least traced back to Anderson (1951). There is much literature after Anderson (1951) about the classic multivariate linear models, see the references in Reinsel and Velu (1998) and Anderson (2004).
With the surge in high dimensional data analysis in the past more than a decade, the multivariate linear models in high dimensional setting are attracting more and more attention than ever before. Many interesting developments in low rank high dimensional multivariate linear models have appeared in literature, see Yuan et al. (2007), Negahban andWainwright (2011), Obozinski et al. (2011), Kong et al. (2017), Bing and Wegkamp (2019), Raskutti et al. (2019), Zheng et al. (2019) and the references therein.
A commonly used approach to deal with the low rank coefficient matrix in a multivariate linear model is based on the idea of decomposing the coefficient matrix, say A with rank r and size p × q, to CDQ, where C and Q are two matrices of size p × r and r × q, respectively, and D is a diagonal matrix of size r, where r is an estimator of r. The estimation of r plays a key role for the success of the approach used. Different approaches may end up with different ways to estimate r, see Yuan et al. (2007) and Bing and Wegkamp (2019). Although the existing approaches for estimating r enjoy nice asymptotic properties, when implementing them, we often come up against a dilemma: we create a new unknown parameter in order to estimate an unknown parameter, this is because we have to select a tuning parameter in the estimation of r. In addition to that, as far as the estimation of A is concerned, which is the ultimate goal for multivariate linear models, even if we knew the rank r, in order to get the estimator of A based on the decomposition, we would have to estimate C, Q and D. Even with the constraints coming with the decomposition, we may have to estimate at least (p + q)r unknown parameters, which is more than the unknown parameters we need to estimate without using the decomposition when r > pq/(p + q). That implies we may end up with a better estimator of A if we simply apply the standard least squares estimation for multivariate linear models when r > pq/(p + q), which clearly shows the limitation of the decomposition based approach.
In this paper, we are going to propose an estimation procedure for the low rank multivariate linear models, in which we only need to estimate r(p + q) -r 2 unknown parameters in order to get the estimator of A. We can easily show r(p + q) -r 2 ≤ pq, because r ≤ min(p, q). Intuitively speaking, the proposed estimation procedure would be more efficient than either of the standard least squares estimation and the decomposition based approach. This conclusion is confirmed to be true by both the asymptotic theory established in Section 4 and simulation studies in Section 5. As part of the proposed estimation procedure, the rank of A is estimated by the BIC, which is free of tuning parameter.
We will show the resulting estimator enjoys good theoretical properties and performs well in simulation studies.
Another advantage of the proposed estimation procedure is it clearly appreciates the high dimensionality by directly imposing a penalty on those entries of A in question, which makes the proposed estimation procedure easily accom-modate the high dimensional cases and enjoy the function of feature selection.
In the context of multi-response data analysis, the proposed estimation procedure also comes with a very nice practical implication, which is the impacts of explanatory variables on some responses are linear combinations of the impacts on certain responses when the matrix coefficient is of low rank, this would be very helpful when it comes to interpreting the results for a given real dataset, and may lead to some interesting findings in the discipline which the dataset comes from.
To implement the proposed estimation procedure, we have also developed an algorithm for the estimation. Our simulation studies show the proposed algorithm is fast and accurate.
The rest of this paper is organised as follows: we begin with a detailed description of the models we are going to address in Section 2. The proposed estimation procedure and associated computational algorithm are described in Section 3. The asymptotic properties of the estimators obtained by the proposed estimation procedure are presented in Section 4. Section 5 is devoted to simulation studies, in which we will examine how well the proposed estimation works. Finally, in Section 6, we apply the low rank multivariate linear models together with the proposed estimation procedure to analyse an environmental data set and predict the concentrations of PM2.5 at the locations concerned.
The results show the proposed method provides more accurate prediction than other methods. We leave the theoretical proofs of all asymptotic properties in the Appendix.
To give a generic description of the models we are going to address, we use Y to denote the vector of all response variables, X the vector of all covariates.
Without any confusion, from now on, we call Y the response variable, X the covariate. We assume Y is of q dimension, X is of p dimension. p and q may tend to ∞ when sample size tends to ∞. The low rank high dimensional multivariate linear models which we are going to address in this paper are
where A is a p × q unknown matrix of unknown rank r, r < q < p, =
( 1 , . . . , q ) is a q dimensional random error, and
Like Yuan et al. (2007), Bing and Wegkamp (2019) and Raskutti et al. (2019),
we assume Σ = σ 2 I q , and σ 2 is unknown.
Suppose we have a sample
model for the sample can be written as
Throughout this paper, for any matrix Ω = (w ij ) of size p × q and any vector
For any integers 1
and λ min (B) and λ max (B) be the smallest and largest eigenvalues of a square matrix B.
Suppose we have an independent and identically distributed sample (Y i , X i ), i = 1, • • • , n, from (Y , X ). A standard penalised least squares estimation would provide us with an estimator à of A, which is the minimiser of
where P λ (•) is a penalty function. However, this estimation does not take into account the information that A is of low rank, which would result in an estimator not as efficient as we could expect. In fact, it is easy to see this estimation is equivalent to componentwise penalised least squares estimation for (2.1).
The proposed estimation procedure will fully make use of the low rank information of A, and the resulting estimator will be more efficient.
The idea, based on which the proposed estimation is constructed, is that each column of A is a linear combination of r linearly independent columns of A.
Based on this idea, we propose the following estimation procedure for A. We start with the case when the rank r of A is known, then propose an estimation for r.
Apply the idea of penalised least squares estimation and minimise
with respect to a j l s, b k s, and
is a penalty function, λ involved is a tuning parameter which can be selected by some criterion, such as BIC.
When r is large, b k s are also likely to have sparsity, in which case, we can add another penalty term into (3.2) to penalise b k s. However, when r is small, which is the case of main interest, there is no need to penalise b k s, this is because for each component of Y i , say the kth component, there are only r b k s, which is not many.
Notice that the minimiser of (3.2) is not unique. We denote a minimiser of (3.2) by
For any j, 1 ≤ j ≤ q, the jth column of A is estimated by
We use Â(r) to denote the estimator of A.
The non-uniqueness of the minimiser of (3.2) is because that there can be more than one ways to choose r independent columns D r so that A Dr is fullrank. Theoretically speaking, as long as A Dr and r can be well estimated, we can recovery the low-rank structure of A regardless of the choice of D r . Our theory shows that the consistency of our proposed estimator holds uniformly in D r ; please refer to Lemmas 1-3 in the Appendix. Hence, this non-unique issue does not affect the performance of the proposed estimation procedure, which is further corroborated via extensive simulations in Section 5.
The estimation of A in section 3.1.1 is built on the assumption that the rank r of A is known, and this assumption is not realistic in reality. In fact, rank r plays a very important role in the estimation of A. If r is underestimated, a substantial bias would creep into the estimation procedure and make the final estimator of A very biased. On the other hand, if r is overestimated, we would have to estimate unnecessarily many unknown parameters, which would make the final estimator of A have big variance. In this paper, we use BIC, which is defined as follows, to estimate r
where h n is a positive diverging sequence, which can be set to be log n. The estimator of r is given by
where r is a pre-specified bound for r. The proposed estimator Â(r) of A is the Â(r), obtained in section 3.1.1, with r being replaced by r.
We will show in Section 4 the proposed BIC estimator r enjoys an excellent asymptotic property, say it tends to identify the true model consistently. If the prediction accuracy is our primary concern, we can consider the multifold crossvalidation (CV), which tends to select the model with the optimal prediction performance (Zhang, 1993). The data are splitted randomly into M groups of equal sizes (assuming that n/M is an integer for simplicity), G m , m = 1, . . . , M .
For each k, 1 ≤ k ≤ q, let Â-m (k) be the estimator of A, obtained by the method in Section 3.1.1 when the rank of A is k, without using the observations of the mth group. The cross-validation sum is defined as
The CV estimator of r is taken to be the minimiser of CV(k).
The minimisation of (3.2) can be difficult. We propose an iterative algorithm to solve this problem. The route of our algorithm is: we first minimise (3.2), for given {j 1 , • • • , j r }, with respect to a j l s, b k s, and denote the resulting
The details of our algorithm are described as follows.
For any given {j 1 , • • • , j r }, we minimise (3.2) with respect to a j l s and b k s by the following iterative approach
(1) We minimise
with respect to a j l s, and denote the minimiser by a (0) j l s. There are many existing methods to do the minimisation in this step, because this is the minimisation for standard penalised least squares estimation.
(2) Minimise
with respect to b k s, and denote the minimiser by b
k enjoys a closed form, therefore, the minimisation in this step is very easy.
(3) Let a k s are the minimser of (3.2), and the minimum of (3.2) is denoted by
However, this approach would have to compute q r F (j 1 , • • • , j r )s, which is computationally too expensive. We shall borrow the idea of forward selection to minimise F (j 1 , • • • , j r ), which is depicted as follows (I) Let F (j 1 ) be F (j 1 , • • • , j r ) when r = 1, and compute F (j 1 ) for each possible j 1 , 1 ≤ j 1 ≤ q. Let ĵ1 be the one which minimises F (j 1 ).
(II) For any k < r, when we have { ĵ1 , • • • , ĵk }, the way to select a j k+1 from
as follows: for each possible j k+1 , we arrange ĵ1 , • • • , ĵk and j k+1 in ascent order, and denote them by j1 < • • • < jk+1 . We compute F ( j1 , • • • , jk+1 ). The selected j k+1 is the one which minimises
). We add the selected j k+1 into the set { ĵ1 , • • • , ĵk }, and sort the elements in the new set in ascent order. With a little bit abuse of notation, we denote the new set by
2), and minimise (3.2) with respect to a ĵl s and b k s. Denote the resulting minimiser by âĵ l s and bk s. We take
as a minimiser of (3.2) with respect to {j 1 , • • • , j r }, a ĵl s and b k s.
In this section we are going to investigate the asymptotic behavior of the proposed estimator of A.
with probability tending to 1. " " and " " are similarly defined.
To make the theoretical derivation more neat, we write the minimisation of (3.2) in matrix form. Specifically, for any given integer k ∈ [1, q), when r = k, the minimisation of (3.2) can be written to the minimisation of the following objective function
Without loss of generality, we use the (adaptive) lasso type penalty
see Zou (2006), where λ j > 0 and u j are the (j, )th element of U.
Throughout this section, we assume that each column of X has been normalized to have L 2 -norm of n. Furthermore, we denote
where a j is the (j, )th entry of A.
In order to establish the asymptotic properties of the proposed methods, we impose the following technical conditions:
Condition 1 There exist positive constants κ and κ such that with probability
Condition 2 For some positive constants C and K,
Condition 3 The elements of A and XA are bounded. The matrix A is s nsparse in the sense that max 1≤ ≤q
where
Remark 1 Condition 1 implies that the predictor matrix has a reasonably good behavior; this is a type of restricted eigenvalue assumption and is commonly used in the literature, e.g., Fan and Peng (2004). Condition 2 requires that each entry of E is sub-Gaussian, which ensures its tail probability decays exponentially. Condition 3 facilitates our derivation but can be much relaxed so that the true signal strength depends on n as well. Condition 4 imposes requirement on the diverging rate in order to obtain consistent estimation when p, r diverges with n. Condition 5 is an identifiability assumption, ensuring that a true low-rank structure can be recognized.
We start with the establishment of the asymptotic property of the minimiser Â(r) = ( Û, V) of (4.1) when k = r. As far as the asymptotic properties are concerned, p, q and r are allowed to depend on n and diverge as n → ∞. The reason for us to suppress the subscript n of p, q and r is to make notations neat.
We define an index set
and its complement set is denoted by M c (D). We have the following theorem:
Theorem 1 Under Conditions 1-5, if γ 0n / rp log p/n → ∞ and rs n γ 2 1n → 0 as n → ∞, there exists, with probability tending to one, a local minimiser { Dr , Â(r)} of L(U, V, D r ; r) satisfying: A Dr is full-rank, and
Theorem 1 implies that we can identify a "correct" Dr in the sense that A Dr is full rank with an overwhelming probability when n is large. The penalised estimators of the zero coefficients are exactly zero under some conditions on γ 0n .
The condition rs n γ 2 1n → 0 together with Condition 4 ensures that the proposed estimator is consistent. From the proof of this theorem in the Appendix, we can see, as a special case, that Â(r) -A = O p ( r(p + q) log p/n) when γ 0n = γ 1n = 0. This is in line with the relevant existing results, see, for example, Negahban and Wainwright (2011).
When the coefficient matrix is sparse, under properly selected tuning parameters, the proposed penalised estimator would enjoy the "oracle property". Specifically, if the adaptive lasso penalty with tuning parameter λ j = λ n ã-1 j is used, where à = (ã jl ) p×q is the standard least-squares estimator of A, we can verify that Â(r) -A = O p ({r(s n + q)(log p/n)} 1/2 ), provided that rqp 2 log p/n → 0 as n → ∞, by setting λ n ∼ log p/n and using the fact that ãj = O p (1) for a j = 0 and ãj = O p ( pq log p/n) for a j = 0.
As mentioned before, the estimation of the rank r of A plays a very important role in the estimation procedure of A. Theorem 2 shows that the proposed BIC estimator r, defined by (3.3), is consistent.
It is very easy to see Theorem 1 together with Theorem 2 imply the proposed estimator of A is consistent.
In this section, we use two simulated examples, low and high dimensional ones, to assess the coefficient matrix estimation, the prediction, and the rank recovery performance of the proposed method. We assess the proposed method and consider the rank being estimated by either 5 fold cross-validation (OUR-CV;
Eq.(3.4)) or BIC (OUR-BIC; Eq.(3.3)). Throughout this section, OUR-CV and OUR-BIC are implemented by the algorithm introduced in Section 3.2 with P λ (•) being the L 1 penalty function. The tuning parameter λ is selected by 5 fold cross-validation.
In the low dimensional example, we compare OUR-CV and OUR-BIC with the Factor Estimation and Selection method (FES) proposed in Yuan et al. (2007), the Rank Selection Criterion (RSC) proposed in Bunea et al. (2011), and the Self-Tuning Rank Selection (STRS) proposed in Bing and Wegkamp (2019). Besides, we also consider the Ordinary Least Squares estimator (OLS) and a Low-rank Matrix Decomposition estimator (LMD) as two benchmarks.
Denote ÂOLS the OLS estimator of A, the LMD estimator of A is obtained from a rank r truncated singular value decomposition of ÂOLS as
where D = diag{σ 1 , . . . , σ r} is a diagonal matrix of r largest positive singular values of ÂOLS , and U = (u 1 , . . . , u r ) and V = (v 1 , . . . , v r) are corresponding left and right singular vectors of ÂOLS , respectively. Further, r is estimated by the eigen-ratio method, see Ahn and Horenstein (2013) and references therein.
In the high dimensional example, we only compare OUR-CV and OUR-BIC with RSC and STRS as FES, OLS and LMD are not applicable.
Consider the multivariate linear regression model (2.2). Similar to Bunea et al.
(2011) and Bing and Wegkamp (2019), we consider a data generating process as follows.
(1) Coefficient matrix A:
and r ≤ min(p, q). The entries of Γ 0 and Γ 1 are independently drawn from N (0, 1). The parameter r controls the rank of A. The parameters b and r together control the signal to noise ratio in (2.2).
(2) Error matrix E: The entries of E are independently drawn from N (0, 1).
The design matrix X is generated with the following two settings to cover the low and high dimensional cases, respectively.
(3a) Design matrix X when n > p (low dimensional case):
are independently drawn from multivariate Normal distribution N p (0, Σ).
For i, j = 1, • • • , p, the (i, j)th entry of Σ is defined as Σ ij = η |i-j| for some η ∈ (0, 1).
(3b) Design matrix X when p > n > q (high dimensional case): Let X = Λ 0 Λ 1 Σ 1/2 , with Λ 0 ∈ IR n×q , and Λ 1 ∈ IR q×p . The entries of Λ 0 and Λ 1 are independently drawn from N (0, 1). The covariance matrix Σ is defined as in (3a).
For each case, we generate a testing sample {Y * , X * } of size n * independent of {Y, X} to assess the prediction performance of each method. To sum up, the parameters that control the data generating process are listed in Table 1.
In the first example, we examine the performance of OUR-CV, OUR-BIC, RSC, STRS, FES, OLS and LMD in the low dimensional case. We set n = 100, n * = 50, p = 25 and q = 20. Then, we vary the rank r = 5, 10, 15, the signal strength parameter b = 0.2, 0.4, and the correlation level η = 0.5, 0.9.
We simulate 200 replications for each scenario.
For each replication, we calculate two re-scaled Frobenius norms as The simulation results of the low dimensional example with b = 0.2 and 0.4 are presented in Tables 2 and3, respectively. In most scenarios, OUR-CV performs slightly better than OUR-BIC but pays a price on the computational cost. When the rank is small (e.g. r = 5), OUR-CV, OUR-BIC and STRS can recovery the correct rank with a rate close to 1 and have small estimation and prediction errors. RSC struggles when η is large and resultes in low correct recovery rates. LMD suffers when the signal to noise ratio is small. When the rank is moderate or large (e.g. r = 10 or 15), the correct rank recovery rates of RSC, STRS and LMD drop. As a result, the estimation and prediction errors of RSC, STRS and LMD are also inflated. Compared with the other methods, OUR-CV and OUR-BIC are less sensitive to rank, correlation level and the signal to noise ratio. In general, FES performs similar as RSC in terms of estimation and prediction, and OLS performs unsatisfactory as it ignores the low rank structure in A.
In the high dimensional example, we examine the performance of OUR-CV and OUR-BIC and compare them with RSC and STRS. We set n = 40, n * = 40, p = 100 and q = 25. Then, we vary the rank r = 10, 20, the signal strength parameter b = 0.2, 0.4, and the correlation level η = 0.5, 0.9. We simulate 200 replications for each scenario. The estimation and prediction errors are still measured by the sample mean and sample standard deviation of the rescaled Frobenius norms defined in (5.1). The estimation accuracy of rank r is measured by the correct rank recovery rate defined in (5.2).
The simulation results of the high dimensional example with b = 0.2 and 0.4 are presented in Tables 4 and5, respectively. Similar to the low dimensional case, OUR-CV outperforms OUR-BIC in most scenarios. When both η and r are large, OUR-CV and OUR-BIC maintains reasonable correct rank recovery rates while RSC and STRS fail to recovery the correct rank in most replications. The columns ∆ and Γ report the sample mean and sample standard deviation (in parentheses) of ∆ k and Γ k , which are defined in (5.1), over 200 replications. The column R reports the rank recovery rate, which is defined in (5.2), over 200 replications.
Particulate matter up to 2.5µm (PM2.5) is a complex mixture of solid particles, chemicals (e.g. sulfates, nitrates) and liquid droplets in the air, which include inhalable particles that are small enough to penetrate the thoracic region of the respiratory system. Short term (days) exposure to inhalable PM2.5 can cause an increase in hospital admissions related to respiratory and cardiovascular morbidity, such as aggravation of asthma, respiratory symptoms, and cardiovascular disorders. Long term (years) exposure to inhalable PM2.5 may The columns ∆ and Γ report the sample mean and sample standard deviation (in parentheses) of ∆ k and Γ k , which are defined in (5.1), over 200 replications. The column R reports the rank recovery rate, which is defined in (5.2), over 200 replications.
lead to an increase in mortality from cardiovascular and respiratory diseases, like lung cancer. The hazardous effects of inhalable PM2.5 on human health have been well-documented, see Riediker et al. (2004), Polichetti et al. (2009), Franck et al. (2011), Xing et al. (2016), Pun et al. (2017) and references therein.
In this section, we investigate the relationship between concentration of PM2.5 and four air pollutants: ozone, sulfur dioxide (SO2), carbon monoxide (CO), and nitrogen dioxide (NO2). The dataset for us to study is available at https://www.epa.gov/outdoor-air-quality-data, The columns ∆ and Γ report the sample mean and sample standard deviation (in parentheses) of ∆ k and Γ k , which are defined in (5.1), over 200 replications. The column R reports the rank recovery rate, which is defined in (5.2), over 200 replications.
it was collected from 37 outdoor monitoring sites across the United States.
Specifically, the concentration of each of the 4 pollutants was measured and collected daily from the 37 sites between January 2017 to April 2019, and it has 729 observations in total. The concentration of PM2.5 was collected in the same manner.
What we are interested in is the association between the concentrations of PM2.5 and the concentrations of the four air pollutants at the 37 monitor sites.
As the concentrations of the four air pollutants at one site may also contribute the concentrations of PM2.5 at other sites, we include the concentrations of the four air pollutants at all 37 sites in the explanatory variables for the concentration of PM2.5 at each site of the 37 sites, this gives us 148 explanatory variables for the concentration of PM2.5 at each site. In Figure 1, we plot the sample means of the concentrations of PM2.5 and of the four air pollutants against the geological locations where they were collected.
We take the first-order difference for each column of the dataset to remove the non-stationarity, and standardize it to make it have mean 0 and variance
) ∈ IR 728×37 be the matrix of the 728 observations of the response variable which contains the concentrations of PM2.5 collected by the 37 sites, and X = (X 1 , • • • , X 148 ) ∈ IR 728×148 be matrix of the 728 observations of predictor which contains the 148 explanatory variables. We apply the multivariate linear regression model (2.2) to fit the dataset, where E ∈ IR 728×37 is the matrix of random errors, and A ∈ IR 148×37 is the coefficient matrix of interest.
We compare the prediction performance of our method with rank estimated by the 5-fold cross-validation (OUR-CV), the Self-Tuning Rank Selection (STRS) proposed in Bing and Wegkamp (2019) where q = 37 and n t = 128, which are the number of columns and the number of rows of Y test , respectively. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.030 0.035 0.040 0.045 OZONE q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 2 3 4 SO2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.2 0.4 0.6 0.8 CO q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 10 20 30 NO2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 15 20 25
Figure 1: Sample means of the concentrations of PM2.5 and of the four air pollutants
We report in Table 6 the prediction error as well as the estimated rank of the coefficient matrix for each method concerned. According to Table 6, OUR-CV achieves the smallest prediction error among the three competitors.
Besides, the rank estimated by OUR-CV is 3 which is more parsimonious than the one estimated by STRS. To justify the rank estimation results, we apply eigen-decomposition to the coefficient matrix estimated by the OLS method, and draw the scree plot with the top 20 eigenvalues in Figure 2. The scree plot shows a clear elbow shape at the third eigenvalue. always dominate the second term with arbitrarily large probability. Consider the third term on the right side of (A.4). By Condition 1 again, we get
and accordingly H -H = O(α n C). By (A.3), we have that
holds uniformly in D, which is of smaller order of ∆ 1 .
For ∆ 3 , observe that which implies that with probability tending to one the case of k > r would not happen as long as k/h n → 0.