Deep Proxy Causal Learning and its Application to Confounded Bandit Policy Evaluation
Year: 2021

Abstract:
Proxy causal learning (PCL) is a method for estimating the causal effect of treatments on outcomes in the presence of unobserved confounding, using proxies (structured side information) for the confounder. This is achieved via two-stage regression: in the first stage, we model relations among the treatment and proxies; in the second stage, we use this model to learn the effect of treatment on the outcome, given the context provided by the proxies. PCL guarantees recovery of the true causal effect, subject to identifiability conditions. We propose a novel method for PCL, the deep feature proxy variable method (DFPV), to address the case where the proxies, treatments, and outcomes are high-dimensional and have nonlinear complex relationships, as represented by deep neural network features. We show that DFPV outperforms recent state-of-the-art PCL methods on challenging synthetic benchmarks, including settings involving high dimensional image data. Furthermore, we show that PCL can be applied to off-policy evaluation for the confounded bandit problem, in which DFPV also exhibits competitive performance.

heading: 1 Introduction
text: In causal learning, we aim to estimate the effect of our actions on the world. For example, we may be interested in measuring the impact of flight ticket prices on sales [2, 34], or the effect of grade retention on cognitive development [4]. We refer to our action as a treatment, which results in a particular outcome. It is often impossible to determine the effect of treatment on outcome from observational data alone, since the observed joint distribution of treatment and outcome can depend on a common confounder which influences both, and which might not be observed. In our example on sales of plane tickets given a price, the two might even be positively correlated in some circumstances, such as the simultaneous increase in sales and prices during the holiday season. This does not mean that raising the price causes sales to increase. In this context, people’s desire to fly is a confounder, since it affects both the number of ticket sales and the prices people are willing to accept. Thus, we need to correct the bias caused by the confounder.
One common assumption to cope with confounding bias is to assume no unobserved confounders exist [8], or more generally, the ignorable treatment assignment assumption [28], which states that the treatment assignment is independent of the potential outcomes caused by the treatment, given the background data available. Although a number of methods are proposed based on this assumption [7, 9, 36], it can be too restrictive, since it is often difficult to determine how the confounder affects treatment assignments and outcomes.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
A less restrictive assumption is that we have access to proxy variables, which contain relevant side information on the confounder. In the flight tickets example, we can use the number of views of the ticket reservation page as a proxy variable, which reflects peoples’ desire for flights. Note that if we can completely recover the confounder from proxy variables, the ignorable treatment assignment assumption can be satisfied. Motivated by this, Lee et al. [14] and Louizos et al. [17] aim to recover the distribution of confounders from proxy variables using modern machine learning techniques such as generative adversarial networks [5] or variational auto-encoders (VAE) [11]. Although these methods exhibit powerful empirical performance, there is little theory that guarantees the correct recovery of the causal effects.
Kuroki and Pearl [13] first considered the necessary conditions on proxy variables to provably recover the underlying causal effect via direct recovery of the hidden confounder. This work was in turn generalized by Miao et al. [20]. In their work, it is shown that two types of proxy variables are sufficient to recover the true causal effects without explicitly recovering the confounder. One is an outcome-inducing proxy, which correlates with confounders and causes the outcome, and the other is a treatment-inducing proxy which is caused by confounders and correlates with the treatment. In the flight ticket example, we can use the number of views of the ticket reservation page as the outcome-inducing proxy, and the cost of fuel as the treatment-inducing proxy. Given these proxy variables, Miao et al. [20] show that the true causal effect can be recovered by solving a Fredholm integral equation, which is referred to as the proxy causal learning (PCL) problem. The PCL problem is known to have an interesting relation to the causal inference with multiple treatments, which uses a subset of treatments as proxy variables in PCL [22, 33].
Although the PCL problem has a solid theoretical grounding, the question of how to estimate the causal effect remains a practical challenge, in particular when we consider nonlinear causal relationships or high-dimensional treatments. In Miao et al. [20], the treatment and outcome are assumed to be categorical variables. In a follow-up study, Miao et al. [21] show that we can learn a linear causal effect by a method of moments. Deaner [3] models the causal effect as a linear combination of nonlinear basis functions, which is learned by solving two stage regression. These two methods are extended by Mastouri and Zhu et al. [18], who learn the causal effect in a predefined reproducing kernel Hilbert space (RKHS). We provide an overview of the PCL problem and the two-stage regression in Section 2. Although these methods enjoy desirable theoretical properties, the flexibility of the model is limited, since all existing work uses pre-specified features.
In this paper, we propose a novel Deep Feature Proxy Variable (DFPV) method, which is the first work to apply neural networks to the PCL problem. The technique we employ builds on earlier work in instrumental variable (IV) regression, which is a related causal inference setting to PCL. A range of deep learning methods has recently been introduced for IV regression [1, 6, 35]. We propose to adopt the Deep Feature Instrumental Variable method [35], which learns deep adaptive features within a two-stage regression framework. Details of DFPV are given in Section 3. In Section 4, we empirically show that DFPV outperforms other PCL methods in several examples. We further apply PCL methods to the off-policy evaluation problem in a confounded bandit setting, which aims to estimate the average reward of a new policy given data with confounding bias. We discuss the setting in Section 3, and show the superiority of DFPV in experiments in Section 4.

heading: 2 Preliminaries
text: In this section, we introduce the proxy causal learning (PCL) problem and describe the existing two-stage regression methods to solve PCL.
Notation. Throughout the paper, a capital letter (e.g. A) denotes a random variable, and we denote the set where a random variable takes values by the corresponding calligraphic letter (e.g. A). The symbol P (·) denotes the probability distribution of a random variable (e.g. P (A)). We use a lowercase letter to denote a realization of a random variable (e.g. a). We denote the expectation over a random variable as E[·] and ‖f‖P(·) as the L2-norm of a function f with respect to P (·); i.e. ‖f‖P(A) = √ EA [f2(A)].
2.1 Problem Setting for Proxy Causal Learning
We begin with a description of the PCL setting. We observe a treatment A ∈ A, whereA ⊂ RdA , and the corresponding outcome Y ∈ Y , where Y ⊂ R. We assume that there exists an unobserved confounder U ∈ U that affects both A and Y . The goal of PCL is to estimate the structural function fstruct(X) defined as
fstruct(a) = EU [EY [Y |A = a,U ]] , which we assume to be continuous. This function is also known as the Average Treatment Effect (ATE). The challenge of estimating fstruct
is that the confounder U is not observable — we cannot estimate the structural function from observations A and Y alone. To deal with this, we assume access to a treatment-inducing proxy variable Z, and an outcome-inducing proxy variable W , which satisfy the following structural assumption and completeness assumption.
Assumption 1 (Structural Assumption [3, 18]). We assume Y ⊥⊥Z|A,U , and W ⊥⊥(A,Z)|U .
Assumption 2 (Completeness Assumption on Confounder [3, 18]). Let l : U → R be any square integrable function ‖l‖P(U) <∞. The following conditions hold:
E [l(U) | A = a,W = w] = 0 ∀(a,w) ∈ A×W ⇔ l(u) = 0 P (U) -a.e. E [l(U) | A = a,Z = z] = 0 ∀(a, z) ∈ A×Z ⇔ l(u) = 0 P (U) -a.e.
Figure 1 shows the causal graph describing these relationships. In our setting, we assume that there is no observable confounder, though this may be easily included [18, 31] as presented in Appendix D. Here, the bidirectional arrows mean that we allow both directions or even a common ancestor variable. Given these assumptions, it is shown that the structural function can be expressed using a solution of an integral equation.
Proposition 1 (Miao et al. [20]). Let Assumptions 1, 2 and Assumptions 4, 5, 6 in Appendix B hold. Then there exists at least one solution to the functional equation
E [Y |A = a,Z = z] = ∫ h(a,W )ρW (w|A = a,Z = z)dw, (1)
which holds for any (a, z) ∈ A×Z . Here, we denote ρW (w|A = a,Z = z) as the density function of the conditional probability P (W |A = a,Z = z). Let h∗ be a solution of (1). The structural function fstruct is given as
fstruct(a) = EW [h∗(a,W )] . (2)
For completeness, we present a proof in Appendix B (Lemma 2 and Corollary 1), which is due to Miao et al. [21] and Deaner [3]. Following Miao et al. [21], we call h∗ the bridge function. From Proposition 1, we can see that the estimation of the structural function reduces to the estimation of the bridge function, since once we obtain the bridge function, the structural function directly follows from (2).

heading: 2.2 Two-stage Regression Scheme for Proximal Causal Learning with Fixed Features
text: To obtain the bridge function, we need to solve the functional equation (1). However, directly solving (1) in a rich function space can be ill-posed (see discussion in Nashed and Wahba [25]). To address this, recent works [3, 18] minimize the following regularized loss LPV to obtain an estimate of the bridge function ĥ:
ĥ = arg min h∈Hh
LPV(h), LPV(h) = EY,A,Z [ (Y −EW |A,Z [h(A,W )])2 ] + Ω(h), (3)
where Hh is an arbitrary space of continuous functions and Ω(h) is a regularizer on h. Note that this loss can be interpreted as the deviation of the r.h.s and the l.h.s in (1) measured in L2-norm with respect to the distribution P (Y,A,Z). Deaner [3] and Mastouri and Zhu et al. [18]solve the minimization problem (3) by successively solving two-stage regression problems. They model the bridge function as
h(a,w) = u>(ψA(a)⊗ψW (w)) (4)
where u is a learnable weight vector, ψA(a),ψW (w) are vectors of fixed basis functions, and ⊗ is a Kronecker product, defined as a⊗ b = vec(ab>) for any finite dimensional vectors a,b.1 An estimate û is obtained by solving the successive regression problems. In Stage 1, we estimate the conditional expectation EW |A=a,Z=z [ψW (W )] as a function of a, z. In Stage 2, we substitute the model (4) into the inner conditional expectation in LPV, EW |A=a,Z=z [h(a,W )] = u>(ψA(a)⊗EW |A=a,Z=z [ψW (W )]), and then minimize LPV with respect to u using the estimate of EW |A=a,Z=z [ψW (W ))] from Stage 1. The above idea can be implemented as follows. We model the conditional expectation as
EW |A=a,Z=z [ψ(W )] = V (φA(a)⊗φZ(z)), where φA(a),φZ(z) are another set of basis functions, and V is a matrix to be learned. Note that we can use different basis functions for φA(a) and ψA(a). In Stage 1, the matrix V is learned by minimizing the following loss,
L1(V ) = EW,A,Z [ ‖ψW (W )−V (φA(A)⊗φZ(Z))‖2 ] + λ1‖V ‖2, (5)
where λ1 > 0 is a regularization parameter. This is a linear ridge regression problem with multiple targets, which can be solved analytically. In Stage 2, given the minimizer V̂ = arg minV L1(V ), we can obtain û by minimizing the loss
L2(u) = EY,A,Z [ ‖Y −u>(ψA(A)⊗ (V̂ (φA(A)⊗φZ(Z))))‖2 ] + λ2‖u‖2, (6)
where λ2 > 0 is another regularization parameter. Stage 2 corresponds to another linear ridge regression from input ψA(A)⊗ (V̂ (φA(A)⊗φZ(Z))) to target Y , and also enjoys a closed-form solution. Given the learned weights û = arg minuL2(u), the estimated structural function is f̂struct(a) = û
>(ψA(a)⊗EW [ψW (W )]). When fixed feature dictionaries are used, this two-stage regression benefits from strong theoretical guarantees [3, 18]. The use of pre-specified feature maps, limits the scope and flexibility of the method, however, especially if the treatment and proxies are high dimensional (e.g. images or text), and the causal relations are nonlinear. To overcome this, we propose to use adaptive features, expressed by neural nets, as described in the next section.

heading: 3 Deep Feature Proxy Causal Learning
text: In this section, we develop the DFPV algorithm, which learns adaptive features modeled by neural nets using a technique similar to Xu et al. [35]. As in Mastouri and Zhu et al. [18], we assume that we do not necessarily have access to observations from the joint distribution of (A,Y,Z,W ). Instead, we are given m observations of (A,Z,W ) for Stage 1 and n observations of (A,Z,Y ) for Stage 2. We denote the stage 1 observations by (ai, zi,wi) and the stage 2 observations by (ãi, z̃i, ỹi). If observations of (A,Y,Z,W ) are given for both stages, we can evaluate the out-of-sample loss of Stage 1 using Stage 2 data and vice versa, and these losses can be used for hyper-parameter tuning of λ1, λ2 (Appendix A). We first introduce two-stage regression with adaptive feature maps and then describe the detailed learning procedure of DFPV.

heading: 3.1 Two-stage regression with adaptive features
text: In DFPV, we consider the following models of the bridge function h(a,w) and conditional feature mean EW |A,Z [ψθW (W )]:
h(a,w) = u>(ψθA(2)(a)⊗ψθW (w)), EW |A=a,Z=z [ψθW (W )] = V (φθA(1)(a)⊗φθZ (z)), where u and V are parameters, and φθA(1) ,φθZ ,ψθA(2) ,ψθW are neural nets parametrized by θA(1), θZ , θA(2), θW , respectively. Again, we may use different neural nets in the treatment features φθA(1) and ψθA(2) .
1Throughout this paper, we assume the number of basis functions to be finite. Mastouri and Zhu et al. [18] consider an infinite number of basis function in a reproducing kernel Hilbert space, and use the definitions of the inner and Kronecker products for that space.
As in the existing work [3, 18], we learn EW |a,z [ψθW (w)] in Stage 1 and h(a,w) in Stage 2, but in addition to the weights u and V , we also learn the parameters of the feature maps. Specifically, in Stage 1, we learn V and parameters θA(1), θZ by minimizing the following empirical loss:
L̂1(V , θA(1), θZ) = 1
m m∑ i=1 ∥∥ψθW (wi)−V (φθA(1)(ai)⊗φθZ (zi))∥∥2 + λ1‖V ‖2. Note that L̂1 is an empirical estimate of L1 in (5) with adaptive feature maps. Although L̂1 depends on θW , at this stage, we do not optimize θW with L̂1 as ψθW (w) is the “target variable” in Stage 1. Given the minimizers (V̂ , θ̂A(1), θ̂Z) = arg min L̂1, we learn weights u and parameters θW , θA(2) by minimizing the empirical stage 2 loss,
L̂2(u, θW , θA(2)) = 1
n n∑ i=1 ( ỹi−u> ( ψθA(2)(ãi)⊗ V̂ ( φθ̂A(1)(ãi)⊗φθ̂Z (z̃i) )))2 + λ2‖u‖2.
Again, the loss L̂2 is an empirical estimate of L2 in (6) with adaptive feature maps. Although the expression of L̂2 does not explicitly contain θW , it implicitly depends on θW as (V̂ , θ̂A(1), θ̂Z) is the solution of a minimization problem involving θW . This implicit dependency makes it challenging to update θW , as we cannot directly obtain its gradient. One possible solution is to use the implicit gradient method [16], but this approach might be too computationally expensive. Instead, we use the method proposed in Xu et al. [35], in which we ignore the dependency of θW on parameters θ̂A(1), θ̂Z , and compute the gradient via the closed-form solution of V̂ .

heading: 3.2 Deep Feature Proxy Variable Method
text: We now describe the learning procedure for DFPV. First, we fix parameters in the adaptive feature maps (θA(1), θZ , θA(2), θW ). Given these parameters, optimal weights V̂ , û can be learned by minimizing the empirical stage 1 loss L̂1 and empirical stage 2 loss L̂2, respectively. These minimizations can be solved analytically, where the solutions are
V̂ (θ) = Ψ>1 Φ1(Φ > 1 Φ1 +mλ1I) −1, û(θ) = ( Φ>2 Φ2 +nλ2I )−1 Φ>2 y2, (7)
where we denote θ = (θA(1), θZ , θA(2), θW ) and define matrices as follows:
Ψ1 = [ψθW (w1), . . . ,ψθW (wm)] > , Φ1 = [v1(a1, z1), . . . ,v1(am, zm)] >, y2 = [ỹ1, . . . , ỹn] >, Φ2 = [v2(ã1, z̃1), . . . ,v2(ãn, z̃n)] >,
v1(a, z) = φθA(1)(a)⊗φθZ (z), v2(a, z) = ψθA(2)(a)⊗ ( V̂ (θ) ( φθA(1)(a)⊗φθZ (z) )) .
Given these weights û(θ), V̂ (θ), we can update feature parameters by a gradient descent method with respect to the residuals of the loss of each stage, while regrading V̂ and û as functions of parameters. Specifically, we take the gradient of the losses
L̂DFPV1 (θ) = 1
m m∑ i=1 ∥∥∥ψθW (wi)− V̂ (θ)(φθA(1)(ai)⊗φθZ (zi))∥∥∥2 + λ1‖V̂ (θ)‖2, L̂DFPV2 (θ) = 1
n n∑ i=1 ( ỹi− û(θ)> ( ψθA(2)(ãi)⊗ V̂ (θ) ( φθA(1)(ãi)⊗φθZ (z̃i) )))2 +λ2‖û(θ)‖2,
where V̂ (θ), û(θ) are given in (7). Given these losses, (θA(1), θZ) are minimized with respect to L̂DFPV1 (θ), and (θA(2), θW ) are minimized with respect to L̂DFPV2 (θ). Finally, we take the empirical mean ofψ
θ (t) W
based on additional output-proxy data SW = {wextrai } nW i=1, which is used for estimating
the structural function. Here, we assume access to SW for proving consistency results, but empirically, we can use stage 1 data to compute this mean.
The complete procedure is presented in Algorithm 1. Note that we may use any sophisticated gradient-based learning method to optimize, such as Adam [10]. As reported in Xu et al. [35], we observe that the learning procedure is stabilized by running several gradient descent steps on the stage 1 parameters (θA(1), θZ) before updating the stage 2 features (θA(2), θW ). Furthermore, we
may use mini-batch updates, which sample subsets of the data at the beginning of each iteration and only use these subsamples to update the parameters.
Algorithm 1 Deep Feature Proxy Causal Learning Input: Stage 1 data S1 = {ai, zi,wi}, Stage 2 data S2 = {ãiz̃i, ỹi}, Additional outcome-
proxy data SW = {wextrai }, Regularization parameters (λ1, λ2). Initial values θ(0) = (θ
(0) A(1), θ (0) Z , θ (0) A(2), θ (0) W ). Learning rate α.
Output: Estimated structural function f̂struct(a) 1: t← 0 2: repeat 3: Compute V̂ (θ(t)), û(θ(t)) in (7) 4: Update parameters in features θ(t+1) ← (θ(t+1)A(1) , θ (t+1) Z , θ (t+1) A(2) , θ (t+1) W ) as follows
θ (t+1) A(1) ← θ (t) A(1)−α∇θA(1)L̂ DFPV 1 (θ)|θ=θ(t) , θ (t+1) Z ← θ (t) Z −α∇θZ L̂ DFPV 1 (θ)|θ=θ(t)
θ (t+1) A(2) ← θ (t) A(2)−α∇θA(2)L̂ DFPV 2 (θ)|θ=θ(t) , θ (t+1) W ← θ (t) W −α∇θW L̂ DFPV 2 (θ)|θ=θ(t)
5: Increment counter t← t+ 1; 6: until convergence 7: Compute û(θ(t)) from (7) 8: Compute mean feature for W using stage 1 dataset: µθW ← 1n ∑ ψ θ (t) W (wextrai )
9: return f̂struct(a) = (û(t))> ( ψ θ̂ (t)
A(2)
(a)⊗µθW )

heading: 3.3 Consistency of DFPV
text: In this section, we show that the DFPV method yields a consistent estimate of the bridge function. Our approach differs from the consistency result of Deaner [3] and Mastouri and Zhu et al. [18], which depend on the assumption that the true bridge function lies in a specific functional space (such as an RKHS [18]) or satisfies certain smoothness properties [3]. Since our features are adaptive, however, we instead build our result on a Rademacher complexity argument [24], which holds for a wide variety of hypothesis spaces. Furthermore, we do not assume that our hypothesis space contains the true bridge function, which makes this theory of independent interest.
To derive the excess risk bound, we need another completeness assumption. Assumption 3 (Completeness Assumption on Outcome-Inducing Proxy [3, 18]). Let l :W → R be any square integrable function ‖l‖P(W ) <∞, We assume the following condition:
E [l(W ) | A = a,Z = z] = 0 ∀(a, z) ∈ A×Z ⇔ l(w) = 0 P (W ) -a.e.
This assumption is not necessary for identification (Proposition 1), but we need it for connecting the two-stage loss and the deviation of the bridge function. Given this assumption, we have the following consistency result. Proposition 2. Let Assumption 3 and Assumptions 7 and 8 in Appendix C hold. Given stage 1 data S1 = {(wi, ai, zi)}mi=1, stage 2 data S2 = {(ỹi, ãi, z̃i)}ni=1, and additional output-proxy data SW = {wextrai } nW i=1, then for the minimizer of L̂2 denoted as (û, θ̂A(2), θ̂Z), we have
‖fstruct− f̂struct‖P(A) ≤ Op (√ κ1 + R̂S1(H1) +
1√ m + κ2+
√ R̂S2(H2) +
1√ n + 1 √ nW
) ,
where H1,H2 are the functional classes that are considered in each stage, comprising both the predictions and their associated losses, which are defined in (12) and (15), respectively; and R̂S1(H1), R̂S2(H2) are their respective empirical Rademacher complexities; the quantities κ1, κ2 measure the misspecification in each regression, and are defined in (13) and (16), respectively.
We formally restate Proposition 2 in Appendix C (Theorem 3), in which the definitions of hypothesis spaces and constants can be found. From this, we can see that if we correctly identify the hypoth-
esis space (κ1, κ2 = 0) and have vanishing Rademacher complexities (R̂S1(H1), R̂S2(H2)→ 0 in probability), then the estimated bridge function converges in probability to the true one. Note that the Rademacher complexity of certain neural networks is known; e.g., ReLU two-layer networks are known to have the complexity of order O( √ 1/n) [26]. However, it might be not straightforward to derive the Rademacher complexities of the specific classes of neural networks we use, because we are employing the outer product of two neural networks in (5) and (6), whose Rademacher complexity remains a challenging open problem. The analysis in Proposition 2 is similar to Xu et al. [35], with two important distinctions: first, the target fstruct is estimated using the bridge function h marginalized over one of its arguments W , but PCL learns the bridge function itself. Second, our result is more general and includes the misspecified case, which was not treated in the earlier work.
One limitation of our analysis is that we leave aside questions of optimization. As we previously discussed, L̂2 does not explicitly include θW , thus it is difficult to develop a consistency result that includes the optimization. We emphasize that Algorithm 1 does not guarantee to converge to the global minimizer, since we ignore the implicit dependency of θW on θA(1), θZ when calculating the gradient. As we will see in the Section 4, however, the proposed method outperforms competing methods with fixed feature dictionaries, which do not have this issue.

heading: 3.4 DFPV for Policy Evaluation
text: PCL methods can estimate not only the structural function, but also other causal quantities, such as Conditional Average Treatment Effect (CATE) [29] and Average Treatment Effect on Treated (ATT) [3, 29]. Beyond the causal context, [32] uses a PCL method to conduct the off-policy evaluation in a partially observable Markov decision process with discrete action-state spaces. Here, we present the application of PCL methods to bandit off-policy evaluation with continuous actions, and propose a method to leverage adaptive feature maps in this context.
In bandit off-policy evaluation, we regard the outcome Y as the reward and aim to evaluate the value function of a given policy. Denote policy π : C → A and define the value function v(π) as
v(π) = EU,C [E [Y |A = π(C),U ]] , where C ∈ C is a variable on which the policy depends. In the flight ticket sales prediction example, a company might be interested in predicting the effect of discount offers. This requires us to solve the policy evaluation with C = A since a new price is determined based on the current price. Alternatively, if the company is planning to introduce a new price policy that depends on the fuel cost, we need to consider the policy evaluation with C = Z, since the fuel cost can be regarded as the treatment-inducing proxy variable. Here, we assume the policy π to be deterministic, but this can be easily generalized to a stochastic policy. Although the value function v contains the expectation over the latent confounder U , it can be estimated through a bridge function.
Proposition 3. Assume that the true bridge function h∗(a,w) : A×W → R is jointly measurable. Suppose C⊥⊥W |U . Under Assumptions 1, 2 and Assumptions 4, 5, 6 in Appendix B, we have
v(π) = EC,W [h∗(π(C),W )] .
The proof can be found in Appendix B (Corollary 2). From Proposition 3, we can obtain the value function by taking the empirical average over the bridge function h(π(C),W ). In DFPV, we estimate the bridge function as h(a,w) = (û)>(ψθ̂A(2)(a)⊗ψθ̂W (w)) where û, θ̂A(2), θ̂W minimize L̂2. Hence, given n′ observation of (C,W ) denoted as {či, w̌i}n ′
i=1, we can estimate the value function v(π) as
v̂(π) = 1
nC n′∑ i=1 û> ( ψθ̂A(2)(π(či))⊗ψθ̂W (w̌i)) ) .
We derive the following consistency result for policy value evaluation in a bandit with confounding bias.
Proposition 4. Let Assumption 3 and Assumptions 7, 8, 9 in Appendix C hold. Denote the minimizers of L̂2 as (û, θ̂A(2), θ̂Z). Given stage 1 data S1 = {(wi, ai, zi)}mi=1, stage 2 data
S2 = {(ỹi, ãi, z̃i)}ni=1, and data for policy evaluation S3 = {(či, w̌i)}n ′ i=1, we have
|v(π)− v̂(π)| ≤ Op (√ κ1 + R̂S1(H1) +
1√ m + κ2+
√ R̂S2(H2) +
1√ n + 1√ n′ ) where H1,H2 are the function classes that are considered in each stage, comprising both the predictions and their associated losses, which are defined in (12) and (15), respectively; and their empirical Rademacher complexities are given by R̂S1(H1), R̂S2(H2); the quantities κ1, κ2 measure the misspecification in each regression, which are defined in (13) and (16), respectively.
The formal statement of Proposition 4 can be found in Appendix C (Theorem 3). As in the structural function case, we can see that estimated value function converges in probability to the true one, provided that we have vanishing Rademacher complexities and there is no misspecification.

heading: 4 Experiments
text: In this section, we report the empirical performance of the DFPV method. First, we present the results of estimating structural functions; we design two experimental settings for low-dimensional treatments and high-dimensional treatments, respectively. Then, we show the result of applying PCL methods to the bandit off-policy evaluation problem with confounding. We include the results for problems considered in prior work in Appendix E. The experiments are implemented using PyTorch [27]. The code is included in the supplemental material. All experiments can be run in a few minutes on Intel(R) Xeon(R) CPU E5-2698 v4 2.20GHz.
Experiments for Structural Function We present two structural function estimation experiments. One is a demand design experiment based on a synthetic dataset introduced by Hartford et al. [6], which is a standard benchmark for the instrumental variable regression. Here, we modify the data generating process to provide a benchmark for PCL methods. We consider the problem of predicting sales Y from ticket price P , where these are confounded by a potential demandD ∈ [0,10]. To correct this confounding bias, we use the fuel price (C1,C2) as the treatment-inducing proxy, which has an impact on price P , and the number of views of the ticket reservation page V as the outcome-inducing proxy. Details of the data generation process can be found in Appendix F.1.
Our second structural function estimation experiment considers high-dimensional treatment variables. We test this using the dSprite dataset [19], which is an image dataset described by five latent parameters (shape, scale, rotation, posX and posY). The images are 64× 64 = 4096-dimensional. Based on this, Xu et al. [35] introduced the causal experiment, where the treatment is each figure, and the confounder is posY. Inspired by this, we consider the PCL setting that learns the same structural functions with nonlinear confounding, which is not possible to handle in the instrumental variable setting. Specifically, the structural function fstruct and outcome Y are defined as
fstruct(a) = ‖Ba‖22− 5000
1000 , Y =
1
12 (posY− 0.5)fstruct(A) + ε, ε ∼ N (0,0.5),
where each element of the matrix B ∈ R10×4096 was generated from Unif(0.0,1.0) and fixed throughout the experiment. We fixed the shape parameter to heart and used other parameters as the treatment-inducing proxy Z. We sampled another image that shared the same posY as treatment A, which is used as output-inducing proxy W . Details of data generation process can be found in Appendix F.2.
We compare the DFPV method to three competing methods, namely KPV [18], PMMR [18], and an autoencoder approach derived from the CEVAE method [17]. In KPV, the bridge function is estimated through the twostage regression as described in Section 2, where feature functions are fixed via their kernel functions. PMMR also models the bridge function using kernel functions, but parameters are learned by moment matching. CEVAE is not a PCL method, however, it represents a state-of-the-art approach in correcting for hidden confounders using observed proxies. The causal graph for CEVAE is shown in Figure 2, and CEVAE uses a VAE [11] to recover the distribution of confounder U from the “proxy” Q. We make two modifications to CEVAE to apply it in our setting. First, we include
both the treatment-inducing proxy Z and output-inducing proxy W as Q in CEVAE (we emphasize that this does not follow the causal graph in Figure 2, since there exist arrows from Q to A,Y ). Second, CEVAE is originally used in the setting where Q is conditioned on a particular value, whereas we marginalize Q. See Appendix F.4 for the choice of the network structure and hyper-parameters. We tuned the regularizers λ1, λ2 as discussed in Appendix A, with the data evenly split for Stage 1 and Stage 2. We varied the dataset size and ran 20 simulations for each setting. Results are summarized in Figure 3.
In both experiments, DFPV consistently outperforms existing methods. This suggests that DFPV is capable of learning complex structural functions by taking the advantage of the flexibility of neural networks. KPV and PMMR perform similarly in all settings, but KPV tends to perform slightly better when the data size is large. This might be because KPV expresses its solution using a larger number of parameters than PMMR. Although CEVAE also learns a flexible model with a neural network, it’s highly unstable in the demand design experiment and underperforms kernel methods in the dSprite experiment. This is because CEVAE does not take advantage of the relations between proxy variables and the structural function.
Experiments for Offline Policy Evaluation We now describe the offline policy evaluation experiment based on the demand design data. We set up synthetic experiments to evaluate two types of policy. In the first case, a policy depends on costs C = (C1,C2), to address the question How much would we gain if we decide the price based on the fuel costs. In the second case, the policy depends on current action C = P , which addresses How much would we gain if we cut the price by 30 percent. The detailed policies can be found in Appendix F.3.
We only compare PCL methods here, as it is not straightforward to apply CEVAE to the off-policy evaluation problem. We evenly split the data for Stage 1, Stage 2, and policy evaluation (i.e we set n = m = n′).We ran 20 simulations for each setting. Results are summarized in Figure 4, in which DFPV performs better than existing PCL methods. This is not surprising since, as demonstrated in the structural function experiments, DFPV can estimate complex bridge functions, which results in a more accurate estimation of the value function.

heading: 5 Conclusion
text: We have proposed a novel approach for proxy causal learning, the Deep Feature Proxy Variable (DFPV) method, which performs two-stage least squares regression on flexible and expressive features. Motivated by the literature on the instrumental variable problem, we showed how to learn these feature maps adaptively with deep neural networks. We also showed that PCL learning can be used for off-policy evaluation in the bandit setting with confounding, and that DFIV performs competitively in this domain. This work thus brings together research from the worlds of deep offline RL and causality from observational data.
In future work, it would be interesting to explore different approaches to learning deep models in the PCL problem. One direction would be to adapt the method of moment matching, which has been studied extensively in the instrumental variable literature [1]. Moreover, Liao et al. [15] recently developed a novel adversarial method to solve a class of functional equations, which would be also a promising approach in the PCL problem. In the RL context, Tennenholtz et al. [32] shows an interesting connection between PCL learning and off-policy evaluation in a partially observable Markov decision process, and we believe that adapting DFPV to this setting will be of value.

Review summaries:
This work extends known Proxy Causal Learning methods that aim at learning the causal effect of treatment in cases where unobserved confounders are assumed yet there is observable information in form of proxy variables to these confounders. In particular, authors extend the theory to allow deep neural nets to be used as building blocks of the two stage regression proposed by Mastouri et al instead of RKHS. The paper also explains how to use the method to perform off-policy evaluation in a bandit setting where arms are confounded. Synthetic experiments are also reported to compare the performance with a few baselines of the same family.

This paper presents a novel method: deep feature proxy variable. It applies a neural network to PCL, allowing non-linear relationships between variables. The authors prove the consistency of the method and demonstrate how this method can be applied to off-policy evaluation. The authors test the methods using two simulation studies, showing that DFPV outperforms existing methods. 

The paper proposes an algorithm that operationalize the proxy variable strategy for causal identification. It builds on existing approaches (Deaner [3] and Mastouri and Zhu et al. [18]) and additionally consider adaptive features as in Deep Feature Instruemtal Variable method [32]. The key idea is to consider features parameterized by neural nets, and alternate between learning $V$ by minimizing the empirical loss and learning the parameters of the feature maps. The paper shows that the additional flexibility of the features helps the proxy variable method achieve better out of sample performance and smaller errors in empirical studies.