Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
46,309 characters · 14 sections · 40 citation commands
Causal Inference in High-Dimensional Generalized Linear Models with Binary Outcomes
\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.
In observational studies, an unconfoundedness assumption—treatment assignment is 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 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) 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.
This example underscores the need to mitigate bias for causal targets when regularization is necessary.
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 belloni2014high. Doubly robust and orthogonal estimators combine outcome and propensity models to obtain $\sqrt n$ inference under first-stage rate conditions, as developed by farrell2015robust, belloni2017program, and chernozhukov2018double, and are related to debiased Lasso, as in van2014asymptotically and 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 graham2012inverse and graham2016efficient, relies on inverse propensities and can be unstable in high dimensions. The ARB estimator of 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 hirshberg2021augmented. Compared with ARB, we target GLM moments, use a different program, and obtain a feasible variance estimator. Riesz-based AML estimators of chernozhukov2022automatic and 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. 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. belloni2016post and 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.
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$.
\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 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$.
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.
Our target parameter is
where $n_t$ is the number of treated units. Under Assumption 2.2,
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$.
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.
Compared to the approach by 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 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 (ref) 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 hirshberg2021augmented. Their analysis relies on Donsker classes; our high-dimensional setting does not, and we instead use standard assumptions for regularized GLMs.
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 (ref). The direct plug-in estimator is
and the plug-in estimator with simple averaging of residuals:
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$:
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 (ref) improves on the plug-in (ref) 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 (ref). In contrast, our estimator (ref) 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 (ref) can be bounded as follows:
In Equation (ref), 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
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 hirano2003efficient. We also considers the Double Machine Learning (DML) estimator of chernozhukov2018double, which combines regression imputation with inverse propensity weighting,
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 chernozhukov2022automatic. It avoids inverse propensity weights by estimating a Riesz representer $\alpha(\cdot)$ and then forming
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.
\setcounter{equation}{0} \setcounter{theorem}{0} \setcounter{assumption}{0} \setcounter{proposition}{0} \setcounter{corollary}{0} \setcounter{lemma}{0} \setcounter{example}{0} \setcounter{remark}{0}
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 van2014asymptotically and javanmard2014confidence. Here we allow $\xi$ to be dense and consider $\theta = \xi^\top \beta_c$. Linked to Algorithm 1, the estimator is given by
where $\boldsymbol{\gamma}$ solves the constrained program
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 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 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 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 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:
While 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. 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 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 negahban2012unified and cai2023statistical.
Lemma (ref) 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 Assumptions 2.3 and 3.2, the condition $\xi^\top \Gamma^{-2} \xi \leq V$ fundamentally serves to constrain the magnitude of $||\xi||_2$.
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 bickel2009simultaneous. However, it is important to note that the sparsity level requirement is slightly stronger than that in athey2018approximate, where $k \ll \frac{\sqrt{n_c}}{\log p}$. This requirement in the condition is necessitated by the GLM framework; the error decomposition (ref) introduces Taylor expansion's remainder terms, and a stronger constraint is required to ensure their negligibility.
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:
The second term is the second-order remainder of the Taylor expansion and negligible under certain conditions. For this target we replace Assumption 3.1 with conditions more common in causal inference:
Theorem (ref) establishes root-$n$ inference for the treatment effect without requiring the consistent estimation of the propensity score, and provides feasible standard errors.
\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 (ref) , IPW by (ref), DML by (ref), AML by (ref), 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.}
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) achieve coverage close to 95% and improve with $n$; DB6 has the shortest intervals.
\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 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$.
Table (ref) 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 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.
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 farrell2015robust, belloni2017program, and 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.