Comment on ‘Review of sparse sufficient dimension reduction’

We congratulate the authors on a very interesting overview of sparse sufficient dimension reduction (SDR). Sparse SDR methods are discussed in both the classical n>p setting as well as the high-dimensional p>n setting. Related topics such as model-free variable selection and variable screening are also discussed in a most logical fashion. Last but not least, two new methodological contributions are made in this review paper. Namely, new variable screening methods are proposed as extensions of Yu et al. (2016), and novel sparse SDR methods are discussed following the sparse SIR in Tan et al. (2020).

While all the methods discussed in this review are in the frequentist domain, our comment will focus on sparse SDR through Bayesian methods. Reich et al. (2011) proposed an SDR approach via Bayesian mixture modelling. Take single-index model as an example. Let $\{({\mathit{x}}_i,Y_i),i=1,\dots ,n\}$ be an i.i.d. sample from $(\mathit{x},Y)$. Let $\mathit{\beta }\in {\mathbb{R}}^p$ be a basis for the central subspace and let ${\lambda }_i={\mathit{\beta }}^{\mathrm{T}}{\mathit{x}}_i$ be the sufficient predictor. Then the conditional distribution of $Y_i$ given ${\mathit{x}}_i$ can be modelled as (1) $p(Y_i|{\lambda }_i)=\sum _{k=1}^K{p}_k({\lambda }_i)\mathcal{N}({\mu }_k,{\sigma }_y^2),$(1) where there are K normal mixture components and $p_k({\lambda }_i)$ denotes the weight of the kth component. By choosing the weights carefully, model (1(1) $p(Y_i|{\lambda }_i)=\sum _{k=1}^K{p}_k({\lambda }_i)\mathcal{N}({\mu }_k,{\sigma }_y^2),$(1) ) can be expressed as $\{\begin{array}{l}{Y}_i\sim \mathcal{N}({\mu }_{g_i},{\sigma }_y^2)\\ g_i=k& \text{if }{\psi }_k<Z_i<{\psi }_{k+1}\\ Z_i\sim \mathcal{N}({\lambda }_i,{\sigma }_z^2).\end{array}$ Here $Z_i$ is a latent continuous variable, $-\mathrm{\infty }={\psi }_1<{\psi }_2<\cdots <{\psi }_{K+1}=\mathrm{\infty }$ are cutpoints. By placing a prior on $\mathit{\beta }$ and the cutpoints, one can compute the conditional distributions and carry out the full Bayesian analysis. For SDR without sparsity, the prior for $\mathit{\beta }=({\beta }_1,\dots ,{\beta }_p{)}^{\mathrm{T}}$ is set as ${\beta }_j\sim \mathcal{N}(0,1)$, $j=1,\dots ,p$. To introduce sparsity, a two-component mixture prior is assumed as $\{\begin{array}{l}{\beta }_j\sim \mathcal{N}(0,{\pi }_j+c^2(1-{\pi }_j))\\ {\pi }_j\sim \mathrm{Ber}(\overline{\pi }),\end{array}$ where 0<c<1 is a fixed constant and $\overline{\pi }$ is the prior inclusion probability. If ${\pi }_j=1$, then the jth predictor is included in the model. Otherwise the jth predictor is removed from the model. Reich et al. (2011) also discussed a similar Bayesian mixture model for the multiple-index model, and the details are omitted.

Motivated by a frequentist SDR method recently proposed by Fang and Yu (2020), Power and Dong (2020) proposed a new Bayesian approach for sparse sufficient dimension reduction. Let $\mathrm{Var}(\mathit{x})=\mathbf{\Sigma }$, $\mathrm{E}(\mathit{x})=\mathit{\mu }$, and denote $\{J_1,\dots ,J_H\}$ as a partition for the support of Y. The classical sliced inverse regression (SIR) (Li, 1991) uses the kernel matrix ${\mathit{M}}_{\mathrm{SIR}}=\sum _{h=1}^H{\mathit{\xi }}_h{\mathit{\xi }}_h^{\mathrm{T}}/p_h$, where ${\mathit{\xi }}_h={\mathbf{\Sigma }}^{-1}\mathrm{E}\{(\mathit{x}-\mathit{\mu }){\delta }_h\}$ with ${\delta }_h=I(Y\in J_h)$ and $p_h=E({\delta }_h)$ for $h=1,\dots ,H$. Note that ${\mathit{\xi }}_h$ can be solved as an optimisation problem (2) ${\mathit{\xi }}_h=\underset{\mathit{\gamma }\in {\mathbb{R}}^p}{\mathrm{argmin}}\mathrm{E}[\{{\delta }_h-p_h-{\mathit{\gamma }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu }){\}}^2].$(2) Fang and Yu (2020) then applied the Mallows model averaging (MMA) of Hansen (2007) to solve the least squares problem (2(2) ${\mathit{\xi }}_h=\underset{\mathit{\gamma }\in {\mathbb{R}}^p}{\mathrm{argmin}}\mathrm{E}[\{{\delta }_h-p_h-{\mathit{\gamma }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu }){\}}^2].$(2) ). Power and Dong (2020) utilised Bayesian model averaging (BMA) (Raftery et al., 1997) to solve (2(2) ${\mathit{\xi }}_h=\underset{\mathit{\gamma }\in {\mathbb{R}}^p}{\mathrm{argmin}}\mathrm{E}[\{{\delta }_h-p_h-{\mathit{\gamma }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu }){\}}^2].$(2) ) instead. Similar to MMA, BMA works well for sparse models and may also adapt to models with dense signals. Furthermore, instead of solving for ${\mathit{\xi }}_h$, $h=1,\dots ,H$, individually, we may solve for them jointly. Let $\mathit{W}=({\mathit{\xi }}_1,\dots ,{\mathit{\xi }}_H)$ and $\mathit{U}=({\delta }_1-p_1,\dots ,{\delta }_H-p_H{)}^{\mathrm{T}}$. Then we have (3) $\begin{array}{r}\mathit{W}=\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times H}}{\mathrm{argmin}}\mathrm{E}[\{\mathit{U}-{\mathbf{\Theta }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu }){\}}^{\mathrm{T}}\{\mathit{U}-{\mathbf{\Theta }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu })\}].\end{array}$(3) To the best of our knowledge, there is no frequentist model averaging approach to solve (3(3) $\begin{array}{r}\mathit{W}=\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times H}}{\mathrm{argmin}}\mathrm{E}[\{\mathit{U}-{\mathbf{\Theta }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu }){\}}^{\mathrm{T}}\{\mathit{U}-{\mathbf{\Theta }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu })\}].\end{array}$(3) ). On the other hand, multi-response BMA (Brown et al., 1998) can be easily adapted to solve (3(3) $\begin{array}{r}\mathit{W}=\underset{\mathbf{\Theta }\in {\mathbb{R}}^{p\times H}}{\mathrm{argmin}}\mathrm{E}[\{\mathit{U}-{\mathbf{\Theta }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu }){\}}^{\mathrm{T}}\{\mathit{U}-{\mathbf{\Theta }}^{\mathrm{T}}(\mathit{x}-\mathit{\mu })\}].\end{array}$(3) ). As shown in Power and Dong (2020), the multi-response BMA outperforms the frequentist MMA for SDR.

We congratulate the authors again for providing a stimulating review of existing sparse SDR techniques, which should motivate further development of new SDR methods. It is our belief that more Bayesian approaches may be applied for this cause.