High-dimensional sliced inverse regression with endogeneity

Sufficient dimension reduction (SDR, Ma and Zhu,, 2013; Li,, 2018) refers to a class of methods that assume the outcome $y\in\mathbb{R}$ depends on a set of covariates $\mathbf{X}\in\mathbb{R}^{p}$ through a small number of linear combinations $\bm{\beta}_{1}^{\top}\mathbf{X},\ldots,\bm{\beta}_{d}^{\top}\mathbf{X}$, with $d\ll p$. These linear combinations are known as sufficient predictors and retain the full regression information between the response and covariates. Since the pioneering work ofLi, (1991), a vast literature has grown on different SDR approaches, ranging from inverse-moment-based and regression-based methods (e.g., Cook and Forzani,, 2008; Bura et al.,, 2022), forward regression (e.g., Xia et al.,, 2002; Dai et al.,, 2024), and semiparametric techniques (e.g., Zhou and Zhu,, 2016; Sheng and Yin,, 2016). Among the various SDR methods, sliced inverse regression (SIR, Li,, 1991) remains the most widely used due to its simplicity and computational efficiency, and as such continues to be the subject of much research (see Tan et al.,, 2018; Xu et al.,, 2022, among many others). When the number of covariates $p$ is large, it is often (further) assumed each linear combination depends only on a few covariates, i.e., the vectors $\bm{\beta}_{1},\ldots,\bm{\beta}_{d}$ are sparse. Recent literature has thus sought to combine SIR and other SDR methods with variable selection to achieve sparse sufficient dimension reduction (e.g., Lin et al.,, 2018, 2019; Nghiem et al., 2023c, ).

Much of the literature dedicated to SIR and many other SDR methods, including those mentioned above, has focused on the setting where all covariates are exogenous, and little has been done to tackle sufficient dimension reduction with endogeneity. While there is some related literature on single index models with endogenity (Hu et al.,, 2015; Gao et al.,, 2020), their approach requires approximation of the link function and is not easy to generalize to multiple index models. More specifically, letting $\mathbf{B}=[\bm{\beta}_{1},\ldots,\bm{\beta}_{d}]$ denote a $p\times d$ matrix, we can formulate SDR as a multiple index model

where $f$ denotes an unknown link function, and $\varepsilon\in\mathbb{R}$ is random noise. In (1), the target of estimation is the column space of $\mathbf{B}$, also known as the central subspace. Almost all current SDR methods assume the random noise $\varepsilon$ is independent of $\mathbf{X}$, implying $y\perp\mathbf{X}\mid\mathbf{B}^{\top}\mathbf{X}$ where ‘$\perp$’ denotes statistical independence. However, this assumption may not be reasonable in many settings where the random noise $\varepsilon$ is dependent upon $\mathbf{X}$ i.e., the covariates are endogenous.

We discuss two examples, corresponding to two data applications in Section 6, where such endogeneity can arise.

• •

  Omitted-variable bias occurs when one or more (truly important) covariate is not included in the model. A classic example relates to studying the effect of education on future earnings. Here, it is well-known future earnings depend on many other variables that are also highly correlated with education level, such as family background or an individual’s overall abilities (Card,, 1999). Such variables are usually not fully known or not collected, so the ordinary least squares fit of a regression model of earnings on education levels and other observed factors produces a biased estimate for the true effect of education (Card,, 1999). Similarly, in Section 6.1, we consider a genomic study that relates a biological outcome (e.g., body weight) to gene expressions of mice. Outside of gene expressions, there can be potentially unmeasured confounding factors, such as those related to environmental conditions, which influence both gene expression levels and the outcome.

• •

  Covariate measurement error arises when some or all covariates are not accurately measured. In the second data application in Section 6.2, we consider a nutritional study relating an individual’s health outcomes to their long–term diet, collected through 24–hour recall interviews that are subject to measurement errors (Carroll et al.,, 2006). Let $\mathbf{X}$ denote the long-term dietary covariates. The observed data can then be modeled by $\mathbf{W}=\mathbf{X}+\mathbf{E}$, with $\mathbf{E}$ denoting an additive measurement error independent of $\mathbf{X}$. In this situation, even when in the true model (1) it holds that $\varepsilon\perp\mathbf{X}$, a naive analysis which replaces $\mathbf{X}$ by $\mathbf{W}$ results in endogeneity, since we can write $y=f(\mathbf{B}^{\top}\mathbf{W}+\mathbf{B}^{\top}\mathbf{E},\varepsilon)$ and $\mathbf{E}$ is correlated with $\mathbf{W}$. This causes naive estimators that ignore measurement errors to be biased and inconsistent.

As an aside, note that while the omitted variables lead to endogeneity because too few covariates are included in the model, Fan and Liao, (2014) argued that endogeneity also arises when a large number of covariates are included in the model i.e., $\mathbf{X}$ is high-dimensional in (1). Specifically, for large $p$ there is a high probability at least one covariate is correlated with the random noise. For instance, genome-wide association studies typically include thousands of gene expression covariates to model a clinical outcome, and at least one gene expression is likely to correlate with random noise reflecting experimental conditions or environmental perturbations (Lin et al.,, 2015).

In this paper, we show under fairly general conditions that endogeneity also invalidates sliced inverse regression, leading to inconsistent estimates of the true column space of $\mathbf{B}$. To address this issue, we propose a two-stage Lasso SIR method to estimate the central space in the presence of endogeneity when valid instrumental variables are available. In the first stage, we construct a sparse multivariate linear model regressing all the covariates against the instruments and then obtain fitted values from this model. In the second stage, we apply a sparse sliced inverse regression between the response and the fitted values from the first stage, to perform simultaneous SDR and variable selection.

Addressing endogeneity via instrumental variables has been extensively studied in the econometrics and statistics literature; see Angrist and Krueger, (2001) for a review. Classical instrumental variable estimators are mostly studied with a linear model for the outcome, i.e., $y=\mathbf{X}^{\top}\bm{\beta}_{0}+\varepsilon$, and often computing a two-stage least squares estimator or a limited information maximum likelihood estimator. Hansen et al., (2008) showed in this setting that these estimators are consistent for $\bm{\beta}_{0}$ when the number of instruments $q$ increases slowly with the sample size. More recently however, a number of methods have combined instrumental variables with regularization techniques on the covariates and/or instrumental variables to improve their performance in high-dimensional settings e.g., the Lasso-type generalized methods of moments (Caner,, 2009; Caner et al.,, 2018), focused generalized methods of moments (Fan and Liao,, 2014), and two-stage penalized least squares Lin et al., (2015); the last of these methods inspires this article. Another line of research has focused on the use of sparsity-inducing penalties to estimate optimal instruments from a large set of a-priori valid instruments (see for instance Belloni et al.,, 2012; Chernozhukov et al.,, 2015), or to screen out invalid instruments (e.g., Windmeijer et al.,, 2019). To the best of our knowledge though, none of these methods allow for the presence of non-linear link functions in the model for the response e.g., the multiple index model (1), and more broadly have not been adapted to the SDR context.

Under a setting where both the number of covariates $p$ and the number of instrumental variables $q$ can grow exponentially with sample size $n$, we derive theoretical bounds on the estimation and selection consistency of the two-stage Lasso SIR estimator. Simulation studies demonstrate the proposed estimator exhibits superior finite sample performance compared with a single stage Lasso SIR estimator when endogeneity is present but ignored, and the two-stage regularized estimator proposed by Lin et al., (2015) when the outcome model is non-linear in nature. We apply the proposed method to two real data applications studying the effects of different nutrients on the cholesterol level (covariate measurement error), and for identifying genes that are relevant to the body weight of mice (omitted variable bias).

The remainder of this article is structured as follows. In Section 2, we briefly review two versions of sliced inverse regression and prove that endogeneity causes SIR to be inconsistent for estimating the target subspace. In Section 3, we formulate the proposed two-stage Lasso estimators, and derive their theoretical properties in Section 4. Section 5 presents simulation studies to compare the performance of the proposed estimator with existing methods, while Section 6 discusses two data applications Finally, Section 7 offers some concluding remarks.

In this section, we first establish the setting of interest, namely a high-dimensional multiple index model for the response along with a high-dimensional linear instrumental variable model for the covariates characterizing endogeneity. We then review SIR and a version augmented with the Lasso penalty (Section 2.1), and then demonstrate inconsistency of SIR under endogeneity in Section 2.2.

Before proceeding, we set out some notation to be used throughout the paper. For any matrix $\mathbf{A}$ with elements $a_{ij}$, let $\|\mathbf{A}\|_{1}=\max_{j}\sum_{i}|a_{ij}|$, $\|\mathbf{A}\|_{\infty}=\max_{i}\sum_{j}|a_{ij}|$, and $\left\lVert\mathbf{A}\right\rVert_{\max}=\max_{i,j}|a_{i,j}|$ denote the matrix $1$-norm, $\infty$-norm, and element-wise maximum norm, respectively. Also, let $\|\mathbf{A}\|_{2}$ denote the $\ell_{2}$-operator norm, i.e. the largest singular value of $\mathbf{A}$. For any vector $\mathbf{v}$ with elements $v_{i}$, let $\|\mathbf{v}\|_{2}=\left(\sum_{i}v_{i}^{2}\right)^{1/2}$ and $\|\mathbf{v}\|_{1}=\sum_{i}|v_{i}|$ and $\|\mathbf{v}\|_{\infty}=\max_{i}|v_{i}|$ denote its Euclidean norm, $\ell_{1}$ norm, and $\ell_{\infty}$ norm, respectively. Finally, for an index set $I$, let $\mathbf{v}_{I}$ denote the subvector of $\mathbf{v}$ formed with the elements of $\mathbf{v}$ whose indices $i\in I$, $|I|$ denote its cardinality, and $I^{c}$ denote its complement.

Consider a set of $n$ independent observations $\{(y_{i},\mathbf{X}_{i}^{\top},\mathbf{Z}_{i}^{\top});i=1,\ldots,n\}$ where, without loss of generality, we assume $E(\mathbf{X}_{i})=E(\mathbf{Z}_{i})=\bm{0}$. Each triplet $(y_{i},\mathbf{X}_{i}^{\top},\mathbf{Z}_{i}^{\top})$ follows the multiple index model (1) coupled with the linear instrumental variable model. That is, the full model is given by

where $\mathbf{B}\in\mathbb{R}^{p\times d}$, $\bm{\Gamma}\in\mathbb{R}^{q\times p}$, $\varepsilon_{i}\in\mathbb{R}$ denotes the random noise for the outcome model, and $\mathbf{E}_{i}\in\mathbb{R}^{p}$ denotes the random noise for the instrumental variable model. The vector $(\mathbf{U}_{i}^{\top},\varepsilon_{i})^{\top}\in\mathbb{R}^{p+1}$ is assumed to have mean zero and finite covariance matrix $\bm{\Sigma}$, whose diagonal elements are denoted by $\sigma_{j}^{2}$ for $j=1,\ldots,p+1$. To characterize endogeneity between $\mathbf{X}$ and $\varepsilon$, the off-diagonal elements in the last row and column of $\bm{\Sigma}$ are generally non-zero, while the instrumental variable $\mathbf{Z}_{i}$ satisfies $\mathbf{Z}_{i}\perp(\mathbf{U}_{i}^{\top},\varepsilon_{i})^{\top}$ for all $i=1,\ldots,n$.

We seek to estimate the column space of $\mathbf{B}$, when the dimensions $p$ and $q$ are growing exponentially with the sample size $n$. A common assumption in such settings is that sparsity holds. Specifically: i) each covariate $X_{j}$ depends upon only a few instrumental variables, such that in the instrumental variable model the $j$th column of $\bm{\Gamma}$ is a sparse vector with at most $r_{j}<q$ non-zero elements, where $r=\max_{j=1,\ldots,p}r_{j}$ is bounded, and ii) the multiple index model is sparse such that $y$ depends only on $s<p$ covariates i.e., at most $s$ rows of $\mathbf{B}$ are non-zero row vectors. Equivalently, if $\mathcal{P}(\mathbf{B})=\mathbf{B}(\mathbf{B}^{\top}\mathbf{B})^{-1}\mathbf{B}^{\top}$ denotes the corresponding projection matrix, then only $s$ diagonal elements of $\mathcal{P}(\mathbf{B})$ are non-zero. Without loss of generality, we assume $S=\{1,\ldots,s\}$ i.e., the first $s$ covariates indexes this set of non-zero elements, or, equivalently, the support set of the column space of $\mathbf{B}$. We let $\mathcal{X},\mathcal{Z}$ and $\mathbf{y}$ denote the $n\times p$ matrix of covariates, $n\times q$ matrix of instrumental variables, and $n\times 1$ vector of responses, whose $i$th rows are given by $\mathbf{X}_{i}^{\top}$, $\mathbf{Z}_{i}^{\top}$ and $y_{i}$, respectively. Let $\mathbf{x}_{j}$ and $\mathbf{z}_{j}$ denote the $j$th column of $\mathcal{X}$ and $\mathcal{Z}$, respectively. For the remainder of this article, we use $y,\mathbf{X}$, and $\mathbf{Z}$ to denote the random variables from which $(y_{i},\mathbf{X}_{i}^{\top},\mathbf{Z}_{i}^{\top})$ is drawn.

To overcome the general inconsistency of the SIR estimator when endogeneity is present, we propose a two-stage approach for SDR and variable selection inspired by the classical two-stage least-squares estimator (Hansen et al.,, 2008). Let $\widehat{\mathbf{X}}=E(\mathbf{X}\mid\mathbf{Z})=\bm{\Gamma}^{\top}\mathbf{Z}$, the conditional expectation of the covariates on the instruments. We first impose a new linearity condition to ensure the conditional expectation $E(\widehat{\mathbf{X}}\mid y)$ behaves similarly to $E(\mathbf{X}\mid y)$ when no endogeneity is present.

Note the covariance matrix $\bm{\Sigma}_{\widehat{\mathbf{X}}}$ need not be of full-rank and invertible: this occurs if the number of endogenous covariates $p$ is greater than the number of instruments $q$. Nevertheless, Lemma 1 motivates us to propose a two-stage Lasso SIR (abbreviated to 2SLSIR) estimator for the column space of $\mathbf{B}$, assuming its entries are sparse, based on the observed vector $\mathbf{y}$, covariate matrix $\mathcal{X}$, and instrument matrix $\mathcal{Z}$. As we shall see, while theoretically we need to understand the behavior of its inverse, computationally $\bm{\Sigma}_{\widehat{\mathbf{X}}}$ does not have to be inverted to construct the 2SLSIR estimator.

The two stages of the 2SLSIR estimator are detailed as follows.

• •

  Stage 1: Identify and estimate the effects of the instruments $\mathbf{Z}$, and subsequently obtain the predicted values of the covariates $\mathbf{X}$. To accomplish this, we fit a multivariate instrumental variable model given by the second equation in (2), with the covariates $\mathbf{X}$ being the $p$-dimensional response and instruments $\mathbf{Z}$ being the $q$-dimensional predictor. We assume each endogenous covariate depends only on a few instrumental variables (Belloni et al.,, 2012; Lin et al.,, 2015), thus motivating a sparse first-stage estimator $\hat{\bm{\Gamma}}$ that can be computed in parallel by computing the estimated $j$th column as

  $\displaystyle\hat{\bm{\gamma}}_{j}=\arg\!\min_{\bm{\gamma}}\dfrac{1}{2n}\|\mathbf{x}_{j}-\mathcal{Z}\bm{\gamma}\|_{F}^{2}+\mu_{1j}\|\bm{\gamma}\|_{1},\;j=1,\ldots,p,$

  where $\mu_{1j}>0$ denotes the tuning parameter for the $j$th covariate in the first stage. After the estimator $\hat{\bm{\Gamma}}=[\hat{\bm{\gamma}}_{1},\ldots,\hat{\bm{\gamma}}_{p}]$ is formed, we construct the first-stage prediction matrix $\widehat{\mathcal{X}}=\mathcal{Z}\widehat{\bm{\Gamma}}$.

• •

  Stage 2: Compute a sparse estimator for the column space of $\mathbf{B}$ by applying the Lasso SIR estimator (3) with the response vector $\mathbf{y}$ and the first-stage prediction matrix $\hat{\mathcal{X}}$, Specifically, let $\widehat{\bm{\Lambda}}=(nc)^{-1}\widehat{\mathcal{X}}^{\top}\mathbf{D}\widehat{\mathcal{X}}=(nc)^{-1}\widehat{\bm{\Gamma}}^{\top}{\mathcal{Z}}^{\top}\mathbf{D}\mathcal{Z}\widehat{\bm{\Gamma}}$, and define the pseudo-response vector $\tilde{\bm{y}}_{k}=(\hat{\lambda}_{k})^{-1}\mathbf{D}\widehat{\mathcal{X}}\hat{\bm{\eta}}_{k}$, where $\hat{\bm{\eta}}_{k}$ denotes the eigenvector associated with the $k$th largest eigenvalue $\hat{\lambda}_{k}$ of $\hat{\bm{\Lambda}}$, and $\mathbf{D}$ is defined as in Section 2.1. Then the 2SLSIR estimator for the $k$th column of $\mathbf{B}$ is given by

  $$\hat{\bm{\beta}}_{k}=\arg\!\min_{\bm{\beta}}\dfrac{1}{2n}\|\tilde{\bm{y}}_{k}-\widehat{\mathcal{X}}\bm{\beta}\|_{2}^{2}+\mu_{2k}\|\bm{\beta}\|_{1},\;k=1,\ldots,d,$$

  (4)

  where $\mu_{2k}>0$ is a tuning parameter for the $k$th pseudo-response in the second stage.

We highlight that, computationally, both stages can be implemented straightforwardly and in parallel using standard coordinate descent algorithms. In this paper, we apply the R package glmnet (Friedman et al.,, 2010) to implement the 2SLSIR estimator, and select the tuning parameters $\mu_{j}$ and $\alpha_{k}$ via default standard ten-fold cross-validation as available in the software.

We study the properties of the 2SLSIR estimator for the column space of $\mathbf{B}$ under the multiple index model (2), when the dimensionality of the covariates $p$ and instruments $q$ can grow exponentially with sample size $n$. Compared with the theoretical results for the linear outcome model in Lin et al., (2015) say, the major challenge here lies in the nonlinear nature of the conditional expectation $E(\widehat{\mathbf{X}}\mid y)$, and subsequently the behavior of the eigenvectors corresponding to the covariance matrix $\bm{\Lambda}=\operatorname{\text{Cov}}\left\{E(\widehat{\mathbf{X}}\mid y)\right\}$ which require more complex technical developments. We will use $C$, $C^{\prime}$, $C^{\prime\prime}$ to denote generic positive constants. Also, without loss of generality, we assume each column of the instrumental matrix $\mathcal{Z}$ is standardized such that $\mathbf{z}_{j}^{\top}\mathbf{z}_{j}=1$ for $j=1,\ldots,q$.

Following the notation in Wainwright, (Chapter 7, 2019), we say an $n\times m$ matrix $\mathbf{A}$ satisfies a restricted eigenvalue (RE) condition over a set $S$ with parameters $(\kappa,\nu)$ if $\|\mathbf{A}^{1/2}\bm{\delta}\|_{2}^{2}\geq\kappa\|\bm{\delta}\|_{2}^{2},$ for all $\bm{\delta}\in\mathcal{C}(S,\nu)=\left\{\delta\in\mathbb{R}^{m}:~{}\|\bm{\delta}_{S^{c}}\|_{1}\leq\nu\|\bm{\delta}_{S}\|_{1}\right\}$. Also, for a vector $\mathbf{v}$ let $\text{supp}(\mathbf{v})=\{j:v_{j}\neq 0\}$ be the support set of $\mathbf{v}$. Let $\sigma_{\max}=\max_{1\leq j\leq p}\sigma_{j}$, where we recall $\sigma^{2}_{j}$ are the diagonal elements of the covariance matrix $\bm{\Sigma}$ corresponding to $(\mathbf{U}_{i}^{\top},\varepsilon_{i})^{\top}\in\mathbb{R}^{p+1}$. We assume each $\sigma_{j}^{2}$ and hence $\sigma_{\max}$ to be a constant with respect to sample size.

The following technical conditions on the distribution of the instruments $\mathbf{Z}$ and the conditional expectation $E(\mathbf{X}\mid\mathbf{Z})=\bm{\Gamma}^{\top}\mathbf{Z}$ are required.

The sub-Gaussian condition is often imposed in high-dimensional data analysis (e.g. Tan et al.,, 2018; Wainwright,, 2019; Nghiem et al., 2023b, ). Note a direct consequence of Condition 4.1 is that the instrumental matrix $\mathcal{Z}$ satisfies the RE condition with parameter $(\kappa,3)$ with high probability when $n>Cr\log q$ for some finite $\kappa$ (Banerjee et al.,, 2014); we will condition on this event in our development below. Furthermore, letting $\widehat{\bm{\Sigma}}_{\mathbf{Z}}=\mathcal{Z}^{\top}\mathcal{Z}/n$ denote the sample covariance of $\mathcal{Z}$ data, then

with probability at least $1-\exp(-C^{\prime}\log q)$ (see Ravikumar et al.,, 2011; Tan et al.,, 2018, for analogous results). Condition 4.2 is also imposed in the theoretical analysis of the two-stage penalized least square methods in Lin et al., (2015), and is necessary to bound the estimation error from the first stage. Finally, in the theoretical analysis below we allow the number of covariates $p$ and instruments $q$, and the sparsity levels $s$ and $r$ to depend on the sample size $n$, subject to the rate of growth of these parameters dictated by Condition 4.3. As a special case, if $r$ and $s$ are bounded then $p$ and $q$ can grow exponentially with $n$.

We conducted a numerical study to evaluate the performance of the proposed 2SLSIR estimator for recovering the (sparsity patterns of the) central subspace. We compared this against two alternatives: the sparse Lasso SIR (LSIR) estimator of Lin et al., (2018, 2019) which ignores endogeneity and does not use the instruments, and, in the case of single index model $(d=1)$, the two-stage Lasso (2SLasso) estimator of Lin et al., (2015) which assumes linearity and additive errors in the response model i.e., $y_{i}=\mathbf{X}_{i}^{\top}\mathbf{B}+\varepsilon_{i}$.

We simulated data $\{(y_{i},\mathbf{X}_{i},\mathbf{Z}_{i});i=1,\ldots,n\}$, from a linear instrumental variable model $\mathbf{X}_{i}=\bm{\Gamma}^{\top}\mathbf{Z}_{i}+\mathbf{U}_{i}$ coupled with one of the following four models for the response:

These single and double index models are similar to those previously used in the literature to assess SDR methods (e.g., Lin et al.,, 2019). The true sparse matrix $\mathbf{B}=\bm{\beta}_{1}$ in the single index models (i) and (ii), and $\mathbf{B}=[\bm{\beta}_{1},\bm{\beta}_{2}]$ in the multiple index models (iii) and (iv), was constructed by randomly selecting $s=5$ nonzero rows of $\mathbf{B}$ and generating each entry in these rows from a uniform distribution on the two disjoint intervals $[-1,-0.5]$ and $[0.5,1]$, which we denote as $U([-1,-0.5]\cup[0.5,1])$. All remaining entries are set to zero. We assumed the number of dimensions $d$ was known throughout the simulation study.

For the instrumental variable model, we followed the same generation mechanism as in Lin et al., (2015). Specifically, to represent different levels of instrument strength, we considered three choices for $\bm{\Gamma}$: i) all instruments are strong, where we sampled $r=5$ nonzero entries for each column from the uniform distribution $U([-b,-a]\cup[a,b])$ with $(a,b)=(0.75,1)$; ii) all instruments are weak, where we sampled the same number of nonzero entries but with $(a,b)=(0.5,0.75)$; iii) instruments are of mixed strength, such that for each column of $\bm{\Gamma}$ we sampled $r=50$ nonzero entries comprising five strong/weak instruments with $(a,b)=(0.5,1)$ and 45 very weak instruments with $(a,b)=(0.05,0.1)$. Next we generated the random noise vector $(\mathbf{U}_{i}^{\top},\varepsilon_{i})^{\top}$ from a $(p+1)$-variate Gaussian distribution with zero mean vector and covariance matrix $\bm{\Sigma}=(\sigma_{ij})$ specified as follows. We set $\sigma_{ij}=(0.2)^{|i-j|}$ for $i,j=$ $1,\ldots,p$, and $\sigma_{p+1,p+1}=1$. For the last column/row of $\bm{\Sigma}$ characterizing endogeneity, the nonzero entries were then given by the nonzero row of $\mathbf{B}$, with five of these entries randomly set to $0.3$. Finally, we generated the instruments $\mathbf{Z}_{i}$ from a $q$-variate Gaussian distribution with zero mean vector and two choices for the covariance matrix $\bm{\Sigma}_{\mathbf{Z}}$, namely the $q\times q$ identity matrix and the AR(1) covariance matrix with autocorrelation $0.5$. In Supplementary Material S4, we conducted additional simulation results where the instruments $\mathbf{Z}$ were generated from a Bernoulli distribution, with largely similar trends to those seen below.

We set $p=q=100$ and varied the sample size $n\in\{500,1000\}$, simulating $100$ datasets for each configuration. For each simulated dataset, all tuning parameters across the three methods (2SLSIR, LSIR, 2SLasso) were selected via ten-fold cross-validation. We assessed performance using the estimation error on the projection matrix $\|\mathcal{P}(\mathbf{B})-\mathcal{P}(\hat{\mathbf{B}})\|_{F}$, along with the F₁ score for variable selection.

Across both single and double index models (Tables 1 and 2, respectively), the proposed 2SLSIR estimator performed best overall in terms of both estimation error and variable selection. Indeed, it is the only estimator among the three whose estimation error consistently decreased when sample size increased from $n=500$ to $n=1000$; this lends empirical support to estimation consistency results established in Section 4. For the single index models (i) and (ii), the two-stage Lasso (2SLasso) estimator performed even worse than the single stage Lasso SIR (LSIR) estimator in terms of estimation error, particularly when the degree of non-linearity in the response model was relatively large.