Latent Variable Modeling for Robust Causal Effect Estimation

TL;DR

This paper introduces a novel framework combining latent variable models with double machine learning to improve causal effect estimation when unobserved confounders are present, demonstrating robustness through extensive experiments.

Contribution

It integrates latent variable modeling into DML, addressing unmeasured confounders in causal inference with a focus on separating representation learning from latent inference.

Findings

Enhanced robustness in causal effect estimation with unobserved confounders.

Effective separation of latent inference from representation learning.

Validated performance on synthetic and real-world datasets.

Abstract

Latent variable models provide a powerful framework for incorporating and inferring unobserved factors in observational data. In causal inference, they help account for hidden factors influencing treatment or outcome, thereby addressing challenges posed by missing or unmeasured covariates. This paper proposes a new framework that integrates latent variable modeling into the double machine learning (DML) paradigm to enable robust causal effect estimation in the presence of such hidden factors. We consider two scenarios: one where a latent variable affects only the outcome, and another where it may influence both treatment and outcome. To ensure tractability, we incorporate latent variables only in the second stage of DML, separating representation learning from latent inference. We demonstrate the robustness and effectiveness of our method through extensive experiments on both synthetic and real-world datasets.

Tables

Table 1. Table 1. Causal effect estimation using the 401(k) dataset with Monte Carlo resampling. Bias values are reported as “bias ± 95% CI width,” with a reference causal effect set at $14.75K based on the DML estimate on full data without missing covariates. All values in the table are in thousands of dollars. Latent DML (Ours) denotes the BIC-selected latent DML method.

Covariates	Method	Bias
Full	DML	0.05 $\pm$ 0.86
Full	Coordinated DML	-0.79 $\pm$ 1.26
Full	Latent DML (Ours)	-0.76 $\pm$ 0.82
w/o e401	DML	-2.96 $\pm$ 0.70
w/o e401	Coordinated DML	-1.84 $\pm$ 0.89
w/o e401	Latent DML (Ours)	-1.64 $\pm$ 1.85
w/o pira	DML	2.66 $\pm$ 0.89
w/o pira	Coordinated DML	1.69 $\pm$ 1.17
w/o pira	Latent DML (Ours)	1.22 $\pm$ 0.78

Table 2. Table 2. Experimental results: causal effect of advertisement timing on ROAS (percentage per 10 days). Latent DML uses BIC to choose between outcome and confounder latent models. Based on this selection, the estimated causal effect of Latent DML is therefore − 60.37 -60.37 .

Method	Estimated Effect	BIC
Ordinary DML	29.45	877.55
Outcome Latent DML	-60.37	839.26
Confounder Latent DML	0.91	858.66

Table 3. Table 3. List of covariates for real data application: causal effect analysis of advertisement timing to ROAS

Variable Name	Description
DaysPostThresholdToAdStart	Days from surpassing a sales threshold to starting advertising
DaysToThresholdForAd	Days to reach a sales threshold for advertisements to be OK
AdStartsWithinEarlySalesPhase	Indicator if advertisement starts within the early sales phase
AdEndsWithinLateSalesPhase	Indicator if advertisement ends within the late sales phase
TotalDaysOfAd	Total days of running advertisements
LogSystemAvgSales	Log of average sales in the system
LogSystemStdDevSales	Log of sales standard deviation in the system
LogSalesInitialPeriod	Log of sales in the initial period
LogSalesRatioMidToInitialPeriod	Log of sales ratio from middle to initial period
LogUnitsSoldInitialPeriod	Log of units sold in the initial period
LogUnitsSoldRatioMidToInitialPeriod	Log of units sold ratio from middle to initial period
LogPreAdAvgSalesRatioOverInitialPeriod	Log of average sales ratio before ad start over the initial period
PreAdStartSalesSlope	Relative sales slope before advertisement start
PreAdEndSalesSlope	Relative sales slope before advertisement end
LogPreAdStartUnitsSoldRatioOverInitialPeriod	Log of average units sold ratio before ad start over the initial period
LogPreAdEndUnitsSoldRatioOverInitialPeriod	Log of average units sold ratio before ad end over the initial period
PreAdStartUnitsSoldSlope	Relative units sold slope before advertisement start
PreAdEndUnitsSoldSlope	Relative units sold slope before advertisement end
HolidayCountDuringProject	Number of holidays during the project period
LogTotalSales	Log of total sales
LogSalesPostInitialPeriod	Log of sales after the initial period
AvgLogReturnAmount	Average of Log of return amount
StdDevLogReturnAmount	Standard deviation of Log of return amount
LogNumberOfReturns	Log of number of returns
ProjectDuration	Duration of the project
NumberOfAdRuns	Number of times advertisements were run
StartPeriod_Q1	Indicator if project started in Q1
StartPeriod_Q2	Indicator if project started in Q2
StartPeriod_Q3	Indicator if project started in Q3
ProjectStartDay_Holiday	Indicator if project started on a holiday
ProjectStartDay_PreHoliday	Indicator if project started on a day before a holiday
Curator_A	Indicator for Curator A
Curator_B	Indicator for Curator B
Curator_C	Indicator for Curator C
Curator_D	Indicator for Curator D
Curator_E	Indicator for Curator E
Curator_F	Indicator for Curator F

Equations

$$Y_i$$

$$D_i$$

$$\mathbb{E}[U_i|X_i]=\mathbb{E}[V_i|X_i]=0$$

$$\mathbb{E}[U_i{V}_i|X_i]=0.$$

$$h(X):=\mathbb{E}[Y|X]=g(X)+\theta m(X).$$

$$\hat{R}=Y-\hat{h}(X),\hat{V}=D-\hat{m}(X).$$

$$R:=Y-h(X)=\theta V+U,$$

$$\hat{\theta }:={\left({\sum }_{i=1}^N{\widehat{V_i}}^2\right)}^{-1}{\sum }_{i=1}^N\widehat{V_i}\widehat{R_i}$$

$$Q(\gamma\,|\,\gamma^{(t)})=\mathbb{E}_{Z\sim p(Z\,|\,X,\gamma^{(t)})}\big{[}\log p(X,Z\,|\,\gamma)\big{]}.$$

$${\gamma }^{(t+1)}=\arg \underset{\gamma }{\max }Q(\gamma |{\gamma }^{(t)}).$$

$$U=f_u(Z)+W_u,V=f_v(Z)+W_v.$$

$$L(\gamma )$$

$$=-\log {\int }_{Z}p(R,V|Z,\gamma )p(Z|\gamma )dZ.$$

$${\sum }_{i=1}^N({\hat{R}}_i^{\mathrm{z}}-\hat{\theta }{\hat{V}}_i){\hat{V}}_i=0,$$

$${\hat{R}}_i^{\mathrm{z}}$$

$$R^z-\theta V=f_u(Z)-\mathbb{E}[f_u(Z)|R,V,{\gamma }_0]+W_u,$$

$$\mathbb{E}[R^z-\theta V|V]$$

$$=0,$$

$$U$$

$$Z\sim \mathcal{E}(\beta )-\beta ,$$

$$p(R,Z|V,\theta ,\beta ,\sigma ):=\mathcal{E}(Z+\beta ;\beta )\mathcal{N}(R;\theta V+Z-\beta ,{\sigma }^2)$$

$$p(R|V,\theta ,\beta ,\sigma )={\int }_{Z}p(R,Z|V,\theta ,\beta ,\sigma )dZ,$$

$$m_i:=\widehat{R_i}-\theta {\hat{V}}_i+\beta -{\sigma }^2/\beta ,$$

$$U=aZ+W_u,V=bZ+W_v,$$

$$Z\sim \mathcal{B}(q)-q,$$

$$p(R,V,Z|\theta ,a,b,q,{\sigma }_u,{\sigma }_v)$$

$$:=\mathcal{B}(Z;q)\mathcal{N}(V;bZ,{\sigma }_v^2)\mathcal{N}(R;\theta V+aZ,{\sigma }_u^2).$$

$${\hat{R}}_i^{\mathrm{z}}:=\widehat{R_i}-a({\pi }_i-q),$$

$${\pi }_i$$

$$\displaystyle\hskip 2.84526pt=\frac{1}{1\!+\!\exp\!\!\bigg{(}\!\!-\frac{(2q-1)a^{\!2}+2a(\hat{R_{i}}-\theta\hat{V}_{i})}{2\sigma_{u}^{2}}-\frac{(2q-1)b^{\!2}+2b\hat{V}_{i}}{2\sigma_{v}^{2}}\!\bigg{)}}.$$

$$\frac{1}{n}\sum _{i=1}^n\psi (W_i,\gamma ,\nu ),$$

$$\Big{(}\mathbb{E}_{X}\Big{[}\big{\|}\hat{\nu}(X)-\nu_{0}(X)\big{\|}^{2}\Big{]}\Big{)}^{1/2}=o_{p}(n^{-1/4}),$$

$$\displaystyle\frac{\partial}{\partial s}\mathbb{E}[\psi(W,\gamma,\nu+s\xi)]\big{|}_{s=0}=0,$$

Full text

latent variable modeling, causal inference, double machine learning

††copyright: acmlicensed††journalyear: 2025††conference: Conference on Information and Knowledge Management; November 10–14, 2025; Coex, Seoul, Korea

1. Introduction

Central to data science are both predictive analytics and causal inference. The primary goal of the latter is to estimate the effect of one variable, known as a treatment, on another, known as an outcome, using observational data (Holland, 1986; Pearl, 2009; Athey, 2017). This task is essential in various domains, where it is required to understand how a treatment impacts an outcome of interest amidst numerous covariates (Imbens and Wooldridge, 2009; Yao et al., 2021). For instance, understanding the efficacy of medical treatments is paramount in healthcare. Researchers often seek to estimate the effect of drug dosage on patient outcomes, such as recovery rates or symptom alleviation, while accounting for patient-specific characteristics and clinical histories (Castro et al., 2020; Frauen et al., 2023). Similarly, in marketing, determining the impact of advertising campaigns on consumer behavior is crucial. Analysts might analyze how different levels of ad exposure influence consumer purchases, taking into account various factors such as demographic data and previous buying behavior (Brodersen et al., 2015; Zhang et al., 2021).

In practical causal inference settings, ensuring the completeness of covariates is often challenging, as highlighted in several studies (Gilchrist and Sands, 2016; Davis and Heller, 2017; er Chen et al., 2021). To address this, we introduce latent variables to account for unobserved factors that may influence the outcome or treatment but are not captured by observed covariates. This approach allows us to incorporate domain knowledge about the data-generating process, such as structured noise, into the model. Even when such knowledge is limited, we propose model selection using information criteria to choose appropriate latent structures.

While latent variable modeling offers a powerful tool for incorporating unobserved factors into causal inference, it often involves iterative procedures, such as alternating between model parameter updates and latent variable inference, that can be computationally intensive. This computational challenge is particularly pronounced when combined with large-scale machine learning models, which are increasingly used in causal inference and already require substantial optimization efforts. To mitigate this, we adopt the double machine learning (DML) framework (Chernozhukov et al., 2018), which decouples the estimation of nuisance components (e.g., outcome and treatment models) from the final causal estimation. Our key idea is to introduce latent variables only into the second stage of DML, enabling efficient estimation of causal effects, while keeping the overall computation tractable. We refer to this framework as latent double machine learning (latent DML).

We consider two distinct modeling scenarios involving latent variables. The first scenario assumes that the latent variable affects only the outcome. For instance, in advertising effectiveness estimation, unobserved social media buzz can significantly amplify the observable impact, such as Return on Advertising Spend (ROAS), without directly influencing the advertising strategies (treatment). The second scenario assumes that the latent variable influences both the treatment and the outcome, acting as an unobserved confounder. One example is a competitor’s unrecorded marketing campaign, which may affect both a store’s marketing decisions (treatment) and its sales figures (outcome). Another example arises in program evaluation studies, where unobserved participant characteristics (e.g., parental income) may influence both the duration of program participation (treatment) and long-term outcomes, potentially biasing the causal effect estimation. Figure 1 illustrates the second scenario by comparing standard DML and our proposed latent DML in the presence of unobserved confounding (see Section 5.1 for details).

Contributions of this paper include:

• •

Integration of latent variable modeling into DML: We propose latent double machine learning (latent DML), a novel framework that integrates latent variable modeling into the DML paradigm to address unobserved factors in causal inference. To the best of our knowledge, this is the first work to develop a general DML-based method tailored to latent variable models.

• •

Algorithmic implementation: We present two implementations of latent DML, including adaptations of the EM algorithm.

• •

Comprehensive validation: We demonstrate the effectiveness of latent DML through extensive empirical validation and provide theoretical guarantees.

2. Preliminary

This section provides an overview of the key concepts in DML and latent variable modeling, which are essential for our paper.

2.1. Double machine learning

Double machine learning (DML) is a framework that integrates machine learning techniques with causal inference, primarily focusing on estimating treatment effects in the presence of high-dimensional covariates (Chernozhukov et al., 2018). This approach builds upon the partially linear model proposed by (Robinson, 1988), enhanced by the incorporation of cross-fitting and the concept of Neyman orthogonality.

A key characteristic of Robinson’s partially linear model is its ability to separate the estimation of causal effects from the nuisance components, which are parts of the model not directly related to the causal effect but potentially correlated with the treatment and outcome. The model is formulated as follows for $N$ i.i.d. data points $\{(X_{i},D_{i},Y_{i})\}_{i=1}^{N}$:

[TABLE]

where $X_{i}\in\mathbb{R}^{d}$ is the covariates, $D_{i}\in\mathbb{R}$ is the treatment variable, $Y_{i}\in\mathbb{R}$ is the outcome, and $\theta\in\mathbb{R}$ is the average causal effect. The unknown functions $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ and $m:\mathbb{R}^{d}\rightarrow\mathbb{R}$ represent the influence of covariates on the outcome and treatment, respectively. The inferential target is the causal effect parameter $\theta$, while $g$ and $m$ are nuisance components, in the sense that their estimation is not of direct interest. The noise terms $U_{i}$ and $V_{i}$ are assumed to be zero mean and independent of the covariates, that is,

[TABLE]

Equation (4) ensures there is no unobservable confounder and is known as the unconfoundedness condition.

In the DML framework following Robinson’s approach, the first step is to estimate the functions $m(X)=\mathbb{E}[D\,|\,X]$ and

[TABLE]

This is achieved by employing machine learning models to approximate these functions. Specifically, a model is used to approximate the expected treatment as $\hat{m}(X)$ and another model to approximate the expected outcome as $\hat{h}(X)$. Once these models are fitted, the estimated functions $\hat{m}(X)$ and $\hat{h}(X)$ are used to compute the residuals:

[TABLE]

It is important to note that $\hat{R}$ represents the residual of the outcome, differing from the original noise term $U_{i}$ in the model. Specifically, the outcome residual is computed as

[TABLE]

reflecting the influence of both the causal effect and the original noise. In contrast, the residual of the treatment corresponds directly to the original noise term $V$. Finally, the causal effect estimate in DML is obtained using the formula:

[TABLE]

This approach effectively isolates the causal effect $\theta$ from the nuisance components represented by $g(X)$ and $m(X)$.

Other feature of DML is the implementation of the cross-fitting procedure. In this process, the data is divided into several folds. Each fold is used to calculate residuals $\hat{U},\hat{V}$ for the estimation of causal effects, while the remaining data serve as the training set for the machine learning models $\hat{h},\hat{m}$. For a detailed description of this procedure, see Algorithm 1. This technique is pivotal in reducing overfitting, particularly in high-dimensional settings, and enhances the convergence rate of the estimators, thereby improving the overall robustness of the causal effect estimation.

2.2. Latent variable modeling

Latent variable modeling is a statistical technique that introduces inductive biases into probabilistic models by incorporating hidden variables, providing additional structure and enhancing model estimation. This approach is often applied in scenarios where direct measurement of certain variables is not feasible, as it helps to account for hidden factors influencing the observed data. Common methods for latent variable models include the Expectation-Maximization (EM) algorithm, variational Bayes, and Markov chain Monte Carlo methods (Bishop, 2006).

We will use the EM algorithm in this study for its effectiveness and simplicity. Here, we assume a parametric stochastic model $p(\{(X_{i},Z_{i})\}_{i=1}^{N}\,|\,\gamma)=\prod_{i=1}^{N}p(X_{i},Z_{i}\,|\,\gamma)$, where $X_{i}$, $Z_{i}$, and $\gamma$ are observed data, latent variable, and model parameter, respectively. For simplicity, we may notate $p(\{(X_{i},Z_{i})\}_{i=1}^{N}\,|\,\gamma)$ as $p(X,Z\,|\,\gamma)$. The EM algorithm iteratively performs two steps until convergence to locally maximize the log likelihood function $\log p(X\,|\,\gamma)$:

• •

Expectation step (E-step): Compute the following expectation with respect to the conditional distribution of the latent variable $Z$ given the data $X$ and the current parameter estimate $\gamma^{(t)}$:

[TABLE]

• •

Maximization step (M-step): Update the parameter estimate by maximizing $Q(\gamma\,|\,\gamma^{(t)})$:

[TABLE]

3. Latent double machine learning

We introduce latent double machine learning (latent DML), a framework that integrates latent variable modeling into the double machine learning (DML) setting for improved causal effect estimation. Section 3.1 presents the overall design, while Sections 3.2 and 3.3 describe two concrete instantiations of the framework. Section 3.4 discusses practical implementation considerations.

3.1. Latent DML framework

We aim to address challenges in causal inference arising from unobserved factors that are not captured by observed covariates but may influence the treatment or outcome. To this end, we propose a latent DML framework that incorporates latent variable modeling into the DML process, combining the flexibility of latent variable models with the computational and statistical advantages of DML.

Since DML inherently follows a two-step procedure, there are multiple ways to incorporate latent variables. However, due to the iterative nature of latent variable modeling, incorporating it into the first step would require repeated training of machine learning models, which is computationally demanding. Therefore, we introduce a latent variable only into the second step of DML, specifically in lines $10$ and $11$ of Algorithm 1. This design allows for efficient latent variable inference without retraining the nuisance components.

In the latent DML framework, we need to specify the distribution of the outcome residual $R=\theta V+U$ in Eq. (6), which consists of the noise terms $U$ and $V$. We consider the following generalized structure for the noise models with parameterized functions $f_{\hskip-0.56905ptu\hskip-0.56905pt}$ and $f_{\hskip-0.56905ptv\hskip-0.56905pt}$ of a latent variable $Z$:

[TABLE]

Here, $Z$, $W_{u}$, and $W_{v}$ are random variables, assumed to be mutually independent and also independent of the covariate $X$. The modeling of these functions and variables, while not trivial, allows for the incorporation of domain knowledge into the model, which could be advantageous when such prior information is available (Gelman et al., 2013; Shams and Beierholm, 2022). Specific instances are discussed in the subsequent subsections, providing two distinct scenarios within our latent DML framework, an outcome-only latent variable scenario in Section 3.2, and an unobservable confounder scenario in Section 3.3.

Most importantly, this framework allows $\mathbb{E}[f_{\hskip-0.56905ptu\hskip-0.56905pt}(Z)f_{\hskip-0.56905ptv\hskip-0.56905pt}(Z)\,|\,X]$ to be nonzero. Thus, it does not require the unconfoundedness condition of Eq. (4) and can handle unobservable confounders when the noise models are well specified.

We denote the set of all parameters involved in the residual modeling as $\gamma$. The set $\gamma$ includes parameters characterizing the functions $f_{\hskip-0.56905ptu\hskip-0.56905pt}$ and $f_{\hskip-0.56905ptv\hskip-0.56905pt}$, parameters for the random variables $Z$, $W_{u}$, and $W_{v}$, and a causal effect parameter $\hat{\theta}$. The objective function for this modeling is the negative log-likelihood,

[TABLE]

We estimate the parameter $\gamma$ by minimizing $L(\gamma)$ using the residual samples $\{(\hat{R}_{i},\hat{V}_{i})\}_{i=1}^{N}$ with a latent variable modeling method such as the EM algorithm, and denote the obtained estimate as $\hat{\gamma}$. If the model space contains the true data-generating process with $\gamma_{0}$, then $\hat{\gamma}$ is expected to be a consistent estimator of $\gamma_{0}$. A more rigorous theoretical analysis of this result is provided in Appendix A.2.

Latent DML estimates the causal effect $\hat{\theta}$ by solving the following score equation, which is derived from moment conditions:

[TABLE]

where $\hat{R}^{\rm z}_{i}$ is an adjusted outcome residual defined as

[TABLE]

$\hat{\theta}$ is an unbiased estimate of the causal effect if the true model is included in the model space and the the nuisance components are estimated as $m=\hat{m}$ and $h=\hat{h}$. The unbiasedness follows from from Eqs. (6), (9), and (12), where the error term $\hat{R}^{\rm z}-\hat{\theta}\hat{V}$ under the true parameter $\gamma_{0}$ satisfies:

[TABLE]

and thus the following condition for the unbiasedness holds:

[TABLE]

since $Z$, $W_{u}$, and $W_{v}$ are assumed to be mutually independent. We summarize the procedure in Algorithm 2.

3.2. Outcome latent model

We present an outcome latent model with $Z\in\mathbb{R}$ for a outcome-only latent variable scenario, where the noise models in Eq. (9) use $f_{\hskip-0.56905ptu\hskip-0.56905pt}(Z)=Z$ and $f_{\hskip-0.56905ptv\hskip-0.56905pt}(Z)=0$, meaning that the treatment residual $V$ is independent of $Z$ and does not need to be modeled explicitly. We model the noise term $U$ as

[TABLE]

where $W_{u}$ follows a normal distribution $\mathcal{N}(0,\sigma^{2})$ of zero mean and standard deviation $\sigma$. Owing to the zero-mean property of $U$ in Eq.(1), the mean of $Z$ must also be zero.

Assuming a normal distribution for $Z$ causes an identifiability issue because it violates the non-singular condition of the Fisher information matrix in Proposition B.2. However, other choices are permissible provided that they maintain a zero mean. As one of the simplest implementations, we model $Z$ as a shifted exponential distribution. Specifically, we assume:

[TABLE]

where $\mathcal{E}(\beta)$ denotes the exponential distribution with mean parameter $\beta$.

The use of an exponential distribution for modeling $Z$ is particularly pertinent in domains where outcomes $Y$ are subject to sudden, substantial increases due to factors not captured by covariates $X$. For instance, in advertising effectiveness analysis, a viral buzz on social media can significantly increase online sales or website traffic. In such cases, the exponential distribution effectively models the potential for upward shifts in these outcomes, providing a realistic representation of these unpredictable events.

The joint distribution is, since the residual of the outcome is $R=\theta V+U$,

[TABLE]

Having defined the stochastic model with latent variables, we now turn our attention to parameter estimation. The optimal parameters are those that maximize the marginal log-likelihood

[TABLE]

and we employ the EM algorithm. Each iteration of the EM algorithm involves updating the parameters in the M-step. While most parameters can be updated in closed form during this step, the optimization of $\beta$ in Eq. 8 poses a unique challenge. This optimization problem cannot be analytically solved, so we employ the natural gradient method for updating $\beta$, using the Fisher information of $\mathcal{E}(\beta)$ as a metric.

Upon convergence of the EM algorithm, the causal effect estimate $\hat{\theta}$ is calculated by solving the score equation (11) with the adjusted outcome residual $\hat{R}^{\rm z}_{i}=\hat{R_{i}}-\mathbb{E}[Z_{i}\,|\,\hat{R_{i}},\hat{V}_{i},\theta,\beta,\sigma]$ of Eq. (12). According to the joint distribution of Eq. (15), $\mathbb{E}[Z_{i}\,|\,\hat{R_{i}},\hat{V}_{i},\theta,\beta,\sigma]$ is calculated as the mean of a truncated distribution of $\mathcal{N}(m_{i},\sigma^{2})$, where

[TABLE]

truncated over the interval $[0,\infty]$.

3.3. Confounder latent model

Unlike the previous subsection, here a latent variable $Z\in\mathbb{R}$ acts as an unobservable confounder, affecting both the treatment and the outcome. We set $f_{\hskip-0.56905ptu\hskip-0.56905pt}(Z)=aZ$ and $f_{\hskip-0.56905ptv\hskip-0.56905pt}(Z)=bZ$ in Eq. (9) as

[TABLE]

where $a\in\mathbb{R}$ and $b\in\mathbb{R}$ are unknown and represent the impact of $Z$. $W_{u}$ and $W_{v}$ follow normal distributions $\mathcal{N}(0,\sigma_{u})$ and $\mathcal{N}(0,\sigma_{v})$, respectively. We model $Z$ as a shifted Bernoulli distribution:

[TABLE]

where $\mathcal{B}(q)$ denotes the Bernoulli distribution with success probability $q\in(0,1)$, ensuring zero mean for $U$ and $V$. The choice of this model for $Z$ is particularly suitable for scenarios involving unobservable, randomly occurring factors, independent of the observed covariates $X$. A relevant example of $Z$ could be a competitor’s unrecorded marketing campaign, influencing both treatment and outcome while remaining unobserved.

Based on this setting, the joint distribution is defined as:

[TABLE]

Similar to the previous section, we employ the EM algorithm for the parameter optimization. In the M-step of this case, while most parameters can be updated in closed form, we update $q$ using the natural gradient method, due to the challenges in obtaining an analytical solution.

Once the EM algorithm have converged, the causal effect estimate $\hat{\theta}$ is calculated from the score equation (11) with the following adjusted outcome residual of Eq. (12),

[TABLE]

where $\pi_{i}$ is the following conditional probability of $Z$,

[TABLE]

This adjustment of the outcome residual $\hat{R_{i}}$ to $\hat{R}^{\rm z}_{i}$ can be regarded as removing the influence of the confounder $Z$ from the noise term $U$ in Eq. (16).

3.4. Practical considerations

As shown in Appendix A.2, under well-specified noise models, the causal effect estimate $\hat{\theta}$ asymptotically converges to the true causal effect even when there are unobservable confounders. This property is particularly beneficial in practice, as it allows for robust causal estimation while keeping residual estimation and latent variable modeling separate. By decoupling these components, our approach reduces computational cost and improves stability, particularly in large-scale models (Farrell et al., 2021). Moreover, it demonstrates that even in the absence of unobservable confounders, latent DML can reduce the variance of the causal effect estimator by leveraging additional latent structure.

Choosing an appropriate model for latent DML can be challenging due to limited domain knowledge. To address this, we propose using model selection criteria, such as the Bayesian Information Criterion (BIC) (Bishop, 2006), to identify a model that best describes the data. The effectiveness of this approach will be empirically evaluated in Section 5.1.2. Moreover, our approach remains valid even under model misspecification by approximating the underlying distribution (Huber, 1967; White, 1982).

4. Related work

Latent variable modeling in causal inference has been developed to address various challenges, such as the impact of outliers and unobserved confounders. Amjad et al. (2018); Harada and Fujisawa (2024) focus on the influence of outliers, aligning with our first scenario where a latent variable affects only the outcome. Other works such as (Louizos et al., 2017; Strobl and Lasko, 2023; Kocacoban and Cussens, 2019; Ashman et al., 2023; Kivva et al., 2023) address unobserved confounders, related to our second scenario where a latent variable influences both treatment and outcome. Notably, our use of discrete latent variables in the second scenario shares similarities with Latent Class Analysis (LCA), which identifies hidden subgroups within a population (Weller et al., 2020; Lanza and Cooper, 2016). Despite these diverse approaches in latent variable modeling, our study is the first to integrate latent variable modeling into the DML framework, allowing for explicit estimation of latent structure within a two-stage causal inference setting.

Addressing unobserved confounders is a central challenge in causal inference. Recent works like RieszNet and ForestRiesz adjust the feature space to improve robustness against hidden confounding, though they do not explicitly model latent variables (Quintas-Martinez et al., 2022; Chernozhukov et al., 2021, 2020). A recent study (Abadie et al., 2025) leverages matrix completion for doubly robust estimation under unobserved confounding, exploiting multiple outcomes to identify low-dimensional latent factors. Feature selection techniques based on orthogonal search or doubly robust methods have also been explored to reduce sensitivity to unobserved confounders (Soleymani et al., 2023; Quinzan et al., 2023). Our latent DML method is complementary to these approaches, as it directly models latent factors rather than merely mitigating their influence. Integration with such methods may further help address unobserved confounding.

Numerous DML extensions have been proposed to handle various practical challenges, including bias reduction (Fingerhut et al., 2022), conditional average treatment effects estimation (Semenova and Chernozhukov, 2021), dynamic treatment effect estimation (Lewis and Syrgkanis, 2021), difference-in-differences (Chang, 2020), off-policy evaluation (Kallus and Uehara, 2020; Narita et al., 2021), kernel-based nonparametric inference (Colangelo and Lee, 2020), and computational efficiency improvements (Kurz, 2021). Our work aligns most closely with Chiang et al. (2021), who address complex data structures such as multiway clustering within the DML setup. In contrast, we introduce latent variable modeling into DML, enabling latent structure to be inferred directly from probabilistic models. Our approach is orthogonal to many existing DML extensions and could potentially enhance their applicability to settings with latent structure.

Our approach is based on the likelihood framework for estimating residual distributions. It remains valid even under model misspecification, as we estimate models with finite-dimensional parameters to approximate the underlying distribution using a certain distance function. This argument has been widely used throughout statistics and econometrics (Beran, 1977; Buja et al., 019a, 019b; Huber, 1967; Rakhlin et al., 2017; Rinaldo and Wasserman, 2010; White, 1980, 1982). The recent literature on causal inference also applies the same argument (Chernozhukov et al., 018b; Cuellar and Kennedy, 2020; Kennedy et al., 2019; Neugebauer and van der Laan, 2007; Semenova and Chernozhukov, 2021; van der Laan and Rubin, 2006).

5. Experiments

In this section, we conduct a series of numerical experiments. Since the true causal effects in real data are generally unknown, we first perform comprehensive experiments using synthetically generated data in Section 5.1. Experiments using real data are also presented in Section 5.2.

5.1. Synthetic data

We first compare and evaluate the latent DML methods against the ordinary DML and regression analysis methods. The focus here is to specifically assess the enhancements achieved by integrating latent variables into the DML framework, rather than conducting an extensive comparison with a wide array of baseline methods. Additionally, experimental results regarding the selection of latent variable models are also presented.

The experimental setup with synthetic data was designed to rigorously evaluate the outcome latent and confounder latent models in the latent DML framework under different conditions. We conducted 100 runs with $300$ samples and $100$ covariates in each evaluation, where covariates, treatment, and outcome are generated randomly with some fixed parameters, ensuring consistent experimental conditions. In the first step of DML, ElasticNet (Hastie et al., 2009) was employed for computing the residuals. We also evaluated ElasticNet as a baseline, where the coefficient corresponding to the treatment variable served as the causal effect estimate. For more details on the experimental setup, see Appendix A.1.1 and B.1.1.

5.1.1. Causal effect estimation under latent variables

In our experiments, we first evaluated the outcome latent model as detailed in Section 3.2, performing tests in an ideal setting without model misspecification. The noise in this setup adhered to the pattern outlined in Eq. (14). As shown in Figure 2, the outcome latent DML method proved effective in estimating causal effects, exhibiting minimal bias and variance.

We then expanded our investigation to scenarios involving unobserved confounders, maintaining the ideal test conditions. This phase included testing cases where unobserved confounders influenced the outcome and treatment in both similar and opposite directions. Figure 3 shows that only the confounder latent DML method was able to accurately estimate the unbiased causal effect in these settings.

5.1.2. Model selection

The robustness of the latent DML framework against model misspecification was a key focus of our investigation. We evaluated five methods: ElasticNet, ordinary DML, outcome latent DML, confounder latent DML, and the BIC-selected latent DML. The BIC-selected latent DML method chooses either the outcome or confounder latent DML variant based on the lowest BIC score.

Data generation encompassed four different scenarios, each with distinct noise structures and parameters. This included conditions without latent variable noise as well as those with Laplace-distributed output noise, the latter representing a case of model misspecification. The results in Figure 4 demonstrate the effectiveness of our approach. They indicate that the latent DML framework is robust, facilitating reliable model selection based on BIC, even in the absence of prior knowledge and amidst the challenges posed by model misspecification.

5.2. Real data applications

We conducted experiments on two real-world datasets: (i) the 401(k) dataset (Abadie, 2003; Chernozhukov et al., 018b), which examines the effect of 401(k) participation on financial wealth (Section 5.2.1), and (ii) an analysis of advertising operations, specifically investigating how the timing of advertisements impacts Return on Advertising Spend (ROAS) (Section 5.2.2).

5.2.1. 401(k) participation and financial wealth

We evaluated our method on the well-studied 401(k) dataset (Abadie, 2003; Chernozhukov et al., 018b; Bach et al., 2023), which contains 9,915 household-level records from the 1991 Survey of Income and Program Participation (SIPP). The outcome variable in our analysis was net_tfa, representing net financial assets, and the treatment variable was p401, indicating participation in a 401(k) plan. The dataset included 10 covariates capturing demographic and financial characteristics, such as age, income, education, marital status, and home ownership.

To assess robustness against missing covariates, we considered two settings where a key variable was unobserved:

• •

e401 – indicating eligibility for a 401(k) plan, highly predictive of treatment but less directly linked to the outcome.

• •

pira – indicating participation in an IRA, which likely correlates with both treatment and outcome via financial planning behavior.

These scenarios emulated latent confounding or selection bias.

We performed Monte Carlo resampling by drawing 2,000-sample subsets from the full dataset and repeated this process 100 times. As a reference, the causal effect of 401(k) participation was estimated at $14.75K using ordinary DML on the full dataset with all covariates. Bias was reported relative to this benchmark.

We compared ordinary DML, Coordinated DML (Fingerhut et al., 2022), and our latent DML method using BIC-based model selection. Table 1 summarizes the results.

As shown in Table 1, latent DML consistently achieved lower bias than standard and coordinated DML when important covariates were omitted. The improvement was particularly notable when pira, a proxy for financial sophistication, was missing—suggesting that latent DML effectively adjusted for unobserved confounders that influenced both treatment and outcome. In contrast, Coordinated DML partially mitigated bias when the omitted variable was only related to treatment (e.g., e401), but fell short under confounding scenarios. These results illustrate the robustness of our method against latent confounding structures.

Further implementation details on Coordinated DML are provided in Appendix A.1.2.

5.2.2. Advertisement timing and ROAS

This experiment investigated the causal effect of advertisement timing on advertising effectiveness using real-world data from a domestic e-commerce platform. The treatment variable was defined as the number of days between when advertising becomes feasible (after sales exceed a threshold) and its actual commencement.

The dataset consisted of 155 cases and 45 covariates, including temporal measures, sales and advertising metrics, project-specific indicators, and curator variables (see Table 3). These covariates were carefully selected to control for potential confounders and ensure a reliable causal analysis. Advertising effectiveness was measured using Return on Advertising Spend (ROAS), a key performance indicator in digital marketing that evaluates how efficiently an advertising campaign converts spending into revenue. ROAS was computed as the ratio of revenue generated from advertisements to the cost of those advertisements.

Table 2 presents the results. The ordinary DML method produced unexpected findings, contrary to domain knowledge. Typically, a shorter delay in advertisement commencement is expected to increase ROAS, implying a negative causal effect. However, ordinary DML estimated a positive causal effect of +29.5% per 100 days, contradicting this expectation. In contrast, Latent DML with model selection based on BIC estimated a negative causal effect of -60.4% per 100 days, which aligns with prior expectations. This finding reinforces the validity of our proposed approach. These results highlight the importance of accounting for latent factors that may influence both the timing of advertising and its effectiveness.

6. Conclusion

This paper introduces a novel framework that integrates latent variable modeling into causal effect estimation, leveraging the structure of double machine learning (DML) for computational and statistical efficiency. The proposed latent DML framework accounts for unobserved factors that may influence treatment or outcome, addressing challenges posed by missing or unmeasured covariates. By introducing latent variables only in the second stage of DML, our method decouples residual estimation from latent variable inference, improving tractability and scalability without sacrificing estimation consistency. We demonstrate the effectiveness of the proposed approach through extensive empirical validation on both synthetic and real-world datasets, supported by theoretical analysis of consistency, underscoring its robustness and practical applicability. Future directions include relaxing noise-model assumptions, analyzing model misspecification, and exploring variational inference methods as an alternative to the EM algorithm.

Appendix A Appendix

Appendix is organized as follows. Section A.1 provides additional details on the simulation study and empirical analysis. Section A.2 presents theoretical results for our estimation procedure. Further details and full proofs are provided in Appendix B.

A.1. Experimental details

This section provides additional details of the experimental setup.

A.1.1. Experimental setup and model training

We describe implementation details of our experiments, including the cross-validation strategy, hyperparameter tuning, and training procedures.

In our experiments, we used ElasticNet both as a baseline method for causal effect estimation and as the ML model for the first step of the DML frameworks. ElasticNet is a linear regression model that combines both L1 and L2 regularization, with two hyperparameters: $\alpha$ and l1_ratio. The hyperparameter $\alpha$ controls the overall strength of the regularization, while l1_ratio determines the relative contribution of L1 and L2 regularization. Specifically, $\texttt{l1\_ratio}=1.0$ corresponds to Lasso regression, and $\texttt{l1\_ratio}=0.0$ corresponds to Ridge regression. For both the baseline and DML-based models, these hyperparameters were selected using 5-fold cross-validation (CV). The search space for $\alpha$ was set to $\{10^{-2},10^{-1},1,10,100\}$, and l1_ratio was varied over $\{0.0,0.25,0.5,0.75,1.0\}$. This approach allows for robust estimation by combining both types of regularization.

For the baseline ElasticNet method, we estimated the causal effect by directly using the coefficient corresponding to the treatment variable in the final fitted ElasticNet model. Unlike DML, this approach does not employ sample splitting or residual estimation, making it a simple but useful benchmark for comparison.

For the DML frameworks, ElasticNet was employed in a slightly different manner. Specifically, we used a 5-fold cross-fitting approach ($K=5$), where the model was trained on four folds, and residuals were computed on the remaining fold in a rotating manner over all five folds, ensuring that each sample was used for residual estimation once. This process is detailed in Algorithm 1 for DML and Algorithm 2 for latent DML.

A.1.2. Implementation details for 401(k) dataset experiment

We used the official implementation available at https://github.com/nitaifingerhut/C-DML for Coordinated DML (Fingerhut et al., 2022). This implementation differed slightly from ours in its treatment of sample splitting. Specifically, it used a two-split strategy, where the model was trained on one subset and residuals were computed on the other. In contrast, our DML and latent DML methods adopted a 5-fold cross-fitting approach, as described in Section A.1.1. While this difference may have caused some variations in estimation performance, we reported the Coordinated DML results without any modifications to maintain consistency with prior work.

A.2. Theoretical analysis

In this section, we present the theoretical properties of our framework, including the consistency and asymptotic normality of the latent DML estimator. The detailed identification analysis and all associated proofs, including those of the theorems stated below, are provided in Appendix B.2.

The population log-likelihood function $L(\gamma)$ is defined in Eq. (10) over the parameter space $\Gamma\subseteq\mathbb{R}^{d_{\gamma}}$, with $d_{\gamma}$ denoting the number of elements of $\gamma$. We denote the gradient by $\nabla_{\gamma}L(\gamma):=\partial L(\gamma)/\partial\gamma$, and let $\gamma_{0}\in\Gamma$ satisfy the first-order condition $\nabla_{\gamma}L(\gamma_{0})=0$.

We now present the consistency and asymptotic normality of the latent DML estimator. To derive these asymptotic properties, we express the score function (the first derivative of the log-likelihood function) by using the moment function $\psi(\cdot)$ as follows:

[TABLE]

where $\{W_{i}\}_{i=1}^{n}$ with $W_{i}:=(D_{i},X_{i},Y_{i})$ is the observed sample, $\gamma$ represents the finite-dimensional parameters, including both causal effect and late model parameters, and $\nu$ denotes the possibly infinite-dimensional nuisance parameters related to the first step in the DML procedure, or $\nu=(h,m)$. In what follows, $\gamma_{0}$ and $\hat{\gamma}$ denote the true finite-dimensional parameter and its estimator, respectively, while $\nu_{0}$ and $\hat{\nu}$ denote the nuisance parameter and its estimator.

The following conditions are imposed on the model space, observations, moment function:

Assumption A.1.

The model space is well specified, meaning that the true parameter $\gamma_{0}$ belongs to $\Gamma$, i.e., $\gamma_{0}\in\Gamma$.

Assumption A.2.

The sample $\{W_{i}\}_{i=1}^{n}$, where $W_{i}:=(D_{i},X_{i},Y_{i})$, is independently and identically distributed (i.i.d.) from $W:=(D,X,Y)$.

Assumption A.3.

The estimator $\hat{\nu}$ of the nuisance function $\nu_{0}$ is $n^{1/4}$-consistent in the squared mean-square-error sense, i.e.

[TABLE]

The nuisance function and its estimator are uniformly bounded by a constant. That is, $\|\nu_{0}(x)\|\leq C$ and $\|\hat{\nu}(x)\|\leq C$ uniformly in $x\in\mathcal{X}$ for some constant $C$ almost surely.

Assumption A.4.

The following conditions hold for the moment function $\psi$:

1. (a)

The moment function satisfies that \mathbb{E}\big{[}\psi\big{(}W,\gamma,\nu_{0}(X)\big{)}\big{]}=0 only when $\gamma=\gamma_{0}$. 2. (b)

The moment function \psi\big{(}W,\gamma,\nu_{0}(X)\big{)} is dominated by some integral function. 3. (c)

The function $\gamma\mapsto\psi(w,\gamma,\nu)$ is continuously differentiable for any $w,\nu$, and its derivative is bounded by an integrable function $c(w)$, i.e., $\|\nabla_{\gamma}\psi(w,\gamma,\nu)\|\leq c(w)$ with $\mathbb{E}[|c(W)|]<\infty$. $\mathbb{E}\left[\nabla_{\gamma}\psi(W,\gamma_{0},\nu_{0})\right]$ is non-singular. 4. (d)

The function $\nu\mapsto\psi(w,\gamma,\nu)$ is twice differentiable for any $w,\gamma$, The spectral norm of $\nabla_{\nu}\psi(w,\gamma,\nu)$ is uniformly bounded by $\sigma$. That is, $\|\nabla_{\nu}\psi(w,\gamma,\nu)\|_{op}\leq\sigma$ for all $w,\gamma,\nu$. the Hessian $\nabla_{\nu\nu}\psi(z,\gamma,\nu)$ has the largest eigenvalue bounded by some constant $\lambda_{\max}$ uniformly for all $\gamma$ and $\nu$.

These regularity conditions above are standard in the literature. See Chernozhukov et al. (2018) for instance.

To establish the consistency and asymptotic normality of the finite-dimensional parameters when the nuisance parameters satisfy Assumption A.3, we employ Neyman orthogonality, which ensures the moment condition is insensitive to perturbations in the nuisance parameters. Formally, this condition is expressed as

[TABLE]

for all $\xi$ in an appropriate space and any real value $s$. The following lemma establishes the Neyman orthogonality for the key moment function of our proposed estimator under the specification considered in this paper. While Neyman orthogonality may not immediately hold for more general specifications of latent distributions, the score function can be reformulated to satisfy this condition, as shown in Morzywolek et al. (2023) and Battey and Reid (2024) among others.

In the lemma below, we consider the specific moment function, given as

[TABLE]

where $\gamma$ represents all finite-dimensional parameters including causal effect parameter $\theta$ as well as the latent variable model. $q_{\gamma}(\cdot)$ and $r_{\gamma}(\cdot)$ are model-specific real-valued functions parametrized by $\gamma$. $h(\cdot)$ and $m(\cdot)$ represent nuisance functions in the double machine learning procedure, i.e., $h,m\in\nu$. In our latent DML scenario in Eq. (11), we set $q_{\gamma}(W)=\mathbb{E}[{f_{\hskip-0.56905ptu\hskip-0.56905pt}}(Z)\,|\,W,\gamma]$ and $r_{\gamma}(W)=0$.

The moment function $\phi(\cdot)$ is assumed to satisfy Assumption A.4. Additionally, the conditional mean-zero restrictions hold:

[TABLE]

Lemma A.5** (Neyman Orthogonality).**

Suppose that Assumption A.4 and the condition in (19)-(20) hold. Then, we have

[TABLE]

for all $\xi$ in an appropriate space with a real value $s$.

We now formally establish the consistency of our estimator under the regularity conditions in the following theorem.

Theorem A.6** (Consistency).**

Under Assumptions A.1-A.4 and the Neyman orthogonality condition in (18), we have, as $n\to\infty$,

[TABLE]

We next turn to the asymptotic distribution of our estimator and establish its limiting behavior. The following theorem shows that our estimator converges to a normal distribution and provides its asymptotic variance.

Theorem A.7** (Asymptotic Normality).**

Under Assumptions A.1-A.4 and the Neyman orthogonality condition in (18), we have, as $n\to\infty$,

[TABLE]

where \Omega:=HVar\big{(}\psi(W,\gamma_{0},\nu_{0})\big{)}H^{\top} with $H:=(\mathbb{E}[\nabla_{\gamma}\psi(W,\gamma_{0},\nu_{0})])^{-1}.$

The asymptotic variance derived in Theorem A.7 demonstrates the potential benefit of controlling for unobserved variables even in the absence of confounding effects. This advantage stems fundamentally from conditioning on additional information, which reduces variability in estimating equations. More precisely, by the law of total variance, the conditioning process necessarily decreases the overall variance of our estimating equations, which translates directly to smaller asymptotic variance and consequently narrower confidence intervals. This efficiency gain occurs despite the absence of endogeneity issue, underscoring that control variables can enhance estimation precision beyond their conventional role in bias reduction.

Usage of generative AI

We used generative AI tools at specific stages of this work. During manuscript preparation, OpenAI ChatGPT (GPT-4) was employed to revise phrasing and improve clarity. GitHub Copilot was used as an assistive tool during code development to suggest code snippets. These tools were not used to generate scientific content, design algorithms, or create experimental results. All core ideas, model development, implementation, and analyses were conducted and validated by the authors.

Appendix B Extended Appendix

We provide additional information to complement the main paper, including experimental settings and theoretical analysis.

B.1. Additional experimental settings

We describe additional experimental settings in Section B.1.1 for the synthetic experiments presented in Section 5.1 of the main paper. Details on the compute resources used are provided in Section B.1.2. We also include the full list of covariates used in the real data analysis of advertisement timing (Table 3).

B.1.1. Synthetic data

Our experiments utilized a synthetic data generation process that closely simulates real-world complexities. We modeled nuisance components $m$ and $g$ as sparse linear models with $100$ covariates, using coefficients drawn from a normal distribution to reflect real-world data sparsity. The noise distributions basically follow Eqs. (13) and (16). The data generator was configured to create datasets with specific attributes such as sample size $300$, number of covariates $100$, true causal effect $1.0$. We used this setup to conduct two sets of experiments: one verifying the outcome latent DML and another focusing on the confounder latent DML method.

For the first experiment, the noise $U$ in the outcome includes an additional noise generated by the exponential distribution with a mean of $5.0$. For the second, we examined cases with both positive and negative correlations between confounders and outcomes, where the confounders are discrete and generated as Bernoulli-distributed variables with a randomly specified success rate $q$. Specifically, in the positive correlation case, the parameters $a$ and $b$ in Eq. (16) were both set to $2.0$, while in the negative correlation case, they were set to $2.0$ and $-2.0$, respectively. In both experiments, the standard deviations of the outcome and treatment noises were set to $1.0$ and $0.5$, respectively. This comprehensive approach allowed for an in-depth evaluation of our latent DML framework, demonstrating its capacity to capture causal effects accurately in a range of simulated environments.

B.1.2. Compute Resources

All experiments were conducted on a standard local machine with 2.3GHz Quad-Core Intel Core i7 CPU and 32 GB Memory. The computational time for each experimental run varied between 2 to 5 hours depending on the complexity of the experiment. Given the relatively modest computational requirements, all experiments were feasible to run on this local machine without the need for specialized hardware or cloud resources.

B.2. Theoretical analysis details

In Section B.2.1, we establish identification conditions for our latent variable models. Section B.2.2 provides the proofs of the theoretical results stated in Section A.2 of the main paper.

B.2.1. Identification of latent variable models

We discuss the conditions necessary for the identification of latent variable models within latent DML. The population log likelihood function $L(\gamma)$is defined in Eq. (10) over the parameter space $\Gamma\subseteq\mathbb{R}^{d_{\gamma}}$, with $d_{\gamma}$ denoting the number of elements of $\gamma$. We denote by the derivative $\nabla_{\gamma}L(\gamma):=\partial L(\gamma)/\partial\gamma$ and let $\gamma_{0}\in\Gamma$ satisfy the first-order condition that $\nabla_{\gamma}L(\gamma_{0})=0$.

We assume the following conditions.

Assumption B.1.

1. (a)

$\Gamma$ is an open, convex set in $\mathbb{R}^{d_{\gamma}}$. 2. (b)

The map $\gamma\mapsto L(\gamma)$ is continuously differentiable over $\Gamma$. 3. (c)

The Fisher information $\mathcal{I}(\gamma):=\mathbb{E}[\nabla_{\gamma}L(\gamma)\nabla_{\gamma}L(\gamma)^{\top}]$ is a matrix with all finite elements and is a continuous function for every $\gamma\in\Gamma$.

The above assumptions are standard in the literature and commonly employed in theoretical analysis. Under Assumption B.1, the proposition below shows the identification condition for $\gamma_{0}$.

Proposition B.2** (Identifiability).**

Under Assumption B.1, $\gamma^{0}$ is locally identified if and only if $\mathcal{I}(\gamma^{\ast})$ is non-singular.

Proof.

Applying the standard identification result for the likelihood estimator, we can establish the desired result. See Theorem 1 of Rothenberg (1971). ∎

The verification of Assumption B.1 and the identifying restriction in the above proposition is straightforward for the models discussed in Section 3.2 and 3.3.

B.2.2. Proofs

Lemma A.5 (Neyman Orthogonality) * Suppose that Assumption A.4 and the condition in (19)-(20) hold. Then, we have*

[TABLE]

*for all $\xi$ in an appropriate space with a real value $s$. *

Proof.

Let $s\in\mathbb{R}$ be an arbitrary scalar value and set $\xi:=(\xi_{h},\xi_{m})^{\top}$ denote arbitrary deviations from $\nu$, so that $h(X)+s\xi_{h}(X)$ and $m(X)+s\xi_{m}(X)$ are well-defined. Under Assumption A.4(b), we can show that $\phi(W,\gamma,\nu+s\xi)$ is dominated by an integrable function. Using the dominated convergence theorem, we can show

[TABLE]

We rewrite the moment function as $\phi(W,\gamma,\nu+s\xi)=\phi_{1}(W,\gamma,\nu+s\xi)\cdot\phi_{2}(W,\gamma,\nu+s\xi),$ where

[TABLE]

Taking the derivative with respect to $s$, we have

[TABLE]

where $\phi_{1}^{\prime}(\cdot)$ and $\phi_{2}^{\prime}(\cdot)$ denote the partial derivatives with respect to $s$, specifically given by

[TABLE]

Also, it follows from (19) that

[TABLE]

Using the law of the iterated expectation, we can obtain the desired result. ∎

Theorem A.6 (Consistency) * Under Assumptions A.1-A.4 and the Neyman orthogonality condition in (18), we have, as $n\to\infty$,*

[TABLE]

Proof.

By applying a first-order Taylor expansion to the empirical moment condition, $n^{-1}\sum_{i=1}^{n}\psi(W_{i},\hat{\gamma},\hat{\nu})=0$, with respect to $\gamma$ around $\gamma_{0}$ with the mean value theorem, we obtain:

[TABLE]

where $\tilde{\gamma}$ represents a convex combination of $\gamma_{0}$ and $\hat{\gamma}$.

For the first term, we apply another first-order Taylor expansion with respect to $\nu$ around $\nu_{0}$:

[TABLE]

where $\tilde{\nu}$ represents a convex combination of $\nu_{0}$ and $\hat{\nu}$. Here, the first term converges to $\mathbb{E}[\psi(W,\gamma_{0},\nu_{0})]=0$ by the law of large numbers, while the second term converges to zero in probability by the Neyman orthogonality condition. It follows from (21) that

[TABLE]

By the regularity conditions, we also have uniform convergence:

[TABLE]

Combined with the identification condition that $\mathbb{E}[\psi(W,\gamma,\nu_{0})]=0$ if and only if $\gamma=\gamma_{0}$, and the non-singularity of $\mathbb{E}[\nabla_{\gamma}\psi(W,\gamma_{0},\nu_{0})]$, the result that \mathbb{E}\big{[}\nabla_{\gamma}\psi(W,\gamma_{0},\nu_{0})\big{]}(\hat{\gamma}-\gamma_{0})=o_{p}(1). leads to the desired conclusion: $\hat{\gamma}\to^{p}\gamma_{0}$. ∎

Theorem A.5 (Asymptotic Normality) * Under Assumptions A.1-A.4 and the Neyman orthogonality condition in (18), we have, as $n\to\infty$,*

[TABLE]

*where \Omega:=HVar\big{(}\psi(W,\gamma_{0},\nu_{0})\big{)}H^{\top} with $H:=(\mathbb{E}[\nabla_{\gamma}\psi(W,\gamma_{0},\nu_{0})])^{-1}.$ *

Proof.

By applying a first-order Taylor expansion to the empirical moment condition with respect to $\nu$ around $\nu_{0}$, coupled with the mean value theorem, we obtain:

[TABLE]

where $\tilde{\gamma}$ represents some convex combination of $\gamma_{0}$ and $\hat{\gamma}$.

We shall first show that the first term converges in probability to some positive-definite matrix and then establish that the second term converges in distribution to a multivariate normal distribution. The desired result follows from Slutsky’s theorem.

Under the regularity conditions on the moments, we obtain a uniform law of large numbers for $(1/n)\sum_{i=1}^{n}\nabla_{\gamma}\psi(W_{i},\gamma,\hat{\nu})$, namely:

[TABLE]

The consistency of $\hat{\gamma}$ implies the consistency of $\tilde{\gamma}$, that is, $\tilde{\gamma}\rightarrow^{p}\gamma_{0}$. Also, $\hat{\nu}$ is a consistent estimator, or $\hat{\nu}\to^{p}\nu_{0}$. Combining these two properties, we obtain that conditional on the data set:

[TABLE]

Because the matrix $\mathbb{E}\left[\nabla_{\nu}\psi(W,\gamma_{0},\nu_{0})\right]$ is non-singular, the continuity of the inverse operator yields:

[TABLE]

Next, we shall show the asymptotic normality. We apply a second-order Taylor expansion around $\nu_{0}$:

[TABLE]

where $\tilde{\nu}$ represents some convex combination of $\nu$ and $\nu_{0}$ by the mean value theorem. In the left hand side of the above equation, observe that the first term represents a normalized sum of $n$ independent and identically distributed random variables, scaled by $\sqrt{n}$. By the Central Limit Theorem, this term converges in distribution to $N(0,V)$, where $V:=Var(\psi(W_{i},\gamma_{0},\nu_{0}))$. It remains to show that both the second and third terms converge in probability to zero.

We consider the second term. It follows from the Neyman orthogonality condition that the conditional mean given the auxiliary information is given by:

[TABLE]

This implies that the mean of the second term is zero by the law of iterated expectations and also that the expectation of the cross-product of elements in the second term is equal to zero under the iid assumption. Thus, the conditional covariance of the second term, given the auxiliary dataset, is bounded above by:

[TABLE]

The consistency of $\hat{\nu}$ ensures that the latter term converges to zero as the sample size goes to infinity.

Since the mean is zero and the variance converges to zero, Chebyshev’s inequality shows that any deviation from the mean must become increasingly rare. This inequality provides a bound on the probability of deviations that is proportional to the variance, which we have shown approaches zero. Therefore, the probability of observing any fixed deviation from zero becomes arbitrarily small as the sample size grows large, establishing convergence in probability to zero.

Now, we consider the third term. Since $\nabla_{\nu\nu}\psi(W,\gamma,\nu)$ has a largest eigenvalue uniformly bounded by $\lambda_{\max}$, the third term is bounded by $(\lambda_{\max}/2)n^{-1/2}\sum_{i=1}^{n}\|\hat{\nu}(X_{i})-\nu_{0}(X_{i})\|^{2}.$ For a fixed auxiliary data set, $n^{-1}\sum_{i=1}^{n}\|\hat{\nu}(X_{i})-\nu_{0}(X_{i})\|^{2}$ converges to $\mathbb{E}[\|\hat{\nu}(X_{i})-\nu_{0}(X_{i})\|^{2}]$. By the $n^{1/4}$-consistency of the estimator of nuisance functions and its regularity, the third term converges to zero in probability. Combining all the results established thus far, we obtain the desired conclusion. ∎