Indirect multivariate response linear regression

We assume that the measured predictor and response pairs $(x_{1},y_{1}),\dots,(x_{n},y_{n})$ are a realization of $n$ independent copies of $(X,Y)$, where $(X^{\prime},Y^{\prime})^{\prime}\sim N_{p+q}(\mu_{*},\Sigma_{*})$. We also assume that $\Sigma_{*}$ positive definite. Define the marginal parameters through the following partitions:

Our goal is to estimate the multivariate regression coefficient matrix $\beta_{*}=\Sigma_{*XX}^{-1}\Sigma_{*XY}$ in the forward regression model

without assuming that $\beta_{*}$ is sparse or that $\|\beta_{*}\|_{F}^{2}$ is small. To do this we will estimate the inverse regression’s coefficient matrix $\eta_{*}=\Sigma_{YY}^{-1}\Sigma_{XY}^{\prime}$ and the inverse regression’s error precision matrix $\Delta_{*}^{-1}$ in the inverse regression model

We connect the parameters of the inverse regression model to $\beta_{*}$ with the following proposition.

We prove Proposition 1 in Appendix A.1. This result leads us to propose a class of estimators of $\beta_{*}$ defined by

where $\hat{\eta}$, $\hat{\Delta}$, and $\hat{\Sigma}_{YY}$ are user-selected estimators of $\eta_{*}$, $\Delta_{*}$, and $\Sigma_{*YY}$. If $n>\max(p,q)$ and the ordinary sample estimators are used for $\hat{\eta}$, $\hat{\Delta}$ and $\hat{\Sigma}_{YY}$, then $\hat{\beta}$ is equivalent to $\hat{\beta}^{({\rm OLS})}$.

We propose to use shrinkage estimators of $\eta_{*}$, $\Delta_{*}^{-1}$, and $\Sigma_{*YY}^{-1}$ in (4). This gives us the potential to indirectly fit an unparsimonious forward regression model by fitting a parsimonious inverse regression model. For example, suppose that $\eta_{*}$ and $\Delta_{*}^{-1}$ are sparse, but $\beta_{*}$ is dense. To fit the inverse regression model, we could use any of the forward regression shrinkage estimators discussed in Section 1.

Lee and Liu, (2012) proposed an estimator of $\beta_{*}$ that also exploits the assumption that $(X^{\prime},Y^{\prime})^{\prime}$ is multivariate Normal; however, unlike our approach that makes no explicit assumptions about $\beta_{*}$, their approach assumes that both $\Sigma_{*}^{-1}$ and $\beta_{*}$ are sparse.

Modeling the inverse regression is a well-known idea in multivariate analysis. For example, when $Y$ is categorical, quadratic discriminant analysis models $(X|Y=y)$ as $p$-variate Normal. There are also many examples of modeling the inverse regression in the sufficient dimension reduction literature (Adragni and Cook,, 2009).

The most closely related work to ours is that by Cook et al., (2013). They proposed indirect estimators of $\beta_{*}$ based on modeling the inverse regression in the special case when the response is univariate, i.e. $q=1$. Under the same multivariate Normal assumption on $(X^{\prime},Y^{\prime})^{\prime}$ that we make, Cook et al., (2013) showed that

They proposed estimators of $\beta_{*}$ by replacing $\Sigma_{*XY}$ and $\Sigma_{*YY}$ in the right hand side of (5) with their usual sample estimators, and by replacing $\Delta_{*}^{-1}$ with a shrinkage estimator. This class of estimators was designed to exploit an abundant signal rate in the forward univariate response regression when $p>n$.

We now describe an estimator of the forward regression coefficient matrix $\beta_{*}$ defined by (4) that exploits zeros in the inverse regression’s coefficient matrix $\eta_{*}$, zeros in the inverse regression’s error precision matrix $\Delta_{*}^{-1}$, and zeros in the precision matrix of the responses $\Sigma_{*YY}^{-1}$. We estimate $\eta_{*}$ with

which separates into $p$ $L_{1}$-penalized least-squares regressions (Tibshirani,, 1996): the first predictor regressed on the response through the $p$th predictor regressed on the response. We select $\lambda_{j}$ with 5-fold cross-validation, minimizing squared prediction error totaled over the folds, in the regression of the $j$th predictor on the response $(j=1,\ldots,p)$. This allows us to estimate the columns of $\eta_{*}$ in parallel.

We estimate $\Delta_{*}^{-1}$ and $\Sigma_{*YY}^{-1}$ with $L_{1}$-penalized Normal likelihood precision matrix estimation (Yuan and Lin,, 2007; Banerjee et al.,, 2008). Let $\hat{\Sigma}_{\gamma,S}^{-1}$ be a generic version of this estimator with tuning parameter $\gamma$ and input $p$ by $p$ sample covariance matrix $S$:

where $\mathbb{S}_{+}^{p}$ is the set of symmetric and positive definite $p$ by $p$ matrices. There are many algorithms that solve (7). Two good choices are the graphical lasso algorithm (Yuan,, 2008; Friedman et al.,, 2008) and the QUIC algorithm (Hsieh et al.,, 2011). We select $\gamma$ with 5-fold cross-validation maximizing a validation likelihood criterion (Huang et al.,, 2006):

where $\mathcal{G}$ is a user-selected finite subset of the non-negative real line, $S_{(-k)}$ is the sample covariance matrix from the observations outside the $k$th fold, and $S_{(k)}$ is the sample covariance matrix from the observations in the $k$th fold centered by the sample mean of the observations outside the $k$th fold. We estimate $\Delta_{*}^{-1}$ using (7) with its tuning parameter selected by (8) and $S=(\mathbb{X}-\mathbb{Y}\hat{\eta}^{{\rm L1}})^{\prime}(\mathbb{X}-\mathbb{Y}\hat{\eta}^{{\rm L1}})/n$. Similarly, we estimate $\Sigma_{*YY}^{-1}$ using (7) with its tuning parameter selected by (8) and $S=\mathbb{Y}^{\prime}\mathbb{Y}/n$.

We propose indirect estimators of $\beta_{*}$ that exploit the assumption that the inverse regression’s coefficient matrix $\eta_{*}$ is rank deficient. We have the following simple proposition that links rank deficiency in $\eta_{*}$ and its estimator to $\beta_{*}$ and its indirect estimator.

We propose the following two example reduced rank indirect estimators of $\beta_{*}$:

1. 1.

   Estimate $\Sigma_{*YY}$ with $\mathbb{Y}^{\prime}\mathbb{Y}/n$ and estimate $(\eta_{*},\Delta_{*}^{-1})$ with Normal likelihood reduced rank inverse regression:

   	$\displaystyle(\hat{\eta}^{(r)},\hat{\Delta}^{-1(r)})=\operatorname*{arg\ min}_{(\eta,\Omega)\in\mathbb{R}^{q\times p}\times\mathbb{S}_{+}^{p}}$	$\displaystyle\left[n^{-1}{\rm tr}\left\{(\mathbb{X}-\mathbb{Y}\eta)^{\prime}(\mathbb{X}-\mathbb{Y}\eta)\Omega\right\}-\log{\rm det}(\Omega)\right]$		(9)
   		$\displaystyle{\rm subject}\ {\rm to}\ {\rm rank}(\eta)=r,$

   where $r$ is selected from $\{0,\ldots,\min(p,q)\}$. The solution to the optimization in (9) is available in closed form (Reinsel and Velu,, 1998).

2. 2.

   Estimate $\eta_{*}$ with $\hat{\eta}^{(r)}$ defined in (9), estimate $\Sigma_{*YY}^{-1}$ with (7) using $S=\mathbb{Y}^{\prime}\mathbb{Y}/n$, and estimate $\Delta_{*}^{-1}$ with (7) using $S=(\mathbb{X}-\mathbb{Y}\hat{\eta}^{(r)})^{\prime}(\mathbb{X}-\mathbb{Y}\hat{\eta}^{(r)})/n$.

Both example indirect reduced rank estimators of $\beta_{*}$ are formed by plugging in the estimators of $\eta_{*},\Delta_{*}^{-1}$, and $\Sigma_{*YY}$ to (4). The first estimator is likelihood-based and the second estimator exploits sparsity in $\Sigma_{*YY}^{-1}$ and $\Delta_{*}^{-1}$. Neither estimator is defined when $\min(p,q)>n$. In this case, which we do not address, a regularized reduced rank estimator of $\eta_{*}$ could be used instead of the estimator defined in (9), e.g. the factor estimation and selection estimator (Yuan et al.,, 2007) or the reduced rank ridge regression estimator (Mukherjee and Zhu,, 2011).

We compared the following indirect estimators of $\beta_{*}$ when the inverse regression’s coefficient matrix $\eta_{*}$ is sparse:

• $I_{L1}$.

  This is the indirect estimator proposed in Section 4.1.

• $I_{S}$.

  This is an indirect estimator defined by (4) with $\hat{\eta}$ defined by (6), $\hat{\Sigma}_{YY}=\mathbb{Y}^{\prime}\mathbb{Y}/n$, and $\hat{\Delta}=(\mathbb{X}-\mathbb{Y}\hat{\eta}^{L1})^{\prime}(\mathbb{X}-\mathbb{Y}\hat{\eta}^{L1})/n$.

• $O_{\Delta}$.

  This is a part oracle indirect estimator defined by (4) with $\hat{\eta}$ defined by (6), $\hat{\Sigma}_{YY}^{-1}$ defined by (7), and $\hat{\Delta}^{-1}=\Delta^{-1}_{*}$.

• $O$.

  This is a part oracle indirect estimator defined by (4) with $\hat{\eta}$ defined by (6), $\hat{\Sigma}_{YY}^{-1}=\Sigma_{*YY}^{-1}$, and $\hat{\Delta}^{-1}=\Delta^{-1}_{*}$.

• $O_{Y}$.

  This is a part oracle indirect estimator defined by (4) with $\hat{\eta}$ defined by (6), $\hat{\Sigma}_{YY}^{-1}=\Sigma_{*YY}^{-1}$, and $\hat{\Delta}^{-1}$ defined by (7).

We also included the following forward regression estimators of $\beta_{*}$:

• OLS/MP.

  This is the ordinary least squares estimator defined by $\operatorname*{arg\ min}_{\beta\in\mathbb{R}^{p\times q}}\|\mathbb{Y}-\mathbb{X}\beta\|_{F}^{2}$. When $n\leq p$, we use the solution $\mathbb{X}^{-}\mathbb{Y}$, where $\mathbb{X}^{-}$ is the Moore-Penrose generalized inverse of $\mathbb{X}$.

• R.

  This is the ridge penalized least squares estimator defined by

  $$\operatorname*{arg\ min}_{\beta\in\mathbb{R}^{p\times q}}\left(\|\mathbb{Y}-\mathbb{X}\beta\|_{F}^{2}+\lambda\|\beta\|_{F}^{2}\right).$$

• $\ell_{2}$.

  This is an alternative ridge penalized least squares estimator defined by

  $$\operatorname*{arg\ min}_{\beta\in\mathbb{R}^{p\times q}}\left(\|\mathbb{Y}-\mathbb{X}\beta\|_{F}^{2}+\sum_{m=1}^{q}\lambda_{m}\sum_{j=1}^{p}\beta_{jm}^{2}\right),$$

  where a separate tuning parameter is used for each response.

We selected the tuning parameters for uses of (6) with 5-fold cross-validation, minimizing validation prediction error on the inverse regression. Tuning parameters for $\ell_{2}$ and R were selected with 5-fold cross-validation, minimizing validation prediction error on the forward regression. We selected tuning parameters for uses of (7) with (8). The candidate set of tuning parameters was $\left\{10^{-8},10^{-7.5},\dots,10^{7.5},10^{8}\right\}$.

For 50 independent replications, we generated a realization of $n$ independent copies of $(X^{\prime},Y^{\prime})^{\prime}$, where $Y\sim N_{q}(0,\Sigma_{*YY})$ and $(X|Y=y)\sim N_{p}(\eta_{*}^{\prime}y,\Delta_{*})$. The $(i,j)$th entry of $\Sigma_{*YY}$ was set to $\rho_{Y}^{|i-j|}$ and the $(i,j)$th entry of $\Delta_{*}$ was set to $\rho_{\Delta}^{|i-j|}$. We set $\eta_{*}=Z\circ A$, where $\circ$ denotes the element-wise product: $Z$ had entries independently drawn from $N(0,1)$ and $A$ had entries independently drawn from the Bernoulli distribution with nonzero probability $s_{*}$. This model is ideal for $I_{L1}$ because $\Delta_{*}^{-1}$ and $\Sigma_{*YY}^{-1}$ are both sparse. Every entry in the corresponding randomly generated $\beta_{*}$ is nonzero with high probability, but the magnitudes of these entries are small. This motivated us to compare our indirect estimators of $\beta_{*}$ to the ridge-penalized least squares forward regression estimators R and $\ell_{2}$.

We evaluated performance with model error (Breiman and Friedman,, 1997; Yuan et al.,, 2007), which is defined by $\|\Sigma_{*XX}^{1/2}(\hat{\beta}-\beta_{*})\|_{F}^{2}$.

We report the average model errors, based on these 50 replications, in Table 1. When $s_{*}=0.1$, the indirect estimators defined by (4) performed well for all choices of $\rho_{Y}$ and $\rho_{\Delta}$. Our proposed estimator $I_{L1}$ was competitive with other indirect estimators also defined by (4), even those that used some oracle information. As $s_{*}$ increased with $\rho_{Y}=0.7$ and $\rho_{\Delta}=0.9$ fixed, the forward regression estimators performed nearly as well as $I_{L1}$.

Similarly, Table 2 shows that when $s_{*}=0.1$, $I_{L1}$ outperforms all three forward regression estimators. However, unlike in the lower dimensional setting illustrated in table 1, when $\eta_{*}$ is not sparse, i.e. $s_{*}\geq.3$, $I_{L1}$ is outperformed by forward regression approaches. The part oracle method $O_{Y}$ that used the knowledge of $\Sigma_{*YY}^{-1}$ outperformed the other two part oracle indirect estimators $O$ and $O_{\Delta}$ when $\rho_{\Delta}=.9$. Also, when $\rho_{\Delta}=.9$, $I_{L1}$ was competitive with the part oracle estimators. Taken together, the results in Tables 1 and 2 suggest that when $\eta_{*}$ is very sparse, our proposed indirect estimator $I_{L1}$ may perform nearly as well as the part oracle indirect estimators and the forward regression estimators.

We compared the performance of the following indirect reduced rank estimators of $\beta_{*}$:

• $I_{ML}^{(r)}$.

  This is the likelihood-based indirect example estimator 1 proposed in Section 4.2.

• $I^{(r)}$.

  This is the indirect example estimator 2 proposed in Section 4.2, which uses sparse estimators of $\Sigma_{*YY}^{-1}$ and $\Delta_{*}^{-1}$ in (4).

• $O^{(r)}$.

  This is a part oracle indirect estimator defined by (4) with $\hat{\eta}$ defined by (9), $\hat{\Delta}^{-1}=\Delta_{*}^{-1}$, and $\hat{\Sigma}_{YY}^{-1}=\Sigma_{*YY}^{-1}$.

• $O^{(r)}_{\Delta}$.

  This is a part oracle indirect estimator defined by (4) with $\hat{\eta}$ defined by (9), $\hat{\Delta}^{-1}$ defined by (7), and $\hat{\Sigma}_{YY}^{-1}=\Sigma_{*YY}^{-1}$.

• $O^{(r)}_{Y}$.

  This is a part oracle indirect estimator defined by (4) with $\hat{\eta}$ defined by (9), $\hat{\Delta}^{-1}=\Delta_{*}^{-1}$, $\hat{\Delta}^{-1}$ defined by (7), and $\hat{\Sigma}_{YY}^{-1}$ defined by (7).

We compared these indirect estimators to the following forward reduced rank regression estimator:

• RR.

  This is the likelihood based reduced rank regression (Izenman,, 1975; Reinsel and Velu,, 1998). The estimator of $\beta_{*}$ and the estimator of the forward regression’s error precision matrix $\Sigma_{*E}^{-1}$ are defined by

  	$\displaystyle(\hat{\beta}^{(r)},\hat{\Sigma}_{E}^{-1(r)})=\operatorname*{arg\ min}_{(\beta,\Omega)\in\mathbb{R}^{p\times q}\times\mathbb{S}_{+}^{q}}$	$\displaystyle\left[n^{-1}{\rm tr}\left\{(\mathbb{Y}-\mathbb{X}\beta)^{\prime}(\mathbb{Y}-\mathbb{X}\beta)\Omega\right\}-\log{\rm det}(\Omega)\right]$
  		$\displaystyle{\rm subject}\ {\rm to}\ {\rm rank}(\beta)=r.$

We selected the rank parameter $r$ for uses of (9) with 5-fold cross-validation, minimizing validation prediction error on the inverse regression. The rank parameter for RR was selected with 5-fold cross-validation, minimizing validation prediction error on the forward regression. We selected tuning parameters for uses of (7) with (8). The candidate set of tuning parameters was $\left\{10^{-8},10^{-7.5},\dots,10^{7.5},10^{8}\right\}$.

For 50 independent replications, we generated a realization of $n$ independent copies of $(X^{\prime},Y^{\prime})^{\prime}$ where $Y\sim N_{q}(0,\Sigma_{*YY})$ and $(X|Y=y)\sim N_{p}(\eta_{*}^{\prime}y,\Delta_{*})$. The $(i,j)$th entry of $\Sigma_{*YY}$ was set to $\rho_{Y}^{|i-j|}$ and the $(i,j)$th entry of $\Delta_{*}$ was set to $\rho_{\Delta}^{|i-j|}$. After specifying $r_{*}\leq\min(p,q)$, we set $\eta_{*}=PQ$, where $P\in\mathbb{R}^{q\times r_{*}}$ and $Q\in\mathbb{R}^{r_{*}\times p}$ had entries independently drawn from $N(0,1)$ so that $r_{*}=\text{rank}(\eta_{*})=\text{rank}(\beta_{*})$. As we did in the simulation in Section 5.1, we measured performance with model error.

We report the model errors, averaged over the 50 independent replications, in Table 3. Under every setting, $I^{(r)}$ outperformed all non-oracle competitors. When $r_{*}\leq 12$, $I^{(r)}$ outperformed both $O_{\Delta}^{(r)}$ and $O_{Y}^{(r)}$, which suggests that shrinkage estimation of $\Delta_{*}^{-1}$ and $\Sigma_{*YY}^{-1}$ was helpful. In each setting, $I_{ML}^{(r)}$ performed similarly to $RR$ even though they are estimating parameters of different condition distributions.

Our simulation studies in the previous sections used inverse regression data generating models. In this section, we compare the estimators from Section 5.2 using a forward regression data generating model.

For 50 independent replications, we generated a realization of $n$ independent copies of $(X^{\prime},Y^{\prime})^{\prime}$ where $X\sim N_{p}(0,\Sigma_{*XX})$ and $(Y|X=x)\sim N_{q}(\beta_{*}^{\prime}x,\Sigma_{*E})$. The $(i,j)$th entry of $\Sigma_{*XX}$ was set to $\rho_{X}^{|i-j|}$ and the $(i,j)$th entry of $\Sigma_{*E}$ was set to $\rho_{E}^{|i-j|}$. After specifying $r_{*}\leq\min(p,q)$, we set $\beta_{*}=ZQ$ where $Z\in\mathbb{R}^{p\times r_{*}}$ had entries independently drawn from $N(0,1)$ and $Q\in\mathbb{R}^{r_{*}\times q}$ had entries independently drawn from $\text{Uniform}(-1/4,1/4)$. In this data generating model, neither $\Delta_{*}^{-1}$ nor $\Sigma_{*YY}^{-1}$ had entries equal to zero.

The model errors, averaged over the 50 replications, are reported in Table 4. Both $I^{(r)}$ and $I_{ML}^{(r)}$ were competitive with RR in most settings. Although neither $\Delta_{*}^{-1}$ nor $\Sigma_{*YY}^{-1}$ were sparse, we again see that $I^{(r)}$ generally outperforms $O_{Y}^{(r)}$ and $O_{\Delta}^{(r)}$, both of which use some oracle information. These results indicate that shrinkage estimators of $\Delta_{*}^{-1}$ and $\Sigma_{*YY}^{-1}$ in (4) are helpful when neither is sparse.