# Supporting Information for ‘Sufficient Dimension Reduction for Populations
with Structured Heterogeneity’ by

Jared D. Huling1, Menggang Yu2

1Division of Biostatistics, University of Minnesota, Minneapolis, Minnesota
2Department of Biostatistics and Medical Informatics,
University of Wisconsin-Madison
huling<EMAIL_ADDRESS>
##  Web Appendix A Additional results

###  Web Appendix A.1 Selection of the bandwidth parameter

Another issue in practice is the selection of the bandwidth in the
nonparametric regressions for the terms in (LABEL:eqn:sample_est_eqn). As is
noted in Ma and Zhu, (2014), the estimator is not critically sensitive to this
choice and can be chosen via cross validation or other standard approaches.
Throughout this paper use a simple but effective subpopulation-specific choice
of the bandwidth,
$c_{\boldsymbol{z}}(3n_{\boldsymbol{z}}/4)^{-1/(d_{\boldsymbol{z}}+4)}$ (Ma
and Zhu,, 2012), where $c_{\boldsymbol{z}}$ is the average standard deviation
of
$(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{\boldsymbol{z}}:i=1,\dots,n,{\boldsymbol{Z}}_{i}={\boldsymbol{z}})$.

###  Web Appendix A.2 Asymptotic results for unconstrained estimator
$\widetilde{\boldsymbol{\beta}}$ in (LABEL:eqn:est_unconstr)

Recall that we let $\mathrm{vecl}(\boldsymbol{A})$ be the operator which
vectorizes the lower $(a-b)\times b$ block of an $a\times b$ matrix
$\boldsymbol{A}$ when $a\geq b$. Denote
$\mathrm{vecl}{(\boldsymbol{\beta})}^{\mathsf{T}}\equiv(\mathrm{vecl}{(\boldsymbol{\beta}_{\boldsymbol{z}})}^{\mathsf{T}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$
as the lower block vectorization of
$\boldsymbol{\beta}\equiv(\boldsymbol{\beta}_{\boldsymbol{z}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$.

######  Theorem S1.

Assume regularity conditions (C1)-(C5) listed in Section Web Appendix B.1 hold
for
$\boldsymbol{\beta}_{\boldsymbol{z}}^{*}:{\boldsymbol{z}}\in\mathcal{Z}^{C}$.
Further assume that ${\boldsymbol{W}}_{n}\xrightarrow{p}\boldsymbol{W}$ as
$n\to\infty$, where $\boldsymbol{W}$ is a symmetric positive semi-definite
matrix. Let

$\displaystyle\boldsymbol{G}_{1}={}$ $\displaystyle
E\big{(}\partial\mathrm{vec}\big{[}\left\\{Y-{E}(Y|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}})\right\\}\left\\{\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})-{E}(\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}})\right\\}\big{]}/\partial\mathrm{vecl}(\boldsymbol{\beta})^{\mathsf{T}}\big{)}.$

Further assume that
$\boldsymbol{\Sigma}_{1}=\\{\boldsymbol{G}_{1}^{\mathsf{T}}\boldsymbol{W}\boldsymbol{G}_{1}\\}^{-1}$
is nonsingular. Then the minimizer $\widetilde{\boldsymbol{\beta}}$ of
(LABEL:eqn:est_unconstr) satisfies

$\sqrt{n}\mathrm{vecl}(\widetilde{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*})\xrightarrow{\mathcal{L}}\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{1}\boldsymbol{\Xi}_{1}\boldsymbol{\Sigma}_{1}),$
(S.1)

where $\xrightarrow{\mathcal{L}}$ is convergence in distribution,
$\boldsymbol{\beta}^{*}=(\boldsymbol{\beta}^{*}_{\boldsymbol{z}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$,
$\boldsymbol{\Xi}_{1}=\boldsymbol{G}_{1}^{\mathsf{T}}\boldsymbol{W}\boldsymbol{V}_{1}\boldsymbol{W}\boldsymbol{G}_{1}$,
and

$\displaystyle\boldsymbol{V}_{1}={}$
$\displaystyle\mbox{cov}\big{(}\mathrm{vec}\big{[}\left\\{Y-{E}(Y|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}})\right\\}\left\\{\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})-{E}(\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}})\right\\}\big{]}\big{)}.$

An important point is that Theorem S1 is for fixed dimensions
$d_{\boldsymbol{z}}$ similar to Ma and Zhu, (2012). In practice these
dimensions need to be estimated. We present an approach for such determination
in the following section, which is a special case of the dimension
determination approach we propose in Section LABEL:sec:method. Based on (S.1),
a sensible choice for ${\boldsymbol{W}}_{n}$ would be
$\widehat{\boldsymbol{V}}_{1n}^{-1}$, where $\widehat{\boldsymbol{V}}_{1n}$ is
any consistent estimator of ${\boldsymbol{V}}_{1}$.

Note that $\boldsymbol{G}_{2}\neq\boldsymbol{G}_{1}$ and
$\boldsymbol{V}_{2}\neq\boldsymbol{V}_{1}$ due to the difference in how the
upper blocks of parameters are imposed in $\widehat{\boldsymbol{\beta}}$
versus $\widetilde{\boldsymbol{\beta}}$, indicating a differing number of free
parameters. However, comparing the variance terms between (S.1) and
(LABEL:asymptotic_distribution_of_beta_hat), we would generally expect
improvement in efficiency due to the product term
$(\boldsymbol{I}-\boldsymbol{P})$ in
(LABEL:asymptotic_distribution_of_beta_hat).

###  Web Appendix A.3 Relationship of $\widetilde{\boldsymbol{\beta}}$ with Li
et al., (2003)

In our estimator, the emphasis is on estimating the central mean subspaces
conditional on subpopulations, as opposed to the central partial mean
subspace, which is the direct sum of the conditional mean subspaces (Li et
al.,, 2003).

We note that the central partial mean subspace of Li et al., (2003) can be
seen as a special case of (LABEL:eqn:est_unconstr) of the main text. To see
this, let $\boldsymbol{\alpha}(\mathbf{x},{\boldsymbol{z}})$ be the vector
that takes all values equal to zero except in the $p$ locations corresponding
to the terms $\boldsymbol{\beta}_{\boldsymbol{z}}$ in $\boldsymbol{\beta}$
where it takes values equal to $\mathbf{x}$. Further, take
${E}(Y|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}})=0$
and
${E}(\mathbf{X}|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}})=0$,
which requires additional assumptions about $\mathbf{X}$ for consistent
estimation. The approach of Li et al., (2003) is based on ordinary least
squares (OLS) regressions for each subpopulation, limiting the dimension of
the central mean subspace for each subpopulation to be 1. Thus, to compare
(LABEL:eqn:est_unconstr) with this approach, assume that each
$d_{\boldsymbol{z}}=1$, and further take
${\boldsymbol{W}}_{n}={\boldsymbol{I}}$. Then (LABEL:eqn:est_unconstr) reduces
to fitting OLS regressions for each subpopulation. Thus, if we use the span of
the resulting subpopulation-specific estimates
$\widetilde{\boldsymbol{\beta}}_{\boldsymbol{z}}$ to recover the central
partial mean subspace, (LABEL:eqn:est_unconstr) relates directly to the
approach of Li et al., (2003). Due to the general form of
(LABEL:eqn:pop_est_eqn), pooled OLS-type estimators are also possible, for one
possibility among many, one could take a common
$\boldsymbol{\alpha}(\mathbf{x},{\boldsymbol{z}})=\text{rep}(\mathbf{x},2^{C})$
for all subpopulations, where $\text{rep}(\cdot)$ forms a vector by repeating
$\mathbf{x}$ for $2^{C}$ times.

###  Web Appendix A.4 Example of construction of the constraint matrix
$\boldsymbol{C}$

Here, we show how the constraint matrix $\boldsymbol{C}$ is constructed to
enforce the hierarchical conditions. For illustrative purposes we consider a
simplified scenario with models stratified by two factors and a specific
choice of structural dimensions, but the principle applies more generally. Let
$\mathcal{Z}^{C}=\\{00,01,10,11\\}$, $d_{00}=1$, $d_{10}=2$, $d_{01}=2$, and
$d_{11}=4$, so that the combined unconstrained parameter matrix
$\boldsymbol{\beta}=(\boldsymbol{\beta}_{{\boldsymbol{z}}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$
can be expressed as

$\displaystyle\boldsymbol{\beta}=(\underbrace{\boldsymbol{\nu}_{00}}_{\boldsymbol{\beta}_{00}},\underbrace{\boldsymbol{\nu}^{1}_{00},\boldsymbol{\nu}_{01}}_{\boldsymbol{\beta}_{01}},\underbrace{\boldsymbol{\nu}^{2}_{00},\boldsymbol{\nu}_{10}}_{\boldsymbol{\beta}_{10}},\underbrace{\boldsymbol{\nu}^{3}_{00},\boldsymbol{\nu}^{1}_{01},\boldsymbol{\nu}^{1}_{10},\boldsymbol{\nu}_{11}}_{\boldsymbol{\beta}_{11}}),$

where it is not necessarily the case that $\boldsymbol{\nu}_{00}$ equals
$\boldsymbol{\nu}^{1}_{00}$, $\boldsymbol{\nu}^{2}_{00}$, or
$\boldsymbol{\nu}^{3}_{00}$ and not necessarily the case that
$\boldsymbol{\nu}_{01}=\boldsymbol{\nu}^{1}_{01}$ or
$\boldsymbol{\nu}_{10}=\boldsymbol{\nu}^{1}_{10}$, but we use this naming
convention as our goal in constructing $\boldsymbol{C}$ is to make
$\boldsymbol{\nu}_{01}=\boldsymbol{\nu}^{1}_{01}$,
$\boldsymbol{\nu}_{10}=\boldsymbol{\nu}^{1}_{10}$, and so on. Here, the
dimension of all $\boldsymbol{\nu}$ terms is 1, so the constraint that the
upper blocks of the $\boldsymbol{\nu}$ terms is identity implies that the free
parameters in each are the lower $p-1$ elements. We denote the free parameters
of $\boldsymbol{\nu}_{00}$ as $\boldsymbol{\nu}_{00,L}$, the free parameters
of $\boldsymbol{\nu}^{1}_{00}$ as $\boldsymbol{\nu}^{1}_{00,L}$, and so on.
Using this notation, here,
$\mathrm{vecl}(\boldsymbol{\beta})=(\boldsymbol{\nu}_{00,L},\boldsymbol{\nu}^{1}_{00,L},\boldsymbol{\nu}_{01,L},\boldsymbol{\nu}^{2}_{00,L},\boldsymbol{\nu}_{10,L},\boldsymbol{\nu}^{3}_{00,L},\boldsymbol{\nu}^{1}_{01,L},\boldsymbol{\nu}^{1}_{10,L},\boldsymbol{\nu}_{11,L})^{T}$,
since the dimension of each $\boldsymbol{\nu}$ term here is 1. To enforce the
hierarchical constraints, we need

$\displaystyle\boldsymbol{C}^{T}\mathrm{vecl}(\boldsymbol{\beta})=\begin{pmatrix}\boldsymbol{\nu}_{00,L}-\boldsymbol{\nu}^{1}_{00,L}\\\
\boldsymbol{\nu}_{00,L}-\boldsymbol{\nu}^{2}_{00,L}\\\
\boldsymbol{\nu}_{00,L}-\boldsymbol{\nu}^{3}_{00,L}\\\
\boldsymbol{\nu}_{01,L}-\boldsymbol{\nu}^{1}_{01,L}\\\
\boldsymbol{\nu}_{10,L}-\boldsymbol{\nu}^{1}_{10,L}\\\
\end{pmatrix}=\boldsymbol{0},$

which is a $5p\times 1$ matrix. The constraint matrix that yields this is

$\boldsymbol{C}^{T}=\begin{pmatrix}\boldsymbol{I}_{p-1}&-\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}&-\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&-\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&-\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{I}_{p-1}&\boldsymbol{0}&\boldsymbol{0}&-\boldsymbol{I}_{p-1}&\boldsymbol{0}\\\
\end{pmatrix},$

where $\boldsymbol{I}_{p-1}$ are $(p-1)\times(p-1)$ identity matrices.

The key principle here is to construct a $\boldsymbol{C}$ that applies the
appropriate contrasts that, when set to zero, enforce the desired constraints.
Of course, if the terms in $\boldsymbol{\beta}$ are switched around,
$\boldsymbol{C}$ will be different, but these differences are unimportant and
is not consequential in any way. When the dimension $k$ of any
$\boldsymbol{\nu}$ term is larger than 1, an additional set of rows of
constraints is necessary, but the dimension of each set of constraints would
be $p-k$ instead of $p-1$.

###  Web Appendix A.5 Determination of structural dimensions

Determining the structural dimension in SDR in the semiparametric framework is
a challenging problem. The challenges inherent in dimension determination are
amplified for the determination of multiple dimensions, especially when the
different dimensions are relevant across multiple subpopulations. In the
setting of this paper, we are interested in estimating the structural
dimensions of the true parameters $\boldsymbol{\nu}^{0}_{\boldsymbol{z}}$ for
${\boldsymbol{z}}\in\mathcal{Z}^{C}$. Due to (LABEL:eqn:equality_condition),
the set of dimensions of $\boldsymbol{\nu}_{\boldsymbol{z}}$ for
${\boldsymbol{z}}\in\mathcal{Z}^{C}$ uniquely determines the set of dimensions
of
$\boldsymbol{\beta}_{\boldsymbol{z}}=(\boldsymbol{\nu}_{{\boldsymbol{z}}^{\prime}}:{\boldsymbol{z}}^{\prime}\preceq{\boldsymbol{z}},{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C})$
for ${\boldsymbol{z}}\in\mathcal{Z}^{C}$; thus we focus on the former.

Ma and Zhang, (2015) proposed a simple, effective, and computationally
tractable information criterion-based approach, called the validated
information criterion (VIC). We extend their approach to our setting with the
imposed hierarchical assumption (LABEL:eqn:hierarchical_assumption_none) and
(LABEL:eqn:none_assumption2).

Denote a given set of candidate dimensions as ${\boldsymbol{k}\
}\equiv(k_{\boldsymbol{z}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$ for the terms
$\boldsymbol{\nu}_{{\boldsymbol{z}}}$ and denote the true dimensions as
${{\boldsymbol{k}\
}\\!\\!}^{0}\\!=(k_{\boldsymbol{z}}^{0}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$
of the parameters ${\boldsymbol{\nu}}_{\boldsymbol{z}}^{0}$. Our proposal
arises by noting that the population estimating equations should be zero when
$\sum_{{\boldsymbol{z}}^{\prime}\preceq{\boldsymbol{z}}}k_{{\boldsymbol{z}}^{\prime}}\equiv{d}_{\boldsymbol{z}}\geq{d}_{\boldsymbol{z}}^{0}\equiv\sum_{{\boldsymbol{z}}^{\prime}\preceq{\boldsymbol{z}}}k_{{\boldsymbol{z}}^{\prime}}^{0}$
for all ${\boldsymbol{z}}$ and should be nonzero when
${d}_{\boldsymbol{z}}<{d}_{\boldsymbol{z}}^{0}$ for any ${\boldsymbol{z}}$.
Denote
${\boldsymbol{d}}=({d}_{\boldsymbol{z}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$,
${\boldsymbol{d}}^{0}=({d}_{\boldsymbol{z}}^{0}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$,
and
$\boldsymbol{\beta}_{(\boldsymbol{d})}=(\boldsymbol{\beta}_{{\boldsymbol{z}}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$
where each $\boldsymbol{\beta}_{{\boldsymbol{z}}}$ is of dimension $p\times
d_{\boldsymbol{z}}$ and defined through (LABEL:eqn:equality_condition). We
estimate $\boldsymbol{\beta}_{(\boldsymbol{d})}$ similar to (S.19) by solving

$\displaystyle\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}=\operatornamewithlimits{argmin}_{\boldsymbol{\beta}_{(\boldsymbol{d})}=(\boldsymbol{\beta}_{{\boldsymbol{z}}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})}{}\biggl{\\{}$
$\displaystyle\frac{1}{2n}\hat{\Psi}_{n}(\boldsymbol{\beta}_{(\boldsymbol{d})})^{\mathsf{T}}\boldsymbol{W}_{n}\hat{\Psi}_{n}(\boldsymbol{\beta}_{(\boldsymbol{d})})\mbox{
s.t.
}\boldsymbol{\boldsymbol{C}}_{(\boldsymbol{d})}^{\mathsf{T}}\mathrm{vecl}(\boldsymbol{\beta}_{(\boldsymbol{d})})=\boldsymbol{0}\biggr{\\}},$

where $\boldsymbol{C}_{(\boldsymbol{d})}$ is constructed similarly to
$\boldsymbol{C}$ but accounting for different dimensions of the terms in
$\boldsymbol{\beta}_{(\boldsymbol{d})}$ from those under true dimensions
${\boldsymbol{d}}^{0}$. Similarly denote the population version as

$\displaystyle\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}=\operatornamewithlimits{argmin}_{\boldsymbol{\beta}_{(\boldsymbol{d})}=(\boldsymbol{\beta}_{{\boldsymbol{z}}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})}{}\biggl{\\{}$
$\displaystyle\frac{1}{2}\Psi_{0}(\boldsymbol{\beta}_{(\boldsymbol{d})})^{\mathsf{T}}\boldsymbol{W}_{n}\Psi_{0}(\boldsymbol{\beta}_{(\boldsymbol{d})})\mbox{
s.t.
}\boldsymbol{\boldsymbol{C}}_{(\boldsymbol{d})}^{\mathsf{T}}\mathrm{vecl}(\boldsymbol{\beta}_{(\boldsymbol{d})})=\boldsymbol{0}\biggr{\\}}.$

Note that under the true structural dimensions $\boldsymbol{d}^{0}$,
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}$ should yield estimated subspaces
that satisfy (LABEL:eqn:hierarchical_assumption_none) and
(LABEL:eqn:none_assumption2) if they are indeed correct, in which case
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}=\boldsymbol{\beta}^{0}$.

Since the estimates $\widehat{\boldsymbol{\beta}}$ are constructed by
minimizing the norm of the estimating equations, evaluation of the estimating
equations at $\widehat{\boldsymbol{\beta}}$ does not necessarily provide a
valid means of assessing the fit of the model. As a workaround, similar to Ma
and Zhang, (2015), the basic idea of the approach is to construct expansions
of the parameters $\boldsymbol{\nu}_{{\boldsymbol{z}}}$ for each population
from dimension $k_{\boldsymbol{z}}$ to $k_{\boldsymbol{z}}+1$. The expanded
parameters are used as a means to assess model fit. In particular, write
${\boldsymbol{\nu}_{{\boldsymbol{z}}}^{0}}^{\mathsf{T}}=(I_{k_{\boldsymbol{z}}},{\boldsymbol{\nu}^{0,U}_{{\boldsymbol{z}}}}^{\mathsf{T}},{\boldsymbol{\nu}^{0,L}_{{\boldsymbol{z}}}}^{\mathsf{T}})$
and
${\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}}^{\mathsf{T}}=(I_{k_{\boldsymbol{z}}},{{}\widehat{\boldsymbol{\nu}}^{U}_{{\boldsymbol{z}}}}^{\mathsf{T}},{{}\widehat{\boldsymbol{\nu}}^{L}_{{\boldsymbol{z}}}}^{\mathsf{T}})$,
where $\boldsymbol{\nu}^{0,U}_{{\boldsymbol{z}}}$ and
$\widehat{\boldsymbol{\nu}}^{U}_{{\boldsymbol{z}}}$ are $1\times
k_{\boldsymbol{z}}$ vectors, $\boldsymbol{\nu}^{0,L}_{{\boldsymbol{z}}}$ and
$\widehat{\boldsymbol{\nu}}^{L}_{{\boldsymbol{z}}}$ are
$(p-1-k_{\boldsymbol{z}})\times k_{\boldsymbol{z}}$ matrices where for
simplicity of presentation we drop explicit reference to the dimension in the
notation. The matrices $\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}(v)$ and
$\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}(v)$ are constructed so that for
any $v$,
$\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}(v)(I_{k_{\boldsymbol{z}}},\boldsymbol{\nu}^{0,U}_{{\boldsymbol{z}}})^{\mathsf{T}}=\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}$
and
$\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}(v)(I_{k_{\boldsymbol{z}}},\widehat{\boldsymbol{\nu}}^{U}_{{\boldsymbol{z}}})^{\mathsf{T}}=\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}$.
Therefore the spaces spanned by the columns of
$\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}(v)$ and
$\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}(v)$ contain the spaces spanned
by the columns of $\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}$ and
$\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}$ respectively. In particular,
write
${\boldsymbol{\nu}_{{\boldsymbol{z}}}^{0}}^{\mathsf{T}}=(I_{k_{\boldsymbol{z}}},{\boldsymbol{\nu}^{0,U}_{{\boldsymbol{z}}}}^{\mathsf{T}},{\boldsymbol{\nu}^{0,L}_{{\boldsymbol{z}}}}^{\mathsf{T}})$
and
${\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}}^{\mathsf{T}}=(I_{k_{\boldsymbol{z}}},{{}\widehat{\boldsymbol{\nu}}^{U}_{{\boldsymbol{z}}}}^{\mathsf{T}},{{}\widehat{\boldsymbol{\nu}}^{L}_{{\boldsymbol{z}}}}^{\mathsf{T}})$,
where $\boldsymbol{\nu}^{0,U}_{{\boldsymbol{z}}}$ and
$\widehat{\boldsymbol{\nu}}^{U}_{{\boldsymbol{z}}}$ are $1\times
k_{\boldsymbol{z}}$ vectors, $\boldsymbol{\nu}^{0,L}_{{\boldsymbol{z}}}$ and
$\widehat{\boldsymbol{\nu}}^{L}_{{\boldsymbol{z}}}$ are
$(p-1-k_{\boldsymbol{z}})\times k_{\boldsymbol{z}}$ matrices where for
simplicity of presentation we drop explicit reference to the dimension in the
notation. Then define

$\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}(v)=\begin{pmatrix}I_{k_{\boldsymbol{z}}}&\boldsymbol{0}_{k_{\boldsymbol{z}}\times
1}\\\ \boldsymbol{0}_{1\times k_{\boldsymbol{z}}}&1\\\
\boldsymbol{\nu}^{0,L}_{{\boldsymbol{z}}}-\boldsymbol{v}\boldsymbol{\nu}^{0,U}_{{\boldsymbol{z}}}&\boldsymbol{v}\end{pmatrix},\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}(v)=\begin{pmatrix}I_{k_{\boldsymbol{z}}}&\boldsymbol{0}_{k_{\boldsymbol{z}}\times
1}\\\ \boldsymbol{0}_{1\times k_{\boldsymbol{z}}}&1\\\
\widehat{\boldsymbol{\nu}}^{L}_{{\boldsymbol{z}}}-\boldsymbol{v}\widehat{\boldsymbol{\nu}}^{U}_{{\boldsymbol{z}}}&\boldsymbol{v}\end{pmatrix},$

where $\boldsymbol{v}$ is a conformable vector with all elements equal to $v$.
It is easy to verify that for any $v$,
$\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}(v)(I_{k_{\boldsymbol{z}}},\boldsymbol{\nu}^{0,U}_{{\boldsymbol{z}}})^{\mathsf{T}}=\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}$
and
$\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}(v)(I_{k_{\boldsymbol{z}}},\widehat{\boldsymbol{\nu}}^{U}_{{\boldsymbol{z}}})^{\mathsf{T}}=\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}$.
The corresponding $\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}}(v)$ is defined
through (LABEL:eqn:equality_condition) with
$\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}(v)$ in place of
$\boldsymbol{\nu}_{{\boldsymbol{z}}}^{0}$ and
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)=(\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}}(v):{\boldsymbol{z}}\in\mathcal{Z}^{C})$.

The corresponding $\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}}(v)$ is defined
through (LABEL:eqn:equality_condition) with
$\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}(v)$ in place of
$\boldsymbol{\nu}_{{\boldsymbol{z}}}^{0}$ and
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)=(\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}}(v):{\boldsymbol{z}}\in\mathcal{Z}^{C})$.
Similarly, $\widehat{\boldsymbol{\beta}}_{{\boldsymbol{z}}}(v)$ is also
defined through (LABEL:eqn:equality_condition) with
$\widehat{\boldsymbol{\nu}}^{0}_{{\boldsymbol{z}}}(v)$ in place of
$\widehat{\boldsymbol{\nu}}_{{\boldsymbol{z}}}^{0}$ and
$\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v)=(\widehat{\boldsymbol{\beta}}_{{\boldsymbol{z}}}(v):{\boldsymbol{z}}\in\mathcal{Z}^{C})$.

Intuitively, if ${d}_{{\boldsymbol{z}}}<{d}_{{\boldsymbol{z}}}^{0}$ for any
${\boldsymbol{z}}\in\mathcal{Z}^{C}$, then
$n^{-1}\Psi_{n}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v))$ is in
general inconsistent for zero. On the other hand, if
${d}_{{\boldsymbol{z}}}\geq{d}_{{\boldsymbol{z}}}^{0}$ for all
${\boldsymbol{z}}\in\mathcal{Z}^{C}$, then
$n^{-1}\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v))$
will still be consistent for zero even if
${d}_{{\boldsymbol{z}}^{\prime}}>{d}_{{\boldsymbol{z}}^{\prime}}^{0}$ for some
${\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C}$. Thus, our proposed criterion
for dimension determination is

$\displaystyle\mbox{VIC}({{\boldsymbol{k}}})={}$
$\displaystyle\frac{1}{rn}\sum_{j=1}^{r}\lVert\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v_{j}))\rVert^{2}+\log(n)p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}{d}_{{\boldsymbol{z}}}.$
(S.2)

We have proposed to use $r$ different values $v_{1},\ldots,v_{r}$ for $v$.
Although a particular $v$ should work but in theory, if
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)=\boldsymbol{\beta}^{0}_{(\boldsymbol{d}+1)}$,
then the VIC criterion for $r=1$ with $v$ will fail to work for dimension
determination. Using $r$ different values for $v_{j}$, it is possible that
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v_{j})=\boldsymbol{\beta}^{0}_{(\boldsymbol{d}+1)}$
for any $r-1$ values $j\in\\{1,\dots,r\\}$. In practice, choosing a moderate
number such as $r=5$ tends to perform well.

To estimate the structural dimensions, one would evaluate (S.2) for all valid
sets of candidate dimensions. The estimates of the structural dimensions for
$\boldsymbol{\nu}^{0}_{{\boldsymbol{z}}}$ for
${\boldsymbol{z}}\in\mathcal{Z}^{C}$, denoted as $\hat{{\boldsymbol{k}\
}}=(\hat{k}_{\boldsymbol{z}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$, is the
minimizer of (S.2). Further, let the estimated structural dimensions for
$\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}}$ for
${\boldsymbol{z}}\in\mathcal{Z}^{C}$ be denoted as $\hat{\boldsymbol{d}\
}=(\sum_{{\boldsymbol{z}}^{\prime}\preceq{\boldsymbol{z}}}\hat{k}_{{\boldsymbol{z}}^{\prime}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$.
It is important to note that
$p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}{d}_{{\boldsymbol{z}}}$, instead of
$p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}k_{{\boldsymbol{z}}}$, should be
used as a penalty in the VIC criterion.

######  Theorem S2.

Define $\boldsymbol{d}\leq{\boldsymbol{d}}^{0}$ to mean
$d_{\boldsymbol{z}}\leq d_{\boldsymbol{z}}^{0}$ for any
${\boldsymbol{z}}\in\mathcal{Z}^{C}$. Assume conditions (C1)-(C5) in the
Supplement hold for
$\boldsymbol{\theta}_{\boldsymbol{z}}\in(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})},\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v_{j}):\boldsymbol{d}\leq{\boldsymbol{d}}^{0},j=1,\dots,r)$
and assume that $E(\Psi_{0}(\boldsymbol{\beta}_{(\boldsymbol{d})}))=0$ has at
most $r$ solutions for $\boldsymbol{d}\leq{\boldsymbol{d}}^{0}$. Then
$\mbox{Pr}(\hat{\boldsymbol{d}\ }=\boldsymbol{d}^{0})\to 1$ as $n\to\infty$.

Thus, minimization of (S.2) leads to consistent determination of the
structural dimensions. Our result shows that despite the constraints on the
parameters, we are still able to develop a VIC-based approach to consistently
estimate the structural dimensions. The presence of these constraints and the
choice of the level of penalty make when using these constraints prevents
direct application of the results of Ma and Zhang, (2015) to our setting for
dimension determination.

###  Web Appendix A.6 Structural dimension determination for
$\widetilde{\boldsymbol{\beta}}$

Similar to $\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}$ from Section Web
Appendix A.5 of the main text, we define

$\displaystyle\widetilde{\boldsymbol{\beta}}_{(\boldsymbol{d})}=\operatornamewithlimits{argmin}_{\boldsymbol{\beta}_{(\boldsymbol{d})}=(\boldsymbol{\beta}_{{\boldsymbol{z}}(d_{\boldsymbol{z}})}:{\boldsymbol{z}}\in\mathcal{Z}^{C})}{}$
$\displaystyle\frac{1}{2n}\hat{\Psi}_{n}(\boldsymbol{\beta}_{(\boldsymbol{d})})^{\mathsf{T}}\boldsymbol{W}_{n}\hat{\Psi}_{n}(\boldsymbol{\beta}_{(\boldsymbol{d})}),$

where $\widetilde{\boldsymbol{\beta}}_{(\boldsymbol{d})}$ can be expressed as
$(\widetilde{\boldsymbol{\beta}}_{{\boldsymbol{z}}(d_{\boldsymbol{z}})}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$.
Using this, we define the non-stochastic expansions
$\widetilde{\boldsymbol{\beta}}_{{\boldsymbol{z}}(d_{\boldsymbol{z}})}(v)$ of
$\widetilde{\boldsymbol{\beta}}_{{\boldsymbol{z}}(d_{\boldsymbol{z}})}$
similarly as for $\widehat{\boldsymbol{\nu}}_{(d_{\boldsymbol{z}})}$.
Similarly as for $\widetilde{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v)$, we
define
$\widetilde{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v)=(\widetilde{\boldsymbol{\beta}}_{{\boldsymbol{z}}(d_{\boldsymbol{z}})}(v):{\boldsymbol{z}}\in\mathcal{Z}^{C})$.
The proposed criterion for structural dimension determination is then

$\displaystyle\mbox{VIC}(\boldsymbol{d})={}$
$\displaystyle\left\\{\frac{1}{rn}\sum_{j=1}^{r}\lVert\hat{\Psi}_{n}(\widetilde{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v_{j}))\rVert^{2}+\log(n)p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}d_{\boldsymbol{z}}\right\\}.$
(S.3)

To estimate the structural dimensions, one would evaluate (S.3) for all valid
sets of candidate dimensions and choose the dimensions that minimize (S.3) as
the estimated dimensions. A selection consistency result like Theorem S2 can
be shown for this minimizer in a similar fashion as the proof of Theorem S2.

###  Web Appendix A.7 Implications of the hierarchy assumption
misspecification

We now consider implications of misspecifying the hierarchical assumption.

######  Corollary S1.

Under the assumptions of Theorem LABEL:thm1, when the hierarchical assumptions
(LABEL:eqn:hierarchical_assumption_none) and (LABEL:eqn:none_assumption2) are
not satisfied, i.e.
$\boldsymbol{C}^{\mathsf{T}}\mathrm{vecl}(\boldsymbol{\beta}^{0})\neq\boldsymbol{0}$,
the solution of (S.19) can be represented as

$\mathrm{vecl}(\widehat{\boldsymbol{\beta}})=(\boldsymbol{I}-\boldsymbol{P})\mathrm{vecl}(\boldsymbol{\beta}^{0})+o_{p}(1)\text{
and}$
$n^{1/2}\big{\\{}\mathrm{vecl}(\widehat{\boldsymbol{\beta}})-(\boldsymbol{I}-\boldsymbol{P})\mathrm{vecl}(\boldsymbol{\beta}^{0})\big{\\}}\xrightarrow{\mathcal{L}}\mathcal{N}(\boldsymbol{0},(\boldsymbol{I}-\boldsymbol{P})\boldsymbol{\Sigma}_{2}\boldsymbol{\Xi}_{2}\boldsymbol{\Sigma}_{2}(\boldsymbol{I}-\boldsymbol{P})),$
(S.4)

where $\boldsymbol{P}$, ${\boldsymbol{\Sigma}_{2}}$, and
$\boldsymbol{\Xi}_{2}$ are as defined in Theorem LABEL:thm1.

Thus, there can be possible bias due to the hierarchical assumption violation.
The bias term
$\boldsymbol{P}\mathrm{vecl}(\boldsymbol{\beta}^{0})={\boldsymbol{\Sigma}_{2}}\boldsymbol{C}(\boldsymbol{C}^{\mathsf{T}}{\boldsymbol{\Sigma}_{2}}\boldsymbol{C})^{-1}\boldsymbol{C}^{\mathsf{T}}\mathrm{vecl}(\boldsymbol{\beta}^{0})$
is closely related to the hierarchical assumption equation,
$\boldsymbol{C}^{\mathsf{T}}\mathrm{vecl}(\boldsymbol{\beta}^{0})$. When the
hierarchical assumption is not correct, but holds approximately, i.e.
$\boldsymbol{C}^{\mathsf{T}}\mathrm{vecl}(\boldsymbol{\beta}^{0})\approx{\boldsymbol{0}}$,
then the potential reduction in variance from assuming subspace hierarchy may
still result in smaller mean squared error than the unconstrained estimator.

Another implication of the corollary is that even when the assumptions
(LABEL:eqn:hierarchical_assumption_none) and (LABEL:eqn:none_assumption2) do
not hold for the true central mean subspaces themselves (i.e. the smallest
dimension reduction subspace), forcing the assumption usually leads to a set
of dimension reduction subspaces that are slightly larger than the true
central mean subspaces. However with those larger subspaces, this assumption
can hold. While such expansion leads to incorrect structural dimensions, its
impact on consequent model building may not be critical as the extra
directions will likely be unimportant. We will explore this further through
numerical studies in Section LABEL:sec:sim.

##  Web Appendix B Proofs and derivations

###  Web Appendix B.1 Regularity conditions

Conditions (C1)-(C5) below are stated using a generic parameter vector
$\boldsymbol{\theta}_{\boldsymbol{z}}$. For Theorem S1, it refers to
$\boldsymbol{\beta}_{\boldsymbol{z}}^{*}$. For Theorem LABEL:thm1 and
Corollary S1 it refers to $\boldsymbol{\beta}_{{\boldsymbol{z}}}^{0}$. For
Theorem S2, it refers to elements in
$\boldsymbol{\theta}_{\boldsymbol{z}}\in\\{\boldsymbol{\beta}^{0}_{\boldsymbol{d}},\boldsymbol{\beta}^{0}_{\boldsymbol{d}}(v_{j}):\boldsymbol{d}\preceq{\boldsymbol{d}}^{0},j=1,\ldots,r\\}$.
Denote $\nu(\boldsymbol{\theta}_{\boldsymbol{z}})$ as the number of columns in
$\boldsymbol{\theta}_{\boldsymbol{z}}$.

  1. (C1)

$n_{\boldsymbol{z}}/n\to p_{\boldsymbol{z}}\in(0,1)$ as $n\to\infty$ for all
${\boldsymbol{z}}\in\mathcal{Z}^{C}$.

  2. (C2)

The univariate kernel function $K(\cdot)$ is Lipschitz and has compact
support. Furthermore for all ${\boldsymbol{z}}\in\mathcal{Z}^{C}$ it satisfies

$\displaystyle\int K(u)\mathrm{d}u=1,\int
u^{i}K(h_{\boldsymbol{z}})\mathrm{d}u=0,1\leq i\leq
m_{\boldsymbol{\theta}_{\boldsymbol{z}}}-1,$ $\displaystyle\mbox{ and
}0\neq\int
u^{m_{\boldsymbol{\theta}_{\boldsymbol{z}}}}K(h_{\boldsymbol{z}})\mathrm{d}u<\infty,$

where $m_{\boldsymbol{\theta}_{\boldsymbol{z}}}$ is an integer related to the
Lipschitz continuity of functions defined in condition (C4). The
$\nu(\boldsymbol{\theta}_{\boldsymbol{z}})$-dimensional kernel function is a
product of $\nu(\boldsymbol{\theta}_{\boldsymbol{z}})$ univariate kernel
functions, i.e.
$K_{h_{\boldsymbol{z}}}(\boldsymbol{u})=K(\boldsymbol{u}/h_{\boldsymbol{z}})/h^{\nu(\boldsymbol{\theta}_{\boldsymbol{z}})}=\prod_{j=1}^{\nu(\boldsymbol{\theta}_{\boldsymbol{z}})}K(u_{j}/h_{\boldsymbol{z}})/h^{\nu(\boldsymbol{\theta}_{\boldsymbol{z}})}$
for
$\boldsymbol{u}=(u_{1},\dots,u_{\nu(\boldsymbol{\theta}_{\boldsymbol{z}})})^{\mathsf{T}}$.
In the above notation, $K$ is used regardless of dimension.

  3. (C3)

The joint density/mass functions of $(\mathbf{X},{\boldsymbol{Z}})$ and
$(\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{Z}},{\boldsymbol{Z}})$,
respectively denoted as
$f_{\mathbf{X},{\boldsymbol{Z}}}(\mathbf{x},{\boldsymbol{z}})=f_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{x})\pi_{{\boldsymbol{Z}}}({\boldsymbol{z}})$
and
$f(\mathbf{x}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{z}},{\boldsymbol{z}})=f_{\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{Z}}|{\boldsymbol{Z}}}(\mathbf{x}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{z}})\pi_{{\boldsymbol{Z}}}({\boldsymbol{z}})$
are bounded from below and above, where
$f_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{x})$ is the density function of
$\mathbf{X}$ given ${\boldsymbol{Z}}$,
$\pi_{{\boldsymbol{Z}}}({\boldsymbol{z}})$ is the mass function of
${\boldsymbol{Z}}$, and
$f_{\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{Z}}|{\boldsymbol{Z}}}(\mathbf{x}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{z}})$
is the density function of
$\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{Z}}$ given
${\boldsymbol{Z}}$. Furthermore,
$E[Y^{2}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{z}}=\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{Z}}={\boldsymbol{z}}]$
and each entry of
$E[\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{z}})\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{z}})^{\mathsf{T}}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{z}}=\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{Z}}={\boldsymbol{z}}]$
are locally Lipschitz-continuous and bounded from above as a function of
$\boldsymbol{\xi}_{\boldsymbol{z}}$ for every
${\boldsymbol{z}}\in\mathcal{Z}^{C}$, where a function
$\boldsymbol{g}(\mathbf{x})$ from $\mathcal{S}\subset\mathbbm{R}^{k}$ to
$\mathbbm{R}^{r}$ is locally Lipschitz-continuous if for every point
$\mathbf{x}\in\mathcal{S}$, there is a constant $L>0$ such that
$\boldsymbol{g}$ is Lipschitz-continuous on the open ball
$\mathcal{B}_{L}(\mathbf{x})$ with center $\mathbf{x}$ and radius $L$. A
function $\boldsymbol{g}$ is Lipschitz continuous on $\mathcal{B}$ if there is
a constant $M>0$ such that
$\|\boldsymbol{g}(\mathbf{x})-\boldsymbol{g}(\mathbf{y})\|\leq
M\|\mathbf{x}-\mathbf{y}\|$ for all $\mathbf{x},\mathbf{y}\in\mathcal{B}$.

  4. (C4)

Let

$\boldsymbol{r}_{1}(\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{z}})=E[\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{z}})|\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{z}}=\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{Z}}={\boldsymbol{z}}]f(\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{z}})$

and

$\boldsymbol{r}_{2}(\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{z}})=E[Y|\mathbf{X}^{\mathsf{T}}\boldsymbol{\theta}_{\boldsymbol{z}}=\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{Z}}={\boldsymbol{z}}]f(\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{z}}).$

The $m_{\boldsymbol{\theta}_{\boldsymbol{z}}}$th derivatives of
$\boldsymbol{r}_{1}(\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{z}})$,
$\boldsymbol{r}_{2}(\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{z}})$, and
$f(\boldsymbol{\xi}_{\boldsymbol{z}},{\boldsymbol{z}})$ are locally Lipschitz-
continuous in $\boldsymbol{\xi}_{\boldsymbol{z}}$ for every
${\boldsymbol{z}}\in\mathcal{Z}^{C}$.

  5. (C5)

For all ${\boldsymbol{z}}\in\mathcal{Z}^{C}$ the bandwidths
$h_{\boldsymbol{z}}=O(n^{-\kappa_{\boldsymbol{\theta}_{\boldsymbol{z}}}})$ for
$1/(4m_{\boldsymbol{\theta}_{\boldsymbol{z}}})<\kappa_{\boldsymbol{\theta}_{\boldsymbol{z}}}<1/(2\nu(\boldsymbol{\theta}_{\boldsymbol{z}}))$
where $\kappa_{\boldsymbol{\theta}_{\boldsymbol{z}}}>0$ is a constant.

Condition (C1) indicates that the sample size for any subpopulation should not
vanish with respect to the overall sample size. Condition (C2) is a standard
requirement in nonparametric kernel regression estimation for multivariate
kernels and merely ensures the kernel used is valid. Conditions (C3) and (C4)
are standard smoothness and boundedness assumptions that are necessary to
ensure order of operations can be exchanged and appropriate terms can be
differentiated. Condition (C5) is a bandwidth assumption necessary to ensure
root-$n$ convergence of the parameters.

###  Web Appendix B.2 Derivation of the orthogonal complement of the nuisance
tangent space for the conditional mean

Recall the model from the main text:

$Y=\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}I({\boldsymbol{Z}}={\boldsymbol{z}})\ell_{\boldsymbol{z}}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{\boldsymbol{z}}^{*})+\epsilon,$
(S.5)

where $\epsilon$ is an error term with
$E(\epsilon|\mathbf{X},{\boldsymbol{Z}})=0$. Note that
$\ell_{{\boldsymbol{z}}}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*})$
is the mean regression function for subpopulation ${\boldsymbol{z}}$. Denote
$\eta_{\mathbf{X}|{\boldsymbol{Z}}}$ as the density of $\mathbf{X}$ given
${\boldsymbol{Z}}$, $\pi_{\boldsymbol{Z}}$ as the probability mass function
(pmf) of ${\boldsymbol{Z}}$, and $\eta_{\epsilon}$ as the density of
$\epsilon$ conditional on $\mathbf{X}$ and ${\boldsymbol{Z}}$ with respect to
some dominating measure. We use the notation $\pi$ for pmfs and $\eta$ for
density functions. We begin by writing down the likelihood of a single
observation:

$\eta_{\epsilon}\left(Y-\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}I({\boldsymbol{Z}}={\boldsymbol{z}})\ell_{{\boldsymbol{z}}}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*})\right)\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{Z}})\pi_{{\boldsymbol{Z}}}({\boldsymbol{Z}})$
(S.6)

In a similar vein as Ma and Zhu, (2012), in order to characterize all
estimation functions associated with (S.5), we seek to identify the orthogonal
complement of the nuisance tangent space in order to derive the space of all
estimating functions for estimation of
$\boldsymbol{\beta}^{*}=(\boldsymbol{\beta}_{\boldsymbol{z}}^{*}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$
in (S.5). Denote as the mean zero Hilbert space
$\mathcal{H}=\\{\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}}):E(\mathbf{f})=0,E(\mathbf{f}^{\mathsf{T}}\mathbf{f})<\infty,\mathbf{f}\text{
is measurable}\\}$ with inner product between any two functions $\mathbf{f}$
and $\mathbf{g}$ defined as $E(\mathbf{f}^{\mathsf{T}}\mathbf{g})$. The
tangent spaces corresponding to the infinite dimensional nuisance parameters
$\eta_{\mathbf{X}|{\boldsymbol{Z}}}$, $\pi_{\boldsymbol{Z}}$, and
$\eta_{\epsilon}$ are

$\displaystyle\Lambda_{\boldsymbol{Z}}=$
$\displaystyle{}\\{\mathbf{f}(\mathbf{X})\in\mathcal{H}:E(\mathbf{f})=0\\}$
$\displaystyle\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}=$
$\displaystyle{}\\{\mathbf{f}(\mathbf{X},{\boldsymbol{Z}})\in\mathcal{H}:E(\mathbf{f}|{\boldsymbol{Z}})=0\\}$
$\displaystyle\Lambda_{\epsilon}=$
$\displaystyle{}\\{\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}})\in\mathcal{H}:E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=0\text{
and }E(\mathbf{f}\epsilon|\mathbf{X},{\boldsymbol{Z}})=0\\}$

The tangent space of the unknown mean functions
$\boldsymbol{\ell}=(\ell_{\boldsymbol{z}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$
is

$\displaystyle\Lambda_{\boldsymbol{\ell}}=$
$\displaystyle\left\\{\frac{\eta_{\epsilon,1}^{\prime}(\epsilon,\mathbf{X},{\boldsymbol{Z}})}{\eta_{\epsilon}(\epsilon,\mathbf{X},{\boldsymbol{Z}})}\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}I({\boldsymbol{Z}}={\boldsymbol{z}})h_{{\boldsymbol{z}}}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*}),\forall
h_{{\boldsymbol{z}}}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*}):{\boldsymbol{z}}\in\mathcal{Z}^{C},h_{{\boldsymbol{z}}}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*}):\mathcal{X}\mapsto\mathbbm{R}\right\\}$

$\Lambda_{\epsilon}$ represented as the intersection of the following two
linear subspaces

$\displaystyle\Lambda_{\epsilon a}=$
$\displaystyle\\{\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}})\in\mathcal{H}:E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=0\\}\mbox{
and}$ (S.7) $\displaystyle\Lambda_{\epsilon b}=$
$\displaystyle\\{\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}})\in\mathcal{H}:E(\mathbf{f}\epsilon|\mathbf{X},{\boldsymbol{Z}})=0\\}.$
(S.8)

$\Lambda_{\boldsymbol{Z}}\perp\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}$ since for
$\mathbf{f}\in\Lambda_{\boldsymbol{Z}}$ and
$\mathbf{g}\in\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}$,
$E(\mathbf{f}\mathbf{g})=E(\mathbf{f}E[\mathbf{g}|{\boldsymbol{Z}}])=0$ and
similarly $\Lambda_{\boldsymbol{Z}}\perp\Lambda_{\epsilon a}$ since for
$\mathbf{f}\in\Lambda_{\boldsymbol{Z}}$ and $\mathbf{g}\in\Lambda_{\epsilon
a}$,
$E(\mathbf{f}\mathbf{g})=E(\mathbf{f}E[\mathbf{g}|\mathbf{X},{\boldsymbol{Z}}])=0$.

We want to express $\Lambda_{\epsilon b}^{\perp}$. Take
$\mathbf{f}\in\Lambda_{\epsilon b}$, then for any
$\alpha(\mathbf{X},{\boldsymbol{Z}})\in\mathcal{H}$,
$E(\mathbf{f}\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})\epsilon)=E(\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})E[\mathbf{f}\epsilon|\mathbf{X},{\boldsymbol{Z}}]))=0$.
Suppose $\mathbf{f}_{\epsilon
b}=\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}})-\epsilon\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})\in\Lambda_{\epsilon
b}$. Thus

$\displaystyle 0=$ $\displaystyle
E[\\{\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}})-\epsilon\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})\\}\epsilon|\mathbf{X},{\boldsymbol{Z}}]$
$\displaystyle=$ $\displaystyle
E(\epsilon\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}})|\mathbf{X},{\boldsymbol{Z}})-\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})E(\epsilon^{2}|\mathbf{X},{\boldsymbol{Z}}).$

Hence $\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})$ is such that
$\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})=E[\epsilon^{2}|\mathbf{X},{\boldsymbol{Z}}]^{-1}E[\epsilon\mathbf{f}(\epsilon,\mathbf{X},{\boldsymbol{Z}})|\mathbf{X},{\boldsymbol{Z}}]$.

Now define
$\mathcal{S}_{1}=\\{\epsilon\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})\mbox{
for all }\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})\text{ that map from
}\mathcal{X}\times\mathcal{Z}^{C}\\}$ and
$\mathcal{S}_{2}=\\{\boldsymbol{\gamma}(\mathbf{X},{\boldsymbol{Z}})\mbox{ for
all }\boldsymbol{\gamma}(\mathbf{X},{\boldsymbol{Z}})\in\mathcal{H}\\}$. Then
$\mathcal{S}_{1}\perp\mathcal{S}_{2}$ because
$E(\epsilon|\mathbf{X},{\boldsymbol{Z}})=0$. Hence we have showed
$\mathcal{S}_{1}=(\Lambda_{\mathbf{X}}+\Lambda_{\boldsymbol{Z}}+\Lambda_{\epsilon})^{\perp}$.
Finally, we need to ensure that
$\epsilon\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})\in\Lambda_{\boldsymbol{\ell}}^{\perp}$
in order to characterize $\Lambda^{\perp}$, since if
$\mathbf{f}\in\mathcal{S}_{1}$ and
$\mathbf{f}\in\Lambda_{\boldsymbol{\ell}}^{\perp}$, then
$\mathbf{f}\in\Lambda^{\perp}$. We use
$\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}}$ to denote
$\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}I({\boldsymbol{Z}}={\boldsymbol{z}})\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*}$.

To construct such an $\boldsymbol{\alpha}$, we note that

$\displaystyle 0=$
$\displaystyle\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}\int\epsilon\boldsymbol{\alpha}(\mathbf{x},{\boldsymbol{z}})\frac{\eta_{\epsilon,1}^{\prime}(\epsilon,\mathbf{x},{\boldsymbol{z}})}{\eta_{\epsilon}(\epsilon,\mathbf{x},{\boldsymbol{z}})}h_{{\boldsymbol{z}}}(\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{x},{\boldsymbol{z}})\pi_{{\boldsymbol{Z}}}({\boldsymbol{z}})\eta_{\epsilon}(\epsilon,\mathbf{x},{\boldsymbol{z}})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{x})\mathrm{d}\mu(\epsilon)$
$\displaystyle=$
$\displaystyle\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}\int\epsilon\boldsymbol{\alpha}(\mathbf{x},{\boldsymbol{z}})h_{{\boldsymbol{z}}}(\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{x},{\boldsymbol{z}})\left\\{\int\epsilon\eta_{\epsilon,1}^{\prime}(\epsilon,\mathbf{x},{\boldsymbol{z}})\mathrm{d}\mu(\epsilon)\right\\}\pi_{\boldsymbol{Z}}({\boldsymbol{z}})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{x})$
$\displaystyle=-E\left[\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}I({\boldsymbol{Z}}={\boldsymbol{z}})h_{\boldsymbol{z}}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*},{\boldsymbol{z}})\right],$

where $\mu_{\boldsymbol{z}}$ is the measure of $\mathbf{X}$ with respect to
population ${\boldsymbol{z}}$. Hence
$E[\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}}]=0$.
Then together, this implies

$\Lambda^{\perp}=\big{\\{}[Y-E(Y|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}})][\,\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})-E[\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}}]\,]\mbox{
for all }\boldsymbol{\alpha}\in\mathcal{H}\big{\\}}.$

###  Web Appendix B.3 Derivation of the orthogonal complement of the nuisance
tangent space for the conditional distribution

If the goal of sufficient dimension reduction is with respect to the full
conditional distribution of the response while considering population
heterogeneity based on ${\boldsymbol{Z}}$, the aim is to estimate matrices
$\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*}$ such that

$F(y|\mathbf{X},{\boldsymbol{Z}}={\boldsymbol{z}})=F(y|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*},{\boldsymbol{Z}}={\boldsymbol{z}})\text{
for }y\in\mathbb{R},$ (S.9)

where $F(y|\mathbf{X},{\boldsymbol{Z}}={\boldsymbol{z}})=Pr(Y\leq
y|\mathbf{X},{\boldsymbol{Z}}={\boldsymbol{z}})$ is the conditional
distribution function of $Y$ given $\mathbf{X}$ and ${\boldsymbol{Z}}$. The
likelihood of one observation is then

$\eta_{Y}\left(Y,\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}I({\boldsymbol{Z}}={\boldsymbol{z}})\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*},{\boldsymbol{Z}}\right)\pi_{{\boldsymbol{Z}}}({\boldsymbol{Z}})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{Z}}),$
(S.10)

where $\eta_{\mathbf{X}|{\boldsymbol{Z}}}$ is the probability mass function of
$\mathbf{X}$ given ${\boldsymbol{Z}}$, probability density function of
$\mathbf{X}$ given ${\boldsymbol{Z}}$, or a mixture distribution,
$\pi_{{\boldsymbol{Z}}}$ is the probability mass function of
${\boldsymbol{Z}}$, and $\eta_{Y}$ is the conditional probability mass/density
function of $Y$ given $\mathbf{X}$ and ${\boldsymbol{Z}}$.

The nuisance tangent spaces of $\eta_{\boldsymbol{Z}}$,
$\pi_{\mathbf{X}|{\boldsymbol{Z}}}$, and $\eta_{Y}$ are denoted as,
respectively, $\Lambda_{\boldsymbol{Z}}$,
$\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}$, and $\Lambda_{Y}$. They are

$\displaystyle\Lambda_{\boldsymbol{Z}}=$
$\displaystyle{}\\{\mathbf{f}(\mathbf{X})\in\mathcal{H}:E(\mathbf{f})=0\\}$
$\displaystyle\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}=$
$\displaystyle{}\\{\mathbf{f}(\mathbf{X},{\boldsymbol{Z}})\in\mathcal{H}:E(\mathbf{f}|{\boldsymbol{Z}})=0\\}$
$\displaystyle\Lambda_{Y}=$
$\displaystyle{}\left\\{\mathbf{f}\left(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}}\right):\text{
for all }\mathbf{f}\text{ such that }\right.$
$\displaystyle\left.E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=E\left(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}}\right)=0\right\\}.$

The following derivation largely follows Ma and Zhu, (2012). It is trivial to
see that $\Lambda_{\boldsymbol{Z}}\perp\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}$,
$\Lambda_{\boldsymbol{Z}}\perp\Lambda_{Y}$, and
$\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}\perp\Lambda_{Y}$ and thus
$\Lambda=\Lambda_{\boldsymbol{Z}}\oplus\Lambda_{\mathbf{X}|{\boldsymbol{Z}}}\oplus\Lambda_{Y}$.
Further, note that
$(\Lambda_{\boldsymbol{Z}}+\Lambda_{\mathbf{X}|{\boldsymbol{Z}}})^{\perp}=\\{\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}}):E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=0\\}\supseteq\Lambda_{Y}$.
We will now demonstrate that
$\Lambda_{Y}^{\perp}=\\{\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}}):E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)\text{
is a function of only
}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})\\}$.

We first show that all elements $\mathbf{f}\in\Lambda_{Y}^{\perp}$ must be
such that
$E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)$
is a function of only
$(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})$.
Consider any $\mathbf{f}\in\Lambda_{Y}^{\perp}$ and denote
$\mathbf{g}=E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)-E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})$.
$E(\mathbf{g}|\mathbf{X},{\boldsymbol{Z}})=E(\mathbf{g}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})=0$,
so $\mathbf{g}\in\Lambda_{Y}$, which indicates by definition that
$E(\mathbf{g}^{\mathsf{T}}\mathbf{f})=0$. Expanding on this, we have

$\displaystyle 0=E(\mathbf{g}^{\mathsf{T}}\mathbf{f})={}$ $\displaystyle
E[\mathbf{g}^{\mathsf{T}}E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)]$
$\displaystyle={}$ $\displaystyle
E(\mathbf{g}^{\mathsf{T}}\mathbf{g})+E[\mathbf{g}^{\mathsf{T}}E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})]$
$\displaystyle={}$ $\displaystyle
E(\mathbf{g}^{\mathsf{T}}\mathbf{g})+E[E(\mathbf{g}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})^{\mathsf{T}}E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})]$
$\displaystyle={}$ $\displaystyle E(\mathbf{g}^{\mathsf{T}}\mathbf{g}),$

which indicates that $\mathbf{g}=0$. The last equality on the first line holds
since $\mathbf{g}$ is a function of only
$\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y$,
the second line’s equality holds due to the definition of $\mathbf{g}$, the
third line’s equality holds by conditioning on
$\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}}$,
and the last equality holds by the tower property of conditional expectation.
Thus, it holds for $\mathbf{f}$ that
$E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)$
is a function of only
$(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})$.

We now aim to show that
$\Lambda^{\perp}=(\Lambda_{\boldsymbol{Z}}+\Lambda_{\mathbf{X}|{\boldsymbol{Z}}})^{\perp}\cap\Lambda_{Y}^{\perp}=\\{\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}})-E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y):E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})\text{
for all }\mathbf{f}\\}$. We denote the set
$\mathcal{A}=\\{\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}})-E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y):E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})\text{
for all }\mathbf{f}\\}$, which is contained in $\Lambda_{Y}^{\perp}$ by the
construction of elements in $\mathcal{A}$. Then to show that
$\mathcal{A}\subseteq(\Lambda_{\boldsymbol{Z}}+\Lambda_{\mathbf{X}|{\boldsymbol{Z}}})^{\perp}$,
we have

$\displaystyle
E[\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}})-E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)|\mathbf{X},{\boldsymbol{Z}}]$
$\displaystyle={}$ $\displaystyle
E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})-E[E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}}]$
$\displaystyle={}$ $\displaystyle
E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})-E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})=0,$

and thus
$\mathcal{A}\subseteq(\Lambda_{\boldsymbol{Z}}+\Lambda_{\mathbf{X}|{\boldsymbol{Z}}})^{\perp}$
and $\mathcal{A}\subseteq\Lambda^{\perp}$. Now it remains to show that
$\Lambda^{\perp}\subseteq\mathcal{A}$. Now take any
$\mathbf{f}\in\Lambda^{\perp}\in\Lambda_{Y}^{\perp}$. We thus have
$E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)=\mathbf{a}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})$
for some function $\mathbf{a}$. The integral representation of this equality
is

$\displaystyle\mathbf{a}(\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{z}},{\boldsymbol{z}})={}$
$\displaystyle\frac{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\mathbf{f}(Y,\mathbf{X},{\boldsymbol{z}}^{\prime})\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\eta_{Y}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\eta_{Y}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}$
$\displaystyle={}$
$\displaystyle\frac{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\mathbf{f}(Y,\mathbf{X},{\boldsymbol{z}}^{\prime})\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}.$

Using this we have

$\displaystyle\mathbf{a}(\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{z}},{\boldsymbol{z}})$
$\displaystyle=\int\mathbf{a}(\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{z}},{\boldsymbol{z}})\eta_{Y}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu(Y)$
$\displaystyle=\int\frac{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\mathbf{f}(Y,\mathbf{X},{\boldsymbol{z}}^{\prime})\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}\eta_{Y}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu(Y)$
$\displaystyle=\frac{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\mathbf{f}(Y,\mathbf{X},{\boldsymbol{z}}^{\prime})\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\eta_{Y}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})\mathrm{d}\mu(Y)}{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}$
$\displaystyle=\frac{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\int\mathbf{f}(Y,\mathbf{X},{\boldsymbol{z}}^{\prime})\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\eta_{Y}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu(Y)\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}{\sum_{{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}={\boldsymbol{z}}}\int_{\mathbf{X}:\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}=\mathbf{x}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{z}}^{\prime}}}\pi_{\boldsymbol{Z}}({\boldsymbol{z}}^{\prime})\eta_{\mathbf{X}|{\boldsymbol{Z}}}(\mathbf{X},{\boldsymbol{z}}^{\prime})\mathrm{d}\mu_{\boldsymbol{z}}(\mathbf{X})}$
$\displaystyle={}$ $\displaystyle 0$

since
$\mathbf{f}\in(\Lambda_{\boldsymbol{Z}}+\Lambda_{\mathbf{X}|{\boldsymbol{Z}}})^{\perp}$
and thus $E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=0$ which means elements in
$\Lambda^{\perp}$ have the form
$\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}})-E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}},Y)$.
Then since
$\mathbf{f}\in(\Lambda_{\boldsymbol{Z}}+\Lambda_{\mathbf{X}|{\boldsymbol{Z}}})^{\perp}$,
$0=E(\mathbf{f}|\mathbf{X},{\boldsymbol{Z}})=E(\mathbf{f}|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})$
by the law of total expectation and thus we have the second condition of
$\mathcal{A}$. Hence $\Lambda^{\perp}=\mathcal{A}$.

Similar to Ma and Zhu, (2012), given any functions
$\mathbf{g}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})$
and $\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})$ we can take
$\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}})$ to be

$[\mathbf{g}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})-E\\{\mathbf{g}(Y,\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}}\\}][\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}})-E\\{\mathbf{f}(Y,\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{\boldsymbol{Z}},{\boldsymbol{Z}}\\}],$

which has the double-robustness property described in Ma and Zhu, (2012).

###  Web Appendix B.4 Proofs

######  Lemma 1.

Assume Conditions (C1)-(C5) hold. Denote the space
$\Omega_{\boldsymbol{\beta}}=\\{(Y,\mathbf{x},{\boldsymbol{z}},\widehat{\boldsymbol{\beta}}):\mathbf{x}\in\mathbb{R}^{p},Y\in\mathbb{R},{\boldsymbol{z}}\in\mathcal{Z}^{C},\mbox{
and }||\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}||\leq
Kn^{-1/2})\\}$, where $||\cdot||$ is the Euclidean norm and $K$ is a constant.
Then there exist bases $\boldsymbol{\beta}_{{\boldsymbol{z}}}^{*}$ of
$\mathcal{S}_{\boldsymbol{z}}$ such that

$\displaystyle\sup_{\Omega_{\boldsymbol{\beta}}}\lvert\widehat{E}\\{\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\widehat{\boldsymbol{\beta}}_{\boldsymbol{Z}},{\boldsymbol{Z}}\\}-\widehat{E}\\{\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}}\\}-{E}\\{\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\widehat{\boldsymbol{\beta}}_{\boldsymbol{Z}},{\boldsymbol{Z}}\\}+$
$\displaystyle{E}\\{\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}}\\}\rvert=O_{p}\left(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}\left[n^{-1/2}h_{\boldsymbol{z}}^{m}+n^{-1}h_{\boldsymbol{z}}^{-(d_{\boldsymbol{z}}+1)}\log
n\right]\right)$ (S.11)

and similarly

$\displaystyle\sup_{\Omega_{\boldsymbol{\beta}}}\lvert\widehat{E}\\{Y|\mathbf{X}^{\mathsf{T}}\widehat{\boldsymbol{\beta}}_{\boldsymbol{Z}},{\boldsymbol{Z}}\\}-\widehat{E}\\{Y|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}}\\}-{E}\\{Y|\mathbf{X}^{\mathsf{T}}\widehat{\boldsymbol{\beta}}_{\boldsymbol{Z}},{\boldsymbol{Z}}\\}+$
$\displaystyle{E}\\{Y|\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}}^{*},{\boldsymbol{Z}}\\}\rvert=O_{p}\left(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}\left[h_{\boldsymbol{z}}^{m}n^{-1/2}+n^{-1}h_{\boldsymbol{z}}^{-(d_{\boldsymbol{z}}+1)}\log
n\right]\right)$ (S.12)

The proof of Lemma 1 is similar to that of Ma and Zhu, (2012) with adjustments
for the fact that there are differing dimensions and bandwidths for each
subpopulation. We prove Lemma 2 in full to demonstrate how this type of
modification carries through.

######  Lemma 2.

Assume Conditions (C1) - (C4) hold. Then

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[Y_{i}-E(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}^{*}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]\left[\widehat{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}^{*},{\boldsymbol{Z}}_{i})-{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}^{*},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle=O_{p}\left(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}\left[1/(nh_{\boldsymbol{z}}^{d_{\boldsymbol{z}}/2})+h_{\boldsymbol{z}}^{m}/n^{1/2}+h_{\boldsymbol{z}}^{2m}+\log^{2}n/(nh_{\boldsymbol{z}}^{d_{\boldsymbol{z}}})\right]\right)$
(S.13) and
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\widehat{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}^{*},{\boldsymbol{Z}}_{i})-E(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}^{*},{\boldsymbol{Z}}_{i})\right]\left[\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})-{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}^{*},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle=O_{p}\left(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}\left[1/(nh_{\boldsymbol{z}}^{d_{\boldsymbol{z}}/2})+h_{\boldsymbol{z}}^{m}/n^{1/2}+h_{\boldsymbol{z}}^{2m}+\log^{2}n/(nh_{\boldsymbol{z}}^{d_{\boldsymbol{z}}})\right]\right)$
(S.14)

###### Proof of Lemma 2.

The result can be shown by similar arguments to those in the proof of Lemma 4
of Ma and Zhu, (2012).

Recall from the Appendix of the main text that
$\mathbf{r}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})=E(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})f(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})$.
For any $i=1,\dots,n$, let
$\varepsilon_{i}=Y_{i}-E(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})$,
$\hat{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})=(n-1)^{-1}\sum_{j\neq
i}I({\boldsymbol{Z}}_{j}={\boldsymbol{Z}}_{i})K_{h_{{\boldsymbol{z}}_{i}}}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}-\mathbf{X}_{j}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{j}})$,
and
$\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})=(n-1)^{-1}\sum_{j\neq
i}I({\boldsymbol{Z}}_{j}={\boldsymbol{Z}}_{i})K_{h_{{\boldsymbol{Z}}_{i}}}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}-\mathbf{X}_{j}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{j}})\boldsymbol{\alpha}(\mathbf{X}_{j},{\boldsymbol{Z}}_{j})$.
Then similar to the proof of Lemma 4 in Ma and Zhu, (2012)

$\displaystyle\frac{1}{n}\varepsilon_{i}\left[\widehat{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle={}$
$\displaystyle\frac{1}{n}\varepsilon_{i}\left[\frac{\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}{\hat{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}-\frac{{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}{{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}\right]$
$\displaystyle={}$
$\displaystyle\frac{1}{n}\varepsilon_{i}\left[\frac{\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}{{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}\right]-\frac{1}{n}\varepsilon_{i}\left[\frac{{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\\{\hat{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\\}}{{f}^{2}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}\right]$
$\displaystyle-\frac{1}{n}\varepsilon_{i}\left[\frac{\\{\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\\}\\{\hat{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\\}}{{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\hat{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}\right]$
$\displaystyle+\frac{1}{n}\varepsilon_{i}\left[\frac{{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\\{\hat{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\\}^{2}}{{f}^{2}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\hat{f}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})}\right]$
(S.15)

Since $n_{\boldsymbol{z}}/n\to c_{\boldsymbol{z}}$ for all
${\boldsymbol{z}}\in\mathcal{Z}^{C}$, then
$\log(n_{\boldsymbol{z}})/\log(n)\to 1$ and by the uniform convergence of
nonparametric regression, the third and fourth terms in (S.15) are
$O_{p}(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}[h_{{\boldsymbol{z}}}^{2m}+\log^{2}(n)/(nh_{\boldsymbol{z}}^{d_{\boldsymbol{z}}})])$.
Due to (C5), it is enough to focus on the convergence rates of the quantity
$n^{-1}\sum_{i=1}^{n}\varepsilon_{i}\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})$,
which we now express as a second-order $U$-statistic.

$\frac{1}{n}\sum_{i=1}^{n}\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\varepsilon_{i}=\frac{1}{n(n-1)}\sum_{i\neq
j}^{n}I({\boldsymbol{Z}}_{i}={\boldsymbol{Z}}_{j})K_{h_{{\boldsymbol{Z}}_{i}}}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}-\mathbf{X}_{j}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{j}})\\{\varepsilon_{i}\boldsymbol{\alpha}(\mathbf{X}_{j},{\boldsymbol{Z}}_{j})+\varepsilon_{j}\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})\\}.$

An application of Lemma 5.2.1.A of Serfling, (2009, page 183) and noticing
that the following difference is a degenerate $U$-statistic yields

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\varepsilon_{i}-\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}E\left\\{K_{h_{{\boldsymbol{Z}}_{i}}}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}-\mathbf{X}_{j}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{j}}){\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i}\right\\}$
$\displaystyle=O_{p}\left(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}1/(nh_{\boldsymbol{z}}^{d_{\boldsymbol{z}}/2})\right).$
(S.16)

By similar arguments as in Ma and Zhu, (2012), we also have

$\displaystyle\frac{1}{n}\sum_{i\neq
j}^{n}\varepsilon_{i}I({\boldsymbol{Z}}_{i}={\boldsymbol{Z}}_{j})\left[E\left\\{K_{h_{{\boldsymbol{Z}}_{i}}}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}}-\mathbf{X}_{j}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{j}}){\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i}\right\\}\right.$
(S.17)
$\displaystyle\left.\vphantom{K_{h_{{\boldsymbol{Z}}_{i}}}}-{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})f(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle=O_{p}\left(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}h_{\boldsymbol{z}}^{m}/n^{1/2}\right).$
(S.18)

Combining (S.16) and (S.18), we have

$\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\left\\{\hat{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-{\mathbf{r}}_{1}(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})f(\mathbf{X}_{i}^{\mathsf{T}}\boldsymbol{\beta}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right\\}=O_{p}\left(\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}(1/(nh_{\boldsymbol{z}}^{d_{\boldsymbol{z}}/2})+h_{\boldsymbol{z}}^{m}/n^{1/2})\right),$

which, combined with (S.15), finishes the proof.

∎

###### Proof of Theorem S1.

For notational convenience, throughout this proof we refer to
$\boldsymbol{\beta}$, $\widetilde{\boldsymbol{\beta}}$, and
$\boldsymbol{\beta}^{*}$ by their vectorized versions and drop the use of
$\mathrm{vec}$ and $\mathrm{vecl}$. We begin by decomposing $\hat{\Psi}_{n}$
as

$\displaystyle n^{-1/2}\hat{\Psi}_{n}(\widetilde{\boldsymbol{\beta}})$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\sum_{i=1}^{n}\left[Y_{i}-\widehat{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]\left[\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})-\widehat{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\sum_{i=1}^{n}\left[Y_{i}-{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]\left[\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})-{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle+n^{-1/2}\sum_{i=1}^{n}\left[Y_{i}-\widehat{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]\left[{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-\widehat{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle+n^{-1/2}\sum_{i=1}^{n}\left[{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-\widehat{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]\left[\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})-{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle+n^{-1/2}\sum_{i=1}^{n}\left[{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-\widehat{E}(Y_{i}|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right]$
$\displaystyle\;\;\;\;\times\left[{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})-\widehat{E}(\boldsymbol{\alpha}(\mathbf{X}_{i},{\boldsymbol{Z}}_{i})|\mathbf{X}_{i}^{\mathsf{T}}{\widetilde{\boldsymbol{\beta}}}_{{\boldsymbol{Z}}_{i}},{\boldsymbol{Z}}_{i})\right].$

By Conditions (C1)-(C5), Lemmas 1 and 2 and similar arguments as in Ma and
Zhu, (2012), the last three terms in the above are $o_{p}(1)$. Thus
$n^{-1/2}\hat{\Psi}_{n}(\widetilde{\boldsymbol{\beta}})=n^{-1/2}{\Psi}_{n}(\widetilde{\boldsymbol{\beta}})+o_{p}(1)$.
Define

$\displaystyle{\boldsymbol{G}}(\boldsymbol{\beta})={}$ $\displaystyle
E\left\\{\frac{\partial}{\partial\boldsymbol{\beta}}\left[Y-{E}(Y|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}},{\boldsymbol{Z}})\right]\left[\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})-{E}(\boldsymbol{\alpha}(\mathbf{X},{\boldsymbol{Z}})|\mathbf{X}^{\mathsf{T}}{\boldsymbol{\beta}}_{{\boldsymbol{Z}}},\mathbf{X})\right]\right\\}$
$\displaystyle\widehat{{\boldsymbol{G}}}(\boldsymbol{\beta})={}$
$\displaystyle
n^{-1}\frac{\partial}{\partial\boldsymbol{\beta}}\hat{\Psi}_{n}(\boldsymbol{\beta}).$

Then by Conditions (C1)-(C5), the first order conditions implied by the
definition of $\widetilde{\boldsymbol{\beta}}$, and Slutsky’s theorem, we have

$\displaystyle 0={}$ $\displaystyle
n^{-1/2}\widehat{{\boldsymbol{G}}}(\widetilde{\boldsymbol{\beta}})^{\mathsf{T}}\boldsymbol{W}_{n}\hat{\Psi}_{n}(\widetilde{\boldsymbol{\beta}})$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{*})^{\mathsf{T}}+o_{p}(1)\\}\boldsymbol{W}\hat{\Psi}_{n}(\widetilde{\boldsymbol{\beta}})$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{*})^{\mathsf{T}}+o_{p}(1)\\}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{*})$
$\displaystyle+\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{*})^{\mathsf{T}}+o_{p}(1)\\}\boldsymbol{W}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{*})+o_{p}(1)\\}n^{1/2}(\widetilde{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*})+o_{p}(1)$
$\displaystyle={}$ $\displaystyle
n^{-1/2}{\boldsymbol{G}}(\boldsymbol{\beta}^{*})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{*})$
$\displaystyle+{\boldsymbol{G}}(\boldsymbol{\beta}^{*})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{*})n^{1/2}(\widetilde{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*})+o_{p}(1).$

Thus

$n^{1/2}(\widetilde{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*})=-n^{-1/2}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{*})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{*})\\}^{-1}{\boldsymbol{G}}(\boldsymbol{\beta}^{*})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{*})+o_{p}(1),$

completing the proof.

∎

###### Proof of Proposition LABEL:thm:achieving_hier_central.

We prove the result for two binary factors. The proof for more than 2 binary
factors can be done in the same fashion, though it is more tedious. Imagine
there exist parameters $\boldsymbol{\beta}^{\ddagger}_{{\boldsymbol{z}}}$ for
${\boldsymbol{z}}\in\\{00,10,01,11\\}$ of smallest dimension such that (S.5)
and (LABEL:eqn:hierarchical_assumption_none) and (LABEL:eqn:none_assumption2)
hold but that (LABEL:eqn:equality_condition) does not hold. Consider the case
that $d_{00}<d_{10},d_{01}<d_{11}$, since the re-parameterizations in
(LABEL:eqn:equality_condition) are trivial for equal dimension cases (e.g.
$d_{00}=d_{10}$). Since (LABEL:eqn:equality_condition) does not hold, there
are no other matrices $\boldsymbol{\beta}^{\dagger}_{00}$ and
$\boldsymbol{\beta}^{\dagger}_{10}$ that can be represented as
$\boldsymbol{\beta}^{\dagger}_{10}=(\boldsymbol{\beta}^{\dagger}_{00},\boldsymbol{\nu}^{\dagger}_{10\dagger})$
such that
$\text{span}(\boldsymbol{\beta}^{\ddagger}_{10})=\text{span}(\boldsymbol{\beta}^{\dagger}_{10})$
and
$\text{span}(\boldsymbol{\beta}^{\ddagger}_{00})=\text{span}(\boldsymbol{\beta}^{\dagger}_{00})$.
But then this implies
$\text{span}(\boldsymbol{\beta}^{\ddagger}_{00})\not\subseteq\text{span}(\boldsymbol{\beta}^{\ddagger}_{10})$,
implying a contradiction. Similar arguments can be made for
$\boldsymbol{\beta}^{\ddagger}_{11}$ and $\boldsymbol{\beta}^{\ddagger}_{01}$.
∎

###### Proof of Theorem LABEL:thm1.

Recall that our proposed estimator is given by the solution
$\widehat{\boldsymbol{\beta}}=({\widehat{\boldsymbol{\beta}}_{\boldsymbol{z}}}:{\boldsymbol{z}}\in\mathcal{Z}^{C})$
to

$\displaystyle\operatornamewithlimits{argmin}_{\boldsymbol{\beta}}{}$
$\displaystyle\frac{1}{2n}\hat{\Psi}_{n}(\boldsymbol{\beta})^{\mathsf{T}}\boldsymbol{W}_{n}\hat{\Psi}_{n}(\boldsymbol{\beta})$
$\displaystyle\mbox{such that
}\boldsymbol{\boldsymbol{C}}^{\mathsf{T}}\mathrm{vecl}(\boldsymbol{\beta})=\boldsymbol{0},$
(S.19)

where $\boldsymbol{W}_{n}$ is a weight matrix possibly chosen for efficiency
improvements.

Similar to the proof of Theorem S1, we drop the use of $\mathrm{vecl}$ and
$\mathrm{vec}$ and henceforth always refer to the vectorized versions of
$\widetilde{\boldsymbol{\beta}}$, $\widehat{\boldsymbol{\beta}}$, and
$\boldsymbol{\beta}^{0}$. The solution $\widehat{\boldsymbol{\beta}}$ of
(S.19) has the following first order conditions:

$\displaystyle\frac{1}{n}\left(\frac{\partial}{\partial\boldsymbol{\beta}}\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}})\right)^{\mathsf{T}}\boldsymbol{W}\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}})-\boldsymbol{C}\widehat{\boldsymbol{\lambda}}$
$\displaystyle{}=\boldsymbol{0}$ $\displaystyle\mbox{and
}\boldsymbol{C}^{\mathsf{T}}\widehat{\boldsymbol{\beta}}$
$\displaystyle{}=\boldsymbol{0}$

Denote ${\boldsymbol{G}}(\boldsymbol{\beta})$ as $\boldsymbol{G}_{2}$ but
evaluated at a given $\boldsymbol{\beta}$ and denote
$\widehat{\boldsymbol{G}}(\boldsymbol{\beta})$ as in the proof of Theorem S1
but accounting for the fact that upper identity block constraints are
different in the estimator under question in Theorem LABEL:thm1.

Similar to the proof of Theorem S1, we take a Taylor expansion of the first
order conditions around $\boldsymbol{\beta}^{0}$, multiply both the first
condition and the second first condition by $n^{1/2}$. By Conditions
(C1)-(C5), Lemmas 1 and 2, we have

$\displaystyle 0={}$ $\displaystyle
n^{-1/2}\widehat{\boldsymbol{G}}(\widehat{\boldsymbol{\beta}})^{\mathsf{T}}\boldsymbol{W}_{n}\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}})-n^{-1/2}\boldsymbol{C}\widehat{\boldsymbol{\lambda}}$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}+o_{p}(1)\\}\boldsymbol{W}\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}})-n^{-1/2}\boldsymbol{C}\widehat{\boldsymbol{\lambda}}$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}+o_{p}(1)\\}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0})$
$\displaystyle+\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}+o_{p}(1)\\}\boldsymbol{W}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})+o_{p}(1)\\}n^{1/2}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})-n^{-1/2}\boldsymbol{C}\widehat{\boldsymbol{\lambda}}+o_{p}(1)$
$\displaystyle={}$ $\displaystyle
n^{-1/2}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0})$
$\displaystyle+{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})n^{1/2}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})-n^{-1/2}\boldsymbol{C}\widehat{\boldsymbol{\lambda}}+o_{p}(1)$
and
$\displaystyle\boldsymbol{C}^{\mathsf{T}}n^{1/2}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})=o_{p}(1).$

Thus

$\displaystyle
n^{1/2}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})={}$
$\displaystyle\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}n^{-1/2}\boldsymbol{C}\widehat{\boldsymbol{\lambda}}$
$\displaystyle-\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}n^{-1/2}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0})+o_{p}(1).$
(S.20)

From
$\boldsymbol{C}^{\mathsf{T}}n^{1/2}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})$,
multiplying (S.20) by $\boldsymbol{C}^{\mathsf{T}}$, we have

$\displaystyle\boldsymbol{C}^{\mathsf{T}}n^{1/2}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{0})=o_{p}(1)={}$
$\displaystyle\boldsymbol{C}^{\mathsf{T}}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}n^{-1/2}\boldsymbol{C}\widehat{\boldsymbol{\lambda}}$
$\displaystyle+\boldsymbol{C}^{\mathsf{T}}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}n^{-1/2}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0})+o_{p}(1)$

and thus

$\displaystyle
n^{-1/2}\widehat{\boldsymbol{\lambda}}=-\left(\boldsymbol{C}^{\mathsf{T}}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}\boldsymbol{C}\right)^{-1}\boldsymbol{C}^{\mathsf{T}}n^{-1/2}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0})+o_{p}(1).$

Combining the above with (S.20) completes the proof.

∎

###### Proof of Corollary S1.

Similar to the proof of Theorem LABEL:thm1, it can be shown that

$\displaystyle\boldsymbol{C}^{\mathsf{T}}\widehat{\boldsymbol{\beta}}=\boldsymbol{0}={}$
$\displaystyle\boldsymbol{C}^{\mathsf{T}}\\{\boldsymbol{\beta}^{0}-n^{-1}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0})\\}$
$\displaystyle+\boldsymbol{C}^{\mathsf{T}}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}\boldsymbol{C}n^{-1}\widehat{\boldsymbol{\lambda}}+o_{p}(1)$

where

$n^{-1}\widehat{\boldsymbol{\lambda}}=\boldsymbol{C}^{\mathsf{T}}\widetilde{\boldsymbol{\beta}}+\\{\boldsymbol{C}^{\mathsf{T}}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}\boldsymbol{C}\\}^{-1}\boldsymbol{C}^{\mathsf{T}}\widetilde{\boldsymbol{\beta}}+o_{p}(1).$

Since
$\widehat{\boldsymbol{\beta}}=\boldsymbol{\beta}^{0}-n^{-1}[\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0})+\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0})\\}^{-1}\boldsymbol{C}\widehat{\boldsymbol{\lambda}}]$,
the bias result then holds by rearranging terms and an application of
Slutsky’s theorem. The proof is completed by plugging in the above value for
$n^{-1/2}\widehat{\boldsymbol{\lambda}}$ into (S.20) and adding
$n^{1/2}\boldsymbol{P}\boldsymbol{\beta}^{0}$ to both sides. ∎

###### Proof of Theorem S2.

Similar to the proof of Theorem S1, it can be shown that

$n^{-1/2}\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v))=n^{-1/2}{\Psi}_{n}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v))+o_{p}(1).$
(S.21)

Following the arguments of the proof of Theorem LABEL:thm1 we have

$\displaystyle
n^{1/2}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}-\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})={}$
$\displaystyle-n^{1/2}\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})\\}^{-1}{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})^{\mathsf{T}}\boldsymbol{W}{\Psi}_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})$
$\displaystyle+o_{p}(1)n^{-1/2}\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})+o_{p}(1),$
(S.22)

where ${\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})$ is defined
as $\boldsymbol{G}_{2}$ in Theorem LABEL:thm1 but allowing for differing
dimensions of $\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}$ from the true
structural dimensions. Unlike in the main text, for clarity we denote the
constituent components of $\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v)$
and $\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)$ as
$\widehat{\boldsymbol{\beta}}_{{\boldsymbol{z}}(\boldsymbol{d})}(v)$ and
$\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}(\boldsymbol{d})}(v)$, respectively,
such that the notation explicitly shows their dimensions. Similar to the proof
of Ma and Zhang, (2015), due to the definitions of
$\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v)$ and
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)$, there exist

$\left\\{(p-1-k_{{\boldsymbol{z}}})(k_{{\boldsymbol{z}}}+1)+\sum_{{\boldsymbol{z}}^{\prime}:{\boldsymbol{z}}^{\prime}\in\mathcal{Z}^{C},{\boldsymbol{z}}^{\prime}\neq{\boldsymbol{z}}}(p-k_{{\boldsymbol{z}}^{\prime}})k_{{\boldsymbol{z}}^{\prime}}\right\\}\times\left\\{\sum_{{\boldsymbol{z}}:{\boldsymbol{z}}\in\mathcal{Z}^{C}}(p-k_{{\boldsymbol{z}}})k_{{\boldsymbol{z}}}\right\\}$

matrices $\boldsymbol{L}_{(\boldsymbol{d}){\boldsymbol{z}}}(v)$ such that
$\mathrm{vecl}(\widehat{\boldsymbol{\beta}}_{{\boldsymbol{z}}(\boldsymbol{d})}(v)-\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}(\boldsymbol{d})}(v))=\boldsymbol{L}_{(\boldsymbol{d}){\boldsymbol{z}}}(v)\mathrm{vecl}(\widehat{\boldsymbol{\beta}}_{{\boldsymbol{z}}(\boldsymbol{d})}-\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}(\boldsymbol{d})})$.
Hence there exists a matrix $\boldsymbol{L}_{(\boldsymbol{d})}(v)$ such that
$\mathrm{vecl}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v)-\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))=\boldsymbol{L}_{(\boldsymbol{d})}(v)\mathrm{vecl}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}-\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})$.
Now by taking a Taylor expansion, using Conditions (C1)-(C5) with
$\boldsymbol{\theta}=\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)$, and (S.22)
above, we can express

$\displaystyle
n^{-1/2}\Psi_{n}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v))={}$
$\displaystyle
n^{-1/2}\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))+n^{1/2}{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))\mathrm{vecl}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v)-\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))+o_{p}(1)$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))+n^{1/2}{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))\boldsymbol{L}_{(\boldsymbol{d})}(v)\mathrm{vecl}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}-\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})+o_{p}(1)$
$\displaystyle={}$ $\displaystyle
n^{-1/2}\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))-n^{-1/2}\boldsymbol{M}_{(\boldsymbol{d})}(v)\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})+o_{p}(1)n^{-1/2}\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})+o_{p}(1),$
(S.23)

where ${\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))$ is
defined in the same way as
${\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})$ and

$\boldsymbol{M}_{(\boldsymbol{d})}(v)={\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)){\boldsymbol{L}}_{(\boldsymbol{d})}(v)\\{{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})^{\mathsf{T}}\boldsymbol{W}{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})\\}^{-1}{\boldsymbol{G}}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})^{\mathsf{T}}\boldsymbol{W}.$

Then gathering (S.21) and (S.23), we can express

$\displaystyle
n^{-1/2}\hat{\Psi}_{n}(\widehat{\boldsymbol{\beta}}_{(\boldsymbol{d})}(v))={}$
$\displaystyle
n^{-1/2}\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))-n^{-1/2}\boldsymbol{M}_{(\boldsymbol{d})}(v)\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})$
$\displaystyle+o_{p}(1)n^{-1/2}\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})+o_{p}(1)$
(S.24) $\displaystyle={}$ $\displaystyle
n^{-1/2}\left\\{\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))+n\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))-\boldsymbol{M}_{(\boldsymbol{d})}(v)\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})\right\\}$
$\displaystyle+o_{p}(1)n^{-1/2}\left\\{\Psi_{n}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})-n\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})\right\\}$
$\displaystyle+n^{1/2}\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))+o_{p}(1)n^{1/2}\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})+o_{p}(1).$
(S.25)

The first term in (S.25) converges to a mean zero normal random variable and
the second term is $o_{p}(1)$. However the behavior of the third and fourth
terms depends on the dimensions $\boldsymbol{d}$. When
$\boldsymbol{d}=\boldsymbol{d}^{0}$,
$\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})})$ is clearly zero and
hence the fourth term is zero. Since the terms
$\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}(\boldsymbol{d})}(v)$ are linear
combinations of $\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}(\boldsymbol{d})}$,
$\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))$ can also be shown to
be zero due to the properties of conditional expectation. Thus
$\text{VIC}(\boldsymbol{d}^{0})=\log(n)p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}d_{\boldsymbol{z}}+O_{p}(1)$.
For vectors $\boldsymbol{v}$ and $\boldsymbol{v}^{\prime}$, we say
$\boldsymbol{v}\prec\boldsymbol{v}^{\prime}$ if $v_{j}<v^{\prime}_{j}$ for all
$j$ and similarly $\boldsymbol{v}\preceq\boldsymbol{v}^{\prime}$ if $v_{j}\leq
v^{\prime}_{j}$ for all $j$. Now, when
$\boldsymbol{d}\prec\boldsymbol{d}^{0}$,
$\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))$ is not in general
$\boldsymbol{0}$ unless by chance
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)=\boldsymbol{\beta}^{0}_{(\boldsymbol{d+1})}$.
When
$\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v)\neq\boldsymbol{\beta}^{0}_{(\boldsymbol{d+1})}$,
we have
$\Psi_{0}(\boldsymbol{\beta}^{0}_{(\boldsymbol{d})}(v))=\boldsymbol{c}(v)\neq
0$ , hence the third term in (S.25) is of order $n^{1/2}$ and the fourth term
is $o_{p}(n^{1/2})$. Thus, when $\boldsymbol{d}\prec\boldsymbol{d}^{0}$,
$\text{VIC}(\boldsymbol{d})=r^{-1}n\\{\sum_{j=1}^{r}{\boldsymbol{c}(v_{j})}^{\mathsf{T}}\boldsymbol{c}(v_{j})\\}+o_{p}(n)+\log(n)p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}d_{\boldsymbol{z}}$,
which clearly dominates $\text{VIC}(\boldsymbol{d}^{0})$ as $n\to\infty$.
Simply due to the magnitude of
$p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}d_{\boldsymbol{z}}$, it is
straightforward to see that for $\boldsymbol{d}$ such that
$\boldsymbol{d}^{0}\prec\boldsymbol{d}$,
$\text{VIC}(\boldsymbol{d})>\log(n)p\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}d^{0}_{\boldsymbol{z}}$
and thus
$\text{Pr}(\text{VIC}(\boldsymbol{d})>\text{VIC}(\boldsymbol{d}^{0}))\to 1$ as
$n\to\infty$ . The result is now shown. ∎

##  Web Appendix C Additional simulation studies and results

###  Web Appendix C.1 Outcome models for simulation in main text

The outcome models used in the simulation studies in the main text are as
follows.

$\displaystyle\begin{matrix}&\textbf{Model 1}&\textbf{Model 2}&\textbf{Model
3}\\\
\ell_{00}=&2(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}&2(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}&2(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}\\\
\ell_{10}=&0.5(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{3}&\dfrac{1}{0.1+0.5(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}}-0.5(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}&0.5(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})(\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{10})^{2}\\\
\ell_{01}=&2\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00}&2(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}&\exp\\{\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00}\\}\\\
\ell_{11}=&\dfrac{\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00}}{0.5+(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00}+1.5)^{2}}&2(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}+(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})(\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{11})&5\exp\\{-(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}\\}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{10})\\\
\\\ \end{matrix}$

The conditions (LABEL:eqn:hierarchical_assumption_none) and
(LABEL:eqn:none_assumption2) are satisfied for all the three data generating
models. Throughout the simulations, all of
$\boldsymbol{\beta}^{0}_{00},\boldsymbol{\nu}^{0}_{10}$, and
$\boldsymbol{\nu}^{0}_{11}$ are $p\times 1$ vectors. Half of the elements of
each of $\boldsymbol{\beta}^{0}_{00},\boldsymbol{\nu}^{0}_{10}$, and
$\boldsymbol{\nu}^{0}_{11}$ are from a uniform distribution on $[-0.25,0.25]$
and the rest from a uniform distribution on $[-0.5,0.5]$. The true dimensions
are
$\boldsymbol{d}^{0}=(d^{0}_{00},d^{0}_{10},d^{0}_{01},d^{0}_{11})=(1,1,1,1)$
for Model 1, $\boldsymbol{d}^{0}=(1,1,1,2)$ for Model 2, and
$\boldsymbol{d}^{0}=(1,2,1,2)$ for Model 3. For $n_{\boldsymbol{z}}$ and $p$,
we consider scenarios with all combinations of
$n_{\boldsymbol{z}}\in\\{200,600,1200\\}$ and $p\in\\{10,20\\}$.

We consider the following data-generating scenario for the covariates. Like
the setup in Ma and Zhu, (2012), for ${\boldsymbol{z}}\in\\{00,10,01,11\\}$,
the first $p/2$ columns of $\mathbf{X}$ are generated from a normal
distribution with mean $0$ and variance-covariance matrix
$(\sigma_{ij})_{(p/2)\times(p/2)}$ where $\sigma_{ij}=0.5^{|i-j|}$. The next
$\lfloor p/4\rfloor$ columns $X_{p/2+i}$, $i=1,\dots,\lfloor p/4\rfloor$ are
generated from Bernoulli distributions with success probabilities
$1/\\{1+\exp(-X_{i})\\}$ where $\lfloor\cdot\rfloor$ is the floor function.
Finally, the last $\lfloor p/4\rfloor$ columns $X_{\lfloor 3p/2\rfloor+i}$,
$i=1,\dots,\lfloor p/4\rfloor$ are generated as $X_{\lfloor
3p/2\rfloor+i}=X_{\lfloor p/4\rfloor+i}^{2}+N(0,1)$. These covariates do not
meet the linearity condition or the constant variance condition required by
many SDR estimation methods.

###  Web Appendix C.2 Structural dimension determination results

To evaluate the performance of the VIC criterion for structural dimension
determination for the “hier sPHD VIC” approach, we considered all combinations
of dimensions such that the largest total dimension is less than 5. Note that
when there are multiple subpopulations, there are many possible choices of
structural dimensions that at least encompass the true underlying structural
dimensions. Hence performance in estimating the subspaces may not be harmed
much by the inclusion of an extra dimension for some subpopulations.

In using (S.2), we used $v_{1}=-1,v_{2}=-1/2,v_{3}=0,v_{4}=1/2$, and
$v_{5}=1$. This value $r=5$ and the choices for $v_{j}$ were somewhat
arbitrarily, as their particular values are not critically-important and $r=5$
works well in practice. Results corresponding to the probability of correct
selection of the structural dimensions and the rank of the VIC value of the
true structural dimensions are displayed in Figure S.1.

Figure S.1: Displayed on the top are the proportion of times over the
simulation runs that the exact set of dimensions was correctly selected.
Displayed on the bottom are average ranks of the VIC values of the true set of
dimensions among the candidate dimension possibilities.

In essentially all scenarios, the probability of correct determination of all
of the structural dimensions increases with sample size and the rank of the
VIC corresponding to the set of true structural dimensions approaches 1.

###  Web Appendix C.3 Simulation results under hierarchy misspecification

In the next set of simulations, data were generated from a modification of
Model 1 of the simulation setup in the main text. The modification is such
that the hierarchical assumptions (LABEL:eqn:hierarchical_assumption_none) and
(LABEL:eqn:none_assumption2) do not hold. We consider three scenarios of
varying degrees of hierarchy misspecification. In all scenarios we allow the
hierarchical assumption to hold approximately by setting
$\boldsymbol{\beta}^{0}_{{\boldsymbol{z}}}=\boldsymbol{\beta}^{0}_{00}+N(0,\tau^{2})$
for ${\boldsymbol{z}}\in\\{10,01,11\\}$. In Scenario 1 we set $\tau=0.05$, in
Scenario 2 we set $\tau=0.1$, and in Scenario 3 we set $\tau=0.25$, thus the
scenarios represent increasingly severe misspecification of hierarchy.

We use the proposed VIC dimension selection approach for selecting the
structural dimensions and evaluate the resulting estimation performance. The
results for the mixed discrete and continuous covariate setting in terms of
the difference norm metric and the angle metric are displayed in Figure S.2.
Since the difference norm
$\frac{1}{2^{C}}\sum_{{\boldsymbol{z}}\in\mathcal{Z}^{C}}||\widehat{\boldsymbol{\beta}}_{\boldsymbol{z}}(\widehat{\boldsymbol{\beta}}_{\boldsymbol{z}}^{\mathsf{T}}\widehat{\boldsymbol{\beta}}_{\boldsymbol{z}})^{-1}\widehat{\boldsymbol{\beta}}_{\boldsymbol{z}}^{\mathsf{T}}-\boldsymbol{\beta}^{0}_{\boldsymbol{z}}({\boldsymbol{\beta}^{0}}_{\boldsymbol{z}}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{\boldsymbol{z}})^{-1}{\boldsymbol{\beta}^{0}}_{\boldsymbol{z}}^{\mathsf{T}}||_{2}.$
and the angle between subspaces do not require the dimensions of
$\boldsymbol{\beta}^{0}_{\boldsymbol{z}}$ and
$\widehat{\boldsymbol{\beta}}_{{\boldsymbol{z}}}$ to be the same, they can be
used to evaluate estimation performance of models with estimated structural
dimensions.

In terms of the angle metric, “hier sPHD VIC” performs the best in all
scenarios even under the strongest misspecification setting, demonstrating its
capability to reconstruct the true underlying central mean subspaces even if
it requires more dimensions than the underlying truth. Despite requiring a
larger dimension, “hier sPHD VIC” is able to recover the true subspace more
accurately and with lower variance. Thus, even if the hierarchical assumption
is violated to a certain degree, “hier sPHD VIC” still performs well. However,
since the difference norm penalizes solutions with incorrect dimensions, “hier
sPHD VIC” performs worse in terms of the difference norm metric than the angle
metric.

Figure S.2: Displayed are the difference norms and angles between true and
estimated subspaces for each method over a variety of simulation settings over
250 datasets when the hierarchical assumption is misspecified. All results are
for the mixed discrete and continuous covariates setting.

Now we consider scenarios with a different form of hierarchy misspecification.
This misspecification is more severe than the previously considered hierarchy
misspecifications. All scenarios are generated under Model 1 of the main text.
In Scenario 1, we set $\boldsymbol{\beta}^{0}_{10}$ and
$\boldsymbol{\beta}^{0}_{11}$ to be equal to $\boldsymbol{\beta}^{0}_{00}$ and
left $\boldsymbol{\beta}^{0}_{01}$ to be completely different. In this
scenario the hierarchical assumption is not satisfied. In Scenario 2, we set
$\boldsymbol{\beta}^{0}_{10}$ to be equal to $\boldsymbol{\beta}^{0}_{00}$ and
left $\boldsymbol{\beta}^{0}_{01}$ and $\boldsymbol{\beta}^{0}_{11}$, to be
completely different, resulting in a more severe hierarchy misspecification
than Scenario 1. While $\boldsymbol{\beta}^{0}_{01}$ and
$\boldsymbol{\beta}^{0}_{11}$ are generated independently of each other, there
may exist some overlap of their column spaces as they are from the same
distribution and not forced to be orthogonal to each other. In Scenario 3 all
parameters are completely independent of each other, thus hierarchy does not
even hold approximately. The results in terms of the angle between estimated
and true subspaces and the difference norms are displayed in Figures S.3 and
S.4, respectively. For Scenario 3, hier sPHD VIC performs better than sPHD
when $p=20$ under the angle metric. Even though the estimated dimensions are
larger than the true structural dimensions, the subspaces are still recovered
well.

Figure S.3: Displayed are the angles between truth and estimate for each
method over a variety of simulation settings over 500 datasets when
(LABEL:eqn:hierarchical_assumption_none) and (LABEL:eqn:none_assumption2) are
misspecified.  Figure S.4: Displayed are the difference norms for each method
over a variety of simulation settings over 500 datasets when
(LABEL:eqn:hierarchical_assumption_none) and (LABEL:eqn:none_assumption2) are
misspecified.

###  Web Appendix C.4 Subpopulation specific simulation results

In this section we provide more detailed information about the simulation from
the main text. In particular, instead of computing the norm distance and the
angle averaged across the subpopulations, we consider subpopulation-specific
difference norms and angles for the simulation settings considered in the main
text with dimension $p=20$. We display the results in terms of the angles
between the estimated and true subpopulation-specific central mean subspaces
in Figure S.5. We display the same results in terms of the difference norm
metric in Figure S.6. Except for the “00” subpopulation, where the PHD
approach works best in several settings, the hier sPHD approaches with known
and estimated dimensions generally have far better estimation performance for
most subpopulations, especially for subpopulation “11” for all three models.
The sPHD approach performs reasonably well for the “01” subpopulation across
all models, but not well for subpopulation “10”, whereas the N&T approach
works reasonably well for “10” but not nearly as well for any of the other
subpopulations.

Figure S.5: Displayed are the angles between estimated and true subpopulation-
specific subspaces for each method over the simulation settings from the main
text with $p=20$ over 250 datasets under Models 1, 2, and 3 when the
covariates are generated under the all continuous setting.  Figure S.6:
Displayed are the difference norms comparing the estimated and true
subpopulation-specific subspaces for each method over the simulation settings
from the main text with $p=20$ over 250 datasets under Models 1, 2, and 3 when
the covariates are generated under the all continuous setting.

###  Web Appendix C.5 Additional simulation results with all continuous
covariates

The results in this section mirror the simulations from the main text except
that covariates are all continuous. Figure S.7 displays the corresponding
estimation results. The results are consistent with these when some covariates
are discrete and thus we omit further discussion.

Figure S.7: Displayed are the difference norms and angles between estimated
and true subspaces for each method over a variety of simulation settings over
250 datasets under Models 1, 2, and 3 when the covariates are generated under
the all continuous setting. All metrics are averages over all subpopulations.

###  Web Appendix C.6 Simulation with more complex heterogeneity and outcome
models

Our simulation setup in this experiment is identical to the setup of the
simulation experiments in the main text, however in this section we explore an
additional outcome model that exhibits additional complexity in the structural
dimensions and outcome models. In this setting, there is more heterogeneity
across subpopulations. The model is displayed below.

$\displaystyle\begin{matrix}&\textbf{Model 4}\\\
\ell_{00}=&2(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}\\\
\ell_{10}=&0.5(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})(\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{10})^{2}+5\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{10}\\\
\ell_{01}=&\exp\\{\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00}\\}+(\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{01})^{2}\\\
\ell_{11}=&5\exp\\{-(\mathbf{X}^{\mathsf{T}}\boldsymbol{\beta}^{0}_{00})^{2}\\}(\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{10})+5\sin(2\pi\mathbf{X}^{\mathsf{T}}\boldsymbol{\nu}^{0}_{01})\\\
\\\ \end{matrix}$

In this model, the true structural dimensions are
$\boldsymbol{d}^{0}=(1,2,2,3)$. The average signal-to-noise ratio under this
model is around 6.5. The results over 250 simulation replications are
displayed below. From Figure S.8, we can see that the structural dimensions
are recovered frequently and that the correct recovery probability tends to 1
with increasing sample size. From Figure S.9, which displays the results in
terms of the angle and difference norm averaged over all the subpopulations,
we can see that the proposed approach works quite well with both true and
estimated structural dimensions and outperforms other approaches. The
subpopulation-specific results are displayed in Figures S.10 and S.11. From
these results, we can see that while the hier sPHD method does not perform
best for the “none” subpopulation, it is substantially better than other
approaches for all sample sizes for the other three subpopulations.

Figure S.8: Displayed on the top are the proportion of times over the
simulation runs that the exact set of dimensions was correctly selected.
Displayed on the bottom are average ranks of the VIC values of the true set of
dimensions among the candidate dimension possibilities.  Figure S.9: Displayed
are the difference norms and angles between estimated and true subspaces for
each method over a variety of simulation settings over 250 datasets under
Model 4 when the covariates are generated under the mixed discrete and
continuous setting. All metrics are averages over all subpopulations. Figure
S.10: Displayed are the angles between estimated and true subpopulation-
specific subspaces for each method over the simulation settings from the main
text with $p=20$ over 250 datasets under Model 4 when the covariates are
generated under the all continuous setting.  Figure S.11: Displayed are the
difference norms comparing the estimated and true subpopulation-specific
subspaces for each method over the simulation settings from the main text with
$p=20$ over 250 datasets under Model 4 when the covariates are generated under
the all continuous setting.

## References

  * Li et al., (2003) Li, B., Cook, R. D., and Chiaromonte, F. (2003).  Dimension reduction for the conditional mean in regressions with categorical predictors.  Annals of Statistics, 31(5):1636–1668.
  * Ma and Zhang, (2015) Ma, Y. and Zhang, X. (2015).  A validated information criterion to determine the structural dimension in dimension reduction models.  Biometrika, 102:409–20.
  * Ma and Zhu, (2012) Ma, Y. and Zhu, L. (2012).  A semiparametric approach to dimension reduction.  Journal of the American Statistical Association, 107(497):168–179.
  * Ma and Zhu, (2014) Ma, Y. and Zhu, L. (2014).  On estimation efficiency of the central mean subspace.  Journal of the Royal Statistical Society. Series B, 76(5):885–901.
  * Serfling, (2009) Serfling, R. J. (2009).  Approximation theorems of mathematical statistics, volume 162.  John Wiley & Sons.