# Sufficient reductions in regression with mixed predictors

Efstathia Bura (Corresponding author<EMAIL_ADDRESS>
Institute of Statistics and Mathematical Methods in Economics
Faculty of Mathematics and Geoinformation
TU Wien
Vienna, 1040, Austria Liliana Forzani<EMAIL_ADDRESS>
Facultad de Ingeniería Química
Universidad Nacional del Litoral
Researcher of CONICET
Santa Fe, Argentina Rodrigo García Arancibia<EMAIL_ADDRESS>
Instituto de Economía Aplicada Litoral-FCE-UNL
Universidad Nacional del Litoral
Researcher of CONICET
Santa Fe, Argentina Pamela Llop<EMAIL_ADDRESS>
Facultad de Ingeniería Química
Universidad Nacional del Litoral
Researcher of CONICET
Santa Fe, Argentina Diego Tomassi<EMAIL_ADDRESS>
Facultad de Ingeniería Química
Universidad Nacional del Litoral
Researcher of CONICET
Santa Fe, Argentina

###### Abstract

Most data sets comprise of measurements on continuous and categorical
variables. In regression and classification Statistics literature, modeling
high-dimensional mixed predictors has received limited attention. In this
paper we study the general regression problem of inferring on a variable of
interest based on high dimensional mixed continuous and binary predictors. The
aim is to find a lower dimensional function of the mixed predictor vector that
contains all the modeling information in the mixed predictors for the
response, which can be either continuous or categorical. The approach we
propose identifies sufficient reductions by reversing the regression and
modeling the mixed predictors conditional on the response. We derive the
maximum likelihood estimator of the sufficient reductions, asymptotic tests
for dimension, and a regularized estimator, which simultaneously achieves
variable (feature) selection and dimension reduction (feature extraction). We
study the performance of the proposed method and compare it with other
approaches through simulations and real data examples.

Keywords: High-dimensional, Multivariate Bernoulli, Regularization, Feature
selection, Feature extraction

##  1 Introduction

Most data sets comprise of measurements on a mixture of categorical and
continuous features. Examples abound in the biomedical and health sciences,
neuro-imaging, genomics, finance, social media, and internet advertising. The
first statistical approach to modeling the dependence structure of mixed data
we found in the literature is the location model of Olkin and Tate (1961). The
location model uses correlation as a measure of dependence and bypasses the
mixed nature of the data by grouping the continuous variables using the
categorical ones and requiring they be normally distributed with different
means but same variance within the groups.

More recently, Markov Networks, or undirected graphical models, that encode
pairwise conditional dependence relationships among random variables have been
used to model multivariate mixed data. With few exceptions (Yang et al.,
2014a, b, 2015; Chen et al., 2014), mixed continuous and categorical data are
modeled with the Gaussian Graphical Model (GGM) in a manner similar to the
location model. Binary variables are used to define the different categories
and GGM requires the continuous variables be conditionally normal and pairwise
conditionally independent within categories. References for GMMs for low-
dimensional mixed data include Lauritzen and Wermuth (1989), Lauritzen (1996),
Yuan and Lin (2007), Wainwright and Jordan (2008), and in the high-dimensional
setting, Cheng et al. (2014, 2017) and Lee and Hastie (2015). In particular,
Cheng et al. (2017) proposed a simplified version of the conditional Gaussian
distribution that reduces the number of parameters significantly while
maintaining flexibility.

Both GGMs and the location model are unsupervised approaches for mixed data
that do not include an output of interest. In the case of a categorical
output, approaches for the treatment of mixed, in particular, binary and
continuous input variables, include methods based on nonparametric density
estimation (Aitchison and Aitken, 1976), the use of logistic discrimination
(Day and Kerridge, 1967), in which the probability of group membership is
assumed to be a logistic function of the observed variates (Anderson, 1972,
1975), and a likelihood ratio classification rule (Krzanowski, 1975) based on
the location model of Olkin and Tate (1961). Krzanowski (1993) surveys and
summarizes the associated developments. More recently, the location model has
been used in multiple imputation [see, e.g., Javaras and van Dyk (2003),
Buuren (2018, Ch. 4, Sec. 4.4)].

In this paper we study the general regression and classification problem with
high-dimensional mixed predictors. Specifically, we consider the conditional
distribution of

$Y\mid{\mathbf{X}},\mathbf{H},$ (1)

where the response $Y$ is either continuous or categorical,
${\mathbf{X}}=(X_{1},X_{2},\dots,X_{p})^{T}$ is a vector of $p$ continuous,
and $\mathbf{H}=(H_{1},H_{2},\dots,H_{q})^{T}$ is a vector of $q$ binary
predictor variables. Our aim is to find a lower dimensional function of the
mixed predictor vector ${\mathbf{Z}}=({\mathbf{X}}^{T},\mathbf{H}^{T})^{T}$
that encapsulates all information the mixed predictors contain for the
response $Y$. Specifically, our target is the identification of a function,
other than the identity, $\mathbf{R}:{}^{p}\times{}^{q}\to{}^{d}$ such that
$F(Y\mid{\mathbf{Z}})=F(Y\mid\mathbf{R}({\mathbf{Z}}))$, where
$F(\cdot|\cdot)$ denotes the conditional cumulative distribution function of
the response given the predictors. Such a function $\mathbf{R}$ is called a
sufficient reduction of the regression of $Y$ on ${\mathbf{Z}}$.

This seemingly ambitious goal turns out to be surprisingly simple using the
inventive tool of inverse regression. When $Y$ and ${\mathbf{Z}}$ are both
random, inverse regression is based on the equivalence of the following two
statements [see Cook (2007), Bura et al. (2016), Bura and Forzani (2015)],

  * (i)

$Y\mid{\mathbf{Z}}\,{\buildrel d\over{=}}\,Y\mid\mathbf{R}({\mathbf{Z}})$

  * (ii)

${\mathbf{Z}}\mid(Y,\mathbf{R}({\mathbf{Z}}))\,{\buildrel
d\over{=}}\,{\mathbf{Z}}\mid\mathbf{R}({\mathbf{Z}})$

where ${\buildrel d\over{=}}$ signifies equal in distribution. Statement (i)
is an alternative definition of a sufficient reduction for the forward
regression in (1) and (ii) is the usual definition of a sufficient statistic
for a parameter $Y$ indexing the distribution of the mixed ${\mathbf{Z}}$. The
equivalence of (i) and (ii) obtains that if one considers $Y$ as a parameter,
the sufficient “statistic” for $Y$ is the sufficient reduction for the
regression of $Y$ on ${\mathbf{Z}}$. In consequence, in order to find a
sufficient reduction for the forward regression of $Y$ on ${\mathbf{Z}}$ in
(1), we can equivalently solve the inverse problem of finding a sufficient
statistic for the regression of ${\mathbf{Z}}$ on $Y$.

Our approach exploits the factorization

$F({\mathbf{X}},\mathbf{H}\mid Y)=F({\mathbf{X}}\mid
Y,\mathbf{H})F(\mathbf{H}\mid Y),$ (2)

by allowing us to treat the continuous and binary predictors separately, while
at the same time we account for their interdependence in their relationship
with $Y$ in Section 2. An advantageous aspect of (2) is that it requires fewer
parameters in order to characterize the distributional structure of the data.

In Section 3 we model $\mathbf{H}\mid Y$ as multivariate Bernoulli, and
${\mathbf{X}}\mid(Y,\mathbf{H})$ as multivariate normal, in analogy to the
Gaussian graphical model and the location model in unsupervised multivariate
analysis of mixed data. We show that the resulting distribution (2) belongs to
the exponential family, and derive sufficient reductions for the regression
$Y\mid({\mathbf{X}},\mathbf{H})$ from the two separate regressions,
${\mathbf{X}}\mid(Y,\mathbf{H})$ and $\mathbf{H}\mid Y$ in Section 3. We
compute the maximum likelihood estimator of the sufficient reduction in
Section 4, its asymptotic distribution in Section 4.3, and an asymptotic test
for the dimension of the sufficient reduction in Section 4.4. We complete our
treatment with a method for simultaneous sufficient dimension reduction and
variable selection in Section 5.

Section 6 contains an extensive simulation study that demonstrates the
competitive performance of our approach. Furthermore, we show the superior
performance of our methods as compared with generalized linear models and a
version of principal component regression that allows for mixed predictors in
the analysis of three data sets in Section 7.

Even though our focus in this paper is the regression of the usually
univariate $Y$ on the mixed ${\mathbf{Z}}$ vector, our development results in
a new multivariate regression method for a mixed continuous and binary
response, on which we comment further as we conclude in Section 8.

##  2 The Model

We start by specifying the notation we use throughout. The $\mathrm{vec}$
operator converts its matrix argument into a column vector. More precisely, if
$\mathbf{G}$ is an $m\times n$ matrix then $\mathrm{vec}(\mathbf{G})$ is an
$mn\times 1$ vector obtained by stacking the columns of $\mathbf{G}$. The
$\mathrm{unvec}$ operator is such that
$\mathrm{unvec}(\mathrm{vec}(\mathbf{G}))=\mathbf{G}$. We let $k_{q}=q(q-1)/2$
and $m_{p}=p(p+1)/2$. The $\mathrm{\,vech}$ operator converts the lower half
of a matrix including the main diagonal to a vector. That is, if $\mathbf{G}$
is a square $q\times q$ matrix then $\mathrm{\,vech}(\mathbf{G})$ is a
$m_{q}\times 1$ vector obtained by stacking the columns of the lower
triangular part of $\mathbf{G}$ including the diagonal. There is a unique
$\mathbf{D}_{q}\in{\mathbb{R}}^{q^{2}\times q(q+1)/2}$ and
$\mathbf{C}_{q}\in{\mathbb{R}}^{q(q+1)/2\times q^{2}}$ such that
$\mathrm{vec}(\mathbf{G})=\mathbf{D}_{q}\mathrm{\,vech}(\mathbf{G})$ and
$\mathrm{\,vech}(\mathbf{G})=\mathbf{C}_{q}\mathrm{vec}(\mathbf{G})$ for any
$\mathbf{G}$ symmetric $q\times q$ matrix.

The matrix $\mathbf{L}_{q}\in{\mathbb{R}}^{q\times q(q+1)/2}$ has entries 1
and 0, so that $\mathbf{L}_{q}\mathbf{C}_{q}$ is equal to $\mathbf{C}_{q}$ but
replacing the values 1/2 by zeros. The matrix
$\mathbf{J}_{q}\in{\mathbb{R}}^{k_{q}\times q(q+1)/2}$ has entries 1 and 0, so
that $\mathbf{J}_{q}\mathbf{C}_{q}$ is equal to $\mathbf{C}_{q}$ but replacing
the ones by zeros. A projection onto to the columns of $b$ is denoted ${\bf
P}_{b}$ and the projection onto the orthogonal complement of $b$ will be
denotes as ${\bf Q}_{b}$.

To regress $({\mathbf{X}},\mathbf{H})$ on $Y$, we model
${\mathbf{X}}\mid(Y,\mathbf{H})$ and $\mathbf{H}\mid Y$ separately and use the
factorization in (2).

###  2.1 The distribution of ${\mathbf{X}}\mid(\mathbf{H},Y)$

We let the $p$-dimensional vector of continuous random variables
${\mathbf{X}}|(\mathbf{H},Y)$ be multivariate normal with

${\mathbf{X}}\mid\mathbf{H},Y\sim
N\left(\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\mathbf{A}(\mathbf{f}_{Y}-\bar{\mathbf{f}}_{Y})+\bm{\beta}(\mathbf{H}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}),\bm{\Delta}\right),$
(3)

where
$\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}=\mathrm{E}_{{\mathbf{X}}}({\mathbf{X}})$,
$\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}=\mathrm{E}_{\scriptstyle{\mathbf{H}}}(\mathbf{H})$,
$\mathbf{f}_{Y}:\real\rightarrow{}^{r}$ is a known function of $Y$,
$\bar{\mathbf{f}}_{Y}=\mathrm{E}_{Y}(\mathbf{f}_{Y})$, $\mathbf{A}:p\times r$,
and $\bm{\beta}:p\times q$, are unconstrained parameter matrices, and
$\bm{\Delta}$ is a $p\times p$ positive definite covariance matrix. For
example, if the response is continuous, $\mathbf{f}_{Y}$ can be a vector of
polynomials of order $r$, or, in order to avoid multicollinearity, of a set of
$r$ orthonormal basis functions. If the response is categorical with values in
one of $h$ categories $C_{k}$, $k=1,\dots,h$, we set $r=h-1$ and let the
$k$-th element of $\mathbf{f}_{Y}$ to be $I(Y\in C_{k})$, where $I$ is the
indicator function. To simplify notation, henceforth $\mathbf{f}_{Y}$ will
signify the centered $\mathbf{f}_{Y}-\bar{\mathbf{f}}_{Y}$.

The probability density function of ${\mathbf{X}}\mid(\mathbf{H},Y)$ in model
(3) is

$\displaystyle f({\mathbf{X}}\mid\mathbf{H},Y=y)$
$\displaystyle=\frac{1}{\sqrt{2\pi}\sqrt{|\bm{\Delta}|}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}({\mathbf{X}}-\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}})-\mathbf{A}\mathbf{f}_{y}-\bm{\beta}(\mathbf{H}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}})\Big{)}^{T}$
$\displaystyle\hskip
113.81102pt\bm{\Delta}^{-1}\Big{(}({\mathbf{X}}-\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}})-\mathbf{A}\mathbf{f}_{y}-\bm{\beta}(\mathbf{H}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}})\Big{)}\bigg{\\}}.$
(4)

###  2.2 The distribution of $\mathbf{H}\mid Y$

The joint distribution of a random vector whose elements are binary random
variables is modelled with the multivariate Bernoulli distribution [see
Whittaker (2009); Dai (2012); Dai et al. (2013)]. Its probability mass
function involves terms representing third and higher order moments of the
random variables. The Ising model (Ising, 1925) is frequently used instead to
alleviate the complexity of modeling as it includes up to second order
interactions among the binary variables. For the multivariate binary
regression $\mathbf{H}|Y$ we use the Ising with covariates model introduced in
Cheng et al. (2014), where covariates are incorporated directly.

Let $\mathcal{H}=$ all possible combinations of $\mathbf{H}\in\\{0,1\\}^{q}$,
$\mathbf{H}_{-j}=(H_{1},\ldots,H_{j-1},H_{j+1},\ldots,H_{q})$,
$\mathbf{H}_{-i,-j}=(H_{1},\ldots,$
$H_{i-1},H_{i+1},\ldots$,$H_{j-1},H_{j+1},\ldots,H_{q}),$ $i,j=1,\ldots,q$.
The joint probability mass function of the $q$-dimensional vector of binary
variables $\mathbf{H}$ conditional on $Y$ is [see Cheng et al. (2014)]

$\displaystyle P(\mathbf{H}\mid Y=y)$
$\displaystyle=\frac{1}{G(\bm{\Gamma}_{y})}\exp\left\\{\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})\mathrm{\,vech}(\bm{\Gamma}_{y})\right\\},$
(5)

where
$G(\bm{\Gamma}_{y})=\sum_{\mathbf{H}\in\mathcal{H}}\exp\left(\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})\mathrm{\,vech}(\bm{\Gamma}_{y})\right),$
and $\bm{\Gamma}_{y}=(\gamma_{ij}^{y})$ is a $q\times q$ symmetric matrix with
elements

$\displaystyle\gamma^{y}_{jj}=\log\left(\frac{\Pr(H_{j}=1\mid\mathbf{H}_{-j}=0,y)}{1-\Pr(H_{j}=1|\mathbf{H}_{-j}=0,y)}\right),$
$\displaystyle\gamma^{y}_{ij}=\log\frac{\Pr(H_{i}=1,H_{j}=1\mid
0,y)\Pr(H_{i}=0,H_{j}=0\mid\mathbf{H}_{-i,-j}=0,y)}{\Pr(H_{i}=1,H_{j}=0\mid\mathbf{H}_{-i,-j}=0,y)\Pr(H_{i}=0,H_{j}=1\mid\mathbf{H}_{-i,-j}=0,y)},$

for $i\not=j$.

A linear model with independent variables $\mathbf{f}_{Y}\in{}^{r}$ is a
natural choice for each $\gamma^{y}_{ij}$,

$\gamma^{y}_{ij}=\tau^{*}_{ij,0}+\bm{\tau}_{ij}^{T}\mathbf{f}_{Y},\hskip
28.45274pti,j=1,\ldots,q,$ (6)

where $\bm{\tau}_{ij}^{T}=(\tau_{ij,1},\ldots,\tau_{ij,r})$ is a vector of
parameters independent of $Y$, and $\tau^{*}_{ij,0}$ is the intercept for each
$(i,j)$. Here again, $\mathbf{f}_{Y}$ is also centered, and can be different
from that in (2.1), even though, as will be seen later, choosing the same
$\mathbf{f}_{Y}$ simplifies the formula for the joint distribution in (8) as
well as the derivation of a sufficient reduction for the regression of $Y$ on
${\mathbf{X}},\mathbf{H}$.

Next we define the $q\times q$ matrices, $\bm{\tau}_{0}$ and $\bm{\tau}_{k}$,
$k=1,\ldots,r$, as $[\bm{\tau}^{*}_{0}]_{ij}=\tau^{*}_{ij,0}$ and
$[\bm{\tau}_{k}]_{ij}=\tau_{ij,k}$ with $i,j=1,\ldots,q$ and $k=1,\ldots,r$.
We let $\bm{\tau}_{0}=\mathrm{\,vech}({\bm{\tau}}^{*}_{0})$, a $q(q+1)/2$
vector, and
$\bm{\tau}=\left(\mathrm{\,vech}({\bm{\tau}}_{1}),\dots,\mathrm{\,vech}({\bm{\tau}}_{r})\right)$,
a $q(q+1)/2\times r$ matrix, so that the $q(q+1)/2$ vector
$\mathrm{\,vech}(\bm{\Gamma}_{y})$ is

$\mathrm{\,vech}(\bm{\Gamma}_{y})=\bm{\tau}_{0}+\bm{\tau}\mathbf{f}_{y}.$

Under (6) the probability mass function of $\mathbf{H}|Y$ in (5) is

$P(\mathbf{H}\mid
Y=y)=\frac{1}{G(\bm{\Gamma}_{y})}\exp\left\\{\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})(\bm{\tau}_{0}+\bm{\tau}\mathbf{f}_{y})\right\\},$
(7)

with
$G(\bm{\Gamma}_{y})=\sum_{\mathbf{H}\in\mathcal{H}}\exp\left(\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})(\bm{\tau}_{0}+\bm{\tau}\mathbf{f}_{y})\right).$

Under (7) and (2.1), the joint distribution of the inverse regression
(${\mathbf{X}},\mathbf{H}\mid Y$) has probability density function

$\displaystyle f({\mathbf{X}},\mathbf{H}\mid Y=y)$
$\displaystyle=f({\mathbf{X}}\mid\mathbf{H},Y=y)f(\mathbf{H}\mid Y=y)$
$\displaystyle=\frac{1}{\sqrt{2\pi}\sqrt{|\bm{\Delta}|}}\exp\bigg{\\{}-\frac{1}{2}\Big{(}({\mathbf{X}}-\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}})-\mathbf{A}\mathbf{f}_{y}-\bm{\beta}(\mathbf{H}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}})\Big{)}^{T}$
$\displaystyle\hskip
128.0374pt\bm{\Delta}^{-1}\Big{(}({\mathbf{X}}-\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}})-\mathbf{A}\mathbf{f}_{y}-\bm{\beta}(\mathbf{H}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}})\Big{)}\bigg{\\}}$
$\displaystyle\hskip
28.45274pt\times\frac{1}{G(\bm{\Gamma}_{y})}\exp\bigg{\\{}\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})\left(\bm{\tau}_{0}+\bm{\tau}\mathbf{f}_{y}\right)\bigg{\\}}.$
(8)

Our regression model for the mixed vector ${\mathbf{Z}}$ is similar to the
regression model of Fitzmaurice and Laird (1997) with the difference that we
do not allow $\bm{\mu}_{\scriptstyle{\mathbf{H}}}$ to vary with $Y$ in (2.1).
This results in different maximum likelihood estimates for the parameters in
(8) in Section 4.1.

##  3 Sufficient Reductions

We focus on the regression problem (1), where we aim to identify a reduction
$\mathbf{R}({\mathbf{Z}})$ such that $Y\mid{\mathbf{Z}}\,{\buildrel
d\over{=}}\,Y\mid\mathbf{R}({\mathbf{Z}})$. Since the latter is equivalent to
${\mathbf{Z}}\mid(Y,\mathbf{R}({\mathbf{Z}}))\,{\buildrel
d\over{=}}\,{\mathbf{Z}}\mid\mathbf{R}({\mathbf{Z}})$, as discussed in the
introduction, we will derive the sufficient reduction
$\mathbf{R}({\mathbf{Z}})$ using (2).

Of central importance to our development is showing that the density of
$({\mathbf{X}},\mathbf{H})\mid Y$ in (8) belongs to the exponential family of
distributions. In Appendix Derivation of Eqn. (9), we express (8) as

$f({\mathbf{X}},\mathbf{H}\mid
Y=y)=h({\mathbf{X}},\mathbf{H})\exp\left(\mathbf{T}^{T}({\mathbf{X}},\mathbf{H})\bm{\eta}_{y}-\psi(\bm{\eta}_{y})\right),$
(9)

which belongs to the natural exponential family of distributions [see, e.g.,
Morris (2006)]. In (9), $h({\mathbf{X}},\mathbf{H})=(2\pi)^{-1/2}$, the
sufficient statistic is

$\mathbf{T}({\mathbf{X}},\mathbf{H})=\left(\begin{array}[]{c}{\mathbf{X}}\\\
\mathbf{H}\\\
-\frac{1}{2}\mathbf{D}_{p}^{T}\mathbf{D}_{p}\mathrm{\,vech}({\mathbf{X}}{\mathbf{X}}^{T})\\\
\mathrm{vec}({\mathbf{X}}\mathbf{H}^{T})\\\
\mathbf{J}_{q}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathrm{\,vech}}(\mathbf{H}\mathbf{H}^{T})\end{array}\right),$
(10)

the natural parameters are

$\displaystyle\bm{\eta}_{y}$
$\displaystyle=\left(\begin{array}[]{c}\bm{\eta}_{y1}\\\ \bm{\eta}_{y2}\\\
\bm{\eta}_{y3}\\\ \bm{\eta}_{y4}\\\
\bm{\eta}_{y5}\end{array}\right)=\left(\begin{array}[]{cccccccc}\mathbf{I}_{p}&\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{p}&0&0&0&0&0&0\\\
0&0&\mathbf{I}_{q}&\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q}&0&0&0&0\\\
0&0&0&0&\mathbf{I}_{m_{p}}&0&0&0\\\ 0&0&0&0&0&\mathbf{I}_{pq}&0&0\\\
0&0&0&0&0&0&\mathbf{I}_{k_{q}}&\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{k_{q}}\end{array}\right)\left(\begin{array}[]{c}\bm{\vartheta}_{1}\\\
\bm{\vartheta}_{2}\\\ \bm{\vartheta}_{3}\\\ \bm{\vartheta}_{4}\\\
\bm{\vartheta}_{5}\\\ \end{array}\right)$ (26)
$\displaystyle\doteq\mathbf{F}_{y}\bm{\vartheta},$ (27)

with
$\bm{\vartheta}^{T}=(\bm{\vartheta}_{1}^{T},\bm{\vartheta}_{2}^{T},\bm{\vartheta}_{3}^{T},\bm{\vartheta}_{4}^{T},\bm{\vartheta}_{5}^{T})^{T}$
where

$\displaystyle\bm{\vartheta}_{1}$
$\displaystyle=\begin{pmatrix}\bm{\vartheta}_{1,0}\\\
\bm{\vartheta}_{1,1}\end{pmatrix}=\begin{pmatrix}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}-\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\mathbf{H}}\\\
\mathrm{vec}(\bm{\Delta}^{-1}\mathbf{A})\end{pmatrix}:(p+pr)\times 1,$
$\displaystyle\bm{\vartheta}_{2}$
$\displaystyle=\begin{pmatrix}\bm{\vartheta}_{2,0}\\\
\bm{\vartheta}_{2,1}\end{pmatrix}=\begin{pmatrix}-\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\mathbf{H}}+\mathbf{L}_{q}\bm{\tau}_{0}-\frac{1}{2}\mathbf{L}_{q}\mathbf{D}_{q}^{T}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})\\\
\mathrm{vec}(\mathbf{L}_{q}\bm{\tau}-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A})\end{pmatrix}:(q+qr)\times
1,$ $\displaystyle{\bm{\vartheta}}_{3}$
$\displaystyle=\bm{\vartheta}_{3,0}=\mathrm{\,vech}(\bm{\Delta}^{-1}):k_{p}\times
1,$ (28) $\displaystyle{\bm{\vartheta}}_{4}$
$\displaystyle=\bm{\vartheta}_{4,0}=\mathrm{vec}(\bm{\Delta}^{-1}\bm{\beta}):pq\times
1,$ $\displaystyle\bm{\vartheta}_{5}$
$\displaystyle=\begin{pmatrix}\bm{\vartheta}_{5,0}\\\
\bm{\vartheta}_{5,1}\end{pmatrix}=\begin{pmatrix}-\frac{1}{2}\mathbf{J}_{q}\mathbf{D}^{T}_{q}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})+\mathbf{J}_{q}\bm{\tau}_{0}\\\
\mathrm{vec}(\mathbf{J}_{q}\bm{\tau})\end{pmatrix},$

and

$\displaystyle\psi(\bm{\eta}_{y})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\log|\mathrm{unvec}(\mathbf{D}_{p}{\bm{\eta}}_{y3})|+\log(G(\bm{\Gamma}_{y}))+\frac{1}{2}{\bm{\eta}}_{y1}^{T}(\mathrm{unvec}(\mathbf{D}_{p}{\bm{\eta}}_{y3}))^{-1}{\bm{\eta}}_{y1}$
$\displaystyle\doteq$
$\displaystyle\psi_{1}(\bm{\eta}_{y})+\psi_{2}(\bm{\eta}_{y})+\psi_{3}(\bm{\eta}_{y}),$

with

$\displaystyle G(\bm{\Gamma}_{y})$ $\displaystyle=$
$\displaystyle\sum_{H}\exp\left[\left(\mathbf{J}_{q}C_{q}\mathrm{vec}(\mathbf{H}\mathbf{H}^{T})\right)^{T}\left({\bm{\eta}}_{y5}+\mathbf{J}_{q}\frac{1}{2}{\mathbf{D}}^{T}_{q}{\mathrm{vec}}(\bar{\bm{\eta}}_{y4}^{T}(\mathrm{unvec}(\mathbf{D}_{p}{\bm{\eta}}_{y3}))^{-1}\bar{\bm{\eta}}_{y4}))\right.\right)$
$\displaystyle\left.+\mathbf{H}^{T}\left({\bm{\eta}}_{y2}+\bar{\bm{\eta}}_{y4}^{T}(\mathrm{unvec}(\mathbf{D}_{p}{\bm{\eta}}_{3}))^{-1}{\bm{\eta}}_{y1}+\frac{1}{2}\mathbf{L}_{q}{\mathbf{D}}_{q}^{T}\mathrm{vec}(\bar{\bm{\eta}}_{y4}^{T}(\mathrm{unvec}(\mathbf{D}_{p}{\bm{\eta}}_{y3}))^{-1}\bar{\bm{\eta}}_{y4}))\right)\right],$

where $\bar{\bm{\eta}}_{y4}=\mathrm{unvec}({\bm{\eta}}_{y4})$.

For any matrix $\mathbf{V}$, let ${\mathcal{S}}_{\mathbf{V}}$ denote the span
of the columns of $\mathbf{V}$; that is,
${\mathcal{S}}_{\mathbf{V}}=\mathrm{span}(\mathbf{V})$. Theorem 1 obtains the
sufficient reduction for the regression of $Y$ on $({\mathbf{X}},\mathbf{H})$
using a result from Bura et al. (2016).

###### Theorem 1

Suppose that $({\mathbf{X}},\mathbf{H})\mid Y$ has density given by (9). The
minimal sufficient reduction for the regression
$Y\mid({\mathbf{X}},\mathbf{H})$ is

$\displaystyle\mathbf{R}({\mathbf{X}},\mathbf{H})$
$\displaystyle=\bm{\alpha}_{\mathbf{a}}^{T}\left(\mathbf{T}({\mathbf{X}},\mathbf{H})-\mathrm{E}(\mathbf{T}({\mathbf{X}},\mathbf{H}))\right),$
(31)

where $\mathbf{T}({\mathbf{X}},\mathbf{H})$ is given by (10) and
$\bm{\alpha}_{\mathbf{a}}$ is a basis for
${\mathcal{S}}_{\bm{\alpha}_{\mathbf{a}}}=\mathrm{span}\\{\bm{\eta}_{Y}-\mathrm{E}({\bm{\eta}_{Y}}),Y\in\mathcal{Y}\\}$,
with $\bm{\eta}_{Y}$ given in (26).

We provide the proof of Theorem 1 in Appendix Proof of Theorem 1, where we see
that the reduction in (31) is characterized by the coefficients of the basis
for
$\mathrm{span}\\{\bm{\eta}_{Y}-\mathrm{E}({\bm{\eta}_{Y}}),Y\in\mathcal{Y}\\}=\mathrm{span}(\mathbf{a})$
with

$\mathbf{a}=\left(\begin{array}[]{c}\bm{\Delta}^{-1}\mathbf{A}\\\
\mathbf{L}_{q}\bm{\tau}-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A}\\\ 0\\\ 0\\\
\mathbf{J}_{q}\bm{\tau}\end{array}\right)=\left(\begin{array}[]{c}\mathrm{unvec}(\bm{\vartheta}_{1,1})\\\
\mathrm{unvec}(\bm{\vartheta}_{2,1})\\\ 0\\\ 0\\\
\mathrm{unvec}(\bm{\vartheta}_{5,1})\end{array}\right).$

Since $\bm{\eta}_{3}\doteq\bm{\eta}_{y3}$ and
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\bm{\eta}_{4}}\doteq\bm{\eta}_{y4}$
do not depend on $y$, Corollary 2 follows.

###### Corollary 2

Suppose the density of $({\mathbf{X}},\mathbf{H})|Y$ is given by (9). A
minimal sufficient dimension reduction for the regression of $Y$ on
$({\mathbf{X}},\mathbf{H})$ is given by

$\displaystyle\mathbf{R}({\mathbf{X}},\mathbf{H})=\bm{\alpha}_{\mathbf{b}}^{T}\left({\mathbf{t}}\left({\mathbf{X}},\mathbf{H}\right)-\mathrm{E}\left({\mathbf{t}}({\mathbf{X}},\mathbf{H})\right)\right),$
(32)

where

$\displaystyle{\mathbf{t}}({\mathbf{X}},\mathbf{H})$ $\displaystyle=$
$\displaystyle\left({\mathbf{X}}^{T},\mathbf{H}^{T},\left[\mathbf{J}_{q}\mathrm{\,vech}(\mathbf{H}\mathbf{H}^{T})\right]^{T}\right)^{T},$
(33)

and $\bm{\alpha}_{\mathbf{b}}$ is a basis for
$\mathcal{S}_{\bm{\alpha}_{\tiny\mathbf{b}}}=\mathrm{span}\\{\mathbf{b}\\}$
with

$\displaystyle\mathbf{b}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}\bm{\Delta}^{-1}\mathbf{A}\\\
\mathbf{L}_{q}\bm{\tau}-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A}\\\
\mathbf{J}_{q}\bm{\tau}\end{array}\right)=\left(\begin{array}[]{c}\mathrm{unvec}(\bm{\vartheta}_{1,1})\\\
\mathrm{unvec}(\bm{\vartheta}_{2,1})\\\
\mathrm{unvec}(\bm{\vartheta}_{5,1})\end{array}\right).$ (40)

As the reduction in (32) is not only sufficient but also minimal, we call it
optimal SDR in the sequel.

###### Corollary 3

When the predictor vector contains only continuous variables; that is, $q=0$
and ${\mathbf{Z}}={\mathbf{X}}$, the sufficient dimension reduction is

$\displaystyle\mathbf{R}({\mathbf{X}})=\bm{\alpha}_{1}^{T}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left({\mathbf{X}}-\mathrm{E}({\mathbf{X}})\right)},$
(41)

where
${\mathcal{S}}_{\bm{\alpha}_{1}}=\mathrm{span}(\bm{\alpha}_{1})=\mathrm{span}(\bm{\Delta}^{-1}\mathbf{A})$,
and $\mathbf{A}:p\times r$ in (3).

The reduction (41) coincides with Principal Fitted Components (PFC) in Cook
and Forzani (2008).

###### Corollary 4

When the predictor vector contains only binary variables; that is, $p=0$ and
${\mathbf{Z}}=\mathbf{H}$, the sufficient dimension reduction is

$\displaystyle\mathbf{R}({\mathbf{X}})$
$\displaystyle=\bm{\alpha}_{2}^{T}({\mathbf{s}}(\mathbf{H})-\mathrm{E}({\mathbf{s}}(\mathbf{H}))),$
(42)

where

${\mathbf{s}}(\mathbf{H})=\left(\mathbf{H}^{T},\left[\mathbf{J}_{q}\mathrm{\,vech}(\mathbf{H}\mathbf{H}^{T})\right]^{T}\right)^{T},$
(43)

and

$\displaystyle{\mathcal{S}}_{\bm{\scriptstyle\alpha}_{2}}$
$\displaystyle=\mathrm{span}(\bm{\alpha}_{2})=\mathrm{span}\left(\begin{array}[]{c}\mathbf{L}_{q}\bm{\tau}\\\
\mathbf{J}_{q}\bm{\tau}\end{array}\right).$ (46)

When the predictors are mixed, we derive a sufficient but not minimal
reduction in Corollary (5), which we call sub-optimal SDR.

###### Corollary 5

Suppose that $({\mathbf{X}},\mathbf{H})|Y$ has density (9). A sufficient
reduction for the regression of $Y$ on $({\mathbf{X}},\mathbf{H})$ is given by

$\displaystyle\mathbf{R}({\mathbf{X}},\mathbf{H})$
$\displaystyle=\bm{\alpha}_{\mathbf{c}}^{T}\left(\mathbf{w}({\mathbf{X}},\mathbf{H})-\mathrm{E}(\mathbf{w}({\mathbf{X}},\mathbf{H}))\right),$
(47)

with

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbf{w}({\mathbf{X}},\mathbf{H})=({\mathbf{X}}^{T},\mathbf{H}^{T},[\mathrm{\,vech}(\mathbf{H}\mathbf{H}^{T})]^{T})^{T}},$
(48)

$\displaystyle\mathcal{S}_{\bm{\alpha}_{\mathbf{c}}}$
$\displaystyle=\mathrm{span}(\bm{\alpha}_{\mathbf{c}})=\mathrm{span}\left(\begin{array}[]{cc}\mathbf{c}_{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{c}_{2}\end{array}\right)=\mathrm{span}\left(\begin{array}[]{cc}\bm{\Delta}^{-1}\mathbf{A}&0\\\
-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A}&0\\\ 0&\bm{\tau}\end{array}\right).$
(54)

If $\mathrm{rank}(\mathbf{c}_{1})=d_{1}\leq\min\\{r,p\\}$ and
$\mathrm{rank}(\mathbf{c}_{2})=d_{2}\leq\min\\{r,q(q+1)/2\\}$, then

$\displaystyle{\mathbf{c}}_{1}=\left(\begin{array}[]{c}{\bm{\Delta}}^{-1}{\mathbf{A}}\\\
-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A}\end{array}\right)=\left(\begin{array}[]{c}{\bm{\alpha}\bm{\xi}}\\\
-\bm{\beta}^{T}\bm{\alpha}\bm{\xi}\end{array}\right),\;\mathbf{c}_{2}=\bm{\kappa}\bm{\iota},$
(59)

where $\mathbf{A}=\bm{\alpha}\bm{\xi}$, $\bm{\alpha}\in{}^{p\times d_{1}}$,
$\bm{\xi}\in{}^{d_{1}\times r}$, $\bm{\kappa}\in{}^{q(q+1)/2\times d_{2}}$ and
$\bm{\iota}\in{}^{d_{2}\times r}$ are full rank matrices. Therefore,

$\displaystyle\mathrm{span}(\mathbf{c}_{1})$
$\displaystyle=\mathrm{span}\\{(\bm{\alpha}^{T},-\bm{\alpha}^{T}\bm{\beta})^{T}\\},$
(60) $\displaystyle\mathrm{span}(\mathbf{c}_{2})$
$\displaystyle=\mathrm{span}(\bm{\kappa}).$ (61)

In Table 1, we summarize the results of this Section and tabulate the
sufficient reductions for mixed normal and binary predictors.

Table 1: Sufficient Reductions in Regressions with Mixed Predictors.

Sufficient Reductions Predictor Distribution Optimal SDR Sub-optimal
$({\mathbf{X}},\mathbf{H})\mid Y$ with density (9)
$\bm{\alpha}_{\mathbf{b}}^{T}({\mathbf{t}}({\mathbf{X}},\mathbf{H})-\mathrm{E}({\mathbf{t}}({\mathbf{X}},\mathbf{H})))$
$\bm{\alpha}_{\mathbf{c}}^{T}({\mathbf{w}}({\mathbf{X}},\mathbf{H})-\mathrm{E}({\mathbf{w}}({\mathbf{X}},\mathbf{H})))$
${\mathbf{t}}({\mathbf{X}},\mathbf{H})$ in (33),
${\mathbf{w}}({\mathbf{X}},\mathbf{H})$ in (48),
$\mathcal{S}_{\bm{\alpha}_{\tiny\mathbf{b}}}=\mathrm{span}\\{\mathbf{b}\\}$,
$\mathbf{b}$ in (40) $\mathcal{S}_{\bm{\alpha}_{\mathbf{c}}}$ in (54), (60),
and (61) ${\mathbf{X}}\mid Y\sim
N\left(\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\mathbf{A}(\mathbf{f}_{Y}-\bar{\mathbf{f}}_{Y}),\bm{\Delta}\right)$
$\bm{\alpha}_{1}^{T}({\mathbf{X}}-\mathrm{E}({\mathbf{X}}))$
${\mathcal{S}}_{\bm{\alpha}_{1}}=\mathrm{span}(\bm{\Delta}^{-1}\mathbf{A})$
$\mathbf{H}\mid Y$ with mass function (7)
$\bm{\alpha}_{2}^{T}({\mathbf{s}}(\mathbf{H})-\mathrm{E}({\mathbf{s}}(\mathbf{H})))$
${\mathbf{s}}(\mathbf{H})$ in (43) ${\mathcal{S}}_{\bm{\alpha}_{2}}$ in (46)

##  4 Reduction Estimators and their Asymptotic Distribution

In this section we derive maximum likelihood estimators for our optimal and
sub-optimal sufficient reductions, the asymptotic normality of the projection
matrix of the optimal SDR, with which we also obtain asymptotic tests of
dimension of both optimal and sub-optimal reductions.

###  4.1 Parameter Estimation via Maximum Likelihood

We assume a random sample $(y_{i},{\mathbf{x}}_{i},\mathbf{h}_{i})$,
$i=1,\ldots,n$, is drawn from the joint distribution of
$(Y,{\mathbf{X}},\mathbf{H})$ and that the conditional distribution models (5)
and (3) hold. Finding the maximum likelihood estimators of the reductions
derived in Section 3 requires first the estimation of the parameters
$\bm{\Delta},\bm{\mu},\bm{\mu}_{\scriptstyle{\mathbf{H}}},\mathbf{A},\bm{\beta},\bm{\tau}_{0},\bm{\tau}$,
in the joint density (8) with log-likelihood

$\sum_{i=1}^{n}\log
f_{{\mathbf{X}},\mathbf{H}}({\mathbf{x}}_{i},\mathbf{h}_{i}|y_{i};\bm{\Delta},\bm{\mu},\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}},\mathbf{A},\bm{\beta},\bm{\tau}_{0},\bm{\tau}).$
(62)

We maximize (62) in two steps. First, we maximize $\sum_{i=1}^{n}\log
f_{{\mathbf{X}}}({\mathbf{x}}_{i}|y_{i},\mathbf{h}_{i};\bm{\Omega})$ to
estimate the parameters
$\bm{\Omega}=\\{\bm{\Delta},\bm{\mu},\bm{\mu}_{\scriptstyle{\mathbf{H}}},\mathbf{A},\bm{\beta}\\}$.
Since ${\mathbf{X}}\mid(\mathbf{H},Y)$ follows a normal distribution, the
maximum likelihood estimator (MLE) of $\bm{\Omega}$ is obtained from fitting a
multivariate normal linear model of ${\mathbf{X}}$ on the centered
$(\mathbf{H},Y)$ via MLE. The MLE of $\mathbf{A}$ and $\bm{\beta}$ are
$(\widehat{\mathbf{A}},\widehat{\bm{\beta}})=\mathbb{X}^{T}\mathbb{L}\left(\mathbb{L}^{T}\mathbb{L}\right)^{-1}$,
where $\mathbb{X}$ denotes the $n\times p$ matrix with rows
$({\mathbf{x}}_{i}-\bar{{\mathbf{x}}})^{T}$, and $\mathbb{L}$ the
$n\times(r+q)$ matrix with rows
$\left((\mathbf{f}_{y_{i}}-\bar{\mathbf{f}}_{y})^{T},(\mathbf{h}_{i}-\bar{\mathbf{h}})^{T}\right)$,
$\bar{{\mathbf{x}}}=\sum_{i=1}^{n}{\mathbf{x}}_{i}/n$,
$\bar{\mathbf{f}}_{y}=\sum_{i=1}^{n}f_{y_{i}}/n$ and
$\bar{\mathbf{h}}=\sum_{i=1}^{n}\mathbf{h}_{i}/n$. The MLE of the covariance
matrix is
$\widehat{\bm{\Delta}}=\left(\mathbb{X}^{T}-(\widehat{\mathbf{A}},\widehat{\bm{\beta}})\mathbb{L}^{T}\right)\left(\mathbb{X}^{T}-(\widehat{\mathbf{A}},\widehat{\bm{\beta}})\mathbb{L}^{T}\right)^{T}/n$.

Next, we estimate $\bm{\Upsilon}=(\bm{\tau}_{0},\bm{\tau})$ maximizing the
conditional log-likelihood function

$\sum_{i=1}^{n}\log f_{\mathbf{H}}(\mathbf{h}_{i}\mid y_{i};\bm{\Upsilon}).$

Using parametrization (6), the joint probability mass function (5) can be
written as

$\displaystyle
P(\mathbf{H}|Y=y)=\exp\bigg{(}\sum_{j=1}^{q}\tau^{*}_{jj0}H_{j}+\sum_{j=1}^{q}\bm{\tau}^{T}_{jj}\mathbf{f}_{y}H_{j}$
$\displaystyle\hskip 85.35826pt+\sum_{1\leq j<j^{\prime}\leq
q}\tau^{*}_{jj^{\prime}0}H_{j}H_{j^{\prime}}+\sum_{1\leq j<j^{\prime}\leq
q}\bm{\tau}^{T}_{jj^{\prime}}\mathbf{f}_{y}H_{j}H_{j^{\prime}}\bigg{)}\frac{1}{G(\bm{\Gamma}_{y})}.$

Following Cheng et al. (2014), we consider a single binary variable $H_{j}$
and condition over the rest
$\mathbf{H}_{-j}=(H_{1},\dots,H_{j-1},H_{j+1},\dots,H_{q})$ to obtain

$\log\frac{P(H_{j}=1\mid\mathbf{H}_{-j},Y)}{P(H_{j}=0\mid\mathbf{H}_{-j},Y)}=\tau^{*}_{jj0}+\bm{\tau}^{T}_{jj}\mathbf{f}_{y}+\sum_{j\neq
j^{\prime}}\tau^{*}_{jj^{\prime}0}H_{j^{\prime}}+\sum_{j<j^{\prime}}\bm{\tau}^{T}_{jj^{\prime}}\mathbf{f}_{y}H_{j^{\prime}}.$
(63)

Thus, the conditional log-odds for a specific binary variable $H_{j}$ is
linear in the parameters. Moreover, the conditional maximum likelihood
estimators for these parameters can be obtained by fitting a logistic
regression of $H_{j}$ on
$(\mathbf{f}_{y},\mathbf{H}_{-j},\mathbf{f}_{y}\mathbf{H}_{-j})$, so that we
obtain estimators for $\bm{\tau}_{0}$ and $\bm{\tau}$ by fitting $q$
univariate logistic regressions. In particular, for the sample points
$(\mathbf{h}_{i}^{T},y_{i})=(h_{i1},\ldots,h_{iq},y_{i})$, for each binary
variable $j$ ($j=1,\ldots,q$), the conditional log-likelihood function is

$\ell_{j}(\bm{\tau}_{0},\bm{\tau};\mathbf{h}_{i},y_{i})=\frac{1}{n}\sum_{i=1}^{n}\log
P(h_{ij}\mid\mathbf{h}_{i,-j},y_{i})=\frac{1}{n}\sum_{i=1}^{n}\left(h_{ij}\epsilon_{ij}-\log(1+\exp(\epsilon_{ij}))\right),$
(64)

where $\mathbf{h}_{i,-j}=(h_{11},\ldots,h_{i,j-1},h_{i,j+1},\ldots,h_{iq})$
and

$\epsilon_{ij}=\log\frac{P(h_{ij}=1\mid\mathbf{h}_{i,-j},y)}{P(h_{ij}=0\mid\mathbf{h}_{i,-j},y_{i})}=\tau^{*}_{jj0}+\bm{\tau}^{T}_{jj}\mathbf{f}_{y_{i}}+\sum_{j\neq
j^{\prime}}\tau^{*}_{jj^{\prime}0}h_{ij^{\prime}}+\sum_{j\neq
j^{\prime}}\bm{\tau}^{T}_{jj^{\prime}}\mathbf{f}_{y_{i}}h_{ij^{\prime}}.$

To estimate $\bm{\Upsilon}$ we use the joint estimation algorithm proposed by
Cheng et al. (2014) that maximizes
$\sum_{j}\ell_{j}(\bm{\tau}_{0},\bm{\tau};\mathbf{h}_{i},y_{i})$.

###  4.2 Maximum Likelihood Estimation of the Reductions

To estimate the optimal SDR $\bm{\alpha}_{\mathbf{b}}$ in Corollary 2 and the
sub-optimal SDR $\bm{\alpha}_{\mathbf{c}}$ in Corollary 5, we need first to
estimate $\mathbf{b}$ in (40) and $\mathbf{c}_{1}$ and $\mathbf{c}_{2}$ in
(54). We use the ML estimators
$(\widehat{\bm{\Delta}},\widehat{\bm{\mu}},\widehat{\bm{\mu}}_{\scriptstyle{\mathbf{H}}},\widehat{\mathbf{A}},\widehat{\bm{\beta}},\widehat{\bm{\tau}}_{0},\widehat{\bm{\tau}})$
of the corresponding parameters in (62) from Section 4.1.

####  4.2.1 Optimal SDR

To estimate the minimal sufficient reduction in (32), or equivalently, derive
a basis estimate of
${\mathcal{S}}_{\bm{\scriptstyle\alpha}_{\mathbf{\scriptstyle b}}}$, we need
to first estimate $\mathbf{b}$ in (40). If
$d=\dim(\mathcal{S}_{\bm{\alpha}_{\tiny\mathbf{b}}})$, with
$d\leq\min\\{r,p+q(q+1)/2\\}$, the rank of $\mathbf{b}$ is also $d$ with
singular value decomposition

${\mathbf{b}}=\mathbf{U}^{T}\left(\begin{array}[]{cc}{\mathbf{K}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{array}\right)\mathbf{R},$ (65)

where $k_{1}\geq\ldots\geq k_{d}>0$ are the singular values of $\mathbf{b}$,
$\mathbf{K}=\mathrm{diag}(k_{1},\ldots,k_{d})$,
$\mathbf{U}^{T}=(\mathbf{U}_{1},\mathbf{U}_{0})$ is a $m\times m$ orthogonal
matrix with $m=p+q(q+1)/2$, $\mathbf{U}_{1}:m\times d$,
$\mathbf{U}_{0}:m\times(m-d)$, and
$\mathbf{R}^{T}=(\mathbf{R}_{1},\mathbf{R}_{0})$ is an $r\times r$ orthogonal
matrix with $\mathbf{R}_{1}:r\times d$, $\mathbf{R}_{0}:r\times(r-d)$. The
submatrices satisfy
$\mathbf{U}_{1}\mathbf{U}_{1}^{T}+\mathbf{U}_{0}\mathbf{U}_{0}^{T}=\mathbf{I}_{m}$,
$\mathbf{U}_{1}^{T}\mathbf{U}_{1}=\mathbf{I}_{d}$,
$\mathbf{U}_{0}^{T}\mathbf{U}_{0}=\mathbf{I}_{m-d}$,
$\mathbf{U}_{0}^{T}\mathbf{U}_{1}=\mathbf{0}$,
$\mathbf{R}_{1}\mathbf{R}_{1}^{T}+\mathbf{R}_{0}\mathbf{R}_{0}^{T}=\mathbf{I}_{r}$,
$\mathbf{R}_{1}^{T}\mathbf{R}_{1}=\mathbf{I}_{d}$,
$\mathbf{R}_{0}^{T}\mathbf{R}_{0}=\mathbf{I}_{r-d}$,
$\mathbf{R}_{0}^{T}\mathbf{R}_{1}=\mathbf{0}$. Then,

$\mathbf{b}=\mathbf{U}_{1}\mathbf{K}\mathbf{R}^{T}_{1},$ (66)

and, as a consequence, $\bm{\alpha}_{\mathbf{b}}$ in Corollary 2 can be set to
$\mathbf{U}_{1}$. Plugging in the ML estimators
$(\widehat{\bm{\Delta}},\widehat{\mathbf{A}},\widehat{\bm{\beta}},\widehat{\bm{\tau}})$
we obtain that the ML estimator of $\mathbf{b}$ is

$\displaystyle\widehat{\mathbf{b}}$
$\displaystyle=\begin{pmatrix}\widehat{\bm{\Delta}}^{-1}\widehat{\mathbf{A}}\\\
\mathbf{L}_{q}\widehat{\bm{\tau}}-\widehat{\bm{\beta}}^{T}\widehat{\bm{\Delta}}^{-1}\widehat{\mathbf{A}}\\\
\mathbf{J}_{q}\widehat{\bm{\tau}}\end{pmatrix}=\begin{pmatrix}\mathrm{unvec}(\widehat{\bm{\vartheta}}_{1,1})\\\
\mathrm{unvec}(\widehat{\bm{\vartheta}}_{2,1})\\\
\mathrm{unvec}(\widehat{\bm{\vartheta}}_{5,1})\end{pmatrix}.$ (67)

The singular value decomposition of the MLE of $\mathbf{b}$ is

$\widehat{\mathbf{b}}=\widehat{\mathbf{U}}^{T}\left(\begin{array}[]{cc}\widehat{\mathbf{K}}_{1}&\mathbf{0}\\\
\mathbf{0}&\widehat{\mathbf{K}}_{0}\end{array}\right)\widehat{\mathbf{R}},$
(68)

where
$\widehat{\mathbf{K}}_{1}=\text{diag}(\widehat{k}_{1},\ldots,\widehat{k}_{d})$,
$\widehat{\mathbf{K}}_{0}=\mathrm{diag}(\widehat{k}_{d+1},\ldots,\widehat{k}_{\min(m,r)})$,
$\widehat{k}_{i}$ are the singular values of $\widehat{\mathbf{b}}$ in
decreasing order, $\widehat{\mathbf{U}}$ is an $m\times m$ orthogonal matrix
whose columns are the left singular vectors of $\widehat{\mathbf{b}}$, and
$\widehat{\mathbf{R}}$ is an $r\times r$ orthogonal matrix, whose columns are
the right-singular vectors of $\widehat{\mathbf{b}}$. Let
$\widehat{\mathbf{U}}_{1}$ be the first $d$ columns of $\widehat{\mathbf{U}}$,
$\widehat{\mathbf{R}}_{1}$ the first $d$ columns of
$\widehat{\mathbf{R}}^{T}$, and
$\widehat{\mathbf{B}}=\widehat{\mathbf{K}}_{1}\widehat{\mathbf{R}}_{1}^{T}$.
An estimator of $\mathbf{b}$ subject to
$d=\dim(\mathcal{S}_{\bm{\alpha}_{\tiny\mathbf{b}}})$ is

$\widehat{\mathbf{b}}^{(d)}=\widehat{\mathbf{U}}_{1}\widehat{\mathbf{K}}_{1}\widehat{\mathbf{R}}_{1}^{T}=\widehat{\mathbf{U}}_{1}\widehat{\mathbf{B}}.$
(69)

and an estimator of the reduction $\bm{\alpha}_{\mathbf{b}}$ in Corollary 2 is

$\widehat{\bm{\alpha}}_{\mathbf{b}}=\widehat{\mathbf{U}}_{1}.$ (70)

####  4.2.2 Sub-optimal SDR: $\widehat{\mathcal{S}}_{\mathbf{c}}$

To obtain an estimator for the space $\mathcal{S}_{\bm{\alpha}_{\mathbf{c}}}$
in (54) that gives the sub-optimal sufficient reduction (47), we set
$\mathbf{c}=(\mathbf{c}_{1},\mathbf{c}_{2})$, where $\mathbf{c}_{1}$ and
$\mathbf{c}_{2}$ are given in (59), with $\hbox{rank}(\mathbf{c}_{1})=d_{1}$
and $\hbox{rank}(\mathbf{c}_{2})=d_{2}$. Plugging in the MLE
$(\widehat{\bm{\Delta}},\widehat{\bm{\mu}},\widehat{\bm{\mu}}_{\scriptstyle{\mathbf{H}}},\widehat{\mathbf{A}},\widehat{\bm{\beta}},\widehat{\bm{\tau}}_{0},\widehat{\bm{\tau}})$
of the corresponding parameters in (62) from Section 4.1, we obtain estimators
of $\mathbf{c}_{1}$ and $\mathbf{c}_{2}$,

$\widehat{\mathbf{c}}_{1}=\left(\begin{array}[]{c}\widehat{\bm{\Delta}}^{-1}\widehat{\mathbf{A}}\\\
-\widehat{\bm{\beta}}^{T}\widehat{\bm{\Delta}}^{-1}\widehat{\mathbf{A}}\end{array}\right),\quad\widehat{\mathbf{c}}_{2}=\widehat{\tau}.$

We then consider their respective SVD decompositions as in Section 4.2.1. Let
$\widehat{\mathbf{U}}_{c_{1}}$ denote the first $d_{1}$ left eigenvectors of
$\widehat{\mathbf{c}}_{1}$ and $\widehat{\mathbf{U}}_{c_{2}}$ the first
$d_{2}$ left eigenvectors of $\widehat{\mathbf{c}}_{2}$. Then, an estimator
for the sub-optimal sufficient reduction in (47) is defined as

$\widehat{\bm{\alpha}}_{\mathbf{c}}=\left(\begin{array}[]{cc}\widehat{\mathbf{U}}_{\mathbf{c}_{1}}&0\\\
0&\widehat{\mathbf{U}}_{\mathbf{c}_{2}}\end{array}\right).$

###  4.3 Asymptotic distribution of the optimal sufficient reduction estimator

In this section we derive the asymptotic distribution of the projection onto
the column space of the estimated optimal sufficient reduction
$\widehat{\bm{\alpha}}_{\mathbf{b}}$ in (70),
$\mathbf{P}_{\widehat{\bm{\alpha}}_{\mathbf{b}}}=\widehat{\bm{\alpha}}_{\mathbf{b}}(\widehat{\bm{\alpha}}_{\mathbf{b}}^{T}\widehat{\bm{\alpha}}_{\mathbf{b}})^{-1}\widehat{\bm{\alpha}}_{\mathbf{b}}^{T}$.
We use this result in the derivation of the asymptotic tests for dimension in
Section 4.4 and for inference about the sufficient dimension reduction.

###### Proposition 6

Suppose that $({\mathbf{X}},\mathbf{H})|Y$ has probability mass function (9)
with the natural parameters $\bm{\eta}_{Y}$ satisfying (26) and that
$\mathbf{b}$ has rank $d$. Then,

$\sqrt{n}\;\mathrm{vec}\left(\mathbf{P}_{\widehat{\bm{\alpha}}_{\mathbf{b}}}-\mathbf{P}_{\bm{\alpha}_{\mathbf{b}}}\right)\xrightarrow{\mathcal{D}}\mathcal{N}\left(\mathbf{0},\mathbf{V}_{\widehat{\bm{\alpha}}_{\mathbf{b}}}\right),$

with

$\mathbf{V}_{\widehat{\bm{\alpha}}_{\mathbf{b}}}=(\mathbf{I}_{m^{2}}\otimes{\mathbf{K}}_{mm})\left(\mathbf{b}^{-}\otimes\mathbf{Q}_{\mathbf{b}}\right)^{T}\mathbf{V}_{rcl}\left(\mathbf{b}^{-}\otimes\mathbf{Q}_{\mathbf{b}}\right)(\mathbf{I}_{m^{2}}\otimes{\mathbf{K}}_{mm}),$
(71)

where $\mathbf{b}^{-}$ is the Moore-Penrose generalized inverse of
$\mathbf{b}$,

$\mathbf{V}_{rcl}=\mathbf{W}\mathbf{M}\mathbf{V}\mathbf{M}^{T}\mathbf{W}^{T},$
(72)

with

$\mathbf{V}^{-1}=\mathrm{E}\left({\mathbf{F}}_{y}^{T}\mathbf{J}{\mathbf{F}}_{y}\right),$
(73)

$\mathbf{F}_{y}$ is defined in (26), $\mathbf{J}$ is the matrix of partial
derivatives given by

$\mathbf{J}=\frac{\partial^{2}\psi(\bm{\eta}_{y})}{\partial\bm{\eta}_{y}\partial\bm{\eta}_{y}^{T}},$
(74) $\displaystyle\mathbf{M}$
$\displaystyle=\begin{pmatrix}\mathbf{0}_{pr\times
p}&\mathbf{I}_{pr}&\mathbf{0}_{pr\times q}&\mathbf{0}_{pr\times
qr}&\mathbf{0}_{pr\times m_{p}}&\mathbf{0}_{pr\times qp}&\mathbf{0}_{pr\times
k_{q}}&\mathbf{0}_{pr\times rk_{q}}\\\ \mathbf{0}_{qr\times
p}&\mathbf{0}_{qr\times pr}&\mathbf{0}_{qr\times
q}&\mathbf{I}_{qr}&\mathbf{0}_{qr\times m_{p}}&\mathbf{0}_{qr\times
qp}&\mathbf{0}_{qr\times k_{q}}&\mathbf{0}_{qr\times rk_{q}}\\\
\mathbf{0}_{rk_{q}\times p}&\mathbf{0}_{rk_{q}\times
pr}&\mathbf{0}_{rk_{q}\times q}&\mathbf{0}_{rk_{q}\times
qr}&\mathbf{0}_{rk_{q}\times m_{p}}&\mathbf{0}_{rk_{q}\times
qp}&\mathbf{0}_{rk_{q}\times k_{q}}&\mathbf{I}_{rk_{q}}\end{pmatrix},$ (75)

and

$\mathbf{W}=\left(\begin{array}[]{c}\mathbf{I}_{r}\otimes\left(\begin{array}[]{l}\mathbf{I}_{p}\\\
\mathbf{0}_{q\times p}\\\ \mathbf{0}_{k_{q}\times
p}\end{array}\right),\mathbf{I}_{r}\otimes\left(\begin{array}[]{l}\mathbf{0}_{p\times
q}\\\ \mathbf{I}_{q}\\\ \mathbf{0}_{k_{q}\times
q}\end{array}\right),\mathbf{I}_{r}\otimes\left(\begin{array}[]{l}\mathbf{0}_{p\times
k_{q}}\\\ \mathbf{0}_{q\times k_{q}}\\\
\mathbf{I}_{k_{q}}\end{array}\right)\end{array}\right).$ (76)

###  4.4 Tests for dimension

We propose two asymptotic tests for the dimension of the sufficient reduction
in optimal SDR. We adapt these tests for the case of sub-optimal SDR, to
estimate the dimension of the continuous predictors separately from the binary
predictors.

The dimension of the sufficient reduction is the rank of $\mathbf{b}$ in (40).
We estimate the rank $d$ of $\mathbf{b}$ by sequentially testing the
hypotheses

$H_{0}:\mathrm{rank}(\mathbf{b})=j\quad\text{vs.}\quad
H_{1}:\mathrm{rank}(\mathbf{b})>j,$ (77)

for $j=0,1,\ldots,\min(r,m)$, where $m=p+q(q+1)/2$. For a fixed level
$\alpha$, the estimated rank is the smallest value of $j$ for which the null
is not rejected.

Bura and Yang (2011) proposed asymptotic tests for the rank of random matrices
in sequential hypothesis testing. To construct the corresponding tests for
dimension, we consider the singular value decomposition of $\mathbf{b}$ in
(66) and $\widehat{\mathbf{b}}$ in (69) with $d=j$.

The first statistic we use to test (77) is
$\Lambda_{1}(j)=n\sum_{i=j+1}^{\min(m,r)}\hat{k}_{i}^{2},$ where
$\widehat{k}_{i}$’s are the singular values of $\widehat{\mathbf{b}}$ in
descending order. Proposition 6 obtains the asymptotic normality of
$\widehat{\mathbf{b}}$ with covariance ${\mathbf{V}}_{rlc}$ in (72). When
$\hbox{rank}(\mathbf{b})=j$,

$\Lambda_{1}(j)\overset{\mathcal{D}}{\longrightarrow}\sum_{i=1}^{s}\omega_{i}X^{2}_{i},$
(78)

where $s=\min(\mathrm{rank}(\mathbf{V}_{rcl}),(r-j)(m-j))$, $X_{i}^{2}$ are
independent chi-squared random variables with 1 degree of freedom, and the
weights are the descending eigenvalues of
$\mathbf{Q}=(\mathbf{R}_{0}^{T}\otimes\mathbf{U}_{0}^{T})\mathbf{V}_{rcl}(\mathbf{R}_{0}\otimes\mathbf{U}_{0})$
[see Bura and Yang (2011)]. In practice, the weights $\omega_{i}$,
$i=1,\ldots,s$ are replaced by
$\widehat{\omega}_{1}\geq\widehat{\omega}_{2}\geq\ldots\geq\widehat{\omega}_{s}$,
the descending eigenvalues of

$\displaystyle\widehat{\mathbf{Q}}$
$\displaystyle=(\widehat{\mathbf{R}}_{0}^{T}\otimes\widehat{\mathbf{U}}_{0}^{T})\widehat{\mathbf{V}}_{rcl}(\widehat{\mathbf{R}}_{0}\otimes\widehat{\mathbf{U}}_{0}),$
(79)

where $\widehat{\mathbf{V}}_{rcl}$ is a consistent estimate of
$\mathbf{V}_{rcl}$. This test rejects $H_{0}$ if $\Lambda_{1}(j)>q_{\alpha}$,
where $q_{\alpha}$ is the $(1-\alpha)$ percentile of the distribution of
$\sum_{i=1}^{s}\widehat{\omega}_{i}X^{2}_{i}$. We estimate $q_{\alpha}$ from
the empirical distribution function of $\Lambda_{1}$, by generating $10000$
realizations of $\sum_{i=1}^{s}\widehat{\omega}_{i}X_{i}^{2}$ and computing
the empirical quantile $\widehat{q}_{\alpha}$.

The second is a Wald test with test statistic,
$\Lambda_{2}(j)=n\mathrm{vec}(\widehat{\mathbf{K}}_{0})^{T}\widehat{\mathbf{Q}}^{{\dagger}}\mathrm{vec}(\widehat{\mathbf{K}}_{0})$,
where $\widehat{\mathbf{K}}_{0}$ is defined in (68) and
$\widehat{\mathbf{Q}}^{{\dagger}}$ is the Moore-Penrose inverse of
$\widehat{\mathbf{Q}}$ in (79).

Following Bura and Yang (2011), since $\widehat{\mathbf{b}}$ is asymptotically
normal, if $j=\mathrm{rank}(\mathbf{b})$, then
$\Lambda_{2}(j)\overset{\mathcal{D}}{\longrightarrow}\chi^{2}(s)$, where the
degrees of freedom are $s=\min(\mathrm{rank}(\mathbf{V}_{rcl}),(r-j)(m-j))$.
The rejection region is $\Lambda_{2}(j)>\chi^{2}_{\alpha}(s)$, where
$\chi^{2}_{\alpha}(s)$ is the $(1-\alpha)$ percentile of the $\chi^{2}(s)$
distribution.

##  5 Variable selection

Identifying variables that are not associated with the outcome is important
for both interpretation and for improving the predictive power of a classifier
or a regression model. We propose a method to simultaneously obtain the
sufficient reduction and carry out variable selection by removing redundant
variables from the reduction. This is obtained jointly with the estimate of
the reduction by introducing structured regularization on a matrix
factorization problem.

In particular, we exploit the factorization of the full rank maximum
likelihood estimate $\widehat{\mathbf{b}}$ into a relevant full-rank factor
$\mathbf{C}\in{\mathbb{R}}^{p+q(q+1)/2\times d}$, which determines the
reduction, and a matrix $\mathbf{B}$ that is immaterial.

The building block of the procedure proposed here is to note that the reduced
rank estimator
$\widehat{\mathbf{b}}^{(d)}=\widehat{\mathbf{U}}_{1}\widehat{\mathbf{B}}$ in
(69) is also the solution to the least squares minimization problem

$\min_{\mathbf{C}\in{\mathbb{R}}^{p+q(q+1)/2\times
d},\mathbf{C}^{T}\mathbf{C}=\mathbf{I},\mathbf{B}\in{\mathbb{R}}^{d\times
r}}(\mathrm{vec}(\widehat{\mathbf{b}})-\mathrm{vec}(\mathbf{C}\mathbf{B}))^{T}(\mathrm{vec}(\widehat{\mathbf{b}})-\mathrm{vec}(\mathbf{C}\mathbf{B})),$
(80)

where $\widehat{\mathbf{b}}$ is the maximum likelihood estimator of
$\mathbf{b}$. The solution can be expressed as
$\widehat{\mathbf{C}}=\widehat{\mathbf{U}}_{1}\mathbf{V}$, for some orthogonal
matrix $\mathbf{V}\in{\mathbb{R}}^{d\times d}$, so that
$\mathrm{span}(\widehat{\mathbf{C}})=\mathrm{span}(\widehat{\mathbf{U}}_{1})$.

All sufficient reductions in Section 4.2.1 are of the form
$\mathbf{R}({\mathbf{X}},\mathbf{H})=\mathbf{U}_{1}^{T}({\mathbf{t}}({\mathbf{X}},\mathbf{H})-\mathrm{E}({\mathbf{t}}({\mathbf{X}},\mathbf{H}))$.
If $t_{j}$ is the $j$th component of ${\mathbf{t}}({\mathbf{X}},\mathbf{H})$,
and $t_{j}$ is not associated with $Y$, the $j$th row of $\mathbf{U}_{1}$ is
zero. Therefore, identifying predictors that are conditionally independent of
$Y$ corresponds to identifying the rows of $\mathbf{C}$ that contain only $0$.
This can be achieved using mixed-norm regularizers that are known to induce
structured sparsity in the estimates (Bach et al., 2012).

The proposed procedure is as follows. For a fixed $d$, once we obtain
$\widehat{\mathbf{b}}^{(d)}=\widehat{\mathbf{U}}_{1}\widehat{\mathbf{B}}$ in
(69), we solve

$\operatorname*{arg\,min}_{\mathbf{C}\in{\mathbb{R}}^{(p+q(q+1)/2)\times
d},\mathbf{C}^{T}\mathbf{C}=\mathbf{I}}\left(\mathrm{vec}(\widehat{\mathbf{b}})-\mathrm{vec}(\mathbf{C}\widehat{\mathbf{B}})\right)^{T}\left(\mathrm{vec}(\widehat{\mathbf{b}})-\mathrm{vec}(\mathbf{C}\widehat{\mathbf{B}})\right)+\lambda\Omega(\mathbf{C}),$
(81)

where $\Omega(\mathbf{C})$ is a mixed-norm regularizer which penalizes the
rows of $\mathbf{C}$ in a similar manner to group-lasso. The specific form of
$\Omega(\mathbf{C})$ depends on the type of predictor variables involved in
the problem, as follows.

  * (a)

When all predictors are continuous (normal), we use the penalty
$\Omega(\mathbf{C})=\sum_{j=1}^{p}||{\mathbf{C}}_{j}||_{2}$, with
$\mathbf{C}_{j}$ the $j$th row of $\mathbf{C}$. In this case the sufficient
reduction contains no interaction terms and each row of $\mathbf{C}$ affects a
single element of ${\mathbf{X}}$. Hence, by shrinking the $j$th row of
$\mathbf{C}$ to $0$, the computed reduction becomes insensitive to the
measured value of $X_{j}$. When all predictors are continuous, under the
assumed model the optimization problem is indeed fairly similar to group lasso
(Yuan and Lin, 2006) as can be seen after rewriting (81) as

$\operatorname*{arg\,min}_{\mathbf{C}\in{\mathbb{R}}^{(p+q(q+1)/2)\times
d},\mathbf{C}^{T}\mathbf{C}=\mathbf{I}}\lVert\mathrm{vec}(\widehat{\mathbf{b}})-(\widehat{\mathbf{B}}^{T}\otimes\mathbf{I})\mathrm{vec}(\mathbf{C})\rVert_{2}^{2}+\lambda\sum_{j=1}^{p}\lVert\mathbf{C}_{j}\rVert_{2}.$

  * (b)

When all predictors are binary, the sufficient reduction includes interaction
effects $H_{i}H_{j}$. Thus, to discard the effect of a given binary variable,
say $H_{j}$, we need to set all the entries in $\mathbf{C}$ related to $H_{j}$
to zero. For a reduction of dimension $d$, there are $d$ such entries related
to the main effects and $d(q-1)$ related to the interaction terms. The
grouping of the entries of $\mathbf{C}$ does not form a partition, since the
entries affecting the interaction terms appear twice. For instance, assume for
simplicity that $d=1$. Parameter $\theta_{13}$ operates on variables $H_{1}$
and $H_{3}$ and then it enters the regularizer in groups
$\\{\eta_{1},\theta_{12},\theta_{13},\dots,\theta_{1q}\\}$ and
$\\{\eta_{3},\theta_{13},\theta_{23},\dots,\theta_{q3}\\}$. Both groups of
parameters overlap at $\theta_{13}$. Thus, the regularizer inducing the
desired sparsity structure is a mixed-norm regularizer with overlapping
groups,
$\Omega(\mathbf{C})=\sum_{g\in\mathcal{G}}\lVert\mathbf{C}_{g}\rVert_{2}$.
Here, $g\subset\\{1,\dots,dq(q+1)/2\\}$ indicates the subset of entries that
affect the binary variable $H_{i}$ and $\mathcal{G}$ is the collection of such
groups. Moreover, each binary variable is associated with two groups, one
derived from the main effects and one from the interaction terms, since they
typically have rather different scales. The obtained regularized problem can
be solved using algorithms for overlapping group lasso, as proposed, for
example, in Liu and Ye (2010).

  * (c)

When the predictors are mixed normal and binary, we combine the regularizers
described in (a) and (b) in a single penalty
$\Omega(\mathbf{C})=\gamma\sum_{j=1}^{p}\lVert{\mathbf{C}}_{j}\rVert_{2}+(1-\gamma)\sum_{g\in\mathcal{G}}\lVert\mathbf{C}_{G_{i}}\rVert_{2}$.
The value of $\gamma$ serves as a tuning weight for the amount of
regularization in the continuous and binary parts, respectively. In sub-
optimal SDR, we carry out variable selection separately for the continuous and
binary variables as described in (a) and (b).

Selection of hyperparameters $(\lambda,\gamma)$ is done using 10-fold cross
validation, with prediction error as the optimization criterion. The procedure
starts by estimating a maximum value $\lambda_{m}$ so that the whole estimate
vanishes for any $\lambda>\lambda_{m}$. We then set a grid of $n_{\lambda}$
candidate values for $\lambda$, uniformly spaced on a logarithmic scale
between $0$ and $\lambda_{m}$. We typically use $n_{\lambda}=100$. For
$\gamma$ we test 11 values uniformly spaced in $[0,1]$. In each fold, an
initial full-rank estimate of the reduction is computed using the training set
and then factorized using truncated SVD to give $\widehat{\mathbf{B}}$ and an
initial estimate for $\mathbf{C}$. Problem (81) is solved for each pair of
candidate values $(\lambda_{k},\gamma_{k})$. The obtained reduction is applied
to both the training and the test sample. With the reduced training set we fit
a prediction model and then we evaluate the prediction error on the reduced
test sample. The average prediction error over the ten cross-validation folds
is then computed for each candidate pair $(\lambda_{k},\gamma_{k})$. We pick
the combination that attains the smallest mean prediction error.

##  6 Simulation Studies

We assess the performance of the proposed methods in estimating the sufficient
reduction and its dimension, out-of-sample prediction, and variable selection
in simulations.

In all our simulations the response is generated from the uniform distribution
on the integers $\\{1,\ldots,r+1\\}$, with $r=5$, and set
$\mathbf{f}_{y}=I(y=j)-n_{j}/n$, where $I$ is the indicator function, $n$
denotes the total sample size and $n_{j}$ the number of observations in
category $j$ for $j=1,\ldots,r$. All reported results are based on sample
sizes $n=100,200,300,500,750$, and 100 repetitions.

###  6.1 Estimation, prediction and dimension tests

We assess the accuracy of estimating $\mathrm{span}(\bm{\alpha})$ with
$\mathrm{span}(\widehat{\bm{\alpha}})$ using
$||\mathbf{P}_{\bm{\scriptstyle\alpha}}-\mathbf{P}_{\widehat{\bm{\scriptstyle\alpha}}}||_{2}$
[see Ye and Lim (2016)]. The prediction error is computed as
$||\mathbf{P}_{\bm{\scriptstyle\alpha}^{T}({\mathbf{X}}_{N},\mathbf{H}_{N})}-\mathbf{P}_{\widehat{\bm{\scriptstyle\alpha}}^{T}({\mathbf{X}}_{N},\mathbf{H}_{N})}||_{2}$,
where $({\mathbf{X}}_{N},\mathbf{H}_{N})$ is a new sample of size $N=2000$
that is independent of the training sample. We estimate the sufficient
reduction using the true $d$.

####  6.1.1 Continuous predictors

We generate $p$-variate continuous predictors as ${\mathbf{X}}\mid
Y=y\sim\mathcal{N}(\bm{\mu}_{y},\bm{\Delta})$ with
$\bm{\mu}_{y}=\mathbf{A}\mathbf{f}_{y}$ for
$\mathbf{A}=\bm{\Delta}\bm{\alpha}\bm{\xi}$, where $\bm{\alpha}\in{}^{p\times
d}$ of $\mathrm{rank}(\bm{\alpha})=d$ and $\bm{\xi}\in{}^{d\times r}$. We let
$p=20$ and $\mathbf{0}_{l}$, $\mathbf{1}_{l}$ denote the $l$-vectors of zeros
and ones, respectively.

  * (a)

For $d=1$, we set $\bm{\xi}=\mathbf{1}_{r}^{T}$
$\bm{\alpha}=(\mathbf{0}_{p/2}^{T},\mathbf{1}_{p/2}^{T})^{T}$,
$\bm{\Delta}=5(\mathbf{I}_{p}+\rho\bm{\alpha}\bm{\alpha}^{T})$ with
$\rho=0.55$.

  * (b)

For $d=2$, we set

$\bm{\xi}=\begin{pmatrix}1&1&1&1&1\\\ 0&0&0&1&1\end{pmatrix},$

and $\bm{\alpha}=(\bm{\alpha}_{1},\bm{\alpha}_{2})$ be an orthonormal basis of
$\mathrm{span}((\mathbf{0}_{p/2}^{T},\mathbf{1}_{p/2}^{T})^{T},(\mathbf{0}_{p/2}^{T},\mathbf{1}_{p/4}^{T},-\mathbf{1}_{p/4}^{T})^{T})$,
$\bm{\Delta}=5(\mathbf{I}_{p}+\rho_{1}\bm{\alpha}_{1}\bm{\alpha}_{1}^{T}+\rho_{2}\bm{\alpha}_{2}\bm{\alpha}_{2}^{T})$
for $\rho_{1}=0.55$ and $\rho_{2}=0.25$.

####  6.1.2 Binary predictors

We generate $q=10$ binary predictors assuming that $\mathbf{H}\mid Y$ follows
an Ising model with parameters $\\{\bm{\tau}_{0},\bm{\tau}\\}$, where
$\bm{\tau}=[\mathrm{\,vech}({\bm{\tau}}_{1}),\dots,\mathrm{\,vech}({\bm{\tau}}_{r})]$,
$\bm{\tau}_{j}$ are $q\times q$ matrices and set $\bm{\tau}_{0}=\mathbf{0}$.

  * (a)

For $d=1$, and $j=1,\dots,r$,
$\bm{\tau}_{j}=3\times{\mathbf{K}_{1}}/\sqrt{\sum_{ij}([\mathbf{K}_{1}]_{ij})}$
with

${\mathbf{K}_{1}}=\left(\begin{array}[]{cccccc|c}1&30&5&0&0&0&\mathbf{0}_{1\times
4}\\\ 30&1&10&0&0&0&\vdots\\\ 5&10&1&30&0&0&\vdots\\\ 0&0&30&1&30&0&\vdots\\\
0&0&0&30&1&30&\vdots\\\ 0&0&0&0&30&1&\vdots\\\
0&0&0&0&0&30&\mathbf{0}_{1\times 4}\\\ \cline{1-6}\cr\mathbf{0}_{3\times
1}&\cdots&\cdots&\cdots&\cdots&\lx@intercol\hfil\mathbf{0}_{3\times
1}\hfil\lx@intercol&\mathbf{0}_{3\times 4}\\\ \end{array}\right).$

  * (b)

For $d=2$,
$\bm{\tau}_{j}=3\times\mathbf{K}_{1}/\sqrt{\sum_{ij}([\mathbf{K}_{1}]_{ij})}$,
for $j=1,3,4,5$, and

$\bm{\tau}_{2}=\frac{12}{\sqrt{6}}\times\left(\begin{array}[]{ll}{\mathbf{I}}_{6}&{\mathbf{0}}_{6\times
4}\\\ {\mathbf{0}}_{4\times 6}&{\mathbf{0}}_{4\times 4}\end{array}\right).$

####  6.1.3 Mixed predictors

  * (a)

For $d=1$, we use the same parameters as for continuous
${\mathbf{X}}\mid(\mathbf{H},Y)$ and binary variables $\mathbf{H}\mid Y$ in
Sections 6.1.1, 6.1.2, respectively. Moreover, we set $\bm{\mu}_{H}=0$ and
$\bm{\beta}=\left(\mathbf{1}_{p\times 6}/10,\mathbf{0}_{p\times
4}\right)\in{}^{p\times q}$. to induce sparsity in the binary predictors.

  * (b)

For $d=2$, we generate $\mathbf{H}\mid Y$ as in Section 6.1.2 with dimension 1
and ${\mathbf{X}}\mid(\mathbf{H},Y)$ as in (a) with dimension 2.

Figure 1: Estimation error and out of sample prediction of optimal SDR with
continuous, binary and mixed predictors for $d=1$ (Left) and $d=2$ (Right).

In Figure 1, we plot the estimation error
$\left\lVert\mathbf{P}_{\widehat{\bm{\alpha}}}-\mathbf{P}_{\bm{\alpha}}\right\rVert_{2}$
and the prediction error
$\left\lVert\mathbf{P}_{\widehat{\bm{\alpha}}^{T}{\mathbf{X}}_{N}}-\mathbf{P}_{\bm{\alpha}^{T}{\mathbf{X}}_{N}}\right\rVert_{2}$
for optimal SDR on the $y$-axis versus the training sample size on the
$x$-axis across all our simulation scenarios. For all types of predictors the
prediction is smaller than the estimation error and both decrease as the
sample size increases. Moreover, both increase as the dimension increases from
1 to 2 in the left and right panels, respectively, across types of predictors.
When comparing types of predictors, continuous predictors exhibit higher
estimation and prediction errors across sample sizes and mixed predictors
result in the highest estimation and prediction errors.

In Figure 2 we plot the estimation and prediction error of sub-optimal SDR,
where the continuous and binary variables are reduced separately. The pattern
of behavior is consistent with that of optimal SDR in Figure 1, with the
continuous variables inducing larger errors of both types across sample sizes
and $d=1,2$. Again, the errors are smaller for dimension 1.

Figure 2: Estimation error and out of sample prediction with mixed predictors.
Suboptimal Reduction for $d=1$ (Left) and $d=2$ (Right).

Under the same simulation settings, we also evaluate the performance of our
simultaneous variable selection and dimension reduction method that is
presented in Section 5. In Table 2 we report the proportion of variables
correctly identified as non-relevant (true positives, TP) and the proportion
of variables erroneously assessed as non-relevant (false negatives, FN).
Between $d=1$ and $d=2$, TP is higher across sample sizes, whereas FN is
lower. Both rates improve substantially as the sample size increases. When all
predictors are continuous both rates are lower across sample sizes. This is
expected since the inclusion of a binary variable results in second order
interaction effects in the reduction. Therefore, to rule out a binary variable
both its own-coefficient and all the coefficients of its interaction terms
must be zero. Overall, our regularized SDR approach achieves high true
positive and small false negative rates for reasonable sample sizes.

Table 2: Accuracy of the regularized estimator in variable selection.

Sample Size Predictors $d$ Rates 100 200 300 500 750 Continuous $1$ TP 0.653
0.751 0.796 0.851 0.889 FN 0.314 0.17 0.095 0.044 0.012 $2$ TP 0.521 0.591
0.629 0.748 0.843 FN 0.165 0.048 0.014 0.004 0.002 Binary $1$ TP 0.188 0.310
0.400 0.55 0.623 FN 0.167 0.117 0.045 0.015 0.018 $2$ TP 0.255 0.300 0.368
0.458 0.528 FN 0.048 0.020 0.012 0.000 0.000 Mixed $1$ TP 0.632 0.592 0.589
0.671 0.674 FN 0.493 0.333 0.196 0.200 0.170 $2$ TP 0.596 0.656 0.639 0.583
0.610 FN 0.451 0.413 0.325 0.163 0.124

In Table 3 we report the proportion of times out of $100$ replications that
the dimension $d$ was correctly estimated based on the sequential tests of
dimension in Section 4.4 for all our simulation settings. The sample size has
a noticeable effect in the accuracy of the estimation of dimension, as
expected since both tests are asymptotic. The weighted $\chi^{2}$ test
accuracy suffers more from increasing the dimension and all binary predictors
as compared to that of the chi-squared test, across sample sizes. For mixed
predictors, as well, the chi-squared test exhibits higher accuracy for both
optimal and sub-optimal SDR across sample sizes.

Table 3: Proportion of correct dimension estimation under the simulation
settings in Section 6.

Predictors Sample Size Dimension Test Method 100 200 300 500 750 Continuous
$d=1$ Weighted $\chi^{2}$ 0.60 0.82 0.83 0.89 0.95 $\chi^{2}$ 0.00 0.99 0.99
1.00 0.98 $d=2$ Weighted $\chi^{2}$ 0.65 0.77 0.86 0.94 0.94 $\chi^{2}$ 0.2
1.00 0.99 0.97 0.95 Binary $d=1$ Weighted $\chi^{2}$ 0 0.80 0.96 0.94 0.94
$\chi^{2}$ 0 0.02 0.94 0.99 0.94 $d=2$ Weighted $\chi^{2}$ 0 0.02 0.30 0.66
0.96 $\chi^{2}$ 0.02 0.20 0.92 0.96 0.94 Mixed $d=1$ Weighted $\chi^{2}$
Optimal 0 0.2 0.45 0.64 0.94 Weighted $\chi^{2}$ Sub-optimal (cts) 0.68 0.84
0.89 0.90 0.95 Weighted $\chi^{2}$ Sub-optimal (bin) 0.5 0.75 0.87 0.94 0.95
$\chi^{2}$ Optimal 0 0.18 1 0.98 0.98 $\chi^{2}$ Sub-optimal (cts) 0 0.95 0.98
0.98 0.95 $\chi^{2}$ Sub-optimal (bin) 0 0.18 0.94 0.98 0.96 $d=2$ Weighted
$\chi^{2}$ Optimal 0 0.06 0.30 0.45 0.90 Weighted $\chi^{2}$ Sub-optimal (cts)
0.60 0.75 0.84 0.95 0.95 Weighted $\chi^{2}$ Sub-optimal (bin) 0 0.08 0.40
0.56 0.96 $\chi^{2}$ Optimal 0.08 0.36 0.64 0.92 0.93 $\chi^{2}$ Sub-optimal
(cts) 0.12 0.98 0.99 0.96 0.95 $\chi^{2}$ Sub-optimal (bin) 0.22 0.30 0.96
0.96 0.95

##  7 Data Analyses

We compare our method with other approaches such as generalized linear models
and principal component regression in two data applications. In particular, we
compare our methods with PCA and PCAmix in Sections 7.1 and 7.2. PCAmix
(Chavent et al., 2012, 2014) is a version of PCA that accommodates mixed
variables and implements PCA with metrics; i.e., Generalized Singular Value
Decomposition (GSVD) of pre-processed data [see Chavent et al. (2014) for
details]. PCAmix is ordinary standard PCA, when all variables are continuous,
and standard multiple correspondence analysis (MCA), when all variables are
categorical (Greenacre and Blasius (2006), Zhu et al. (2011), Camiz and Gomes
(2013)).

###  7.1 Krzanowski Data Sets

Krzanowski (1975) studied the problem of discriminating between two groups in
the presence of both binary and continuous explanatory variables. Krzanowski
(1975) modeled the mixed predictors using the location model (Olkin and Tate,
1961) and proposed an allocation rule to two groups similar to Fisher’s
discriminant function. The location model transforms the $q$ binary variables
$H_{1},\ldots,H_{q}$ to the corresponding $2^{q}$-category multinomial vector
and requires the continuous variables be conditionally normal in each of the
$2^{q}$ categories with different means and same variance-covariance matrix.
He showed that the simple linear discriminant function often gives
satisfactory results, except when there is interaction between the mixed
variables.

We analyze four of the five data sets in Krzanowski’s paper which contains
continuous and binary predictors and a binary response.

  1. 1.

Data Set 1: Ten variables recorded on 40 patients who were surgically treated
for renal hypertension. Seven of the variables were continuous and three
binary. After one year, 20 patients were classified as improved and 20 as
unimproved.

  2. 2.

Data Set 2: Seven variables recorded on 93 patients suffering from jaundice.
Four of the variables were continuous and three binary. The two groups were
patients requiring medical and surgical treatment.

  3. 3.

Data Set 3: Twelve variables recorded on 62 patients suffering from jaundice.
Eight of the variables were continuous and four binary. The two groups were
patients requiring medical and surgical treatment.

  4. 4.

Data Set 4: Eleven variables recorded on 186 patients who underwent ablative
surgery for advanced breast cancer between 1958 and 1965 at Guy’s Hospital,
London. Six of the variables were continuous and three binary. The two groups
were patients for which the treatment was deemed to be successful and failure.

Some of the continuous variables were transformed to normality across all data
sets. Since the response is binary, $\mathbf{f}_{y}$ in (8) is a vector of
frequencies with $r=1$, so that the dimension either SDR method can detect
cannot exceed 1. We reduced the mixed predictors using our two methods, SDR
Optimal and SDR Suboptimal, and also PCA and PCAmix setting $d=1$. In order to
assess the classification accuracy of each method, the reduced predictors
serve as independent variables in a logistic regression model. For comparison,
we also fit an unreduced logistic regression model with all the original
predictors, which we refer to as Full.

In Table 4 we report the leave-one-out misclassification rates and the area
under the receiver operator characteristics curve, AUC (Pepe, 2003, p. 67),
with the smallest and largest values, respectively, in boldface. Sub-optimal
SDR emerges as the best method to summarize the mixed predictors with respect
to misclassification error, followed by SDR Optimal that has better
performance for data set 1. With respect to AUC, SDR Suboptimal is always the
best.

In Table 4, we also provide the leave-one-out misclassification rates of
Fisher’s LDA, logistic regression and Krzanowski’s allocation rule based on
the location model, as reported in Krzanowski (1975, Tab. 3). Sub-optimal SDR
exhibits better performance than Krzanowski’s location model across data sets.
Optimal SDR performs the best in all data sets except for data set 2 where it
is on par with Fisher’s linear discriminant analysis. Moreover, the Optimal
and Sub-optimal SDR misclassification rates are smaller than all other methods
in Krzanowski (1975), as well as mixed nonparametric kernel methods
(Vlachonikolis and Marriott, 1982). Taken all together, our SDR methods for
mixed predictors consistently produce targeted data reductions that provide
better fit and prediction.

Table 4: Leave-one-out misclassification rates and AUC values for four data
sets in Krzanowski (1975).

Set Optimal SubOpt. Full PCA PCAmix Location Fisher Logistic 1 MR 0.250 0.300
0.375 0.325 0.425 0.350 0.325 0.325 AUC 0.918 0.918 0.885 0.675 0.575 - - - 2
MR 0.280 0.204 0.258 0.387 0.290 0.290 0.280 0.301 AUC 0.857 0.858 0.837 0.513
0.469 - - - 3 MR 0.161 0.145 0.226 0.484 0.500 0.226 0.177 0.222 AUC 0.949
0.951 0.944 0.623 0.646 - - - 4 MR 0.296 0.290 0.392 0.457 0.430 0.328 0.382
0.371 AUC 0.784 0.785 0.738 0.544 0.572 - - -

###  7.2 Governance index application

Considerable social science and economics research is devoted to the
construction of indexes for descriptive and predictive purposes (Vyas and
Kumaranayake, 2006; Kolenikov and Angeles, 2009; Filmer and Scott, 2012;
Merola and Baulch, 2014; Forzani et al., 2018). An index is a statistical
summary measure of change in a representative group of individual data points.
It usually synthesizes the information contained in a set of $p$ variables
${\mathbf{X}}\in{}^{p}$ via a linear combination,
$\mathbf{R}({\mathbf{X}})=\bm{\omega}^{T}{\mathbf{X}}\in\real$, where
$\bm{\omega}$ is the vector of weights of the composite index.

In this example, we study the impact of governance on economic growth in the
twelve South American countries as measured by per capita Gross Domestic
Product (GPD) using the World Bank Governance Indicators.111Governance
Indicators and per capita GDP data can be downloaded from Worldwide Governance
Indicators and The World Bank Data, respectively. The World Bank considers the
following six aggregate indicators of governance that combine the views of a
large number of enterprise, citizen and expert survey respondents: control of
corruption ($X_{1}$); rule of law ($X_{2}$); regulatory quality ($X_{3}$);
government effectiveness ($X_{4}$); political stability ($X_{5}$); voice and
accountability ($X_{6}$). They are standardized to have mean zero and standard
deviation one, with values from -2.5 to 2.5, approximately, where higher
values correspond to better governance. All six are highly positively
correlated, and are all positively correlated with the per capita GDP; i.e.,
economic growth is positively associated with better governance indicators.

Our aim is to build a Composite Governance index (CG) to predict $Y$, the
logarithm of per capita Gross Domestic Product (GPD), measured in 2010 US
dollars, over the period 1996 to 2018. Using the set of governance indicator
variables, we start by constructing the CG index via standard Principal
Component Analysis (PCA) and Principal Fitted Components (PFC) [see Corollary
3] setting $d=1$ and $\mathbf{f}_{y}=\log(GDP)$ in (8).

In the left panel of Figure 3, we plot $\log(GDP)$ versus the CG indexes based
on PCA, which is the standard approach in such index construction (Mazziotta
and Pareto, 2019). In the right panel of Figure 3, the response is plotted
versus the index based on PFC. Both plots indicate dependence of the response
on the indexes but the nature of relationship is the data pattern is hard to
understand. A linear trend appears stronger in the right panel, which is
reflected in the better fit of the linear regression model (black) with
$R^{2}=0.27$ versus 0.17 for PCA. However, the PCA-based index in
nonparametric kernel regression (blue) results in better fit. Using the np R
package, the value of the nonparametric version of $R^{2}$ is 0.32 for the
PFC-based CG index, which is much lower than 0.54, the value for the PCA-based
index.

Figure 3: Log of per capita GDP versus Standard PCA and PFC based composite
governance indexes.

In Figure 4, we plot $\log(GDP)$ versus the PCA and PFC composite governance
indexes by country. The plots indicate that the PFC index gives a much better
visualization of the relationship of $\log(GDP)$ within each country,
suggesting that adjusting the index by country could improve its predictive
performance.

Figure 4:  Log of per capita GDP versus standard PCA and standard PFC
composite governance indexes by country.

We add country effect by introducing eleven binary variables $\mathbf{H}$. In
Figure 5 we plot the log of GDP versus the CG index constructed by PCA for
mixed variables (PCAmix) in the left panel and by our mixed optimal SDR
approach in the right panel. Hardly any difference between the plots in the
left panels of Figures 3 and 5 is noticeable. The PCAmix based CG index is
very similar to the conventional PCA based CG that does not include country
effect, with $R^{2}$ equal to 0.17 and 0.61 for the linear and nonparametric
models, respectively. Moreover, neither PCA based CG index exhibits an easy to
understand or model relationship with the response.

In contrast, a very clear and simple pattern appears in the right panel of
Figure 5, where the response is plotted versus our optimal SDR based index.
The pattern suggests modeling $\log(GDP)$ as a linear function of the GC
index. This is a distinct improvement over PCA and PCAmix (left panels of
Figures 3 and 5) but also the SDR method PFC, which does not account for
country effect (right panel of Figure 3). As a result, both the linear (black)
and the kernel (blue) regression models for the regression of the log per
capita GDP on the optimal SDR for mixed predictors based CG index have
excellent fit with respective $R^{2}$ values of 0.91 and 0.93.

Figure 5: Log of per capita GDP versus Composite Governance index with country
effect.

The average of the leave-one-out mean square prediction errors of the linear
and kernel regression models in Table 5, provides an unbiased measure of
predictive performance. The logarithm of the per capita GDP is regressed on
the unsupervised CG indexes, constructed by PCA using only continuous
predictors ($\text{PCA}({\mathbf{X}})$) and its extension for mixed variables
($\text{PCAmix}({\mathbf{X}},\mathbf{H})$), and the supervised CG Indexes,
constructed by PFC only on continuous predictors ($\text{PFC}({\mathbf{X}}$))
and our mixed predictor SDR methods, $\text{\sc
Optimal}({\mathbf{X}},\mathbf{H})$ and $\text{\sc
SubOptimal}({\mathbf{X}},\mathbf{H})$.

The leave-one-out mean squared prediction errors of the supervised PFC based
CG index are smaller than both PCA and PCAmix for the linear model, even
though PFC does not account for country effect. Nevertheless, when the kernel
regression model is fitted, the PCA based index exhibits better performance
than PFC. The dramatic drop in prediction error results from using optimal and
sub-optimal SDR, as it is between 5 to 9 times smaller than the PCA, PCAmix
and PFC errors for both the linear regression and the kernel regression
models.

Table 5: Leave-one-out mean squared prediction errors for the per capita log
GDP in South-American countries.

Predictive Model Index Type Method Linear Non-Parametric Unsupervised
$\text{PCA}({\mathbf{X}})$ 0.319 0.189
$\text{PCAmix}({\mathbf{X}},\mathbf{H})$ 0.320 0.209 Supervised
$\text{PFC}({\mathbf{X}})$ 0.292 0.282 $\text{\sc
SDROptimal}({\mathbf{X}},\mathbf{H})$ 0.029 0.028 $\text{\sc
SDRSubOptimal}({\mathbf{X}},\mathbf{H})$ 0.028 0.022

The regularized estimation of the SDROptimal reduction selects all five
continuous predictors except for rule of law. Political stability and voice
and accountability have the highest weights in the CG index. Rule of law is
the most correlated with four of the other variables, with correlation
coefficient values over 0.80. We stipulate that our method drops it as its
relationship with GDP is mostly absorbed by the other four. The binary
variables are all selected. That is, our method finds a significant country
effect on GDP.

##  8 Discussion

Our approach falls within model-based inverse regression for sufficient
dimension reduction (SDR) (Cook (2007); Cook and Forzani (2008); Bura and
Forzani (2015); Bura et al. (2016)). Model-based SDR requires knowledge of the
family of distributions of the inverse predictors in contrast to moment-based
SDR, such as SIR Li (1991), SAVE Cook and Weisberg (1991), or DR Li and Wang
(2007), that impose conditions on the moments of the marginal distribution of
the predictors. Because of this, our approach provides exhaustive
identification and statistically efficient estimation of sufficient reductions
for the conditional distribution of an output given mixed variables that
contain all information in the mixed predictors for the output $Y$.

Furthermore, outside the context of dimension reduction for the forward
regression problem of $Y$ on mixed predictors ${\mathbf{Z}}$, the modeling we
use to accommodate the factorization in (2) in developing our SDR methods, is
a new multivariate modeling approach for response vectors comprised of mixed
variables. That is, if one were to only consider the multivariate regression
of the mixed vector ${\mathbf{Z}}=({\mathbf{X}}^{T},\mathbf{H}^{T})^{T}$ on
some other variables, say $\mathbf{F}$, the models we use for the continuous
and binary elements of ${\mathbf{Z}}$ in our development provides a new
regression tool for mixed responses. Specifically, since the joint
distribution of ${\mathbf{Z}}\mid\mathbf{F}$ belongs to the exponential family
(9), our approach yields sufficient statistics for the unknown natural
parameters $\bm{\vartheta}$ in (26), as well as optimal (efficient) maximum
likelihood estimators, in a similar manner to generalized linear modeling for
univariate responses.

Acknowledgments

EB would like to acknowledge support for this project from the Austrian
Science Fund (FWF P 30690-N35) and the Vienna Science and Technology Fund
(WWTF ICT19-018).

## Appendix A. Proofs and Derivations for Section 3

### Derivation of Eqn. (9)

From Eqn. (8), the density $f({\mathbf{X}},\mathbf{H}|Y=y)$, up to the
constant $1/\sqrt{2\pi}$, equals

$\displaystyle\exp\bigg{\\{}-\frac{1}{2}\Big{(}({\mathbf{X}}-\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}})-\mathbf{A}\mathbf{f}_{y}-\bm{\beta}(\mathbf{H}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}})\Big{)}^{T}\bm{\Delta}^{-1}\Big{(}({\mathbf{X}}-\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}})-\mathbf{A}\mathbf{f}_{y}-\bm{\beta}(\mathbf{H}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}})\Big{)}$
$\displaystyle\hskip
113.81102pt+\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})\left(\bm{\tau}_{0}+\bm{\tau}\mathbf{f}_{y}\right)+\frac{1}{2}\log(|\bm{\Delta}|^{-1})-\log(G(\bm{\Gamma}_{y}))\bigg{\\}}.$

After some algebra and rearrangement of terms we obtain

$f({\mathbf{X}},\mathbf{H}|Y=y)=h({\mathbf{X}},\mathbf{H})\exp\left(\mathbf{T}^{T}({\mathbf{X}},\mathbf{H})\bm{\eta}_{y}-\psi(\bm{\eta}_{y})\right),$

with $h({\mathbf{X}},\mathbf{H})=(2\pi)^{-1/2}$,

$\displaystyle\mathbf{T}^{T}({\mathbf{X}},\mathbf{H})\bm{\eta}_{y}$
$\displaystyle={\mathbf{X}}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}-{\mathbf{X}}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}+{\mathbf{X}}^{T}\bm{\Delta}^{-1}\mathbf{A}\mathbf{f}_{y}$
$\displaystyle\hskip
28.45274pt-\mathbf{H}^{T}\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\mathbf{H}^{T}\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}-\mathbf{H}^{T}\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A}\mathbf{f}_{y}$
$\displaystyle\hskip
28.45274pt-\frac{1}{2}{\mathbf{X}}^{T}\bm{\Delta}^{-1}{\mathbf{X}}+{\mathbf{X}}^{T}\bm{\Delta}^{-1}\bm{\beta}\mathbf{H}$
$\displaystyle\hskip
28.45274pt-\frac{1}{2}\mathbf{H}^{T}\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\mathbf{H}+\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})\bm{\tau}_{0}+\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})\bm{\tau}\mathbf{f}_{y},$
(82)

and

$\displaystyle\psi(\bm{\eta}_{y})$
$\displaystyle=\frac{1}{2}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\frac{1}{2}\mathbf{f}_{y}^{T}\mathbf{A}^{T}\bm{\Delta}^{-1}\mathbf{A}\mathbf{f}_{y}+\frac{1}{2}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}^{T}\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}$
$\displaystyle\hskip
28.45274pt+\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}^{T}\bm{\Delta}^{-1}\mathbf{A}\mathbf{f}_{y}-\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}-\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}^{T}\beta^{T}\bm{\Delta}^{-1}\mathbf{A}\mathbf{f}_{y}$
(83) $\displaystyle\hskip
28.45274pt-\frac{1}{2}\log(|\bm{\Delta}|^{-1})+\log(G(\bm{\Gamma}_{y})).$

Since
$\mathrm{\,tr}(\mathbf{A}^{T}\mathbf{B})=\mathrm{vec}(\mathbf{A})^{T}\mathrm{vec}(\mathbf{B})$,
$\mathrm{vec}(\mathbf{A}\mathbf{B}\mathbf{C})=(\mathbf{C}^{T}\otimes\mathbf{A})\mathrm{vec}(\mathbf{B})$
and $\mathbf{D}_{q}$ in Section 2 is such that
$\mathrm{vec}(\mathbf{A})=\mathbf{D}_{q}\mathrm{\,vech}(\mathbf{A})$,
(Derivation of Eqn. (9)) becomes

$\displaystyle\mathbf{T}^{T}({\mathbf{X}},\mathbf{H})\bm{\eta}_{y}$
$\displaystyle={\mathbf{X}}^{T}\left(\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}-\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{p})\mathrm{vec}(\bm{\Delta}^{-1}\mathbf{A})\right)$
$\displaystyle\hskip
8.5359pt+\mathbf{H}^{T}\left(-\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}-(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q})\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A})\right)$
$\displaystyle\hskip
8.5359pt-\frac{1}{2}(\mathbf{D}_{p}\mathbf{D}_{p}^{T}\mathrm{\,vech}({\mathbf{X}}{\mathbf{X}}^{T}))^{T}\mathrm{\,vech}(\bm{\Delta}^{-1})+\mathrm{vec}({\mathbf{X}}\mathbf{H}^{T})^{T}\mathrm{vec}(\bm{\Delta}^{-1}\bm{\beta})$
$\displaystyle\hskip
8.5359pt+\mathrm{\,vech}(\mathbf{H}\mathbf{H}^{T})^{T}\left(-\frac{1}{2}\mathbf{D}_{q}^{T}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})+\bm{\tau}_{0}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q(q+1)/2})\mathrm{vec}(\bm{\tau})\right).$

Finally, using the matrices $\mathbf{J}_{q}$ and $\mathbf{L}_{q}$ defined in
Section 2, we obtain Eqns (10) and (26) from

$\displaystyle\mathbf{T}^{T}({\mathbf{X}},\mathbf{H})\bm{\eta}_{y}$
$\displaystyle={\mathbf{X}}^{T}\left(\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}-\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{p})\mathrm{vec}(\bm{\Delta}^{-1}\mathbf{A})\right)$
$\displaystyle\hskip
8.5359pt+\mathbf{H}^{T}\left(-\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\textsc{\tiny$\mathbf{H}$}}-(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q})\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A})\right)$
$\displaystyle\quad-\frac{1}{2}(\mathbf{D}_{p}\mathbf{D}_{p}^{T}\mathrm{\,vech}({\mathbf{X}}{\mathbf{X}}^{T}))^{T}\mathrm{\,vech}(\bm{\Delta}^{-1})+\mathrm{vec}({\mathbf{X}}\mathbf{H}^{T})^{T}\mathrm{vec}(\bm{\Delta}^{-1}\bm{\beta})$
$\displaystyle\quad+\mathbf{H}^{T}\left(-\frac{1}{2}\mathbf{L}_{q}\mathbf{D}_{q}^{T}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})+\mathbf{L}_{q}\bm{\tau}_{0}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q})\mathrm{vec}(\mathbf{L}_{q}\bm{\tau})\right)$
$\displaystyle\quad+(\mathbf{J}_{q}\mathrm{\,vech}(\mathbf{H}\mathbf{H}^{T}))^{T}\left(-\frac{1}{2}\mathbf{J}_{q}\mathbf{D}_{q}^{T}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})+\mathbf{J}_{q}\bm{\tau}_{0}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{k_{q}})\mathrm{vec}(\mathbf{J}_{q}\bm{\tau})\right)$
$\displaystyle={\mathbf{X}}^{T}\bm{\eta}_{y1}+\mathbf{H}^{T}\bm{\eta}_{y2}-\frac{1}{2}(\mathbf{D}_{p}^{T}\mathbf{D}_{p}\mathrm{vec}({\mathbf{X}}{\mathbf{X}}^{T}))^{T}\bm{\eta}_{3}+\mathrm{vec}({\mathbf{X}}\mathbf{H}^{T})^{T}\bm{\eta}_{4}$
$\displaystyle\qquad+(\mathbf{J}_{q}\mathrm{\,vech}(\mathbf{H}\mathbf{H}^{T}))^{T}\bm{\eta}_{y5},$

where $\mathbf{T}({\mathbf{X}},\mathbf{H})$ is defined in (10) and

$\bm{\eta}_{y1}=\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}-\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\mathbf{H}}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{p})\mathrm{vec}(\bm{\Delta}^{-1}\mathbf{A})=\mathbf{F}_{y1}\bm{\vartheta}_{1},$

with
$\mathbf{F}_{y1}=(\mathbf{I}_{p},\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{p})$,
$\bm{\vartheta}_{1}=(\bm{\vartheta}_{10}^{T},\bm{\vartheta}_{11}^{T})^{T}$,
$\bm{\vartheta}_{10}=\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}-\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\mathbf{H}}$,
$\bm{\vartheta}_{11}=\mathrm{vec}(\bm{\Delta}^{-1}\mathbf{A})$,

$\displaystyle\bm{\eta}_{y2}$
$\displaystyle=-\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\mathbf{H}}-(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q})\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A})$
$\displaystyle\qquad-\frac{1}{2}\mathbf{L}_{q}{{\mathbf{D}}}^{T}_{q}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})+\mathbf{L}_{q}\bm{\tau}_{0}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q})\mathrm{vec}(\mathbf{L}_{q}\bm{\tau})$
$\displaystyle=\mathbf{F}_{y2}\bm{\vartheta}_{2},$

with
$\mathbf{F}_{y2}=(\mathbf{I}_{q},\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q})$,
$\bm{\vartheta}_{2}=(\bm{\vartheta}_{20}^{T},\bm{\vartheta}_{21}^{T})^{T}$,
$\bm{\vartheta}_{20}=-\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\mu}_{\textsc{\tiny${\mathbf{X}}$}}+\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta}\bm{\mu}_{\mathbf{H}}+\mathbf{L}_{q}\bm{\tau}_{0}-\frac{1}{2}\mathbf{L}_{q}{{\mathbf{D}_{q}^{T}}}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})$,
$\bm{\vartheta}_{21}=\mathrm{vec}(\mathbf{L}_{q}\bm{\tau}-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A})$,

$\bm{\eta}_{3}=\bm{\eta}_{3y}=\mathrm{\,vech}(\bm{\Delta}^{-1}),$
$\bm{\eta}_{4}=\bm{\eta}_{4y}=\mathrm{vec}(\bm{\Delta}^{-1}\bm{\beta}),$

and

$\displaystyle\bm{\eta}_{y5}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\mathbf{J}_{q}\mathbf{D}_{q}^{T}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})+\mathbf{J}_{q}\bm{\tau}_{0}+(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{k_{q}})\mathrm{vec}(\mathbf{J}_{q}\bm{\tau})$
$\displaystyle=$ $\displaystyle\mathbf{F}_{y5}\bm{\vartheta}_{5},$

with
$\mathbf{F}_{y5}=(\mathbf{I}_{k_{q}},\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{k_{q}})$,
$\bm{\vartheta}_{5}=(\bm{\vartheta}_{50}^{T},\bm{\vartheta}_{51}^{T})^{T}$,
$\bm{\vartheta}_{50}=-{{\frac{1}{2}\mathbf{J}_{q}\mathbf{D}^{T}_{q}}}\mathrm{vec}(\bm{\beta}^{T}\bm{\Delta}^{-1}\bm{\beta})+\mathbf{J}_{q}\bm{\tau}_{0}$
and $\bm{\vartheta}_{51}=\mathrm{vec}(\mathbf{J}_{q}\bm{\tau})$.

By Eqn. (7),

$G(\bm{\Gamma}_{y})=\sum_{H}\exp\big{[}\mathrm{\,vech}^{T}(\mathbf{H}\mathbf{H}^{T})\left(\bm{\tau}_{0}+\bm{\tau}\mathbf{f}_{y}\right)\big{]}.$

Plugging in matrices $\mathbf{J}_{q}$, $\mathbf{L}_{q}$, $\mathbf{D}_{q}$ and
$\mathbf{C}_{q}$, defined in Section 2, and letting
$\bar{\bm{\eta}}_{4}=\mathrm{unvec}({\bm{\eta}}_{4})$, we obtain

$\displaystyle G(\bm{\Gamma}_{y})$ $\displaystyle=$
$\displaystyle\sum_{H}\exp\left[(\mathbf{J}_{q}\mathbf{C}_{q}\mathrm{vec}(\mathbf{H}\mathbf{H}^{T}))^{T}\left(\bm{\eta}_{y5}+\frac{1}{2}\mathbf{J}_{q}{{\mathbf{D}}}^{T}_{q}{\mathrm{vec}}(\bar{\bm{\eta}}_{4}^{T}(\mathrm{unvec}(\mathbf{D}_{p}\bm{\eta}_{3}))^{-1}\bar{\bm{\eta}}_{4}))\right.\right)$
$\displaystyle\left.+\mathbf{H}^{T}\left(\bm{\eta}_{y2}+\bar{\bm{\eta}}_{4}^{T}(\mathrm{unvec}(\mathbf{D}_{p}\bm{\eta}_{3}))^{-1}\bm{\eta}_{y1}+\frac{1}{2}\mathbf{L}_{q}{{\mathbf{D}}}_{q}^{T}\mathrm{vec}(\bar{\bm{\eta}}_{4}^{T}(\mathrm{unvec}(\mathbf{D}_{p}\bm{\eta}_{3}))^{-1}\bar{\bm{\eta}}_{4}))\right)\right].$

Finally, using the matrix $\mathbf{D}_{p}$ defined in Section 2, Eqn.
(Derivation of Eqn. (9)) yields

$\displaystyle\psi(\bm{\eta}_{y})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bm{\eta}_{y1}^{T}(\mathrm{unvec}(\mathbf{D}_{p}\bm{\eta}_{3}))^{-1}\bm{\eta}_{y1}+\log
G(\bm{\Gamma}_{y})-\frac{1}{2}\log|\mathrm{unvec}(\mathbf{D}_{p}\bm{\eta}_{3})|$
(85) $\displaystyle=$
$\displaystyle\psi_{1}(\bm{\eta}_{y})+\psi_{2}(\bm{\eta}_{y})+\psi_{3}(\bm{\eta}_{y}),$

with $G(\bm{\Gamma}_{y})$ given in (Derivation of Eqn. (9)).

### Proof of Theorem 1

Since the density of ${\mathbf{X}},\mathbf{H}\mid Y$ belongs to the full rank
exponential family (Eqn. (9)), the minimal sufficient reduction for the
regression $Y\mid({\mathbf{X}},\mathbf{H})$ is given by

$\displaystyle\mathbf{R}({\mathbf{X}},\mathbf{H})$
$\displaystyle=\bm{\alpha}_{\mathbf{a}}^{T}\left(\mathbf{T}({\mathbf{X}},\mathbf{H})-\mathrm{E}(\mathbf{T}({\mathbf{X}},\mathbf{H}))\right),$

where $\bm{\alpha}_{\mathbf{a}}$ is a basis for
${\mathcal{S}}_{\bm{\alpha}_{\mathbf{a}}}=\mathrm{span}\\{\bm{\eta}_{Y}-\mathrm{E}({\bm{\eta}_{Y}}),Y\in\mathcal{Y}\\}$,
with $\bm{\eta}_{Y}$ given in (26) [see Bura et al. (2016, Thm 1)]. Since
$\mathrm{E}(\mathbf{f}_{Y})=\bm{0}$, applying Eqns. (26) and (28) obtains

$\bm{\eta}_{y}-\mathrm{E}(\bm{\eta}_{y})=\left(\begin{array}[]{c}(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{p})\mathrm{vec}(\bm{\Delta}^{-1}\mathbf{A})\\\
(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{q})\mathrm{vec}(\mathbf{L}_{q}\bm{\tau}-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A})\\\
0\\\ 0\\\
(\mathbf{f}_{y}^{T}\otimes\mathbf{I}_{k_{q}})\mathrm{vec}(\mathbf{J}_{q}\bm{\tau})\end{array}\right)=\left(\begin{array}[]{c}\bm{\Delta}^{-1}\mathbf{A}\mathbf{f}_{y}\\\
(\mathbf{L}_{q}\bm{\tau}-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A})\mathbf{f}_{y}\\\
0\\\ 0\\\ \mathbf{J}_{q}\bm{\tau}\mathbf{f}_{y}\end{array}\right).$

Then,
$\mathrm{span}(\mathbf{a})=\mathrm{span}\left(\bm{\eta}_{y}-\mathrm{E}({\bm{\eta}_{Y}}),y\in\mathcal{Y}\right)$
with

$\mathbf{a}=\left(\begin{array}[]{c}\bm{\Delta}^{-1}\mathbf{A}\\\
\mathbf{L}_{q}\bm{\tau}-\bm{\beta}^{T}\bm{\Delta}^{-1}A\\\ 0\\\ 0\\\
\mathbf{J}_{q}\bm{\tau}\end{array}\right).$

### Proof of Corollary 3

If follows from Corollary 2 since, in this case, $\bm{\vartheta}_{2,1}=0$ and
$\bm{\vartheta}_{5,1}=0$.

### Proof of Corollary 4

If follows from Corollary 2 since in this case $\bm{\vartheta}_{1,1}=0$ and
$\bm{\vartheta}_{2,1}=\mathbf{L}_{q}\bm{\tau}$.

### Proof of Corollary 5

It suffices to show that
$\mathrm{span}({\mathbf{b}})\subset\mathrm{span}(\bm{\alpha}_{\mathbf{c}})$.
We can write b as

${\mathbf{b}}=\left(\begin{array}[]{cc}\bm{\Delta}\mathbf{A}&\mathbf{0}\\\
-\bm{\beta}^{T}\bm{\Delta}^{-1}\mathbf{A}&\mathbf{L}_{1}\bm{\tau}\\\
\mathbf{0}&\mathbf{J}_{q}\bm{\tau}\end{array}\right)\left(\begin{array}[]{c}\mathbf{I}_{r}\\\
\mathbf{I}_{r}\end{array}\right)=\mathbf{\widetilde{b}}\left(\begin{array}[]{c}\mathbf{I}_{r}\\\
\mathbf{I}_{r}\end{array}\right),$

with
$\mathrm{span}(\mathbf{\widetilde{b}})=\mathrm{span}(\bm{\alpha}_{\mathbf{c}}).$
As a consequence,
$\mathrm{span}({\mathbf{b}})\subset\mathrm{span}(\bm{\alpha}_{\mathbf{c}})$,
and therefore $\mathbf{R}({\mathbf{X}},\mathbf{H})$ in Eqn. (47) is a
sufficient dimension reduction, not necessary minimal. The rest of the
corollary immediately follows.

## Appendix B: Proof of Proposition 6

We first derive the asymptotic distribution of $\widehat{\mathbf{b}}$ in (67)
and prove auxiliary lemmas in Section in order to prove Proposition 6.

### Asymptotic distribution of $\widehat{\mathbf{b}}$

###### Proposition 7

$\sqrt{n}\mathrm{vec}(\widehat{\mathbf{b}}-\mathbf{b})\xrightarrow{\mathcal{D}}\mathcal{N}(\mathbf{0},\mathbf{V}_{rcl})$
with

$\mathbf{V}_{rcl}=\mathbf{W}\mathbf{M}\mathbf{V}\mathbf{M}^{T}\mathbf{W}^{T},$

as in equation (72), where $\mathbf{M}$, $\mathbf{W}$ and $\mathbf{V}$ are
defined in Eqns. (75), (76) and (73), respectively.

Proof  We rewrite

$\displaystyle\mathbf{b}$
$\displaystyle=\begin{pmatrix}{\bm{\Delta}}^{-1}{\mathbf{A}}\\\
\mathbf{L}_{q}{\bm{\tau}}-{\bm{\beta}}^{T}{\bm{\Delta}}^{-1}{\mathbf{A}}\\\
\mathbf{J}_{q}{\bm{\tau}}\end{pmatrix}=\begin{pmatrix}\mathrm{unvec}({\bm{\vartheta}}_{1,1})\\\
\mathrm{unvec}({\bm{\vartheta}}_{2,1})\\\
\mathrm{unvec}({\bm{\vartheta}}_{5,1})\end{pmatrix},$

as follows. Let
$\widetilde{\mathbf{b}}=\left(\bm{\vartheta}_{1,1}^{T},\bm{\vartheta}_{2,1}^{T},\bm{\vartheta}_{5,1}^{T}\right)^{T}$.
Then, $\widetilde{\mathbf{b}}=\mathbf{M}\bm{\vartheta}$, with $\mathbf{M}$
given in (75), so that

$\mathrm{vec}(\mathbf{b})=\mathbf{W}\mathbf{\widetilde{b}}=\mathbf{W}\mathbf{M}\bm{\vartheta},$
(86)

with $\mathbf{W}$ defined on (76). Then,

$\mathrm{vec}(\widehat{\mathbf{b}})=\mathbf{W}\mathbf{M}\widehat{\bm{\vartheta}}.$
(87)

The asymptotic normality of $\widehat{\mathbf{b}}$ follows from the asymptotic
normality of $\widehat{\bm{\vartheta}}$, which is derived in Lemma 8, with

$\displaystyle\mathrm{avar}(\sqrt{n}\widehat{\mathbf{b}})$
$\displaystyle=\mathbf{W}\mathbf{M}\mathrm{avar}\left(\sqrt{n}\widehat{\bm{\vartheta}}\right)\mathbf{M}^{T}\mathbf{W}^{T}$
$\displaystyle=\mathbf{W}\mathbf{M}\mathbf{V}\mathbf{M}^{T}\mathbf{W}^{T}.$

###### Lemma 8

If $\mathrm{avar}\left(\sqrt{n}\widehat{\bm{\vartheta}}\right)=\mathbf{V}$,
then

$\mathbf{V}^{-1}=E\left[{\mathbf{F}}_{y}^{T}\mathbf{J}{\mathbf{F}}_{y}\right],$

as in (73), where $\mathbf{F}_{y}$ is defined in (26) and

$\mathbf{J}=\frac{\partial^{2}\psi(\bm{\eta}_{y})}{\partial\bm{\eta}_{y}\partial\bm{\eta}_{y}^{T}}$

in (74).

Proof  Since $\widehat{\bm{\vartheta}}$ is the maximum likehood estimator,

$\mathbf{V}=\mathrm{avar}\left(\sqrt{n}\widehat{\bm{\vartheta}}\right)=-\left(E\left[\frac{\partial^{2}\log
f({\mathbf{X}},\mathbf{H}\mid
Y=y)}{\partial\bm{\theta}\partial\bm{\theta}^{T}}\right]\right)^{-1}$

We plug in $\bm{\eta}_{y}=\mathbf{F}_{y}\bm{\vartheta}$ (from Eqn. (26)) in
Eqn. (9) to obtain

$\displaystyle\log f({\mathbf{X}},\mathbf{H}\mid Y=y)$ $\displaystyle=\log
h({\mathbf{X}},\mathbf{H})+{\mathbf{T}}^{T}({\mathbf{X}},\mathbf{H})\bm{\eta}_{y}-\psi(\bm{\eta}_{y})$
$\displaystyle=\log
h({\mathbf{X}},\mathbf{H})+{\mathbf{T}}^{T}({\mathbf{X}},\mathbf{H}){\mathbf{F}}_{y}\bm{\vartheta}-\psi({\mathbf{F}}_{y}\bm{\vartheta}).$

Then,

$\displaystyle\frac{\partial\log f({\mathbf{X}},\mathbf{H}\mid
Y=y)}{\partial\mathrm{vec}^{T}(\bm{\theta})}$
$\displaystyle={\mathbf{T}}^{T}({\mathbf{X}},\mathbf{H}){\mathbf{F}}_{y}-\frac{\partial\psi(\bm{\eta}_{y})}{\partial\bm{\eta}_{y}^{T}}{\mathbf{F}}_{y},$
$\displaystyle\frac{\partial^{2}\log f({\mathbf{X}},\mathbf{H}\mid
Y=y)}{\partial\mathrm{vec}(\bm{\theta})\mathrm{vec}^{T}(\partial\bm{\theta})}$
$\displaystyle=-{\mathbf{F}}_{y}^{T}\frac{\partial^{2}\psi(\bm{\eta}_{y})}{\partial\bm{\eta}_{y}\partial\bm{\eta}_{y}^{T}}{\mathbf{F}}_{y}=-{\mathbf{F}}_{y}^{T}{\mathbf{J}}{\mathbf{F}}_{y}.$

Therefore
$\mathbf{V}^{-1}=\mathrm{E}\left[\mathbf{F}_{y}^{T}\mathbf{J}\mathbf{F}_{y}\right]$
from which Proposition 8 follows.

In order to compute $\mathbf{J}$, the first and second derivatives of
$\psi(\bm{\eta}_{y})$ with respect to $\bm{\eta}_{y}$ are required. The
computation is carried out in Section LABEL:supplement (Supplementary
Material).

### Auxiliary lemmas for Proposition 6

###### Lemma 9

Let
$\widehat{\mathbf{H}}=\widehat{\mathbf{U}}_{1}\widehat{{\mathbf{K}}}_{1}\widehat{\mathbf{R}}_{1}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1}$.
Then,

$\sqrt{n}\mathrm{vec}(\widehat{\mathbf{H}}-\mathbf{U}_{1})\rightarrow\mathcal{N}(0,({\mathbf{K}}^{-1}\mathbf{R}_{1}^{T}\otimes\mathbf{I}_{m}){\mathbf{V}}_{rlc}(\mathbf{R}_{1}{\mathbf{K}}^{-1}\otimes\mathbf{I}_{m})).$

where $\mathbf{V}_{rlc}$ is defined in Eqn. (72),
$\widehat{\mathbf{U}}_{1},\widehat{\mathbf{K}}_{1}$ and
$\widehat{\mathbf{R}}_{1}$ in Eqn. (40), and ${\mathbf{U}}_{1},{\mathbf{K}}$
and $\mathbf{R}_{1}$ in Eqn. (65).

Proof  By Eqn. (65),
$\mathbf{b}={\mathbf{U}}_{1}\mathbf{K}{\mathbf{R}}_{1}^{T}$ and by Eqn. (40),
$\widehat{\mathbf{b}}=\widehat{\mathbf{U}}_{1}\widehat{{\mathbf{K}}}_{1}\widehat{\mathbf{R}}_{1}^{T}+\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}$.
Then,

$\displaystyle\widehat{\mathbf{H}}-\mathbf{U}_{1}$
$\displaystyle=\widehat{\mathbf{U}}_{1}\widehat{{\mathbf{K}}}_{1}\widehat{\mathbf{R}}_{1}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1}-\mathbf{U}_{1}$
$\displaystyle=\widehat{\mathbf{b}}\mathbf{R}_{1}{\mathbf{K}}^{-1}-\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1}-\mathbf{U}_{1}$
$\displaystyle=(\widehat{\mathbf{b}}-\mathbf{b})\mathbf{R}_{1}{\mathbf{K}}^{-1}-\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1}.$

Thus,

$\displaystyle\sqrt{n}\mathrm{vec}(\widehat{\mathbf{H}}-\mathbf{U}_{1})$
$\displaystyle=\sqrt{n}\mathrm{vec}((\widehat{\mathbf{b}}-\mathbf{b})\mathbf{R}_{1}{\mathbf{K}}^{-1}-\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1})$
$\displaystyle=\sqrt{n}(\mathbf{K}^{-1}\mathbf{R}_{1}^{T}\otimes\mathbf{I}_{m})\mathrm{vec}(\widehat{\mathbf{b}}-\mathbf{b})-\mathrm{vec}(\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1}).$
(88)

From Proposition 7 we have

$\sqrt{n}\mathrm{vec}(\widehat{\mathbf{b}}-\mathbf{b})\xrightarrow{\mathcal{D}}\mathcal{N}(\mathbf{0},\mathbf{V}_{rcl}),$

so that

$\sqrt{n}(\mathbf{K}^{-1}\mathbf{R}_{1}^{T}\otimes\mathbf{I}_{m})\mathrm{vec}(\widehat{\mathbf{b}}-\mathbf{b})\xrightarrow{\mathcal{D}}\mathcal{N}(\mathbf{0},\bm{\Sigma}_{\mathbf{U}}),$
(89)

with
$\bm{\Sigma}_{\mathbf{U}}=({\mathbf{K}}^{-1}\mathbf{R}_{1}^{T}\otimes\mathbf{I}_{k})\mathbf{V}_{rcl}({\mathbf{K}}^{-1}\mathbf{R}_{1}^{T}\otimes\mathbf{I}_{k})^{T}=({\mathbf{K}}^{-1}\mathbf{R}_{1}^{T}\otimes\mathbf{I}_{k})\mathbf{V}_{rcl}(\mathbf{R}_{1}{\mathbf{K}}^{-1}\otimes\mathbf{I}_{k})$.

Since
$\sqrt{n}\left(\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\right)=O_{p}(1)$
and
$\mathbf{P}_{\mathbf{R}_{1}}=\left(\mathbf{P}_{\widehat{\mathbf{R}}_{1}}+O_{p}(n^{-1/2})\right)$,
we get

$\displaystyle\sqrt{n}(\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1})$
$\displaystyle=\sqrt{n}(\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T})\mathbf{P}_{\mathbf{R}_{1}}\mathbf{R}_{1}{\mathbf{K}}^{-1}$
$\displaystyle=\sqrt{n}\left(\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\right)\left(\mathbf{P}_{\widehat{\mathbf{R}}_{1}}+O_{p}(n^{-1/2})\right)\mathbf{R}_{1}{\mathbf{K}}^{-1}$
$\displaystyle=\sqrt{n}\left(\widehat{\mathbf{U}}_{0}\widehat{{\mathbf{K}}}_{0}\widehat{\mathbf{R}}_{0}^{T}\right)O_{p}(n^{-1/2})\mathbf{R}_{1}{\mathbf{K}}^{-1}$
$\displaystyle=O_{p}(n^{-1/2}),$

where we use
$\widehat{\mathbf{R}}_{0}^{T}\widehat{\mathbf{R}}_{1}=\mathbf{0}$. As a
consequence,
$\sqrt{n}\mathrm{vec}(\widehat{\mathbf{U}}_{0}\widehat{\mathbf{K}}_{0}\widehat{\mathbf{R}}_{0}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1})\rightarrow\mathbf{0}$
in probability, which, together with (89) in (Auxiliary lemmas for Proposition
6), obtain the result.

###### Lemma 10

Let $\bm{\Gamma}$ be a matrix of dimension $p\times d$ of full rank $d$ with
$d\leq p$, and let $\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}$ be the orthogonal
projection onto the column space of $\bm{\Gamma}$ and
$\mathbf{Q}_{\bm{{\scriptstyle\Gamma}}}=\mathbf{I}-\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}$.
Also, let $\mathbf{K}_{pm}\in{\mathbb{R}}^{pm\times pm}$ be the unique matrix
such that, for any symmetric $p\times m$ matrix $\mathbf{A}$,
$\mathrm{vec}(\mathbf{A}^{T})=\mathbf{K}_{pm}\mathrm{vec}(\mathbf{A})$. Then,

$\frac{\partial{\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}}}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}=(\mathbf{I}_{p^{2}}+\mathbf{K}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{Q}_{\bm{{\scriptstyle\Gamma}}}).$
(90)

Proof  We will use the following two identities.

  * (i)

Let ${\mathbf{X}}$ be a matrix and $\mathbf{F}({\mathbf{X}}):m\times p$ and
$\mathbf{G}({\mathbf{X}}):p\times q$ differentiable matrix valued functions of
${\mathbf{X}}$. Then,

$\frac{\partial\mathrm{vec}[\mathbf{F}({\mathbf{X}})\mathbf{G}({\mathbf{X}})]}{\partial\mathrm{vec}^{T}({\mathbf{X}})}=(\mathbf{G}^{T}\otimes\mathbf{I}_{m})\frac{\partial\mathrm{vec}[\mathbf{F}({\mathbf{X}})]}{\partial\mathrm{vec}^{T}({\mathbf{X}})}+(\mathbf{I}_{q}\otimes\mathbf{F})\frac{\partial\mathrm{vec}[\mathbf{G}({\mathbf{X}})]}{\partial\mathrm{vec}^{T}({\mathbf{X}})}.$
(91)

  * (ii)

Let $\mathbf{F}({\mathbf{X}})={\mathbf{X}}^{T}$ and
$\mathbf{G}({\mathbf{X}})={\mathbf{X}}$ with ${\mathbf{X}}:p\times q$. By
(91),

$\displaystyle\frac{\partial\mathrm{vec}({\mathbf{X}}^{T}{\mathbf{X}})}{\partial^{T}\mathrm{vec}({\mathbf{X}})}$
$\displaystyle=$
$\displaystyle(\mathbf{I}_{q^{2}}+\mathbf{K}_{qq})(\mathbf{I}_{q}\otimes{\mathbf{X}}^{T}).$
(92)
$\displaystyle\frac{\partial\mathrm{vec}({\mathbf{X}}^{T}{\mathbf{X}})^{-1}}{\partial^{T}\mathrm{vec}({\mathbf{X}})}$
$\displaystyle=$
$\displaystyle-(({\mathbf{X}}^{T}{\mathbf{X}})^{-1}\otimes({\mathbf{X}}^{T}{\mathbf{X}})^{-1}))\frac{\partial\mathrm{vec}({\mathbf{X}}^{T}{\mathbf{X}})}{\partial^{T}\mathrm{vec}({\mathbf{X}})}.$

Applying (91) yields

$\displaystyle\frac{\partial\mathrm{vec}\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}$
$\displaystyle=\frac{\partial\mathrm{vec}(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T})}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}$
$\displaystyle=(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{I}_{p})\frac{\partial\mathrm{vec}(\bm{\Gamma})}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}+(\mathbf{I}_{p}\otimes\bm{\Gamma})\frac{\partial\mathrm{vec}((\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T})}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}$
$\displaystyle=(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{I}_{p})+(\mathbf{I}_{p}\otimes\bm{\Gamma})\frac{\partial\mathrm{vec}((\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T})}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}.$

Let

$\mathbf{H}=\frac{\partial\mathrm{vec}((\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T})}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}.$

Using (91), (92) and ((ii)), we get

$\displaystyle\mathbf{H}$ $\displaystyle=$
$\displaystyle(\bm{\Gamma}\otimes\mathbf{I}_{d})\frac{\partial\mathrm{vec}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}+(\mathbf{I}_{p}\otimes(\bm{\Gamma}^{T}\bm{\Gamma})^{-1})\mathbf{K}_{pd}$
$\displaystyle=$
$\displaystyle-(\bm{\Gamma}\otimes\mathbf{I}_{d})((\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}))(\mathbf{I}_{d^{2}}+\mathbf{K}_{dd})(\mathbf{I}_{d}\otimes\Gamma^{T})+(\mathbf{I}_{p}\otimes(\bm{\Gamma}^{T}\bm{\Gamma})^{-1})\mathbf{K}_{pd}.$

Then,

$\displaystyle\frac{\partial\mathrm{vec}\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}$
$\displaystyle=$
$\displaystyle(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{I}_{p})+(\mathbf{I}_{p}\otimes\bm{\Gamma})\times$
$\displaystyle\left[-(\bm{\Gamma}\otimes\mathbf{I}_{d})((\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes(\bm{\Gamma}^{T}\bm{\Gamma})^{-1})(\mathbf{I}_{d^{2}}+\mathbf{K}_{dd})(\mathbf{I}_{d}\otimes\bm{\Gamma}^{T})+(\mathbf{I}_{p}\otimes(\bm{\Gamma}^{T}\bm{\Gamma})^{-1})\mathbf{K}_{pd}\right]$
$\displaystyle=$
$\displaystyle(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{I}_{p})+(\mathbf{I}_{p}\otimes\bm{\Gamma}(\bm{\Gamma}\bm{\Gamma})^{-1}){\mathbf{K}}_{pd}-(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T})$
$\displaystyle-(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}){\mathbf{K}}_{dd}(\mathbf{I}_{d}\otimes\bm{\Gamma}^{T})$
$\displaystyle=$
$\displaystyle(\mathbf{I}_{p^{2}}+{\mathbf{K}}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{I}_{p})-(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T})$
$\displaystyle-(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1})(\bm{\Gamma}^{T}\otimes\mathbf{I}_{d}){\mathbf{K}}_{pd}$
$\displaystyle=$
$\displaystyle(\mathbf{I}_{p^{2}}+{\mathbf{K}}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{I}_{p})-(\mathbf{I}_{p^{2}}+{\mathbf{K}}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{P}_{\bm{{\scriptstyle\Gamma}}})$
$\displaystyle=$
$\displaystyle(\mathbf{I}_{p^{2}}+{\mathbf{K}}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{I}_{p}-\mathbf{P}_{\bm{{\scriptstyle\Gamma}}})$
$\displaystyle=$
$\displaystyle(\mathbf{I}_{p^{2}}+{\mathbf{K}}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{Q}_{\bm{{\scriptstyle\Gamma}}}).$

###### Lemma 11

Suppose the two matrices $\widehat{\bm{\Gamma}}$ and $\bm{\Gamma}$ are of
order $p\times d$ with $d\leq p$ with $\bm{\Gamma}$ of full rank $d$. Assume
that $\widehat{\bm{\Gamma}}$ is asymptotically normal with

$\sqrt{n}\mathrm{vec}(\widehat{\bm{\Gamma}}-\bm{\Gamma})\xrightarrow{\mathcal{D}}\mathcal{N}(0,\mathbf{V}).$

Then,
$\sqrt{n}\mathrm{vec}(\mathbf{P}_{\widehat{\bm{{\scriptstyle\Gamma}}}}-\mathbf{P}_{\bm{{\scriptstyle\Gamma}}})$
is asymptotically normal with mean $\mathbf{0}$ and variance

$(\mathbf{I}_{p^{2}}+\mathbf{K}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{Q}_{\bm{{\scriptstyle\Gamma}}})\mathbf{V}((\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T}\otimes\mathbf{Q}_{\bm{{\scriptstyle\Gamma}}})(\mathbf{I}_{p^{2}}+{\mathbf{K}}_{pp}).$

Proof  Let
$\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}=\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T}$
be the orthogonal projection onto the column space of $\bm{\Gamma}$ and let
$g$ be a function defined in the subspace of the $p\times d$ matrices of full
rank $d$ such that
$g(\bm{\Gamma})=\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\bm{\Gamma}^{T}=\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}$.
From Lemma 10 we have that

$\nabla
g(\bm{\Gamma})=\frac{\partial{\mathbf{P}_{\bm{{\scriptstyle\Gamma}}}}}{\partial\mathrm{vec}^{T}(\bm{\Gamma})}=(\mathbf{I}_{p^{2}}+\mathbf{K}_{pp})(\bm{\Gamma}(\bm{\Gamma}^{T}\bm{\Gamma})^{-1}\otimes\mathbf{Q}_{\bm{{\scriptstyle\Gamma}}}).$

By the Delta method,

$\sqrt{n}\left(g(\widehat{\bm{\Gamma}})-g(\bm{\Gamma})\right)\rightarrow\mathcal{N}\left(0,\nabla
g(\bm{\Gamma})\mathbf{V}\nabla^{T}g(\bm{\Gamma})\right),$

which completes the proof.

### Proof of Proposition 6

From (70), $\widehat{\bm{\alpha}}_{\mathbf{b}}=\widehat{\mathbf{U}}_{1}$ and
therefore $\bm{\alpha}_{\mathbf{b}}=\mathbf{U}_{1}$ and
$\mathrm{span}(\widehat{\mathbf{U}}_{1})=\mathrm{span}(\widehat{\mathbf{H}})$
with
$\widehat{\mathbf{H}}=\widehat{\mathbf{U}}_{1}\widehat{{\mathbf{K}}}_{1}\widehat{\mathbf{R}}_{1}^{T}\mathbf{R}_{1}{\mathbf{K}}^{-1}$
defined in Lemma 9, which also provides the asymptotic distribution of
$\widehat{\mathbf{H}}$. Applying Lemma 11 with
$\widehat{\bm{\Gamma}}=\widehat{\mathbf{H}}$ and
$\bm{\Gamma}=\mathbf{H}=\mathbf{U}_{1}$ we obtain the asymptotic distribution
with asymptotic variance

$\displaystyle(\mathbf{I}_{p^{2}}+\mathbf{K}_{pp})(\mathbf{U}_{1}{\mathbf{K}}^{-1}\mathbf{R}_{1}^{T}\otimes\mathbf{Q}_{\mathbf{U}_{1}}){\mathbf{V}}_{rlc}(\mathbf{R}_{1}{\mathbf{K}}^{-1}\mathbf{U}_{1}^{T}\otimes\mathbf{Q}_{\mathbf{U}_{1}}))(\mathbf{I}_{p^{2}}+{\mathbf{K}}_{pp})$

since $\mathbf{U}_{1}^{T}\mathbf{U}_{1}=\mathbf{I}_{d}$. By (66),
$\mathbf{b}=\mathbf{U}_{1}\mathbf{K}\mathbf{R}^{T}_{1}$, therefore
$\mathbf{b}^{-}=\mathbf{R}_{1}\mathbf{K}^{-1}\mathbf{U}_{1}^{T}$ and the
result follows.

## References

  * Aitchison and Aitken (1976) J. Aitchison and C. G. G. Aitken.  Multivariate binary discrimination by the kernel method.  _Biometrika_ , 63(3):413–420, 1976.  ISSN 00063444.  URL http://www.jstor.org/stable/2335719.
  * Anderson (1972) J. A. Anderson.  Separate sample logistic discrimination.  _Biometrika_ , 59(1):19–35, 1972.  ISSN 00063444.  URL http://www.jstor.org/stable/2334611.
  * Anderson (1975) J. A. Anderson.  Quadratic logistic discrimination.  _Biometrika_ , 62(1):149–154, 1975.  ISSN 00063444.  URL http://www.jstor.org/stable/2334497.
  * Bach et al. (2012) Francis Bach, Rodolphe Jenatton, Julien Mairal, and Guillaume Obozinski.  Structured sparsity through convex optimization.  _Statist. Sci._ , 27(4):450–468, 11 2012.  doi: 10.1214/12-STS394.  URL https://doi.org/10.1214/12-STS394.
  * Bura and Yang (2011) E. Bura and J. Yang.  Dimension estimation in sufficient dimension reduction: A unifying approach.  _Journal of Multivariate Analysis_ , 102(1):130 – 142, 2011.  ISSN 0047-259X.  doi: https://doi.org/10.1016/j.jmva.2010.08.007.  URL http://www.sciencedirect.com/science/article/pii/S0047259X10001661.
  * Bura et al. (2016) E. Bura, S. Duarte, and L. Forzani.  Sufficient reductions in regressions with exponential family inverse predictors.  _Journal of the American Statistical Association_ , 111(515):1313–1329, 2016.
  * Bura and Forzani (2015) Efstathia Bura and Liliana Forzani.  Sufficient reductions in regressions with elliptically contoured inverse predictors.  _Journal of the American Statistical Association_ , 110(509):420–434, 2015.  doi: 10.1080/01621459.2014.914440.  URL https://doi.org/10.1080/01621459.2014.914440.
  * Buuren (2018) Stef van Buuren.  _Flexible imputation of missing data_.  CRC Press, 2nd edition, 2018.
  * Camiz and Gomes (2013) S. Camiz and G.C. Gomes.  Joint correspondence analysis versus multiple correspondence analysis: a solution to an undetected problem.  In _Classification and data mining_ , Stud. Classification Data Anal. Knowledge Organ., pages 11–18. Springer, Heidelberg, 2013.
  * Chavent et al. (2012) M. Chavent, V. Kuentz-Simonet, B. Liquet, and J. Saracco.  Orthogonal rotation in pcamix.  _Advances in Data Analysis and Classification_ , 6:131–146, 2012.
  * Chavent et al. (2014) Marie Chavent, Vanessa Kuentz-Simonet, Amaury Labenne, and Jérôme Saracco.  Multivariate analysis of mixed data: The r package pcamixdata, 2014.
  * Chen et al. (2014) Shizhe Chen, Daniela M. Witten, and Ali Shojaie.  Selection and estimation for mixed graphical models.  _Biometrika_ , 102(1):47–64, 12 2014.  ISSN 0006-3444.  doi: 10.1093/biomet/asu051.  URL https://doi.org/10.1093/biomet/asu051.
  * Cheng et al. (2014) Jie Cheng, Elizaveta Levina, Pei Wang, and Ji Zhu.  A sparse Ising model with covariates.  _Biometrics_ , 70(4):943–953, 2014.  ISSN 0006-341X.  doi: 10.1111/biom.12202.  URL https://doi.org/10.1111/biom.12202.
  * Cheng et al. (2017) Jie Cheng, Tianxi Li, Elizaveta Levina, and Ji Zhu.  High-dimensional mixed graphical models.  _Journal of Computational and Graphical Statistics_ , 26(2):367–378, 2017.  doi: 10.1080/10618600.2016.1237362.  URL https://doi.org/10.1080/10618600.2016.1237362.
  * Cook (2007) R.D. Cook.  Fisher lecture: Dimension reduction in regression (with discussion).  _Statistical Science_ , 22:1–26, 2007.
  * Cook and Forzani (2008) R.D. Cook and L. Forzani.  Principal fitted components for dimension reduction in regression.  _Statistical Science_ , 23:485–501, 2008.
  * Cook and Weisberg (1991) R.D. Cook and S. Weisberg.  Discussion of sliced inverse regression for dimension reduction.  _Journal of the American Statistical Association_ , 86:328–332, 1991.
  * Dai (2012) Bin Dai.  Multivariate bernoulli distribution models.  Technical report, Dept. Statistics, Univ. Wisconsin, Madison, WI 53706, July 2012.
  * Dai et al. (2013) Bin Dai, Shilin Ding, and Grace Wahba.  Multivariate bernoulli distribution.  _Bernoulli_ , 19(4):1465–1483, 09 2013.  doi: 10.3150/12-BEJSP10.  URL https://doi.org/10.3150/12-BEJSP10.
  * Day and Kerridge (1967) N. E. Day and D. F. Kerridge.  A general maximum likelihood discriminant.  _Biometrics_ , 23(2):313–323, 1967.  ISSN 0006341X, 15410420.  URL http://www.jstor.org/stable/2528164.
  * Filmer and Scott (2012) D. Filmer and K. Scott.  Assessing Asset Indices.  _Demography_ , 49:359–392, 2012.
  * Fitzmaurice and Laird (1997) Garrett M. Fitzmaurice and Nan M. Laird.  Regression models for mixed discrete and continuous responses with potentially missing values.  _Biometrics_ , 53(1):110–122, 1997.  ISSN 0006341X, 15410420.  URL http://www.jstor.org/stable/2533101.
  * Forzani et al. (2018) L. Forzani, R. García-Arancibia, P. Llop, and D. Tomassi.  Supervised dimension reduction for ordinal predictors.  _Computational Statistics and Data Analysis_ , 125, 2018.
  * Greenacre and Blasius (2006) M. Greenacre and J. Blasius, editors.  _Multiple correspondence analysis and related methods_.  Statistics in the Social and Behavioral Sciences Series. Chapman & Hall/CRC, Boca Raton, FL, 2006.  ISBN 978-1-58488-628-0; 1-58488-628-5.  doi: 10.1201/9781420011319.
  * Ising (1925) Ernst Ising.  Beitrag zur theorie des ferromagnetismus.  _Zeitschrift für Physik_ , 31(1):253–258, Feb 1925.  ISSN 0044-3328.  doi: 10.1007/BF02980577.  URL https://doi.org/10.1007/BF02980577.
  * Javaras and van Dyk (2003) Kristin N. Javaras and David A. van Dyk.  Multiple imputation for incomplete data with semicontinuous variables.  _Journal of the American Statistical Association_ , 98(463):703–715, 2003.  ISSN 01621459.  URL http://www.jstor.org/stable/30045298.
  * Kolenikov and Angeles (2009) S. Kolenikov and G. Angeles.  Socioeconomic status measurement with discrete proxy variables: Is principal component analysis a reliable answer?  _The Review of Income and Wealth_ , 55(1):128–165, 2009.
  * Krzanowski (1993) W. J. Krzanowski.  The location model for mixtures of categorical and continuous variables.  _Journal of Classification_ , 10(1):25–49, Jan 1993.  ISSN 1432-1343.  doi: 10.1007/BF02638452.  URL https://doi.org/10.1007/BF02638452.
  * Krzanowski (1975) W.J. Krzanowski.  Discrimination and classification using both binary and continuous variables.  _Journal of the American Statistical Association_ , 70(352):782–790, 1975.
  * Lauritzen (1996) S. L. Lauritzen.  _Graphical Models_.  Oxford University Press, Oxford, 1996.
  * Lauritzen and Wermuth (1989) S. L. Lauritzen and N. Wermuth.  Graphical models for associations between variables, some of which are qualitative and some quantitative.  _Ann. Statist._ , 17(1):31–57, 03 1989.  doi: 10.1214/aos/1176347003.  URL https://doi.org/10.1214/aos/1176347003.
  * Lee and Hastie (2015) Jason D. Lee and Trevor J. Hastie.  Learning the structure of mixed graphical models.  _Journal of Computational and Graphical Statistics_ , 24(1):230–253, 2015.  doi: 10.1080/10618600.2014.900500.  URL https://doi.org/10.1080/10618600.2014.900500.  PMID: 26085782.
  * Li and Wang (2007) B. Li and S. Wang.  On directional regression for dimension reduction.  _Journal of the American Statistical Association_ , 102(479):997–1008, 2007.
  * Li (1991) K. C. Li.  Sliced inverse regression for dimension reduction (with discussion).  _Journal of the American Statistical Association_ , 86:316–342, 1991.
  * Liu and Ye (2010) Jun Liu and Jieping Ye.  Fast overlapping group lasso.  _arXiv:1009.0306v1_ , 2010.
  * Mazziotta and Pareto (2019) Matteo Mazziotta and Adriano Pareto.  Use and misuse of pca for measuring well-being.  _Social Indicators Research_ , 142(2):451–476, Apr 2019.  ISSN 1573-0921.  doi: 10.1007/s11205-018-1933-0.
  * Merola and Baulch (2014) G. Merola and B. Baulch.  Using sparse categorical principal components to estimate asset indices new methods with an application to rural south east asia.  2014\.
  * Morris (2006) Carl N. Morris.  _Natural Exponential Families_.  American Cancer Society, 2006.  ISBN 9780471667193.  doi: 10.1002/0471667196.ess1759.pub2.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/0471667196.ess1759.pub2.
  * Olkin and Tate (1961) I. Olkin and R. F. Tate.  Multivariate correlation models with mixed discrete and continuous variables.  _Ann. Math. Statist._ , 32(2):448–465, 06 1961\.  doi: 10.1214/aoms/1177705052.  URL https://doi.org/10.1214/aoms/1177705052.
  * Pepe (2003) M.S. Pepe.  _The Statistical Evaluation of Medical Tests for Classification and Prediction_.  Oxford University Press, New York, 2003.
  * Vlachonikolis and Marriott (1982) I. G. Vlachonikolis and F. H. C. Marriott.  Discrimination with mixed binary and continuous data.  _Journal of the Royal Statistical Society. Series C (Applied Statistics)_ , 31(1):23–31, 1982.  ISSN 00359254, 14679876.  URL http://www.jstor.org/stable/2347071.
  * Vyas and Kumaranayake (2006) S. Vyas and L. Kumaranayake.  Constructing socio-economic status indices: How to use principal components analysis.  _Health Policy and Planning_ , 21(6):459–468, 2006.
  * Wainwright and Jordan (2008) Martin J. Wainwright and Michael I. Jordan.  Graphical models, exponential families, and variational inference.  _Foundations and Trends® in Machine Learning_ , 1(1–2):1–305, 2008.  ISSN 1935-8237.  doi: 10.1561/2200000001.  URL http://dx.doi.org/10.1561/2200000001.
  * Whittaker (2009) Joe Whittaker.  _Graphical Models in Applied Multivariate Statistics_.  Wiley Publishing, 2009.  ISBN 0470743662, 9780470743669.
  * Yang et al. (2014a) Eunho Yang, Yulia Baker, Pradeep Ravikumar, Genevera Allen, and Zhandong Liu.  Mixed Graphical Models via Exponential Families.  In Samuel Kaski and Jukka Corander, editors, _Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics_ , volume 33 of _Proceedings of Machine Learning Research_ , pages 1042–1050, Reykjavik, Iceland, 22–25 Apr 2014a. PMLR.  URL http://proceedings.mlr.press/v33/yang14a.html.
  * Yang et al. (2014b) Eunho Yang, Pradeep Ravikumar, Genevera I. Allen, Yulia Baker, Ying-Wooi Wan, and Zhandong Liu.  A general framework for mixed graphical models, 2014b.
  * Yang et al. (2015) Eunho Yang, Pradeep Ravikumar, Genevera I. Allen, Zh, and ong Liu.  Graphical models via univariate exponential family distributions.  _Journal of Machine Learning Research_ , 16(115):3813–3847, 2015.  URL http://jmlr.org/papers/v16/yang15a.html.
  * Ye and Lim (2016) Ke Ye and Lek-Heng Lim.  Schubert varieties and distances between subspaces of different dimensions.  _SIAM Journal on Matrix Analysis and Applications_ , 37(3):1176–1197, 2016.  doi: 10.1137/15M1054201.  URL https://doi.org/10.1137/15M1054201.
  * Yuan and Lin (2006) Ming Yuan and Yi Lin.  Model selection and estimation in regression with grouped variables.  _Journal of the Royal Statistical Society. Series B, statistical methodology_ , 2006.  ISSN 1369-7412.
  * Yuan and Lin (2007) Ming Yuan and Yi Lin.  Model selection and estimation in the gaussian graphical model.  _Biometrika_ , 94(1):19–35, 2007.  ISSN 00063444.  URL http://www.jstor.org/stable/20441351.
  * Zhu et al. (2011) Q Zhu, L. Lin, M.-L. Shyu, and S.-C. Chen.  Effective supervised discretization for classification based on correlation maximization.  _IEEE International Conference on Information Reuse & Integration_, pages 390–395, 2011.