[table]capposition=top

# Deep Partial Least Squares for Instrumental Variable Regression

Maria Nareklishvili
Booth School of Business
University of Chicago
Nicholas Polson111Email<EMAIL_ADDRESS>
Booth School of Business
University of Chicago
Vadim Sokolov
Department of Systems Engineering
and Operations Research
George Mason University

###### Abstract

In this paper, we propose deep partial least squares for the estimation of
high-dimensional nonlinear instrumental variable regression. As a precursor to
a flexible deep neural network architecture, our methodology uses partial
least squares for dimension reduction and feature selection from the set of
instruments and covariates. A central theoretical result, due to Brillinger,
(2012), shows that the feature selection provided by partial least squares is
consistent and the weights are estimated up to a proportionality constant. We
illustrate our methodology with synthetic datasets with a sparse and
correlated network structure and draw applications to the effect of
childbearing on the mother’s labor supply based on classic data of Angrist and
Evans, (1996). The results on synthetic data as well as applications show that
the deep partial least squares method significantly outperforms other related
methods. Finally, we conclude with directions for future research.

Keywords: Dimensionality Reduction, Deep Learning, Instrumental Variables,
Partial Least Squares.

## 1\. Introduction

Nonlinear instrumental variable (IV) regression is a vital tool for estimating
the causal effect of exposure on certain outcomes. In fact, virtually all
legitimate techniques for causal inference can be seen as manifestations of
instrumental variables, including but not limited to randomized clinical
trials (considered a perfect instrument), intention-to-treat analysis
(stemming from random incentive allocation), natural experiments, and
regression discontinuity (Heckman,, 1995). A valid instrumental variable meets
the following assumptions: it has a significant effect on the treatment (i.e.,
relevance). It does not influence the outcome directly, through channels other
than the treatment (i.e., the exclusion restriction), and a valid instrument
is not associated with the unobserved characteristics that affect the outcome
(i.e., exogeneity). Under these assumptions, the IV approach recovers
consistent coefficients of interest. Traditionally, the techniques to estimate
the causal parameter in the IV approach are well-suited for low and middle-
scale data (White,, 1982; Terza et al.,, 2008). These techniques most often
fail with high-dimensional instruments (McCullagh and Polson,, 2018) and
struggle to extract meaningful insights when instruments sparsely relate to
each other. Several papers investigate the approach with many, possibly
correlated and sparse instrumental variables (Belloni et al.,, 2011;
Chernozhukov et al., 2015b, ; Chernozhukov et al., 2015a, ).

For predictive tasks, contemporary research has predominantly centered on
narrow deep neural networks (Hartford et al.,, 2016; Polson and Sokolov,,
2017; Liu et al.,, 2020; Beise et al.,, 2021), which are distinguished by
their “self-featurizing” basis functions. In other words, feature extraction
and dimension reduction are integrated into the pattern-matching algorithm.
However, there are limitations to this approach, including a lack of
theoretical understanding, difficulty in quantifying uncertainty, and limited
capacity for probabilistic reasoning. To enhance the efficiency of deep neural
networks, some authors introduce linear methods as a precursor to implementing
a deep learning model. Debiolles et al., (2004); Jia et al., (2016) use
partial least squares (PLS) to extract features before learning the output
with deep neural networks. PLS is a linear method that selects relevant
dimensions that are particularly relevant for predicting the outcome.

To mitigate the issues inherent to high-dimensional instrumental variables, we
propose deep partial least squares for IV regression. One can view our
approach as merging two cultures of machine learning tools and statistical
inference (Brillinger,, 2012; Sarstedt et al.,, 2022; Monis et al.,, 2022).
Specifically, we employ linear techniques, such as partial least squares to
identify meaningful instruments and form predictors that are ensemble averages
of deep learners. We show that, under several data-generating process
assumptions, the method exhibits desirable large sample properties.

A traditional approach to estimating the IV regression is a two-stage least-
squares (2SLS) technique. In the first stage, valid instrumental variables
predict the treatment. In the second stage, we estimate the effect of the
predicted treatment (and possibly other control variables) on the outcome of
interest (Mogstad et al.,, 2021). When the relation of the treatment and
instruments is nonlinear, Terza et al., (2008) show that the 2SLS method is
inconsistent. Iwata, (2001) uses a re-centered and re-scaled outcome variable
and proves the consistency of the method for censored and truncated data.
Terza et al., (2008) and Adkins, (2012) use predicted treatment residuals in
addition to the treatment and other exogenous variables in the second-stage
regression. Controlling for residuals in the outcome regression has a long
history (Florens et al.,, 2002; Heckman and Navarro-Lozano,, 2004; Navarro,,
2010). The methods introduced by Terza et al., (2008); Adkins, (2012) can also
be identified as control function methods.

Our goal in this article is to extend the deep partial least squares (DPLS)
method to extract relevant instruments in the first stage. More importantly,
when the observed outcome variables consist of errors, we show that the method
is consistent under a re-centered and re-scaled outcome described by Iwata,
(2001) or after controlling for the predicted residuals in the second stage
regression as illustrated by Heckman and Navarro-Lozano, (2004). The first
layer of DPLS consists of a partial least squares method to extract features
via hyperplanes in a high-dimensional setting. The subsequent layers
efficiently post-process them based on ReLU networks as a deep learner
(Yarotsky,, 2017). PLS is an attractive method for feature extraction.
However, as an alternative, the generalized method of moments (GMM, Hansen and
Singleton,, 1982) can also be applied to our feature selection stage to
provide variance stabilized estimators. This results in a statistical
improvement on traditional stochastic descent estimators that are commonplace
in machine learning (see, for example, Iwata and Ghahramani,, 2017).

Several related articles apply instrumental variables to address challenges
and research questions in business and industry. Abadie et al., (2002) examine
the effect of job training program participation on earnings. By comparison,
Chernozhukov and Hansen, (2004) investigate the effect of 401(k) participation
on wealth. Moreno and Terwiesch, (2014) mine text sentiment from user-
generated comments on an online service platform, and subsequently estimate
the impact of (predicted) sentiment on buyers’ purchasing decisions. Hartford
et al., (2017) estimate the effect of the airline ticket price on demand.
Another related article (Cengiz et al.,, 2022) builds a boosted tree model to
identify minimum wage workers based on their demographics, and then examines
the effect of minimum wage policies on labor market outcomes for these
workers. In this article, we revisit the effect of childbearing on women’s
labor supply (Angrist and Evans,, 1996). The instruments consist of the second
and third kids being twins and the interactions with parental characteristics
that are relevant for predicting the mother’s labor supply.

In addition, we conduct two simulation experiments to demonstrate the
predictive performance of DPLS in IV regression. The aim of these experiments
is to mimic high-dimensional and sparse data, most often prevalent in
business, finance, and economics. In the first experiment, we consider a high-
dimensional IV space, where some of the instruments are redundant for
predicting the treatment. The instruments are uncorrelated with each other. In
the second experiment, we extend the first one by introducing a network
structure among the instruments, as proposed by Jeong et al., (2003). In this
setup, the instruments are sparsely correlated, and some of them are not
related to the treatment. Our aim is to illustrate whether DPLS can identify
relevant instruments and improve prediction accuracy in both scenarios.

Statistical properties of deep neural networks are sparse but growing. Polson
and Ročková, (2018) discuss theoretical foundations of sparse rectified-
linear-unit (ReLU) networks and the advantage of using Spike-and-Slab prior as
an alternative to Dropout. They show that the resulting posterior prediction
of ReLU networks with Spike-and-Slab regularization converge to a true
function at a rate of $log^{\delta}(n)/n^{-K}$ for $\delta>1$ and a positive
constant $K$ (with $n$ number of observations). Polson and Sun, (2019) provide
a theoretical connection between the Spike-and-Slab priors and $L_{0}$ norm
regularization. They demonstrate that the regularized estimators can result in
improved out-of-sample prediction performance. To emphasize the advantages of
regularization in deep neural networks, Bhadra et al., (2020) discuss
theoretical and empirical justifications (as well as challenges) for the
horseshoe prior.

By comparison, Hanin and Rolnick, (2019) derive the sharp upper bound for the
number of activation function regions in ReLU neural networks. They find that
this number in practice is far from the maximum possible and depends on the
number of neurons in the network, rather than the depth. A recent paper by
Farrell et al., (2021) provides groundbreaking results on the asymptotic
theory of deep neural networks. Under the assumption that the number of hidden
layers, i.e., the depth of the network, grows with the sample size, they
provide a high probability convergence rate for ReLU neural networks.

Our paper sheds light on the dimension reduction importance in sparse ReLU
networks. Specifically, DPLS provides a simple, interpretable framework for
modeling instrumental variables when the policy (treatment) and outcome are
measured with errors. DPLS consists of a system of equations that can
theoretically be viewed as an infinite sequential generalization of 2SLS.
Beyond the consistency of this method, we demonstrate that the shrinkage in
the first layer can substantially improve prediction performance in an
instrumental variable regression with many instruments.

The results based on simulated experiments and the application show that DPLS
significantly outperforms other related methods. Specifically, the first stage
prediction performance improvement of DPLS relative to OLS, LASSO
(Tibshirani,, 2011), and DeepIV (Hartford et al.,, 2016) is more than 79%, 57%
and 22%, respectively. We find that a deep, multi-layered structure of DPLS
significantly increases predictive performance and representation learning
ability. To incorporate the uncertainty of the model coefficients, we also
consider the extension of DPLS to a Bayesian framework and discuss
implications.

One area for future research is the study of full uncertainty quantification.
It is well known that posterior distributions for IV regression require
careful assessment dating back to Zellner, (1975) (see also Hoogerheide et
al.,, 2007; Lopes and Polson,, 2010). Puelz et al., (2017) discuss parameter
uncertainty and variable selection in linear factor models. Hahn et al.,
(2020) provide a fully Bayesian model for posterior uncertainty for causal
inference using Bayesian Additive Regression Trees (BART) to model the
nonlinearity in the outcome equation. This provides a gold standard for
comparison to a Bayesian DPLS method.

The rest of the paper is outlined as follows. Section 2 describes our general
nonlinear IV model and specific Tobit variation. In section 3 we discuss
dimensionality reduction with partial least squares. In Section 4 we introduce
deep partial least squares and examine the asymptotic theory. Section 5
illustrates the applications of the method. Finally, Section 6 concludes with
directions for future research.

## 2\. Nonlinear IV Model

One goal in this article is to predict the outcome $y^{\star}$ that possibly
nonlinearly depends on a policy (treatment) $p^{\star}$ and predictors $x$ for
each observation $i=1\dots,n$:

$\displaystyle y^{\star}=p^{\star}\beta+x\beta_{x}+u.$ (2.1)

We assume, neither $y^{\star}$ nor $p^{\star}$ are directly observable.
Instead, they consist of errors. Define the following variables

$\displaystyle(y^{\star},y)\in\mathbb{R}=$ $\displaystyle\;\;{\rm
the\;potential\;and\;observed\;outcome\;variables},$
$\displaystyle(p^{\star},p)\in\mathbb{R}=$ $\displaystyle\;\;{\rm
the\;potential\;and\;observed\;policy\;variables},$ $\displaystyle
x\in\mathbb{R}^{k}=$ $\displaystyle\;\;{\rm observable\;features},$
$\displaystyle z\in\mathbb{R}^{m}=$ $\displaystyle\;\;{\rm
instrumental\;variables},$ $\displaystyle(u,w,v,\varepsilon)\in\mathbb{R}=$
$\displaystyle\;\;{\rm
latent/error\;variables\;that\;affect\;}(y^{\star},p^{\star},p,y),\;{\rm
respectively,}$ $\displaystyle m+k<$ $\displaystyle n.$

Instead of (2.1), we have access to the following structural equation model:

$\displaystyle p=p^{\star}+v,$ $\displaystyle
p^{\star}=g(z\alpha+x\alpha_{z})+w,$ (2.2) $\displaystyle
y=\tau(y^{\star})+\varepsilon.$

$g(\cdot)$ and $\tau(\cdot)$ are potentially non-linear continuous functions,
possibly deep learners (Hartford et al., 2016). $\tau(\cdot)$ is a known
transformation. In this study, we define the Tobit model
$\tau(y^{\star})=1(y^{\star}>0)\cdot y^{\star}$. The errors $v$ and
$\varepsilon$ are assumed to be uncorrelated with any other latent error.
However, policy $p^{\star}$ is allowed to correlate with $u$. Covariates $x$
are typically independent of $w$, $v$ and $u$. $\beta$ and $\beta_{x}$
represent the effects of $p^{\star}$ and $x$ on $y$, while $\alpha$ and
$\alpha_{z}$ are the effects of instruments z and covariates $x$ on the policy
$p^{\star}$, correspondingly. The errors $w$ and $u$ are correlated. For
example, consider, $y$ is a customer’s decision to buy an airline ticket
(observed for a particular group of customers), and $p^{\star}$ is the price
of this ticket. In that case, policy $p^{\star}$ is said to be endogenous
when, conditional on $x$, $(p^{\star},u)$ correlate. i.e,
$\mathbb{E}(u|x,p^{\star})\neq 0$. For example, ticket prices might increase
during conferences that are unobservable to a researcher, and in that case,
$\mathbb{E}(p^{\star}u|x)\neq 0$. Classical estimation methods, such as OLS,
will lead to a spurious positive relation between prices and sales.

Define the joint error $\eta=w+v$. Moreover, from (2) we see that
$p^{\star}=p-v$. In that case, (2.1) becomes

$\displaystyle y^{\star}=p\beta+x\beta_{x}+\xi,$ (2.3)

where $\xi=u-v\beta$. Then by combining (2.3) and (2), we get:

$\displaystyle y$ $\displaystyle=\tau(p\beta+x\beta_{x}+\xi)+\varepsilon,$
(2.4) $\displaystyle p$ $\displaystyle=g(z\alpha+x\alpha_{z})+\eta,$

where the latent errors $\xi$ and $\eta$ are correlated through $w$. The
structural equation model in (2.4) defines the nonlinear IV framework. In a
traditional nonlinear IV setup (when $\tau(\cdot)$ is an identity function),
the presence of valid instruments $z$ that satisfy Assumptions 1-3 allows us
to predict the unbiased mean outcome.

######  Assumption 1 (Relevance).

Instruments $z$ strongly relate to policy $p$, i.e., the density of $p$,
$F(p|z,x)$, is not constant in $z$.

######  Assumption 2 (Exclusion Restriction).

$z$ is conditionally orthogonal to outcome $y$:

$\displaystyle z\perp y|\ (x,p,\xi).$

######  Assumption 3 (Exogeneity).

$z$ is conditionally orthogonal to the latent error term $\xi$:

$\displaystyle z\perp\xi|x.$

Specifically, when $y=p\beta+x\beta_{x}+\epsilon$ with
$\epsilon=\xi+\varepsilon$, valid instruments efficiently separate information
in $p$ that is unrelated to $\xi$, and result in the consistent estimate of
the effect of the policy on the outcome ($\beta$). For example, if the fuel
cost represents an instrument for ticket prices, and is unrelated to
conferences, it can recover the exogenous variation in ticket prices and
correctly estimate the negative relation between prices and ticket sales.

Under Assumptions 1-3, a standard approach to estimating $\beta$ is a 2SLS
method. The method entails predicting $\hat{p}$ in the first stage (“treatment
network”). Then the predicted policy $\hat{p}$ replaces $p$ in the outcome
equation, and we estimate a second stage regression (“outcome network”) with
another consistent method (Angrist and Evans,, 1996; Hartford et al.,, 2016;
Mogstad et al.,, 2021).

Define the predicted mean outcome:

$\displaystyle\mathbb{E}(y|x,z)=\mathbb{E}[f(p,x)|x,z]+\mathbb{E}[\xi|x].$
(2.5)

Then we can evaluate the effect of a marginal change in policy (e.g., prices)
from $p_{0}$ to $p_{1}$ on the outcome of interest (treatment effect):

$\displaystyle\mathbb{E}(y|p_{1},x)-\mathbb{E}(y|p_{0},x)=\mathbb{E}[f(p_{1},x)|x,z]-\mathbb{E}[f(p_{0},x)|x,z].$

The issue is that $\tau(.)$ typically is not an identity function. In such
nonlinear regression models, latent errors are no longer additively separable
from the true regressors ($\hat{p}$ and $u$ are still related), and hence, the
true relationship breaks down with errors in variables. The orthogonality
condition of the instruments and the outcome error is violated:
$\mathbb{E}\big{(}z(y-p\beta)|x\big{)}\neq 0$. As a result, the 2SLS estimator
fails to be consistent for nonlinear errors-in-variables models (Amemiya,,
1985). Moreover, when instruments and/or covariates are close to the number of
observations and nonlinearly relate to the policy and outcomes, OLS is no
longer an efficient solution of (2.4).

The goal of this paper is to predict outcome $y^{\star}$ with the policy that
depends on many instruments $z$ and regressors $x$. To do so, we extend deep
partial least squares with a recentered and rescaled outcome (Iwata,, 2001)
and additionally illustrate it with a control function approach (Terza et
al.,, 2008).

### A. Predicton in Low-Dimensional Data

In this section, our attention is directed towards predicting the outcome
using a Tobit model, wherein $\tau(y^{\star})=1(y^{\star}>0)y^{\star}$. We
demonstrate that by recentering and rescaling the outcome, two-stage least
squares (2SLS) can yield a consistent estimator of $\beta$. Alternatively, the
control function approach can provide a valid prediction of the outcome.

####  2SLS with the Tobit Model

To establish consistency, we impose an essential assumption regarding the
underlying process responsible for generating the observed data.

######  Assumption 4 (Elliptical Distribution).

$(y^{\star},z)$ have a joint elliptical (or Gaussian in the simplest setting)
distribution.

Assumption 4 requires that the instrumental variables and the outcome are
jointly elliptically distributed. The assumption can be strong in many
settings. Though Brillinger, (2012) and Iwata, (2001) show that the results
are robust to significant deviations from the assumption. This assumption can
also be relaxed and generalized to other distributions (see e.g. Adcock,,
2007).

Define the recentered and rescaled output variable that mimics the unobserved
outcome $y^{\star}$,

$\tilde{y}=\psi_{1}^{-1}(y-\psi_{2}),\;\;{\rm
where}\;\;\psi_{1}=\mathrm{cov}(y,y^{\star})/\mathrm{var}(y^{\star}),\psi_{2}=\mathbb{E}(y)-\psi_{1}\mathbb{E}(y^{\star}).$

Specifically, $\psi_{1}$ is the ratio of the covariance between $y$ and
$y^{\star}$ to the variance of $y^{\star}$, while $\psi_{2}$ is the intercept
term when we consider a linear projection of $y$ on $y^{\star}$. Next,
consider a linear projection of z on $y^{\star}$

$\displaystyle\mathbb{E}(z|y^{\star})=\mathbb{E}(z)+\frac{\mathrm{cov}(z,y^{\star})}{\mathrm{Var}(y^{\star})}\big{[}y^{\star}-\mathbb{E}(y^{\star})].$

Hence,

$\displaystyle\mathbb{E}(z\mid
y^{\star})-\mathbb{E}(z)=\gamma(y^{\star}-\mathbb{E}(y^{\star})),$ (2.6)

where $\gamma=\mathrm{cov}(z,y^{\star})/\mathrm{var}(y^{\star})$ 222Elliptical
contours allow for discrete outcomes.. By Stein’s lemma (Landsman and
Nešlehová,, 2008), we can calculate the covariance of instruments and rescaled
output as

$\displaystyle\mathrm{cov}(z,\tilde{y})=\psi_{1}^{-1}cov(z,y)=\psi_{1}^{-1}\mathbb{E}_{y^{\star}}\left(\mathbb{E}\big{[}(z-\mathbb{E}(z))y\mid
y^{\star}\big{]}\right).$ (2.7)

Conditional expectations are linear under elliptical contours. Additionally,
$\mathbb{E}(y\mid y^{\star})=\tau(y^{\star})$,
$\gamma=\mathrm{cov}(z,y^{\star})/var(y^{\star})$, and replacing the last
equality from (2.7) with (2.6) gives

$\displaystyle\mathrm{cov}(z,\tilde{y})$
$\displaystyle=\psi_{1}^{-1}\mathbb{E}_{y^{\star}}\left(\mathbb{E}(\gamma(y^{\star}-\mathbb{E}(y^{\star}))y\mid
y^{\star}\right)$
$\displaystyle=\gamma\psi_{1}^{-1}\mathbb{E}((y^{\star}-\mathbb{E}(y^{\star}))y)=\gamma\psi_{1}^{-1}\mathrm{cov}(y^{\star},y)$
$\displaystyle=\psi_{1}^{-1}\frac{cov(z,y^{\star})}{var(y^{\star})}\mathrm{cov}(y,y^{\star})=\psi_{1}^{-1}cov(z,y^{\star})\psi_{1}=cov(z,y^{\star}).$

For additional details, see Iwata, (2001). 2SLS first regresses $p$ on $z$ and
$x$ to get the predicted policy $\hat{p}$. Define $\bar{Z}=[z,x]$ and
$\hat{p}=P_{z}\bar{Z}$ where
$P_{z}=\bar{Z}(\bar{Z}^{T}\bar{Z})^{-1}\bar{Z}^{T}$ is the projection matrix
onto the instrumental variable space (including covariates). In addition,
define $\bar{P}=[\hat{p},x]$ and let $e=\tilde{y}-p\beta$ be the residual from
the re-scaled regression, then

$\mathrm{cov}(z,e|x)=\mathrm{cov}(z,\tilde{y}-p\beta|x)=\mathrm{cov}(z,y^{\star}-p\beta|x)=\mathrm{cov}(z,u|x)=0.$

The re-scaled instrumental variable estimator is given by

$\displaystyle\hat{\beta}^{GMM}=(\bar{Z}^{T}P_{z}\bar{Z})^{-1}P_{z}\bar{Z}^{T}\tilde{y}=(\bar{P}^{T}\bar{P})^{-1}\bar{P}^{T}\tilde{y}$
(2.8)

From the above, we then have $\mathbb{E}(\hat{\beta}^{GMM})=\beta_{p}$ with
$\beta_{p}=[\beta,\beta_{x}]$ and $\mathrm{E}(\tilde{y})=\bar{P}\beta_{p}$.
Note that $\hat{\beta}^{GMM}$ is a $(m+k)\times 1$ vector, where the first
element is the policy effect on outcome.

####  Estimating the Proportionality Constants

In the Tobit case, the constants $\phi_{1},\phi_{2}$ can be calculated
theoretically from the model as

$\phi_{1}=\Phi\left(\delta\right)\;\;{\rm
and}\;\;\phi_{2}=\sigma^{\star}_{y}\phi(\delta),$

where $\Phi\left(\delta\right)$ and $\phi(\delta)$ are the cumulative and
probability density functions of a normally distributed random variable,
respectively, and

$\delta=\mu_{X}^{T}\beta/\sigma^{\star}_{y}\;\;{\rm
and}\;\;(\sigma^{\star}_{y})^{2}=\sigma^{2}_{y^{\star}}+\beta^{T}\Sigma_{xx}\beta.$

Iwata, (2001) shows that, with the truncated outcome variable, we can obtain
their estimators:

$\displaystyle\hat{\psi}_{1}=\hat{\Phi}=\dfrac{1}{n}\sum_{i=1}^{n}\mathds{1}(y>0),$
$\displaystyle\hat{\psi}_{2}=\hat{\sigma}^{\star}_{y}\hat{\phi}\;\;{\rm
where}\;\hat{\phi}=\phi\big{(}\Phi^{-1}(\hat{\psi}_{1})\big{)}.$

The estimator of the variance of the outcome variable is given as

$\hat{\sigma}_{y}^{\star
2}=\dfrac{1}{n\;c(\hat{k})}\sum_{i=1}^{n}(y_{i}-\bar{y})^{2},$

where

$c(\hat{k})=\hat{\Phi}-(\hat{\phi}-\Phi^{-1}(\hat{\Phi})(1-\hat{\Phi}))(\hat{\phi}+\Phi^{-1}(\hat{\Phi})\hat{\Phi}),$

and $\hat{\Phi}=\hat{\psi}_{1}$ and $\phi$ are the cumulative and probability
density functions of a normally distributed random variable, respectively.
$\mathds{1}(y>0)$ is a binary variable and equals one if the outcome is
positive. Iwata, (2001) shows that the scaling constant $\widehat{\psi}_{1}$
and the estimator $\hat{\beta}^{GMM}$ defined in (2.8) are asymptotically
jointly normally distributed, with

$\displaystyle\sqrt{n}(\hat{\psi}_{1}-\psi_{1})\rightarrow\mathcal{N}\big{(}0,\hat{\psi}_{1}(1-\hat{\psi}_{1})\big{)},$
(2.9)
$\displaystyle\sqrt{n}(\hat{\beta}^{GMM}-\beta)\rightarrow\mathcal{N}\big{(}0,\Sigma_{\star}-\Phi(1-\Phi)\beta\beta^{T}\big{)}.$
(2.10)

Based on Hansen and Singleton, (1982), $\Sigma_{\star}$ is the variance of the
standard GMM estimator and can be computed as follows 333See also Baum et al.,
(2003) for the discussion of a GMM estimator for the IV regression. :

$\displaystyle\widehat{\Sigma}_{\star}$
$\displaystyle=n(\bar{P}^{T}\bar{Z}\widehat{G}\bar{Z}^{T}\bar{P})^{-1}\bar{P}^{T}\bar{Z}\widehat{G}(n\widehat{A})\widehat{G}\bar{Z}^{T}\bar{P}(\bar{P}^{T}\bar{Z}\widehat{G}\bar{Z}^{T}\bar{P})^{-1},$
$\displaystyle\widehat{A}$
$\displaystyle=\frac{1}{n}\sum_{i}\widehat{e}^{2}_{i}\bar{Z}_{i}\bar{Z}_{i}^{T},\
\ \text{with}\ \
\widehat{e}^{2}=\widehat{\tilde{y}}-\bar{P}\hat{\beta}^{GMM},$
$\displaystyle\widehat{G}$
$\displaystyle=\widehat{H}\big{(}\widehat{H}^{-1}-\bar{Z}^{T}\bar{Z}/\bar{Z}^{T}\widehat{H}\bar{Z}\big{)}\widehat{H},$
$\displaystyle\widehat{H}$ $\displaystyle=\widehat{A}^{-1}.$

### B. Control Function Approach

An alternative approach to estimate the parameters is a control function
approach (Terza et al.,, 2008). Instead of substituting the predicted policy
in the second stage regression, we control for the predicted residuals of the
first stage:

$\displaystyle\widehat{\eta}$ $\displaystyle=p-\widehat{p},$ (2.11)
$\displaystyle y$
$\displaystyle=\tau(p\beta+\widehat{\eta}\beta_{\eta})+\varepsilon^{2SRI},$
(2.12)

where the residuals of the treatment network $\widehat{\eta}$ can be predicted
by any consistent method. Note that $\varepsilon^{2SRI}$ is not identical to
$\varepsilon$ due to the substitution of $\xi$ with $\widehat{\eta}$. The
control function approach can be viewed as a special case of 2SLS.
Specifically, the inclusion of $\widehat{\eta}$ in the outcome equation of
(2.4) allows us to control for the correlation of $\eta$ and $\xi$, and
predict the outcome.

## 3\. Dimensionality Reduction

The estimation of $\hat{\beta}^{GMM}$ in (2.8) rests on the assumption that
the number of instruments (as well as independent variables) is strictly lower
than the number of observations. Throughout the proof, we maintain this
assumption.

######  Assumption 5.

The number of covariates is strictly smaller than the number of observations,
$m+k<n$, where $m$ is the dimension of the instruments and $k$ is the number
of covariates.

Nevertheless, when the dimension of the independent variables is close to $N$,
$\hat{\beta}^{GMM}$ captures the unwanted variation reflected in the predicted
policy $\hat{p}$ (Gagnon-Bartsch et al.,, 2013). Since $\bar{Z}$ are
independently and identically distributed elliptical random variables with
covariance $\Sigma_{zz}$, Brillinger, (2012) shows that

$\displaystyle\text{cov}(\bar{Z},p)=\text{cov}(\bar{Z},f(U))=\text{cov}(\bar{Z},U)\text{cov}(f(U),U)/\text{var}(U)=k\Sigma_{zz}\alpha_{z},$
(3.1)

where $\bar{Z}=[z,x]$, $U=\bar{Z}\alpha_{z}$ with
$\alpha_{z}=[\alpha,\alpha_{x}]$, and the constant
$k=\text{cov}(f(U),U)/\text{var}(U)$. Based on this result, Brillinger, (2012)
shows that the OLS coefficient is a consistent estimate of $\alpha$ up to a
proportionality constant. In this article, we show that the deep partial least
squares method has the same property.

### C. Partial Least Squares

Partial least squares is a dimensionality reduction method that generalizes
and combines features from the principal component analysis and multiple
regression (Abdi and Williams,, 2010).

Consider the augmented instrumental variable $\bar{Z}=[z,x]$ with the
dimension $n\times(m+k)$, where $m$ and $k$ are the dimensions of instruments
and covariates, respectively. The policy $p=P$ is a $n\times 1$ vector as
before. PLS can be summarized by the following relationship:

$\displaystyle\bar{Z}=TV+F,$ (3.2) $\displaystyle P=UQ+E,$

where $T$ and $U$ are $n\times L$ projections (scores) of $\bar{Z}$ and $P$,
respectively. $V$ and $Q$ are orthogonal projection matrices (loadings).
Maximizing the covariance between the augmented instruments and the policy
leads to the first PLS projection pair $(v_{1},q_{1})$:

$\displaystyle\max_{v,q}\ (\bar{Z}v_{1})^{T}(Pq_{1})$
$\displaystyle\text{subject to}\ ||v_{1}||=||q_{1}||=1.$

The corresponding scores are $t_{1}=\bar{Z}v_{1}$ and $u_{1}=Pq_{1}$. It is
clear that the directions (loadings) for the policy $P$ and the augmented
instruments $\bar{Z}$ are the right and left singular vectors of
$\bar{Z}^{T}P$, respectively. PLS in the next step performs an ordinary
regression of $U$ on $T$, namely $U=T\beta$. Then the next projection pair
$(v_{2},q_{2})$ is found by calculating the singular vectors of the residual
matrix $(\bar{Z}-t_{1}v_{1}^{T})^{T}(P-T\beta q^{T})$. Lastly, the final
regression of interest is $U=T\beta$ (the tutorial Geladi and Kowalski, (1986)
contains further details).

The key property of PLS is that it is consistent for estimating parameters of
interest even in the presence of nonlinearity via a sequence of covariance
calculations. Fisher, (1922) first observed this in the Probit regression.
Naik and Tsai, (2000) show the consistency of PLS based on the result of
Brillinger, (2012). One desirable property of PLS is that it has a closed-form
solution. The PLS estimator of $\alpha$ is given as follows (Helland,, 1990;
Stone and Brooks,, 1990):

$\displaystyle\widehat{\alpha}^{PLS}=\hat{R}(\hat{R}^{T}S_{zz}\hat{R})^{-1}\hat{R}^{T}s_{zp},$
(3.3)

where $\hat{R}=(s_{zp},S_{zz}s_{zp},\dots,S_{zz}^{q-1}s_{zp})$ is the
$(m+k)\times q$ matrix of the Krylov sequence with a $(m+k)\times(m+k)$ matrix
$S_{zz}$ and a $(m+k)\times 1$ vector $s_{zp}$ defined as follows:

$\displaystyle S_{zz}=\frac{\bar{Z}^{T}(I-11^{T}/n)\bar{Z}}{n-1},$
$\displaystyle
s_{zp}=\frac{(\bar{Z}-\mathbb{E}(\bar{Z}))^{T}(P-\mathbb{E}(P))}{n-1},$

where $I$ is an identity matrix and $1$ is a matrix of ones. Intuitively, PLS
searches for factors that capture the highest variability in $\bar{Z}$, and at
the same time maximizes the covariance between $\bar{Z}$ and $P$. If the
number of factors equals the dimension of instruments, $q=m+k$, the method is
equivalent to OLS (Helland,, 1990).

## 4\. Deep Partial Least Squares for IV Regression

A useful generalization of PLS is to consider a deep-layered feed-forward
neural network structure:

$\displaystyle\hat{p}^{(1)}$ $\displaystyle=f(\bar{Z}\hat{\alpha}^{PLS}),$
$\displaystyle\hat{p}^{(2)}$
$\displaystyle=f(\hat{p}^{(1)}\hat{\alpha}^{(2)}),$ $\displaystyle\vdots$
(4.1) $\displaystyle\hat{p}^{(L)}$
$\displaystyle=f(\hat{p}^{(L-1)}\hat{\alpha}^{(L)}),$

where $f(\cdot)=\max(\cdot,0)$ is a rectified linear unit (ReLU) activation
function. Note that $\bar{Z}=[z,x]$ as before and $\hat{\alpha}^{PLS}$ is a
$q\times 1$ vector where $q$ is the number of PLS factors in the treatment
network. To reduce the dimensionality of the instruments, we predict the
treatment in the first layer by PLS. The parameters $\alpha^{(\ell)}$ in
subsequent layers $\ell=2,\dots,L$ can be identified by OLS or PLS. The
following proposition shows that the parameters are consistent in each layer
up to a proportionality constant.

######  Proposition 1.

Let $S_{zz}$ and $s_{zp}$ converge in probability to $\Sigma_{zz}$ (the
population variance of $z$) and $\sigma_{zp}$ (the population covariance of
$z$ and $p$) when $n\rightarrow\infty$. Moreover, let there exist a pair of
eigenvectors and eigenvalues $(v_{j},\lambda_{j})$ for which
$\sigma_{zp}=\sum_{j=1}^{M}\gamma_{j}v_{j}$ (with $\gamma_{j}$ non-zero for
each $j=1,\dots,M$). Assume also $\mathbb{E}(|g(U)|)<\infty$ and
$\mathbb{E}(U|g(U)|)<\infty$ and q = M. Then
$\hat{\alpha}=\\{\hat{\alpha}^{PLS},\hat{\alpha}^{(2)},\dots,\hat{\alpha}^{(L)}\\}$
are consistent up to a proportionality constant.

###### Proof.

We follow the approach by Naik and Tsai, (2000). Let
$\alpha^{\star}=\Sigma_{zz}^{-1}\sigma_{zp}$. Define
$R=(\sigma_{zp},\Sigma_{zz}\sigma_{zp},\dots,\Sigma_{zz}^{q-1}\sigma_{zp})$.
Then, based on the assumption that $S_{zz}\rightarrow\Sigma_{zz}$ and
$s_{zp}\rightarrow\sigma_{zp}$ when n $\rightarrow\infty$, we have:

$\displaystyle\hat{\alpha}^{PLS}\rightarrow
R(R^{T}\Sigma_{zz}R)^{-1}\Sigma_{zz}\alpha^{\star}\ \text{ in probability
when}\ n\rightarrow\infty.$

The assumptions $q=M$ and $\sigma_{zp}=\sum_{j=1}^{M}\gamma_{j}v_{j}$ imply
that $\alpha^{\star}$ is contained in the space spanned by $R$. Consequently,
$\Sigma_{zz}^{1/2}\alpha^{\star}$ is contained in the space spanned by
$R^{\star}=\Sigma_{zz}^{1/2}R$. Therefore,

$\displaystyle
R^{\star}(R^{\star^{T}}R^{\star})^{-1}R^{\star^{T}}\Sigma_{zz}^{1/2}\alpha^{\star}=\Sigma_{zz}^{1/2}\alpha^{\star},$

and

$\displaystyle
R(R^{\star^{T}}R^{\star})^{-1}R^{T}\Sigma_{zz}\alpha^{\star}=\alpha^{\star}.$

Hence, $\hat{\alpha}^{PLS}\rightarrow\alpha^{\star}$. Equation (3.1) implies
that $\sigma_{zp}=k\Sigma_{zz}\alpha$, and therefore,
$\hat{\alpha}^{PLS}\rightarrow\Sigma_{zz}^{-1}\sigma_{zp}=k\alpha$. $k$ is
defined in (3.1). This proves that the PLS in the first layer is consistent.

Now, consider the second layer $\ell=2$. If $\hat{\alpha}^{(2)}$ is estimated
by either OLS or PLS, (3.1) directly imply that $\hat{\alpha}^{(2)}=k\alpha$.
Next, to run OLS (or PLS) of $p$ on $\hat{p}^{(1)}$, we construct

$\displaystyle\hat{p}^{(1)}=\max(z\hat{\alpha}^{PLS},0)=\max(zk\alpha^{(1)},0)=k\max(z\alpha^{(1)},0)=kz^{(1)}.$

The model then becomes

$\displaystyle\hat{p}^{(2)}=max(kz^{(1)}k_{2}\alpha^{(2)},0)=k\cdot
k_{2}max(z^{(1)}\alpha^{(2)},0)=k^{(2)}z^{(2)}.$ (4.2)

By induction, we can consistently estimate $\alpha$ in each subsequent layer,
up to $L$. See Polson et al., (2021) for the detailed discussion of such a
deep learning structure.

Lastly, we note that the consistent estimator of the effect of the treatment
on the outcome is defined as

$\displaystyle\hat{\beta}^{GMM}=\big{(}(\bar{P}^{(L)})^{T}\bar{P}^{(L)}\big{)}^{-1}(\bar{P}^{(L)})^{T}\tilde{y},$
(4.3)

where $\bar{P}^{(L)}$ is the augmented treatment $[\hat{p},x]$ predicted by
layer $L$. Moreover, the predicted outcome is
$\hat{y}=\bar{P}\hat{\beta}^{GMM}$.

Define $\hat{e}^{L}$ as the predicted residuals from the treatment network in
layer $L$ and consider, $\tilde{p}=[p,\hat{e}^{L},x]$. Then, similarly for
2SRI:

$\displaystyle\hat{\beta}^{GMM}=(\tilde{p}^{T}\tilde{p})^{-1}\tilde{p}^{T}y.$
(4.4)

and $\hat{y}=\tilde{p}\hat{\beta}^{GMM}$.

∎

### D. Prediction and Bayesian Shrinkage

One possible extension of deep partial least squares is a quantification of
uncertainty in the density of the outcome. Our probabilistic model takes the
form:

$\displaystyle\tilde{y}\mid f,\bar{P}$ $\displaystyle\sim p(\tilde{y}\mid
f,\bar{P}),$ $\displaystyle f$
$\displaystyle=g(\bar{P}\beta_{p})+\varepsilon,$

where $\tilde{y}$ is an $n\times 1$ rescaled and recentered outcome variable
as before, $\bar{P}=[\hat{p},x]$ is an $n\times(1+k)$ augmented (predicted)
policy. $\beta_{p}=[\beta,\beta_{x}]$ is a $(1+k)\times 1$ vector of
coefficients in the outcome network. Here $g$ is a deep partial least squares
method. To estimate parameters in the first layer, the method uses the SIMPLS
algorithm (De Jong,, 1993). Subsequent layers use the stochastic gradient
descent (SGD) method (Bottou,, 2012) for optimizing and training the
parameters.

The key result, due to Brillinger, (2012) and Naik and Tsai, (2000), is that
$\beta_{p}$ can be estimated consistently, up to a constant of proportionality
using PLS, irrespective of the nonlinearity of $g$. Given a specification of
$g$, the constant of proportionality can also be estimated consistently with
$\sqrt{n}$-asymptotics. It is worth noting that typically, standard SGD
methods will not yield asymptotically normally distributed parameters.
However, they can substantially increase the precision of the coefficients of
interest.

Suppose that we wish to predict the outcome at a new level $\bar{P}^{\star}$.
Then, we can use the predictive distribution to make a forecast as well as
provide uncertainty bounds:

$y_{\star}\sim
p\left(y\;\mid\;g(\bar{P}^{\star}\hat{\beta}_{p}^{DPLS})\right).$

The advantage of modeling a probabilistic model is the flexibility and the
possibility to incorporate uncertainty in the parameters of interest. We
approximate the posterior distribution of $\hat{\beta}^{GMM}$ with its
asymptotic distribution based on the Bernstein-Von Mises theorem (see e.g.,
Van der Vaart,, 2000; Bhadra et al.,, 2019). Define
$\text{Data}=(\tilde{y},p,z,x)$, then the densities
$P(\hat{\beta}^{GMM}|\text{Data})$ and $P(\hat{\psi}_{1})$ come from a normal
distribution with the mean and variance depicted in (2.9). To shrink the
effect of redundant instruments in the treatment network, we consider a Ridge
regression estimator (Marquardt and Snee,, 1975):

$\displaystyle\hat{\alpha}^{Ridge}=(\bar{Z}^{T}\bar{Z}+\lambda
I_{m})^{-1}\bar{Z}^{T}p.$ (4.5)

However, as pointed out by a referee, this still underestimates the
uncertainty due to the estimation of
$(\hat{\beta}_{p}^{DPLS},\hat{\alpha}^{Ridge})$. A fully Bayesian model with a
non-uniform prior density of the coefficients can increase the precision of
uncertainty bounds.

## 5\. Applications

In this section, we evaluate the predictive performance of the DPLS-IV method
relative to benchmark methods. We generate synthetic data to mimic the high-
dimensional and sparse nature of instruments as described by Kong and Yu,
(2018). To illustrate how well the method captures the complex nonlinear and
sparse relation of covariates, our design mimics their setup. In addition to
synthetic designs, this section illustrates findings on data provided by
Angrist and Evans, (1996). Throughout the experiments, we split data into two
partitions. To find optimal parameters in each method, we use one sub-sample
(further partitioned into train and validation data) and illustrate prediction
performance measures on the other. We implement the experiments in R and the
code to replicate our findings is available upon request.

### E. Sparse Uncorrelated Instruments

The data-generating process is summarized by the following structural equation
model:

$\displaystyle y$ $\displaystyle=f(p\beta+x\beta_{x}+\xi)+\varepsilon,$ (5.1)
$\displaystyle p$
$\displaystyle=g(z\alpha+\text{sigmoid}(z^{2})\gamma+x\alpha_{x})+w,$
$\displaystyle w;u\sim\mathcal{N}(0,\Sigma),$
$\displaystyle\varepsilon\sim\mathcal{N}(0,\sigma_{\varepsilon}^{2}),$
$\displaystyle z$ $\displaystyle\sim\mathcal{N}(0,\Sigma_{z}),\ x\sim
N(0,\Sigma_{x}),$
$\displaystyle\alpha;\alpha_{x},\gamma,\beta,\beta_{x}\sim\mathcal{N}(0,1).$

(5.1) holds for each observation $i=1,\dots n$, where $n=1000$. $w$ and $\xi$
are correlated and jointly normally distributed with a mean vector $0$ and a
variance-covariance matrix $\Sigma=\begin{pmatrix}3.000&-0.087\\\
-0.087&0.010\end{pmatrix}$. In this simulation setup, $g$ and $f$ are modeled
by ReLU (or leakyReLu) (Dubey and Jain,, 2019). In addition to a $n\times 50$
matrix of instruments $z$, the treatment $p$ contains nonlinear
transformations of $z$ based on the sigmoid function. To introduce sparsity in
the design, out of fifty instruments, ten are redundant. In particular, only
forty instruments are relevant for predicting the treatment $p$. $x$ is an
$n\times 25$ matrix of covariates, and 20 of them have no influence on outcome
$y$. $\Sigma_{z}$ and $\Sigma_{x}$ represent covariance matrices of $z$ and
$x$, respectively. In this setting, the covariance between the instruments is
small (0.001) and the features $x$ are uncorrelated.

Figure 1 visualizes predictions of the policy $p$ in the treatment network
(Figure 10 in Appendix H shows predicted outcome $y$). In each case, $g$ and
$f$ represent a rectified linear unit function. A visual inspection of Figures
1 and 10 verifies that DPLS-IV results in more accurate predictions relative
to the other methods. In this setting, DeepIV (Hartford et al.,, 2016) is the
second-best alternative. According to Figures 1 and 10, predictions of the
treatment, $\widehat{p}$, appear to have a higher variability compared to the
outcome predictions. This is not surprising, as $p$ contains instruments in
addition to $x$. Note that, even though prediction performance measures of
DeepIV are close to those of DPLS-IV, Table V in Appendix H shows that,
compared to DeepIV, DPLS-IV is more robust to changes in the activation
function.

To evaluate the effectiveness of different methods, we examine their out-of-
sample $R^{2}$ and root mean squared error (RMSE) as we increase the parameter
values of $\Sigma$. The results, as depicted in Figure 2, demonstrate that
DPLS-IV exhibits robustness to increasing errors and is capable of accounting
for the endogeneity of $p$ reflected in the covariance of error terms $w$ and
$\xi$. Furthermore, our analysis, as shown in Table I, indicates that DPLS-IV
outperforms OLS, PLS, LASSO, and DeepIV methods by a significant margin. These
findings suggest that DPLS-IV holds considerable promise as a powerful and
reliable tool for predicting outcomes in the presence of endogeneity and
measurement error.

(a) OLS (b) PLS

(c) DeepIV (d) LASSO

(e) DPLS-IV

Figure 1: First stage prediction performance. The X-axis depicts predicted
treatment $\widehat{p}$, and the Y-axis represents true values of $p$. DPLS-IV
denotes the method introduced in this study. We use test data for evaluating
the methods.

To illustrate the predictive power of DPLS-IV, Figure 3 compares coefficients
estimated by the OLS and PLS methods in the treatment network. Parameters
estimated by PLS are closer to their true values relative to OLS.
Additionally, Figure 4 shows the absolute bias of these parameters. According
to Figure 4, the cumulative distribution function of the absolute bias of the
parameters recovered by PLS stochastically dominates the ones based on OLS and
LASSO. Specifically, smaller values of the absolute bias of the coefficients
are more likely under PLS, relative to OLS and LASSO. The sum of the absolute
bias is the smallest under PLS (14.406), followed by LASSO (15.233). OLS leads
to the highest value of the sum of the absolute bias (21.610).

Treatment network |  |  |  |  |  |
---|---|---|---|---|---|---
Measures | PLS | OLS | DeepIV | LASSO | DPLS-IV |
$R^{2}$ | 0.751 | 0.688 | 0.932 | 0.753 | 0.956 |
RMSE | 10.603 | 21.540 | 5.789 | 10.497 | 4.508 |
Outcome network |  |  |  |  |  |
$R^{2}$ | 0.932 | 0.933 | 0.889 | 0.933 | 0.938 |
RMSE | 1.668 | 1.718 | 4.573 | 1.670 | 1.622 |

Table I: Prediction performance of DPLS-IV relative to other methods. We
present out-of-sample $R^{2}$ and RMSE. To present the maximum prediction
performance of OLS, PLS and LASSO in the outcome network, we use residuals
predicted by DPLS-IV in the first stage. Figure 2: $\text{R}^{2}$ and RMSE for
increasing values of $\text{cov}(v,\varepsilon)$. The black circles represent
the estimated values which are combined with the dotted lines. Each color
corresponds to the corresponding $\text{R}^{2}$ and RMSE in the treatment and
outcome networks. Specifically, $\text{R}^{2}_{1}$ and $\text{RMSE}_{1}$
denote prediction performance measures in a treatment network. $R^{2}_{2}$ and
$\text{RMSE}_{2}$ are the prediction performance measures in the outcome
network.

(a) OLS (b) PLS

Figure 3: Estimated coefficients and their corresponding true values in the
treatment network. The X-axis shows coefficients predicted by OLS (left) and
PLS (right), and the Y-axis depicts the corresponding true values.  (a)
Figure 4: Empirical cumulative distribution functions (CDF) of the absolute
bias of the parameters estimated by PLS, OLS and LASSO.

To introduce uncertainty in the parameters of interest, we extend DPLS-IV to a
Bayesian setup and compare it to a classical Bayesian IV approach.

Outcome network |  |
---|---|---
Measures | Bayesian IV | Bayesian DPLS-IV
$R^{2}$ | 0.778 | 0.838
RMSE | 4.375 | 3.092
Table II: Prediction performance of Bayesian DPLS-IV and IV methods. The
measures are computed based on the mean prediction out of 10,000 predicted
outcome variables.

Table II shows that Bayesian DPLS-IV (with ReLU as an activation function and
two hidden layers) significantly outperforms its’ linear counterpart. Figure 5
verifies that the Bayesian DPLS-IV closely replicates the density of the
original outcome variable.

(a)
Figure 5: The density of the original and predicted outcome variables (on the
test data).

### F. Instruments with the Network Structure

The second simulation closely follows the network data structure described by
Kong and Yu, (2018). The data-generating process is the same as shown in
(5.1). However, the instrumental variable network comes from the preferential
attachment algorithm (Jeong et al.,, 2003). Each node of the network
represents one feature. The resulting network follows a power-law degree
distribution, and thus, is scale-free. That means, only a few instruments in
the network have a relatively large number of “neighbors”. The distance
between two instruments is the shortest path between them in the network. We
calculate a $p\times p$ ($p=50$) pairwise distance matrix $D$. Next, this
distance matrix is transformed into a covariance matrix
$\Sigma_{z,(i,j)}=0.7^{D_{(i,j)}}$, where $(i,j)$ represents the element in
each row $i$ and column $j$ of a matrix $D$ ($i,j=1,\dots,p$). Figure 6 shows
an example of such a network with 100 nodes.

(a)
Figure 6: A network with 100 nodes.  (a)
Figure 7: Left - CDF of the outcome variable; Right - CDF of the absolute
value of the bias of the parameters (estimated by the corresponding method) in
the treatment network.

For comparability, DPLS-IV and DeepIV consist of the same number of layers and
neurons in each layer. In particular, the input layer consists of 200 neurons;
the hidden layers consist of 200, 100, and 50 neurons, respectively. We use
ReLU as an activation function. Figure 7 shows cumulative distribution
functions of the predicted outcome variable and the absolute bias of the
estimated coefficients (right) in the treatment network. Based on the results
in Figure 7, the outcome predicted by DPLS-IV is closer to the CDF of the
simulated outcome variable. Moreover, DPLS-IV yields a CDF of the absolute
bias of the parameters that stochastically dominates the CDF of the other
benchmark methods.

We also investigate prediction performance in the outcome network. Table III
shows $R^{2}$ and RMSE of DPLS-IV relative to other methods. In this setting,
DPLS-IV outperforms other benchmark methods.

Outcome network |  |  |  |  |  |
---|---|---|---|---|---|---
Measures | PLS | OLS | DeepIV | LASSO | DPLS-IV |
$R^{2}$ | 0.940 | 0.930 | 0.912 | 0.942 | 0.944 |
RMSE | 1.631 | 1.846 | 2.108 | 1.609 | 1.585 |
Table III: Prediction performance of DPLS-IV relative to other methods. We
present out-of-sample $R^{2}$ and RMSE. To present the maximum prediction
performance of OLS, PLS, and LASSO in the outcome network, we use residuals
predicted by DPLS-IV in the first stage. (a)
Figure 8: Predicted outcome $y$ and the corresponding 95% confidence
intervals.

Additionally, we compare a Bayesian DPLS-IV to a Bayesian IV approach.
Bayesian DPLS-IV consists of two hidden layers and a ReLU activation function.
Note that, the RMSE based on the mean predicted outcome out of 10,000
predictions increases (2.262) relative to DPLS-IV. However, $R^{2}$ is
unchanged (0.941). Figure 8 shows predicted outcome values against their true
counterparts. Bayesian DPLS-IV leads to outcome predictions that are closer to
the true values.

### G. Labor Supply of Women

Angrist and Evans, (1996) examine the effect of childbearing on women’s labor
supply. They use a mixed sibling-sex composition and twins as instruments for
the size of the family. To illustrate the method, our analysis uses 1980 U.S.
Census data that include all women with two or more children. The model that
we are going to estimate is defined by the following structural equation
model:

$\displaystyle y$
$\displaystyle=f(\text{kids}\cdot\beta+x\beta_{x}+\xi)+\varepsilon.$ (5.2)
kids $\displaystyle=g(\text{twins}\cdot\alpha+\text{twins}\cdot
x\gamma+x\alpha_{x})+w,$

where $y$ is the outcome variable and measures the logarithm of hours worked
per week by the mother. The outcome is observed when the age of the mother is
more than the average age of the mothers in the population. In particular,

$\displaystyle y=\text{log(hourswm)}\times\mathbb{I}(agem>\mathbb{E}(agem)),$
(5.3)

where $\mathbb{I}(agem>\mathbb{E}(agem))$ is an indicator variable and equals
one if the age of the mother is more than the population mean, and zero
otherwise. The treatment, kids, is the number of total kids in a family. twins
represents an instrumental variable and equals one if the second and third
children are twins, otherwise zero. Additionally, we use interactions of the
instrument with covariates, $\text{twins}\cdot x$ as instruments for the
number of kids. The covariates include the gender and age of the first and
second child, the mother’s age, marital status, race, education, and the age
of the mother when she first gave birth. See Angrist and Evans, (1996) for a
detailed summary of the variables.

The prediction performance advantage of DPLS-IV is also clearly evident in
Table IV. The first stage prediction performance of DPLS-IV is similar in each
method, however, $R^{2}$ (RMSE) is considerably high (low) in the outcome
network. It is worth noting that the measures become substantially worse in
the outcome network compared to the ones in the treatment network. One of the
potential reasons is that the outcome distribution is bimodal. Figure 11
presents the original density of the outcome variable and the density of the
outcome predicted by Bayesian DPLS-IV. Figure 11 shows that Bayesian DPLS-IV
successfully captures the bimodal nature of the outcome.

Treatment network |  |  |  |  |
---|---|---|---|---|---
Measures | PLS | OLS | DeepIV | LASSO | DPLS-IV
$R^{2}$ | 0.204 | 0.205 | 0.226 | 0.205 | 0.233
RMSE | 0.655 | 0.655 | 0.647 | 0.647 | 0.644
Outcome network |  |  |  |  |
$R^{2}$ | 0.532 | 0.536 | 0.771 | 0.536 | 0.772
RMSE | 12.838 | 12.782 | 8.980 | 12.780 | 8.960
Table IV: Prediction performance of DPLS-IV relative to other methods. We
present out-of-sample $R^{2}$ and RMSE. To present the maximum prediction
performance of OLS, PLS and LASSO in the outcome network, we use residuals
predicted by DPLS-IV in the first stage.

We also investigate the layer-by-layer transformation of the proposed method.
Figure 9 illustrates the feature representation of the original outcome $y$
(on the test data), the score matrix ($T$), and the final feed-forward neural
network layer in DPLS-IV. We consider the projection of features on two
distinct classes of the outcome variable. In particular, when $y>0$, we label
it as a class 1, and when $y=0$, it represents a class 0. Figure 9 shows that
the intermediate layers in DPLS-IV significantly improve the representation of
the original covariate space. The borders of the two classes become highly
evident in the last layer of DPLS-IV.

(a) Original (b) PLS Scores T

(c) Final Layer

Figure 9: Feature representation of DPLS-IV based on the t-SNE algorithm
(Wattenberg et al.,, 2016), with the perplexity parameter equal to 50.
Features are reduced to two distinct classes of the outcome $y$, with a class
1 when $y>0$, and 0 otherwise. The final layer represents the predicted $y$ of
individual neurons in the last layer of DPLS-IV. For computational benefits,
we randomly sample 1000 observations from the original data.

## 6\. Discussion

In this article, we propose deep partial least squares for reducing the
dimension of the instrumental variable space. The deep partial least squares
method efficiently extracts features based on partial least squares and
further processes the input with a feed-forward deep learner. The method is
well-tailored for correlated instruments with sparse and nonlinear structures.
More importantly, deep partial least squares are consistent, up to a
proportionality constant.

The applications on synthetic data as well as the application to the effect of
childbearing on the mother’s labor supply (Angrist and Evans,, 1996) show that
the deep partial least squares method outperforms other related methods.
Moreover, a flexible number of layers allows us to efficiently capture
nonlinearities embedded in the instrumental variable network.

An interesting extension of this work is to consider a Bayesian model with
various priors on the coefficients of interest. Datta and Ghosh, (2013)
consider asymptotic properties of Bayes risk with the horseshoe prior. We
believe, investigating different prior beliefs can result in the increased
predictive performance of deep partial least squares. Another useful extension
is to draw applications to eminent domain with judge characteristics as
instruments (Belloni et al.,, 2012). In addition to prediction problems, the
consistency of DPLS allows us to address the estimation of the treatment
effect. In the future, we plan to demonstrate the precision of the treatment
effect and compare it to other benchmark algorithms.

## References

  * Abadie et al., (2002) Abadie, A., Angrist, J., and Imbens, G. (2002).  Instrumental variables estimates of the effect of subsidized training on the quantiles of trainee earnings.  Econometrica, 70(1):91–117.
  * Abdi and Williams, (2010) Abdi, H. and Williams, L. J. (2010).  Principal component analysis.  Wiley Interdisciplinary Reviews: Computational Statistics, 2(4):433–459.
  * Adcock, (2007) Adcock, C. (2007).  Extensions of stein’s lemma for the skew-normal distribution.  Communications in Statistics—Theory and Methods, 36(9):1661–1671.
  * Adkins, (2012) Adkins, L. C. (2012).  Testing parameter significance in instrumental variables probit estimators: some simulation.  Journal of Statistical Computation and Simulation, 82(10):1415–1436.
  * Amemiya, (1985) Amemiya, Y. (1985).  Instrumental variable estimator for the nonlinear errors-in-variables model.  Journal of Econometrics, 28(3):273–289.
  * Angrist and Evans, (1996) Angrist, J. and Evans, W. N. (1996).  Children and their parents’ labor supply: Evidence from exogenous variation in family size.  National Bureau of Economic Research.
  * Baum et al., (2003) Baum, C. F., Schaffer, M. E., and Stillman, S. (2003).  Instrumental variables and gmm: Estimation and testing.  The Stata Journal, 3(1):1–31.
  * Beise et al., (2021) Beise, H.-P., Da Cruz, S. D., and Schröder, U. (2021).  On decision regions of narrow deep neural networks.  Neural Networks, 140:121–129.
  * Belloni et al., (2012) Belloni, A., Chen, D., Chernozhukov, V., and Hansen, C. (2012).  Sparse models and methods for optimal instruments with an application to eminent domain.  Econometrica, 80(6):2369–2429.
  * Belloni et al., (2011) Belloni, A., Chernozhukov, V., and Hansen, C. (2011).  Inference for high-dimensional sparse econometric models.  arXiv preprint arXiv:1201.0220.
  * Bhadra et al., (2020) Bhadra, A., Datta, J., Li, Y., and Polson, N. (2020).  Horseshoe regularisation for machine learning in complex and deep models 1.  International Statistical Review, 88(2):302–320.
  * Bhadra et al., (2019) Bhadra, A., Datta, J., Polson, N. G., and Willard, B. (2019).  Lasso meets horseshoe: A survey.  Statistical Science, 34(3):405–427.
  * Bottou, (2012) Bottou, L. (2012).  Stochastic gradient descent tricks.  In Neural Networks: Tricks of the Trade, pages 421–436. Springer.
  * Brillinger, (2012) Brillinger, D. R. (2012).  A generalized linear model with “gaussian” regressor variables.  In Selected Works of David Brillinger, pages 589–606. Springer.
  * Cengiz et al., (2022) Cengiz, D., Dube, A., Lindner, A., and Zentler-Munro, D. (2022).  Seeing beyond the trees: Using machine learning to estimate the impact of minimum wages on labor market outcomes.  Journal of Labor Economics, 40(S1):S203–S247.
  * Chernozhukov and Hansen, (2004) Chernozhukov, V. and Hansen, C. (2004).  The effects of 401 (k) participation on the wealth distribution: an instrumental quantile regression analysis.  Review of Economics and statistics, 86(3):735–751.
  * (17) Chernozhukov, V., Hansen, C., and Spindler, M. (2015a).  Post-selection and post-regularization inference in linear models with many controls and instruments.  American Economic Review, 105(5):486–490.
  * (18) Chernozhukov, V., Hansen, C., and Spindler, M. (2015b).  Valid post-selection and post-regularization inference: An elementary, general approach.  Annu. Rev. Econ., 7(1):649–688.
  * Datta and Ghosh, (2013) Datta, J. and Ghosh, J. K. (2013).  Asymptotic properties of bayes risk for the horseshoe prior.  Bayesian Analysis, 8(1):111–132.
  * De Jong, (1993) De Jong, S. (1993).  Simpls: an alternative approach to partial least squares regression.  Chemometrics and Intelligent Laboratory Systems, 18(3):251–263.
  * Debiolles et al., (2004) Debiolles, A., Oukhellou, L., and Aknin, P. (2004).  Combined use of partial least squares regression and neural network for diagnosis tasks.  In Proceedings of the 17th International Conference on Pattern Recognition, 2004. ICPR 2004., volume 4, pages 573–576. IEEE.
  * Dubey and Jain, (2019) Dubey, A. K. and Jain, V. (2019).  Comparative study of convolution neural network’s relu and leaky-relu activation functions.  In Applications of Computing, Automation and Wireless Systems in Electrical Engineering, pages 873–880. Springer.
  * Farrell et al., (2021) Farrell, M. H., Liang, T., and Misra, S. (2021).  Deep neural networks for estimation and inference.  Econometrica, 89(1):181–213.
  * Fisher, (1922) Fisher, R. A. (1922).  On the mathematical foundations of theoretical statistics.  Philosophical Transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character, 222(594-604):309–368.
  * Florens et al., (2002) Florens, J.-P., Heckman, J., Meghir, C., and Vytlacil, E. (2002).  Instrumental variables, local instrumental variables and control functions.  Technical report, Cemmap Working Paper.
  * Gagnon-Bartsch et al., (2013) Gagnon-Bartsch, J. A., Jacob, L., and Speed, T. P. (2013).  Removing unwanted variation from high dimensional data with negative controls.  Berkeley: Tech Reports from Dep Stat Univ California, pages 1–112.
  * Geladi and Kowalski, (1986) Geladi, P. and Kowalski, B. R. (1986).  Partial least-squares regression: a tutorial.  Analytica chimica acta, 185:1–17.
  * Hahn et al., (2020) Hahn, P. R., Murray, J. S., and Carvalho, C. M. (2020).  Bayesian regression tree models for causal inference: Regularization, confounding, and heterogeneous effects (with discussion).  Bayesian Analysis, 15(3):965–1056.
  * Hanin and Rolnick, (2019) Hanin, B. and Rolnick, D. (2019).  Deep relu networks have surprisingly few activation patterns.  Advances in Neural Information Processing Systems, 32.
  * Hansen and Singleton, (1982) Hansen, L. P. and Singleton, K. J. (1982).  Generalized instrumental variables estimation of nonlinear rational expectations models.  Econometrica: Journal of the Econometric Society, pages 1269–1286.
  * Hartford et al., (2016) Hartford, J., Lewis, G., Leyton-Brown, K., and Taddy, M. (2016).  Counterfactual Prediction with Deep Instrumental Variables Networks.  arXiv:1612.09596 [cs, stat].
  * Hartford et al., (2017) Hartford, J., Lewis, G., Leyton-Brown, K., and Taddy, M. (2017).  Deep iv: A flexible approach for counterfactual prediction.  In International Conference on Machine Learning, pages 1414–1423. PMLR.
  * Heckman and Navarro-Lozano, (2004) Heckman, J. and Navarro-Lozano, S. (2004).  Using matching, instrumental variables, and control functions to estimate economic choice models.  Review of Economics and Statistics, 86(1):30–57.
  * Heckman, (1995) Heckman, J. J. (1995).  Randomization as an instrumental variable.  National Bureau of Economic Research.
  * Helland, (1990) Helland, I. S. (1990).  Partial least squares regression and statistical models.  Scandinavian Journal of Statistics, pages 97–114.
  * Hoogerheide et al., (2007) Hoogerheide, L., Kleibergen, F., and van Dijk, H. K. (2007).  Natural conjugate priors for the instrumental variables regression model applied to the angrist–krueger data.  Journal of Econometrics, 138(1):63–103.
  * Iwata, (2001) Iwata, S. (2001).  Recentered and Rescaled Instrumental Variable Estimation of Tobit and Probit Models with Errors in Variables.  Econometric Reviews, 20(3):319–335.
  * Iwata and Ghahramani, (2017) Iwata, T. and Ghahramani, Z. (2017).  Improving output uncertainty estimation and generalization in deep learning via neural network gaussian processes.  arXiv preprint arXiv:1707.05922.
  * Jeong et al., (2003) Jeong, H., Néda, Z., and Barabási, A.-L. (2003).  Measuring preferential attachment in evolving networks.  EPL (Europhysics Letters), 61(4):567.
  * Jia et al., (2016) Jia, W., Zhao, D., and Ding, L. (2016).  An optimized rbf neural network algorithm based on partial least squares and genetic algorithm for classification of small sample.  Applied Soft Computing, 48:373–384.
  * Kong and Yu, (2018) Kong, Y. and Yu, T. (2018).  A deep neural network model using random forest to extract feature representation for gene expression data classification.  Scientific Reports, 8(1):1–9.
  * Landsman and Nešlehová, (2008) Landsman, Z. and Nešlehová, J. (2008).  Stein’s lemma for elliptical random vectors.  Journal of Multivariate Analysis, 99(5):912–927.
  * Liu et al., (2020) Liu, R., Shang, Z., and Cheng, G. (2020).  On deep instrumental variables estimate.  arXiv preprint arXiv:2004.14954.
  * Lopes and Polson, (2010) Lopes, H. F. and Polson, N. G. (2010).  Extracting sp500 and nasdaq volatility: The credit crisis of 2007-2008.  Handbook of Applied Bayesian Analysis, pages 319–342.
  * Marquardt and Snee, (1975) Marquardt, D. W. and Snee, R. D. (1975).  Ridge regression in practice.  The American Statistician, 29(1):3–20.
  * McCullagh and Polson, (2018) McCullagh, P. and Polson, N. G. (2018).  Statistical sparsity.  Biometrika, 105(4):797–814.
  * Mogstad et al., (2021) Mogstad, M., Torgovitsky, A., and Walters, C. R. (2021).  The causal interpretation of two-stage least squares with multiple instrumental variables.  American Economic Review, 111(11):3663–98.
  * Monis et al., (2022) Monis, J. B., Sarkar, R., Nagavarun, S., and Bhadra, J. (2022).  Efficient net: Identification of crop insects using convolutional neural networks.  In 2022 International Conference on Advances in Computing, Communication and Applied Informatics (ACCAI), pages 1–7. IEEE.
  * Moreno and Terwiesch, (2014) Moreno, A. and Terwiesch, C. (2014).  Doing business with strangers: Reputation in online service marketplaces.  Information Systems Research, 25(4):865–886.
  * Naik and Tsai, (2000) Naik, P. and Tsai, C.-L. (2000).  Partial least squares estimator for single-index models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(4):763–771.
  * Navarro, (2010) Navarro, S. (2010).  Control functions.  In Microeconometrics, pages 20–28. Springer.
  * Polson et al., (2021) Polson, N., Sokolov, V., and Xu, J. (2021).  Deep learning partial least squares.  arXiv preprint arXiv:2106.14085.
  * Polson and Ročková, (2018) Polson, N. G. and Ročková, V. (2018).  Posterior concentration for sparse deep learning.  Advances in Neural Information Processing Systems, 31.
  * Polson and Sokolov, (2017) Polson, N. G. and Sokolov, V. O. (2017).  Deep learning for short-term traffic flow prediction.  Transportation Research Part C: Emerging Technologies, 79:1–17.
  * Polson and Sun, (2019) Polson, N. G. and Sun, L. (2019).  Bayesian l 0-regularized least squares.  Applied Stochastic Models in Business and Industry, 35(3):717–731.
  * Puelz et al., (2017) Puelz, D., Hahn, P. R., and Carvalho, C. M. (2017).  Variable selection in seemingly unrelated regressions with random predictors.  Bayesian Analysis, 12(4):969–989.
  * Sarstedt et al., (2022) Sarstedt, M., Hair, J. F., Pick, M., Liengaard, B. D., Radomir, L., and Ringle, C. M. (2022).  Progress in partial least squares structural equation modeling use in marketing research in the last decade.  Psychology & Marketing, 39(5):1035–1064.
  * Stone and Brooks, (1990) Stone, M. and Brooks, R. J. (1990).  Continuum regression: cross-validated sequentially constructed prediction embracing ordinary least squares, partial least squares and principal components regression.  Journal of the Royal Statistical Society: Series B (Methodological), 52(2):237–258.
  * Terza et al., (2008) Terza, J. V., Basu, A., and Rathouz, P. J. (2008).  Two-stage residual inclusion estimation: addressing endogeneity in health econometric modeling.  Journal of Health Economics, 27(3):531–543.
  * Tibshirani, (2011) Tibshirani, R. (2011).  Regression shrinkage and selection via the lasso: a retrospective.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(3):273–282.
  * Van der Vaart, (2000) Van der Vaart, A. W. (2000).  Asymptotic statistics, volume 3.  Cambridge university press.
  * Wattenberg et al., (2016) Wattenberg, M., Viégas, F., and Johnson, I. (2016).  How to use t-sne effectively.  Distill, 1(10):e2.
  * White, (1982) White, H. (1982).  Instrumental variables regression with independent observations.  Econometrica: Journal of the Econometric Society, pages 483–499.
  * Yarotsky, (2017) Yarotsky, D. (2017).  Error bounds for approximations with deep relu networks.  Neural Networks, 94:103–114.
  * Zellner, (1975) Zellner, A. (1975).  Bayesian analysis of regression error terms.  Journal of the American Statistical Association, 70(349):138–144.

## Appendix

### H. Prediction Performance Measures with LeakyReLu

Treatment network |  |  |  |  |
---|---|---|---|---|---
Measures | PLS | OLS | LASSO | DeepIV | DPLS-IV
$R^{2}$ | 0.757 | 0.696 | 0.760 | 0.939 | 0.957
RMSE | 10.503 | 21.404 | 10.398 | 5.498 | 4.449
Outcome network |  |  |  |  |
$R^{2}$ | 0.933 | 0.933 | 0.934 | 0.878 | 0.938
RMSE | 1.666 | 1.720 | 1.667 | 2.224 | 1.623
Table V: Prediction performance of DPLS-IV relative to other methods. We
present out-of-sample $R^{2}$ and RMSE. To present the maximum prediction
performance of OLS, PLS and Tobit in the outcome network, we use residuals
predicted by DPLS-IV in the first stage. We use LeakyReLU to simulate the
outcome variables. DPLS-IV consists of an input layer with 50 neurons, and a
hidden layer with 30 neurons. DeepIV consists of 100 neurons in the input
layer, followed by three hidden layers with 100, 100 and 30 neurons,
respectively. We use ReLU to activate neurons in DeepIV and DPLS-IV.

(a) OLS (b) PLS

(c) DeepIV (d) LASSO

(e) DPLS-IV

Figure 10: The second stage prediction performance. The X-axis depicts
predicted outcome $\widehat{y}$, and the Y-axis represents true values of $y$.
DPLS-IV denotes the method introduced in this study. We use test data for
evaluating the methods.  (a)
Figure 11: The original distribution of the outcome variable and the one
predicted by Bayesian DPLS-IV. Bayesian DPLS-IV consists of two hidden layers
and a ReLU activation function.