The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
46,309 characters
Causal Inference in High-Dimensional Generalized Linear Models with Binary Outcomes
\begin{abstract}
This paper proposes a debiased estimator for causal effects in high-dimensional generalized linear models with binary outcomes and general link functions. The estimator augments a regularized regression plug-in with weights computed from a convex optimization problem that approximately balances link-derivative–weighted covariates and controls variance; it does not rely on estimated propensity scores. Under standard conditions, the estimator is $\sqrt{n}$-consistent and asymptotically normal for dense linear contrasts and causal parameters. Simulation results show the superior performance of our approach in comparison to alternatives such as inverse propensity score estimators and double machine learning estimators in finite samples. In an application to the National Supported Work training data, our estimates and confidence intervals are close to the experimental benchmark.
\keywords{Causal Parameters, High-dimensional Models, Link Functions, Weighting}
\end{abstract}
\section{Introduction}\label{sec:intro}
\setcounter{equation}{0}
\setcounter{theorem}{0}
\setcounter{assumption}{0}
\setcounter{proposition}{0}
\setcounter{corollary}{0}
\setcounter{lemma}{0}
\setcounter{example}{0}
\setcounter{remark}{0}
Causal inference from observational data is central across economics, epidemiology, and related fields. Binary outcomes are common, and generalized linear models (GLMs) provide a convenient framework. This paper proposes a method to estimate causal parameters in high-dimensional GLMs with binary outcomes.
\subsection{Motivation}
In observational studies, an unconfoundedness assumption—treatment assignment is \emph{as good as random} conditional on covariates—is often invoked. To make this plausible, researchers include many covariates or flexible functions of them. When the number of covariates is large relative to the sample size, regularization method such as Lasso by \citet{tibshirani1996regression} is used, but shrinkage introduces bias and the sampling distribution of the original Lasso estimator is not generally tractable for inference.
For illustration, consider a single-covariate two-stage model:
\[Y_i \sim \text{Bernoulli}(g(\beta_{D}D_i+\beta_{X}X_i)), \quad D_i \sim \text{Bernoulli}(g^*(X_i)),\]
where $X_i \sim N(0,1)$ independently drawn across $i \in \{1,...,n\}$ for $n = 100$, and parameters $\beta_{D} = \beta_{X} = 1$, $g(x) = g^{*}(x) = \exp(x)/(1+\exp(x))$.
The respective regression plug-in estimates for a causal parameter: average treatment effect on the treated (ATET), using standard maximum likelihood estimates (MLE) and penalized MLE with Lasso penalty, are:
\[\hat{\tau}_{MLE} = \frac{1}{n_t}\sum_{\{i:D_i = 1\}} g(\hat{\beta}_{D,MLE}+\hat{\beta}_{X, MLE}X_i)-g(\hat{\beta}_{X, MLE}X_i),\]
\[\hat{\tau}_{Lasso} = \frac{1}{n_t}\sum_{\{i:D_i = 1\}} g(\hat{\beta}_{D,Lasso}+\hat{\beta}_{X, Lasso}X_i)-g(\hat{\beta}_{X, Lasso}X_i).\]
Figure \ref{fig:introduction} displays histograms of the estimated ATET using MLE and penalized MLE with Lasso each based on 500 replications; the dashed black line marks the true value. The two panels show that the MLE estimator is centred close to the truth and approximates a normal distribution, whereas the Lasso-based estimator exhibits bias and noticeable deviations from normality. This pattern is consistent with the fact that the regression plug-in with MLE performs well when the outcome model is correctly specified.
\begin{figure}
\begin{center}
\begin{minipage}[b]{0.45\textwidth}
\includegraphics[scale=0.17]{Intro_MLE.eps}
\end{minipage}
\hspace{5mm}
\begin{minipage}[b]{0.45\textwidth}
\includegraphics[scale=0.17]{Intro_PMLE.eps}
\end{minipage}
\end{center}
\caption{Histograms of the estimated ATET using MLE and Penalized MLE (Lasso)}
\label{fig:introduction}
\end{figure}
This example underscores the need to mitigate bias for causal targets when regularization is necessary.
\subsection{Contributions and Related Work}
This paper proposes a debiased estimator for causal effects in high-dimensional generalized linear models (GLMs) with binary outcomes that does not rely on estimated propensity score. The estimator starts from a regularized GLM fit for the outcome model and then solves a single convex optimization to construct weights. The optimization program approximately balances covariate functions implied by the link derivative and simultaneously control the variance of the final estimator.
The key principle is a bias–variance trade-off. The GLM link function determines, through its derivative, which functions must be matched between treated and controls to remove the leading first-stage bias. The program therefore enforces approximate balance in these derivative-weighted covariates. At the same time, GLMs imply heteroskedastic residuals, so the objective includes an explicit penalty on the variance. In the special case of the identity link, corresponding to the linear model, the derivative is constant and the requirement reduces to equalizing treated and control covariate means, which coincides with standard covariate balance.
We establish two sets of asymptotic results: (i) for linear functionals $\xi^\top\beta_c$ (which may be dense), and (ii) for the causal parameter itself (a sample-average target rather than a fixed linear contrast). In both cases, the balance constraint makes the plug-in bias vanish faster than root-$n$, so the estimator is $\sqrt n$-consistent and asymptotically normal. We also provide a feasible variance estimator from the same optimization.
In simulations (main text and supplement), the estimator performs better than inverse probability weighting (IPW), double/debiased machine learning (DML), automatic debiased machine learning (AML), and linear approximate residual balancing (ARB). The improvements are largest when the propensity model is complex. Confidence interval coverage approaches the nominal level as $n$ grows. The full-sample implementation has the lowest mean-squared error and the shortest intervals among our variants. In an empirical re-analysis of the National Supported Work (NSW) dataset, our estimate is closest to the experimental benchmark, and its 95\% confidence interval almost coincides with the randomized-experiment interval.
This paper contributes to three strands. First, it relates to doubly robust and orthogonal estimation in high dimensions. Regularized outcome regression can induce first-order bias in ATEs, as shown by \citet{belloni2014high}. Doubly robust and orthogonal estimators combine outcome and propensity models to obtain $\sqrt n$ inference under first-stage rate conditions, as developed by \citet{farrell2015robust}, \citet{belloni2017program}, and \citet{chernozhukov2018double}, and are related to debiased Lasso, as in \citet{van2014asymptotically} and \citet{zhang2014confidence}. We debias without inverting an estimated propensity score by imposing GLM-specific balance and minimizing variance in a single optimization problem.
Second, it connects to balancing-based debiasing and minimax linear estimation. Exact balancing via inverse-probability tilting or augmented scores, as in \citet{graham2012inverse} and \citet{graham2016efficient}, relies on inverse propensities and can be unstable in high dimensions. The ARB estimator of \citet{athey2018approximate} addresses linear models by balancing residualized covariates to remove regularization bias. We extend this idea to GLMs, where the relevant moments are link-derivative-weighted and the variance is heteroskedastic. Our optimization has a minimax-MSE interpretation and is related to augmented minimax linear estimation in \citet{hirshberg2021augmented}. Compared with ARB, we target GLM moments, use a different program, and obtain a feasible variance estimator. Riesz-based AML estimators of \citet{chernozhukov2022automatic} and \citet{chernozhukov2024automatic} avoid explicit propensities by learning the Riesz representer; our approach instead solves directly for weights that achieve debiasing in the GLM metric and explicitly optimizes variance, improving finite-sample stability.
Finally, the paper relates to inference in high-dimensional GLMs and to post-selection inference. \citet{cai2023statistical} provide debiased inference for single coefficients under general links using a two-step weighted correction. We allow dense linear contrasts and supply a single weighting scheme that controls both bias and variance. \citet{belloni2016post} and \citet{belloni2018uniformly} develop post-selection and uniform inference for linear models, including dense contrasts. Their procedures estimate both the outcome and an auxiliary propensity-type regression (double selection); our debiasing focuses on the outcome regression and an optimization that enforces balance.
\subsection{Organization and Notations}
The remainder of the paper is organized as follows. Section 2 formulates the problem and describes the estimator and algorithm. Section 3 develops the asymptotic theory for dense linear contrasts and causal parameters in a high-dimensional framework, deriving convergence rates and limit distributions. Section 4 reports simulation results. Section 5 presents an empirical application. Section 6 concludes. Proofs and additional simulation results are provided in the supplementary materials.
Throughout, for a vector $\boldsymbol{b} = (b_1,...,b_n)^\top \in \mathbb{R}^n$, we define the $\ell_p$ norm by $||\boldsymbol{b}||_p = \left(\sum_{i=1}^n|b_i|^p\right)^{1/p}$, and the $\ell_\infty$ norm by $||\boldsymbol{b}||_\infty = \max_{1\leq i \leq n} |b_i|$. For a matrix $A \in \mathbb{R}^{p\times q}$, $\sigma_i$ stands for the $i$-th largest singular value of $A$ and $\sigma_{\max}(A) = \sigma_1(A)$, $\sigma_{\min}(A) = \sigma_{\min\{p,q\}}(A)$.
For a smooth function $f(x)$ defined on $\mathbb{R}$, we define $f'(x) = df(x)/dx$ and $f''(x) = d^2f(x)/dx^2$. We define $\phi(x)$ and $\Phi(x)$ as the density function and cdf of the standard Gaussian random variable, respectively. We denote $\stackrel{p}{\rightarrow}$, $\stackrel{d}{\rightarrow}$ as convergence in probability and in distribution respectively.
For positive sequences $\{a_n\}$ and $\{b_n\}$, we write $a_n = o(b_n)$, $a_n \ll b_n$ or $b_n \gg a_n$ if $\lim_n a_n/b_n = 0$, and write $a_n = \mathcal{O}(b_n)$, $a_n \lesssim b_n$ or $b_n \gtrsim a_n$ if there exists a constant $C$ such that $a_n \leq Cb_n$ for all $n$. We write $a_n \simeq b_n$ if $a_n \lesssim b_n$ and $a_n \gtrsim b_n$.
For random variable sequences $\{X_n\}$ and positive sequences $\{a_n\}$, we write $X_n = o_p(a_n)$ if $X_n/a_n \stackrel{p}{\rightarrow} 0$, and write $X_n = O_p(a_n)$ if $\limsup_n P[|X_n/a_n| \geq M] \to 0$.
\section{Causal inference under high-dimensional GLM}
\setcounter{equation}{0}
\setcounter{theorem}{0}
\setcounter{assumption}{0}
\setcounter{proposition}{0}
\setcounter{corollary}{0}
\setcounter{lemma}{0}
\setcounter{example}{0}
\setcounter{remark}{0}
\subsection{Problem Formulation}
We estimate the average treatment effect on the treated (ATET) under the potential outcomes framework.\footnote{The average treatment effect (ATE) and the average treatment effect on the control (ATEC) can be handled analogously; they require modelling the missing potential outcome for the corresponding group.} For units $i=1,\ldots,n$, let $Y_i(1),Y_i(0)\in\{0,1\}$ denote potential outcomes, $D_i\in\{0,1\}$ the treatment indicator, $Y_i^{\text{obs}}=Y_i(D_i)$ the observed outcome, and $X_i\in\mathbb{R}^p$ a (possibly high-dimensional) covariate vector. We assume unconfoundedness for the control outcome and a GLM for $Y_i(0)\mid X_i$.
\begin{assumption}
$D_i \perp\!\!\!\perp Y_i(0)|X_i$ for all $i \in \{1,...,n\}$.
\end{assumption}
\begin{assumption}
The response function conditional on covariates satisfies
\begin{equation}
Y_i(0)|X_i \sim \text{Bernoulli}(g(X_i^\top \beta_c)),
\end{equation}
where $g: \mathbb{R} \to (0,1)$ is a known link function and $\beta_c \in \mathbb{R}^p$. In other words, we have $Y_i(0) = g(X_i^\top \beta_c)+\varepsilon_i$ where $\mathbb{E}(\varepsilon_i|X_i)=0$ and $var(\varepsilon_i|X_i)=g(X_i^\top \beta_c)(1-g(X_i^\top \beta_c))$ for all $i \in \{1,...,n\}$.
\end{assumption}
\begin{assumption}
The general link function $g$ satisfies the following regularity conditions: (a) The link function $g$ is twice differentiable, monotonic increasing, Lipschitz on $\mathbb{R}$, and concave on $\mathbb{R}_{+}$; and for any $x\in \mathbb{R}$, it holds that $g(x)+g(-x)=1$; (b) There exist some constants $C_1$, $C_2>0$ such that, for all $x \geq 0$, $g(x) \leq \Phi(C_1x)$ where $\Phi(x)$ is the standard Gaussian cdf, and $\max\{g'(x)/(1-g(x)), x^2g'(x)\} < C_2$; (c) There exist some constants $c>0$ such that $\sup_{x\in \mathbb{R}}|g''(x)| \leq c$; (d) For $\ell_g(\beta)$ defined by
\begin{equation}
\ell_g(\beta) = -\frac{1}{n}\sum_{i=1}^n y_i\log\left[\frac{g(X_i^\top\beta)}{1-g(X_i^\top\beta)}\right]-\frac{1}{n}\sum_{i=1}^n\log(1-g(X_i^\top\beta)),
\label{eq:lg}
\end{equation}
there exists some constant $C>1$ such that the Hessian matrix $\ell''_g(\beta)$ can be expressed as $\ell''_g(\beta) = \sum_{i=1}^n h(\beta;y_i,X_i)X_iX_i^\top/n$ for some $h(\beta;y_i,X_i)>0$ satisfying
\begin{equation*}
\max_{1\leq i\leq n}\left|\log h(\beta+b;y_i,X_i) - \log h(\beta;y_i,X_i)\right| \leq C(|X_i^\top \beta|^2+|X_i^\top b|^2+|X_i^\top b|).
\end{equation*}
\end{assumption}
\textit{Assumption 2.3} allows many links used in practice, including but not limited to canonical links. The supplement verifies these conditions for the examples below.
\begin{example}
\textnormal{The Logistic regression model prescribes the following link function
\[g(x) = \frac{\exp(x)}{1+\exp(x)},\]
which has a characteristic, $g'(x) = g(x)(1-g(x)) = \exp(x)/(1+\exp(x))^2$.}
\end{example}
\begin{example}
\textnormal{Consider a latent variable model:
\[y_i^* = X_i^\top \beta + \epsilon_i, \quad y_i = \boldsymbol{I}(y_i^* \geq 0),\]
where $g(\cdot)$ is the cdf of $-\epsilon_i$. Beside the logistic link function, examples include:
\setcounter{bean}{0}
\begin{list}
{(\alph{bean})}{\usecounter{bean}}
\item Probit link function: $g(x) = \Phi (x)$ where $\Phi(x)$ denoting the standard Gaussian cdf and $g'(x) = \phi(x)$ being the standard Gaussian pdf.
\item Student's $t_v$-distributions with $\nu \in \mathbb{N} $: $g'(x) = \frac{\Gamma_d(\frac{\nu+1}{2})}{\sqrt{\nu \pi}\Gamma_d(\frac{\nu}{2})}\left(1+\frac{x^2}{\nu}\right)^{-(\nu+1)/2}$ where $\nu$ is the degree of freedom and $\Gamma_d$ is the gamma function.
\end{list}}
\end{example}
\begin{remark}
\textnormal{We treat the link as known (\textit{Assumption 2.2}). When it is unknown, \citet{liu2021singleindex} propose an average-derivative estimator for fixed linear combinations of coefficients. Their estimator could in principle be combined with our balancing step; we leave this extension to future work.}
\end{remark}
Our target parameter is
\begin{equation}
\tau := \frac{1}{n_t} \sum_{\{i:D_i = 1\}} \mathbb{E}[Y_i(1) - Y_i(0)|X_i],
\end{equation}
where $n_t$ is the number of treated units. Under \textit{Assumption 2.2},
\begin{equation}
\tau = \mu_t - \mu_c,\quad \quad \mu_c = \frac{1}{n_t}\sum_{\{i:D_i = 1\}} g(X_i^\top \beta_c)
\end{equation}
We can estimate $\mu_t$ unbiasedly by a simple average of $\{Y_i^{\text{obs}}, D_i = 1\}$. The main challenge is high-dimensional estimation of $\mu_c$.
\subsection{Debiased Estimation Procedure}
We estimate $\mu_c$ by combining a regularized GLM plug-in with a residual reweighting step. The plug-in can suffer from regularization bias, especially when treated and control covariates differ and extrapolation is required. We therefore reweight control residuals using a convex program that (i) balances link-derivative–weighted covariates between treated and controls and (ii) controls variance.
\begin{algorithm}
\textnormal{
Input: design matrix $\boldsymbol{X} \in \mathbb{R}^{n \times p}$, treatment indicator $\boldsymbol{D} \in \{0,1\}^n$, outcomes $\boldsymbol{Y} \in \{0,1\}^n$.}
\setcounter{bean}{0}
\begin{center}
\begin{list}
{\textnormal{\textsc{Step} \arabic{bean}.}}{\usecounter{bean}}
\item \textnormal{Randomly split the control sample into two equal parts, $\mathcal{S}_1$ and $\mathcal{S}_2$, each containing $n_c$ observations.}
\item \textnormal{Fit a penalized GLM on $\mathcal{S}_1$:
\begin{equation}
\hat{\beta}_{c1} = \operatorname{argmin}_{\beta} \{\ell_g(\beta)+\lambda_n ||\beta||_1\}
\end{equation}
where $\ell_g(\beta)$ is defined by \eqref{eq:lg}.}\footnote{Throughout this paper, we use the lasso method to compute the plug-in regression estimate \(\hat{\beta}_{c1}\). However, other methods that offer similar \(L_1\) estimation error bounds can also be used.}
\item \textnormal{Compute weights $\gamma_i$ by solving
\normalsize{
\begin{equation*}
\begin{aligned}
\min_{\boldsymbol{\gamma}} &\left[\ \ \quad \quad (1-\zeta)\sum_{\{i \in \mathcal{S}_2\}} \gamma_i^2 g(X_i^\top \hat{\beta}_{c1})(1 - g(X_i^\top \hat{\beta}_{c1})) +\right. \\
&\left. \quad \zeta \left\Vert \frac{1}{n_t}\sum_{\{i:D_i=1\}} g'(X_i^\top\hat \beta_{c1})X_i^\top - \sum_{\{i \in \mathcal{S}_2\}} \gamma_i g'(X_i^\top\hat \beta_{c1})X_i^\top \right \Vert_\infty^2 \quad \right]
\end{aligned}
\end{equation*}
}
\normalsize{
\begin{equation}
s.t. \ \ \sum_{i \in \mathcal{S}_2} \gamma_i = 1, \quad 0\leq \gamma_i \leq \log(n_c)/n_c.
\label{eq:optm_initial}
\end{equation}
}}
\item \textnormal{Form
\begin{equation}
\hat{\mu}_c = \frac{1}{n_t} \sum_{\{i:D_i = 1\}}g(X_i^\top \hat{\beta}_{c1}) + \sum_{i \in \mathcal{S}_2} \gamma_i(Y_i - g(X_i^\top \hat{\beta}_{c1})).
\label{eq:mu_chat}
\end{equation}
and estimate ATET by:
\begin{equation}
\hat{\tau} = \frac{1}{n_t}\sum_{\{i:D_i = 1\}} Y_i^{\text{obs}} - \hat{\mu}_c
\label{eq:tau_hat}
\end{equation}}
\end{list}
\end{center}
\end{algorithm}
Compared to the approach by \citet{athey2018approximate},
the key differences are: our objective accounts for GLM heteroskedasticity when penalizing variance, and the balance conditions use link-derivative–weighted covariates rather than residualized linear covariates. Relative to \citet{cai2023statistical}, who target single coefficients, our program targets $\mu_c$ which depends on all entries of $\beta_c$ through $g(X^\top\beta_c)$ and the treated covariates.
The optimization function in \eqref{eq:optm_initial} has a minimax interpretation. The objective combines a variance term and a squared worst-case design-conditional bias term; $\zeta\in[0,1]$ governs the trade-off. At $\zeta=1/2$ the optimization problem chooses weights that minimize a maximal mean squared error (MSE). This links our estimator to the augmented minimax linear estimator of \citet{hirshberg2021augmented}. Their analysis relies on Donsker classes; our high-dimensional setting does not, and we instead use standard assumptions for regularized GLMs.
\begin{remark}
\textnormal{Note that the sample-splitting process is primarily used to simplify theoretical analysis and does not pose a limitation for the applications.
In simulations of Section 4, cross-fitting and full-sample implementations perform better by using all observations.
This mirrors the strategy used by \citet{cai2023statistical}, who, while demonstrating the favourable properties of their estimator through sample-splitting, also effectively showed its practical efficacy using full sample in simulations.}
\end{remark}
\begin{remark}
\textnormal{While our primary focus is on estimating the ATET, analogous estimators can be constructed, such as the ATE \(\sum_{i=1}^n \mathbb{E}[Y_i(1)-Y_i(0)|X_i]/n\),
and the ATEC \(\sum_{\{i:D_i=0\}} \mathbb{E}[Y_i(1)-Y_i(0)|X_i]/(n-n_t)\), by changing the target sample over which the missing potential outcomes are imputed.}
\end{remark}
\subsection{Comparison with Different Methods}
In this subsection, we focus on comparing different approaches for estimating $\mu_c$.\footnote{For simplicity, the intuitive comparison is based on no-splitting implementations for all estimators (direct plug-in and plug-in with different weightings).} First, let us consider a direct plug-in estimator, a plug-in with equal residual weights, and our plug-in plus optimized weights in~\eqref{eq:mu_chat}. The direct plug-in estimator is
\begin{equation}
\hat{\mu}_{c,plug} = \frac{1}{n_t} \sum_{\{i:D_i=1\}} g(X_i^\top \hat \beta_c)
\label{eq:mu_chat_plugin}
\end{equation}
and the plug-in estimator with simple averaging of residuals:
\begin{equation}
\hat{\mu}_{c,reg} = \frac{1}{n_t} \sum_{\{i:D_i=1\}} g(X_i^\top \hat \beta_c) + \frac{1}{n_c} \sum_{\{i: D_i=0\}} (Y_i^{obs}-g(X_i^\top \hat \beta_c))
\label{eq:mu_chat_reg}
\end{equation}
where $\hat{\beta}_c$ is obtained by Lasso.
Examine error analyses by contrasting the discrepancies among estimators $\hat{\mu}_{c,plug}$, $\hat{\mu}_{c,reg}$, $\hat{\mu}_{c}$ and the true value $\mu_c$:
\begin{equation}
\hat \mu_{c,plug} -\mu_c = \left[\frac{1}{n_t}\sum_{\{i:D_i=1\}} g'(X_i^\top\hat \beta_c)X_i^\top \right](\hat{\beta}_c - \beta_c)
+ \frac{1}{n_t}\sum_{\{i:D_i=1\}} \Delta_i^t,
\label{eq:mu_chat_plugin_error}
\end{equation}
\begin{equation}
\begin{aligned}
\hat \mu_{c,reg} -\mu_c &= \left[\frac{1}{n_t}\sum_{\{i:D_i=1\}} g'(X_i^\top\hat \beta_c)X_i^\top - \frac{1}{n_c}\sum_{\{i:D_i=0\}} g'(X_i^\top\hat \beta_c)X_i^\top \right](\hat{\beta}_c - \beta_c) \\
&+ \frac{1}{n_c} \boldsymbol{1}_{n_c}^\top\boldsymbol{\varepsilon} + \frac{1}{n_t}\sum_{\{i:D_i=1\}} \Delta_i^t + \frac{1}{n_c}\sum_{\{i:D_i=0\}} \Delta_i^c,
\label{eq:mu_chat_reg_error}
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\hat \mu_{c} -\mu_c &= \left[\frac{1}{n_t}\sum_{\{i:D_i=1\}} g'(X_i^\top\hat \beta_c)X_i^\top - \sum_{\{i:D_i=0\}} \gamma_i g'(X_i^\top\hat \beta_c)X_i^\top \right](\hat{\beta}_c - \beta_c) \\
&+ \boldsymbol{\gamma}^\top\boldsymbol{\varepsilon} + \frac{1}{n_t}\sum_{\{i:D_i=1\}} \Delta_i^t + \sum_{\{i:D_i=0\}} \gamma_i \Delta_i^c.
\label{eq:mu_chat_error}
\end{aligned}
\end{equation}
where $\Delta_i^{(\cdot)} = g''(X_i^\top\hat \beta_c + t_{i}^{(\cdot)}X_i^\top(\beta-\hat \beta_c))[X_i^\top(\hat \beta_c-\beta)]^2$ for some $t_{i}^{(\cdot)} \in (0,1)$.
In low dimensions, the direct plug-in estimator with MLE is $\sqrt{n_c}$-consistent and asymptotically normal. In high dimensions, plugging in regularized estimates need not retain these properties due to shrinkage bias. When treated and control groups are (approximately) balanced, the regression–imputation estimator in \eqref{eq:mu_chat_reg_error} improves on the plug-in \eqref{eq:mu_chat_plugin_error} by averaging residuals and can be valid under additional conditions. However, with noticeable imbalance, the simple residual average does not control the first linear term in \eqref{eq:mu_chat_reg_error}. In contrast, our estimator \eqref{eq:mu_chat_error} are flexible to be able to chooses weights to balance the derivative-weighted covariates and let the first term to be as small as possible.
By applying H\"older's inequality, the absolute value of \eqref{eq:mu_chat_error} can be bounded as follows:
\begin{equation}
\begin{aligned}
|\hat \mu_{c} -\mu_c| & \leq
\underbrace{\left\Vert \frac{1}{n_t}\sum_{\{i:D_i=1\}} g'(X_i^\top\hat \beta_c)X_i^\top - \sum_{\{i:D_i=0\}} \gamma_i g'(X_i^\top\hat \beta_c)X_i^\top \right \Vert_\infty ||\hat{\beta}_c - \beta_c||_1}_{\text{(a) Main term of bias}} \\
&+ \underbrace{\left|\sum_{\{i:D_i=0\}} \gamma_i \varepsilon_i \right|}_{\text{(b) Main term of variance}} + \underbrace{\left|\frac{1}{n_t}\sum_{\{i:D_i=1\}} \Delta_i^t \right| + \left|\sum_{\{i:D_i=0\}} \gamma_i \Delta_i^c\right|}_{\text{(c) Negligible terms}}.
\label{eq:mu_chat_error_byH}
\end{aligned}
\end{equation}
In Equation \eqref{eq:mu_chat_error_byH}, term (a) is the primary bias, which is our main concern, term (b) denotes the variance, anticipated to be of order $O_p(1/\sqrt{n}_c)$, and term (c) consists of negligible components.
The bias term (a) equals the imbalance in derivative-weighted covariates times the $L_1$ error of the first stage estimates.
Under standard conditions, this error scales as $O_p(k\sqrt{\log(p)/n_c})$. Therefore, our optimization should drive the imbalance to $O_p(\sqrt{\log(p)/n_c})$.
to ensure that the entire bias term (a) is negligible relative to the variance term (b) when $k \ll \sqrt{n_c}/\log p$.
This key insight is the main foundation of our proofs of theorems presented in Section 3.
We next consider benchmark our estimator against methods that use inverse propensity weighting. For example,
a pure weighting method $\hat{\mu}_{c,weight} = \sum_{\{i:D_i =0\}} \gamma_i Y_i^{\text{obs}}$ with
\begin{equation}
\gamma_i = \frac{\frac{p(X_i)}{1-p(X_i)}}{\sum_{\{i:D_i = 0\}}\frac{p(X_i)}{1-p(X_i)}}.
\label{eq:mu_chat_weig}
\end{equation}
where $p(x) := \mathbb{P}(D=1|X=x)$ is the propensity score.
With nonparametric $p(\cdot)$ and a fixed number of covariates, this estimator is asymptotically linear and attains the semiparametric efficiency bound, as shown by \citet{hirano2003efficient}.
We also considers the Double Machine Learning (DML) estimator of \citet{chernozhukov2018double}, which combines regression imputation with inverse propensity weighting,
\begin{equation}
\hat{\mu}_{c,DL} = \frac{1}{n_t} \sum_{\{i:D_i = 1\}} g(X_i^\top \hat{\beta}_c) + \sum_{\{i:D_i = 0\}} \frac{\frac{\hat{p}(X_i)}{1-\hat{p}(X_i)}(Y_i^{\text{obs}}-g(X_i^\top \hat{\beta}_c))}{\sum_{\{i:D_i=0\}}\frac{\hat{p}(X_i)}{1-\hat{p}(X_i)}}
\label{eq:mu_chat_DL}
\end{equation}
With sample-splitting or cross-fitting, DML exhibits desirable asymptotic performance, as long as the estimations of both the outcome and propensity score model are moderately well.
However, in high dimensions and finite sample, inverse weights can be unstable. Specifically, a minor perturbation in the estimated value of $1-p(X_i)$ can translate into a significant effect, especially when it is close to 0, due to the involvement of its inverse in the weighting. This challenge is indicative of a broader issue that often emerges when inverting estimated values.
Our estimator differs in how the weights are formed: we strategically relies on the GLM structure to choose weights through a single convex program that balances link-derivative–weighted covariates and penalizes variance. This avoids the instability associated with inverse weights. Conversely, while our method is confined to the GLM setting, inverse propensity score weights debiasing accommodates any general form of outcome models, highlighting a flexibility not inherent to our approach.
Finally, we consider AML, as developed by \citet{chernozhukov2022automatic}. It avoids inverse propensity weights by estimating a Riesz representer $\alpha(\cdot)$ and then forming
\begin{equation}
\hat{\mu}_{c,AL} = \frac{1}{n_t} \sum_{\{i:D_i = 1\}} g(X_i^\top \hat{\beta}_c) + \sum_{\{i:D_i = 0\}} \hat{\alpha}(X_i)(Y_i^{\text{obs}}-g(X_i^\top \hat{\beta}_c))
\label{eq:mu_chat_AML}
\end{equation}
where \(\hat{\alpha}(\cdot)\) can be learned by a minimum-distance program (e.g., Lasso).
In contrast to our method, which controls both bias and variance, the AML approach mainly focuses on bias reduction, even under heteroskedasticity case, through its goal of Riesz representer estimation.
\section{Theoretical Properties}
\setcounter{equation}{0}
\setcounter{theorem}{0}
\setcounter{assumption}{0}
\setcounter{proposition}{0}
\setcounter{corollary}{0}
\setcounter{lemma}{0}
\setcounter{example}{0}
\setcounter{remark}{0}
\subsection{Debiased Estimation for Dense Contrasts}
In the preceding section, the focus was on the causal parameter estimation. To formally establish its asymptotic properties, it is convenient to begin with a linear contrast $\xi^\top \beta_c$ where $\xi$ is a $p \times 1$ vector.
For instance, if $\xi$ is a basis vector like $(1,0,...,0)$, i.e., our target of interest is the first element of $\beta_c$, the problem becomes akin to the widely-studied debiased lasso estimator such as \citet{van2014asymptotically} and \citet{javanmard2014confidence}.
Here we allow $\xi$ to be dense and consider $\theta = \xi^\top \beta_c$. Linked to Algorithm 1, the estimator is given by
\begin{equation}
\hat{\theta} = \xi^\top \hat{\beta}_{c1} + \sum_{i \in \mathcal{S}_2} \gamma_i(Y_i^{obs} - g(X_i^\top \hat{\beta}_{c1})),
\label{eq:est_dense}
\end{equation}
where $\boldsymbol{\gamma}$ solves the constrained program
\begin{equation*}
\min_{\boldsymbol{\gamma}} \sum_{\{i\in \mathcal{S}_2\}} \gamma_i^2 g(X_i^\top \hat \beta_{c1})(1-g(X_i^\top \hat \beta_{c1}))
\end{equation*}
\begin{equation}
\begin{aligned}
s.t.\quad \quad \left\Vert \xi - \sum_{\{i \in \mathcal{S}_2\}} \gamma_i g'(X_i^\top\hat \beta_{c1})X_i \right \Vert_\infty & \leq C_3 \sqrt{\log p/n}, \\
\max_{i \in \mathcal{S}_2} |\gamma_i| & \leq C_4 \frac{\sqrt{\log n_c}}{n_c}.
\label{eq:optm_dense}
\end{aligned}
\end{equation}
Note that the optimization problem here is represented in a constraint form rather than the Lagrange form in Algorithm 1, facilitating theoretical proofs.
Intuitively, the goal of this problem is to find the weights that minimize the variance, while simultaneously ensuring that the bias remains constrained to a low level.
This formulation bridges important cases. For $\xi=e_1$ it aligns with \citet{cai2023statistical}. For $\xi=\frac{1}{n_t}\sum_{D_i=1} g'(X_i^\top\hat\beta_{c1})X_i$, $\hat\theta$ is equivalent to $\hat\mu_c$.
It should be emphasized that applying the method of \citet{cai2023statistical} directly within our framework is not suitable.
Specifically, Theorem 1 in their work indicates that their debiased estimator \(\tilde{\beta}\) satisfies \(\sqrt{n_c}(\tilde{\beta}_{c,i} - \beta_{c,i}) = A_{n,i} + B_{n,i}\), where \(A_{n_i}\) is asymptotically normally distributed and \(B_{n_i}\) is negligible per coordinate. But the aggregated error term \(\sum_{i=1}^p \xi_i B_{n,i}\) need not be negligible, as it could increase polynomially with the dimension of covariates $p$.
In addition, the simultaneous inference approach suggested by Section 2.3 of \citet{cai2023statistical} is conservative for our purposes, given our interest in a linear combination of coefficients rather than a subset of coefficients.
Therefore, our approach extends \citet{cai2023statistical}, accommodating various \(\xi\) values and thereby generalizing from the debiasing of individual coefficients to encompass any combination of coefficients.
We impose two additional assumptions:
\begin{assumption}
$\{X_i\}_{1\leq i\leq 2n_c}$ satisfies $X_i = Q_i \Sigma^{1/2}$, where $\mathbb{E}(Q_{ij}) = 0$, var$(Q_{ij})=1$, for all $i$ and $j$, and all individual entries $Q_{ij}$ are sub-Gaussian and independent.
\end{assumption}
\begin{assumption}
Let that $\Sigma = \mathbb{E}[X_iX_i^\top] \in \mathbb{R}^{p\times p}$ and
\[\Theta(k) = \{(\beta_c,\Sigma):||\beta||_0 \leq k, ||\beta||_2 \leq M_1, M_2^{-1} \leq \lambda_{\min}(\Sigma)\leq \lambda_{\max}(\Sigma) \leq M_2\}\]
for constants $M_1>0$ and $M_2>1$ independent of $n_c$ and $p$. Suppose $(\beta_c, \Sigma) \in \Theta(k)$.
\end{assumption}
While \textit{Assumption 3.1} may impose restrictions in practical applications, it is an important technical condition to ensure the feasibility of the minimization problem.
In causal estimation, this assumption can be relaxed under the overlap assumption.
\textit{Assumption 3.2}, on the other hand, is conventional within the literature under the exact sparsity and GLM framework.\footnote{The exact sparsity can be weakened to approximate sparsity, see \citet{belloni2016post}, as long as the conditions on the estimation error bound hold, since our theoretical results depend only on the error bound on coefficients estimation. For simplicity, we restrict our attention to the exact sparsity case for theoretical properties. However, in our simulation study (Section 4), we include a decaying-coefficient design, illustrating performance when the approximate sparsity holds.}
The positive definiteness of the covariance matrix $\Sigma$ ensures the restricted strong convexity of the loss function, playing a pivotal role in delineating the Lasso risk bounds, such as in \citet{negahban2012unified} and \citet{cai2023statistical}.
\begin{lemma}
Assuming \textit{Assumptions 2.2-2.3} and \textit{3.1-3.2} hold and given that $\xi^\top \Gamma^{-2} \xi \leq V$ for some constant $V>0$, where
\begin{equation*}
\Gamma = \mathbb{E} [g'(X_i^\top\hat \beta_{c1})X_i X_i^\top],
\end{equation*}
it follows that the optimization problem defined by \eqref{eq:optm_dense} is feasible with probability tending to 1.
Specifically, the constraints are met by setting $\gamma_i^* = \frac{1}{n_c}\xi^\top \Gamma^{-1}X_i$.
\label{lemma:optm}
\end{lemma}
Lemma \ref{lemma:optm} plays a pivotal role in facilitating subsequent theorems related to the asymptotic analysis.
It provides a feasible solution, $\gamma_i^*$, that forms an upper bound for the asymptotic variance of the debiased estimator.
Given that $\Gamma$ is invertible for any given $\hat{\beta}_{c1}$ by \textit{Assumptions 2.3} and \textit{3.2}, the condition $\xi^\top \Gamma^{-2} \xi \leq V$ fundamentally serves to constrain the magnitude of $||\xi||_2$.
\begin{theorem}
Assuming the conditions of Lemma \ref{lemma:optm} and the following constraints: a minimum estimand size of
$||\xi||_{\infty} \geq \kappa >0$,
a lasso penalty parameter is $\lambda_n = O(\sqrt{\frac{\log p}{n_c}})$, and a sparsity level $k \ll \frac{\sqrt{n_c}}{\log p \log n_c}$,
the estimator defined by \eqref{eq:est_dense} is asymptotically normal. Specifically,
\[(\hat{\theta} - \theta)/\sqrt{v} \stackrel{d}{\rightarrow} N(0,1),\]
with $v:=\sum_{i\in \mathcal{S}_2} \gamma_i^2 g(X_i^\top \beta_c)(1-g(X_i^\top \beta_c))$ and $v = O_p(1/n_c)$. In addition, the variance $v$ can be consistently estimated by $\hat{v}:=\sum_{i \in \mathcal{S}_2} \gamma_i^2 g(X_i^\top\hat{\beta}_{c1})(1-g(X_i^\top\hat{\beta}_{c1}))$, i.e., $\hat{v} \stackrel{p}{\rightarrow} v$.
\label{thm:dense}
\end{theorem}
The lower bound on $\|\xi\|_\infty$ excludes near-zero targets such as $\xi = (1/\sqrt{p},...,1/\sqrt{p})$.
The selection of the lasso tuning parameter's order aligns with established standards found in the literature, as exemplified by \citet{bickel2009simultaneous}.
However, it is important to note that the sparsity level requirement is slightly stronger than that in \citet{athey2018approximate}, where $k \ll \frac{\sqrt{n_c}}{\log p}$.
This requirement in the condition is necessitated by the GLM framework; the error decomposition \eqref{eq:mu_chat_error_byH} introduces Taylor expansion's remainder terms, and a stronger constraint is required to ensure their negligibility.
\subsection{Debiased Estimation for Causal Parameters}
We now analyse ATET.
The main distinction lies in the vector $\xi$: it is now considered as a random vector, whereas in the previous analysis, it was treated as a deterministic vector.
More precisely, we have $\xi = \frac{1}{n_t}\sum_{{i:D_i=1}} g'(X_i^\top \hat{\beta}_{c1})X_i$, leading to:
\begin{equation}
\hat{\mu}_c - \mu_c = \hat{\theta} - \theta +\frac{1}{n_t} \sum_{i:D_i = 1} \Delta_i^t.
\end{equation}
The second term is the second-order remainder of the Taylor expansion and negligible under certain conditions.
For this target we replace \textit{Assumption 3.1} with conditions more common in causal inference:
\begin{assumption}
$\{X_i\}_{1\leq i\leq n}$ conditional on $D_i=d$, where $d\in\{0,1\}$, are independent and identically distributed sub-Gaussian random vectors, that is, there exists a constant $c \in \mathbb{R}$ satisfying $\mathbb{E}[\exp\{v^\top X\}|D=d]\leq \exp\{||v||_2^2 c^2/2\}$ for all $v \in \mathbb{R}^p$ after re-centring.
\end{assumption}
\begin{assumption}
There is a constant $\epsilon>0$ such that $\epsilon \leq p(x)\leq 1-\epsilon$ for all $x\in \mathbb{R}$.
\end{assumption}
\begin{theorem}
Under \textit{Assumptions 2.1-2.3} and \textit{3.2-3.4} and assuming that the overall odds of receiving treatment, denoted by $\mathbb{P}(D=1)/\mathbb{P}(D=0)$, converge to a limit $ \rho \in (0,\infty)$,
along with the constraints on the sparsity level $k \ll \frac{\sqrt{n_c}}{\log p \log n_c}$, the lasso penalty parameter $\lambda_n = O(\sqrt{\log p/n_c})$ and $\ell_\infty$ norm of $||\mathbb{E}[g'(X^\top \hat{\beta}_{c1})X|D=1]||_\infty \geq \kappa>0$,
the estimator defined by \eqref{eq:mu_chat}, with weights as given by the constrained form of \eqref{eq:optm_initial}, is asymptotically normal:
\[(\hat{\mu}_c - \mu_c)/\sqrt{v_c} \stackrel{d}{\rightarrow} N(0,1),\]
where $v_c :=\sum_{i\in \mathcal{S}_2} \gamma_i^2 g(X_i^\top \beta_c)(1-g(X_i^\top \beta_c))$, $v_c= O_p(1/n_c)$ and the variance estimate $\hat{v}_c:=\sum_{i \in \mathcal{S}_2} \gamma_i^2 g(X_i^\top\hat{\beta}_{c1})(1-g(X_i^\top\hat{\beta}_{c1}))$ converges in probability to $v_c$ at a rate faster than $1/n_c$.
Consequently, the following results hold:
\begin{list}{(\alph{bean})}{\usecounter{bean}}
\item For the $\mu_c$ estimation:
\[(\hat{\mu}_c - \mu_c)/\sqrt{\hat{v}_c} \stackrel{d}{\rightarrow} N(0,1).\]
\item For the causal parameter estimation, given by $\hat{\tau} = \bar{Y}_t - \hat{\mu}_c$, the following result holds:
\begin{equation}
\frac{\hat{\tau} - \tau}{\sqrt{\hat{v}_\tau}} \stackrel{d}{\rightarrow} N(0,1)
\end{equation}
with $\hat{v}_\tau = \hat{v}_c+\hat{v}_t$ and $\hat{v}_t := \frac{1}{n_t^2}\sum_{\{i:D_i = 1\}}(Y_i - \bar{Y}_t)^2$, where $\hat{v}_\tau = O_p(1/n)$.
\end{list}
\label{thm:formalCF}
\end{theorem}
Theorem \ref{thm:formalCF} establishes root-$n$ inference for the treatment effect without requiring the consistent estimation of the propensity score, and provides feasible standard errors.
\section{Simulation Study}
\setcounter{equation}{0}
\setcounter{theorem}{0}
\setcounter{assumption}{0}
\setcounter{proposition}{0}
\setcounter{corollary}{0}
\setcounter{lemma}{0}
\setcounter{example}{0}
\setcounter{remark}{0}
We study a two-stage model with $Y\sim\mathrm{Bernoulli}\{g(X^\top\beta_Y+D)\}$ and $D\sim\mathrm{Bernoulli}\{g(X^\top\beta_D)\}$, where $g(x)=\exp(x)/(1+\exp(x))$. Covariates are $X\sim N(0,\Sigma)$ with $\Sigma_{ij}=\rho^{|i-j|}$.
The outcome coefficients are sparse, $(\beta_Y)_j\propto 1/j^2$ and the propensity coefficients are either sparse $(\beta_D)_j\propto 1/j^2$ or dense $(\beta_D)_j\propto 1/\sqrt{j}$.
We compare our estimator with sample splitting (DB1), cross-fitting (DB2), and no splitting (DB6) to Na\"{i}ve, Regression with simple averaged residuals by \eqref{eq:mu_chat_reg} , IPW by \eqref{eq:mu_chat_weig}, DML by \eqref{eq:mu_chat_DL}, AML by \eqref{eq:mu_chat_AML}, and ARB. Performance is summarized by relative mean-squared error (MSE).\footnote{More designs, full tuning details, and additional metrics are reported in the Online Supplement.}
\begin{table}
\label{tab:RMSE_1}
\caption{Mean‐squared error of different estimators}
\begin{center}
\begin{tabular*}{\textwidth}{@{}lcccccccc@{}}
\hline\hline
\multirow{3}{*}{Method}
& \multicolumn{4}{c}{Sparse $\beta_D$}
& \multicolumn{4}{c}{Dense $\beta_D$} \\
\cline{2-5} \cline{6-9}
& \multicolumn{2}{c}{$\|\beta_D\|_2=1$}
& \multicolumn{2}{c}{$4$}
& \multicolumn{2}{c}{$\|\beta_D\|_2 = 1$}
& \multicolumn{2}{c}{$4$} \\
\cline{2-3}\cline{4-5}\cline{6-7}\cline{8-9}
& $\|\beta_Y\|_2=1$ & $4$
& $1$ & $4$
& $\|\beta_Y\|_2=1$ & $4$
& $1$ & $4$ \\
\hline
Naïve & 1.176 & 4.565 & 2.214 & 8.785 & 0.430 & 1.539 & 0.646 & 2.338 \\
Regression & 0.162 & 0.224 & 0.262 & 0.411 & 0.112 & 0.165 & 0.144 & 0.210 \\
IPW & 0.196 & 0.630 & 0.373 & 1.222 & 0.162 & 0.487 & 0.227 & 0.711 \\
DML1 & 0.164 & 0.255 & 0.330 & 0.431 & 0.236 & 0.305 & 0.494 & 0.698 \\
DML2 & 0.113 & 0.153 & 0.290 & 0.339 & 0.123 & 0.201 & 0.631 & 0.369 \\
AML & 0.115 & 0.169 & 0.170 & 0.361 & 0.098 & 0.156 & 0.112 & 0.168 \\
ARB & 0.094 & 0.141 & 0.210 & 0.315 & 0.095 & 0.158 & 0.113 & 0.191 \\
DB1 & 0.126 & 0.162 & 0.202 & 0.291 & 0.125 & 0.186 & 0.146 & 0.204 \\
DB2 & 0.087 & 0.123 & 0.144 & 0.221 & 0.098 & 0.146 & 0.110 & 0.164 \\
DB6 & 0.067 & 0.089 & 0.116 & 0.161 & 0.070 & 0.099 & 0.078 & 0.108 \\
\hline\hline
\end{tabular*}
\end{center}
\footnotesize
\textbf{Note:}
$n=500$, $p=800$, $\rho=0.5$, $n_{\mathrm{sim}}=1000$. DML1 = 2-fold sample-split; DML2 = 2-fold cross-fit. AML = 2-fold cross-fit. DB1 = proposed (2-fold sample-split), DB2 = proposed (2-fold cross-fit), DB6 = proposed (no split).
\end{table}
According to Table 1, DB6 attains the lowest MSE in all cells; DB2 is second. When the propensity model is complex (dense $\beta_D$), IPW and DML underperform the simple Regression imputation, consistent with instability from inverse weights. Estimators that avoid inverse weights such as ARB, AML, and ours, are more stable. However, when the linear approximation is poor, as shown in the supplementary materials, ARB approach deteriorates due to misspecification. DB2 and DB6 outperform AML, reflecting our method’s explicit variance control in addition to bias reduction.
In addition, according to the coverage table in the supplementary materials, confidence intervals based on Theorem~\ref{thm:formalCF} achieve coverage close to 95\% and improve with $n$; DB6 has the shortest intervals.
\section{Empirical Illustration}
\setcounter{equation}{0}
\setcounter{theorem}{0}
\setcounter{assumption}{0}
\setcounter{proposition}{0}
\setcounter{corollary}{0}
\setcounter{lemma}{0}
\setcounter{example}{0}
\setcounter{remark}{0}
We reanalyse the NSW dataset of \citet{dehejia1999causal}. The treated group comprises 185 men randomized to job training between December 1975 and January 1978; the control group is a non-experimental PSID comparison sample of 2,490 men. We define the binary outcome $Y_i$ to indicate whether an individual’s 1978 earnings exceed zero (i.e.\ employment status). To adjust for confounding, we construct a covariate vector $X_i$ of 60 features, including age, years of education, indicators for Black and Hispanic, an indicator for marital status, and earnings in 1974 and 1975, together with their polynomial and interaction terms. We assume unconfoundedness conditional on $X_i$.
\begin{table}
\caption{Point estimates of ATET in NSW application}
\begin{center}
\begin{tabular*}{\textwidth}{@{}lccccccc@{}}
\hline
\hline
Method & Na\"{i}ve & Regression & IPW & DML2 & AML & ARB & DB6 \\
\hline
Estimates & -0.1585 & 0.2523 & 0.1174 & -0.5635 & 0.0587 & 0.1507 & 0.1122 \\
\hline \hline
\end{tabular*}
\end{center}
\label{tab:EmpIll}
\footnotesize
\textbf{Note:} See note in Table 1.
\end{table}
Table \ref{tab:EmpIll} reports the estimates of ATET using different methods.
For reference, the simple difference‐in‐means estimate within the randomized experiment is
0.1106. Our estimator without splitting (DB6) is 0.1122, closest to this benchmark; IPW is 0.1174 and also close. The other approaches deviate more substantially from the experimental result. IPW performs well here because the logistic propensity appears well specified, as discussed by \citet{dehejia1999causal}, and trimming extreme weights has little effect on the estimate.
Moreover, the 95 \% confidence interval for our DB6 estimator is
[0.0301,0.1944], nearly identical to the randomized‐experiment interval
[0.0253,0.1959]. This suggests that our method delivers both accurate point estimates and valid inference when applied to non‐experimental data.
\section{Conclusion}
This paper proposes a debiased estimator for causal effects in high-dimensional GLMs with binary outcomes. The estimator augments a regularized GLM plug-in with weights from a single convex program that balances link-derivative–weighted covariates and controls variance, avoiding the inversion of estimated propensity scores. We establish $\sqrt n$-consistency, asymptotic normality, and a feasible variance estimator.
Our analysis focuses on ATET under unconfoundedness for $Y(0)$ and a GLM for $Y(0)\mid X$. Analogous estimators for ATE and ATEC are immediate, but identification then requires $D\!\perp\! Y(1)\mid X$ and a GLM for $Y(1)\mid X$. Sample splitting is used for theory; simulations indicate that the no-split implementation performs best, with cross-fitting close behind.
Compared to doubly robust strategies such as \citet{farrell2015robust}, \citet{belloni2017program}, and \citet{chernozhukov2018double} that utilize the estimation of both outcome and propensity score models, our method leans on the structural assumption of generalized linearity in the outcome model.
While their methods offer robustness to misspecifications, our approach could serve as a viable alternative when generalized linearity is believed, which is not a overly restrictive assumption under the high-dimensional context, as it circumvents the instability of inverse propensity score estimates.
Two directions could be useful for future work: a data-driven choice of the bias–variance tuning parameter and extensions that relax the known-link assumption (e.g., single-index links) while retaining the balancing-and-variance optimization.
\bibliographystyle{chicago}
\bibliography{ref.bib}