EconBase
← Back to paper

Heterogeneous Treatment Effect Bounds under Sample Selection with an Application to the Effects of Social Media on Political Polarization

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.

104,581 characters · 22 sections · 131 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Heterogeneous Treatment Effect Bounds under Sample Selection with an Application to the Effects of Social Media on Political Polarization

titlepage\thispagestyle{empty} \begin{abstract} \singlespacing We propose a method for estimation and inference for bounds for heterogeneous causal effect parameters in general sample selection models where the treatment can affect whether an outcome is observed and no exclusion restrictions are available. The method provides conditional effect bounds as functions of policy relevant pre-treatment variables. It allows for conducting valid statistical inference on the unidentified conditional effects. We use a flexible debiased/double machine learning approach that can accommodate non-linear functional forms and high-dimensional confounders. Easily verifiable high-level conditions for estimation, misspecification robust confidence intervals, and uniform confidence bands are provided as well. We re-analyze data from a large scale field experiment on Facebook on counter-attitudinal news subscription with attrition. Our method yields substantially tighter effect bounds compared to conventional methods and suggests depolarization effects for younger users. \end{abstract} Keywords: Affective polarization; Debiased/double machine learning; Effect bounds; Facebook; Partial identification \\ JEL classification: C14, C21, D72, L82

\setlength\abovedisplayskip{3pt} {3pt}

\pagenumbering{arabic}

Introduction

In this paper, we propose a novel method for estimation and inference for bounds of heterogeneous causal effects when outcome data is only selectively observed and no exclusion restrictions or instruments are available. In particular, we are concerned with the case when the treatment of interest itself can affect the selection process and when effects are heterogeneous along both observable and unobservable dimensions. The bounds are derived from a conditional monotonicity assumption in the selection equation. They can be used to study the effects of interventions on always-taker units, e.g. the effects of active labor market policies on earnings on the population that is working regardless of whether they were subject to the intervention or not.\footnote{Note that in contrast to the typical setup and nomenclature in the instrumental variables literature, the principal strata are defined with respect to potential selection state caused by treatment, not potential treatment states caused by an instrument. The always-taker stratum is also sometimes referred to as inframarginal or always-observed.} They can also be applied to obtain credible bounds in experimental studies where the original treatment can affect selection.

When applying established partial identification approaches for similar sample selection problems in practice, unconditional or subgroup specific effect bounds horowitz2000nonparametric,zhang2003estimation,lee2009training,semenova2023generalized are often wide and thus too uninformative for assisting policy. Narrower bounds that also exploit covariate information can be used for a better targeting of interventions under weaker, i.e. more credible, conditions compared to restrictive point-identified methods that require exclusion restrictions and/or distributional assumptions. The nonparametric heterogeneity based approach in this paper helps to tighten bounds along policy relevant pre-treatment variables. This is due to the fact that the severity of the identification problem, i.e. the width of the identified set, can vary substantially along the confounding dimensions most associated with the heterogeneity variables of interest. In addition, the procedure has significant advantages over calculating bounds within discrete partitions of the data: It can accommodate continuous variables and, as it extracts signals for the bounds before conditioning on heterogeneity variables, exploits larger samples as well as potential group patterns/restrictions for modeling selection probabilities and other relevant nuisance functions.

The method can also incorporate a high-dimensional number of confounders building on debiased machine learning (DML) methodology chernozhukov2018double,semenova2021debiased. We derive explicit high-level conditions regarding the quality of the nuisance quantity estimators that can be verified in a variety of settings for popular nonparametric or machine learning estimators such as high-dimensional sparse regression, deep neural networks, or random forests. We provide analytical confidence intervals for heterogeneous effects that are robust against different types of model misspecification as well as uniform confidence bands using a multiplier bootstrap.

figure[figure omitted — 1,072 chars of source]

Figure (ref) illustrates the proposed method for a one-dimensional heterogeneity analysis. It plots identified sets, confidence intervals, and confidence band for the causal effect in dependence of a pre-treatment variable $z$. These could e.g. be bounds on the effect of job training on earnings as a function of pre-training earnings. We can see that unconditional analysis cannot rule out a zero effect while the confidence intervals of the heterogeneous bounds clearly suggest significant negative ($0.5 < z < 0.8$) and positive effects ($z > 0.95$) leading to different policy recommendations. Note that these conclusions are achieved not just by heterogeneous locations of the bounds but also by narrower widths for certain $z$ values. Thus, heterogeneous bounds can reduce uncertainty stemming from weaker identification assumptions for certain sub-populations. Monte Carlo simulations suggest that the presented confidence intervals perform well in finite samples.

Our application is concerned with the effect of social media news consumption on political polarization. In 2022, over 70% of US adults consumed news on social media.\footnote{Pew Research Center, Survey of U.S. adults, July 18-Aug 21, 2022.} The consequences of social media and online news consumption on political polarization are of major importance: High partisan attitudes threaten the functioning of society and democracy as well as trust in public and private institutions phillips2022affective. Many democratic societies have been experiencing significant changes in partisan attitudes since the broad roll-out of social networks such as Facebook and Twitter (now X). However, the exact contributions are under debate (haidtONGOINGsocialmediaPolitical, ongoing).

We study the effect of Facebook news subscription on affective polarization. Facebook is the most dominant social media site for news consumption among US adults (31% regularly get news on this site). Affective polarization measures relative attitudes toward opposing partisans in terms of (dis)like and (dis)trust. We re-examine data collected by levy2021social who employs a large scale field experiment on Facebook where units are nudged towards subscribing to popular media outlets with clear partisan ideology such as FoxNews or MSNBC. In particular, we measure the effects of the counter-attitudinal treatment in terms of political leaning on affective polarization after two months. The outcome measure suffers from large differential attrition rates between treatment and control groups and within political ideology (over 50% total).

Overall, the findings do not contradict levy2021social who suggests a decrease in affective polarization by $-0.06$ standard deviations not corrected for attrition. The relative size of the identified set benefits from the inclusion of covariates even under the randomized treatment assignment. Independently of nuisance parameter specifications, our identified sets are around $[-0.08,0.01]$ and 58% tighter compared to conventional monotonicity bounds that cannot rule out a moderate increase in affective polarization ($[-0.16,0.06]$). Looking at subgroup heterogeneity, we are getting closer to point identification for some groups. For example, for conservatives and 18-year olds, the identified sets are $[-0.07,-0.04]$ and $[-0.09,-0.04]$ respectively. Additionally accounting for the statistical uncertainty, there is weak evidence in favor of depolarization effect for the younger users.

The paper is structured as follows: Section (ref) discusses the methodological literature. Section (ref) contains model, estimator, and confidence intervals. Section (ref) provides technical assumptions, confidence bands, and large sample properties. Section (ref) presents Monte Carlo simulations. Section (ref) contains the empirical study. Section (ref) concludes. All proofs and extensions are in the Appendix. The \verb|R| package \verb|HeterogeneousBounds| and replication notebook can be found on the author's personal web-page.

Methodological Literature

Bounds for causal effects under weak assumptions have been considered in a series of papers by Charles Manski and others, see molinari2020microeconometrics for a comprehensive overview. horowitz2000nonparametric develop nonparametric bounds for treatment effects in selected samples. zhang2003estimation consider bounds for always-taker units under a monotonicity assumption regarding the effect of treatment on selection, commonly imposed in generic sample selection models, and/or stochastic dominance assumptions on the potential outcomes. imai2008sharp demonstrates sharpness of these bounds. huber2015sharp consider similar sharp bounds for other principal strata. lee2009training provides asymptotic theory for zhang2003estimation bounds using (conditional) monotonicity. Without further assumptions, these bounds are only applicable unconditionally or for low-dimensional discrete partitions of the covariate space and now commonly referred to as “Lee bounds”. semenova2023generalized provides “generalized Lee bounds” under a conditional monotonicity assumption. Our paper uses similar identification assumptions. semenova2023generalized also allows for high-dimensional and continuous confounders and generalizes the approach to multiple outcomes but does neither consider heterogeneity analysis nor misspecification-robust inference. bartalotti2021identifying also propose identification of bounds for always-takers within a marginal treatment effect framework. Using monotonicity and stochastic dominance, they tighten effect bounds based on underlying treatment propensities. However, they do neither address heterogeneity beyond the propensity score, asymptotic properties, inference, nor potential misspecification. Moreover, their method is not suitable for many confounding variables without imposing additional parametric assumptions.

Our work is also directly related to the literature on robust or (Neyman-)orthogonal moment functions and DML. For point-identified parameters, there are now many approaches that use machine learning and orthogonal moments in both experimental and observational studies, see e.g. belloni2014inference, farrell2015robust, chernozhukov2018double, and wager2018estimation. chernozhukov2018double develop a canonical framework that can be used for inference on low-dimensional target parameters such as the average treatment effect. In the context of heterogeneity analysis, orthogonal moments have been exploited in point-identified problems by using them as pseudo-outcomes in (nonparametric) regression models to obtain predictive causal summary parameters, see e.g. lee2017doubly, fan2020estimation, semenova2021debiased, and heiler2021effect or knaus2022double for an overview. Our localization approach is closest to semenova2021debiased who estimate heterogeneity parameters via nonparametric projections using least squares series methods belloni2015some. Employing a series approach is particularly useful for studying potential misspecification since the heterogeneity step can then be reformulated as a relaxed moment inequality problem with deviation parameters that are proportional to the difference between effect bounds and their linear predictors. It also allows for construction of uniform confidence bands under correct specification.

semenova2023generalized constructs DML estimators for unconditional Lee-type bounds. semenova2023debiased considers partially identified parameters for linear moment functions with DML. A crucial point in both of these papers is that, while the effect of interest might not be point-identified, bounds themselves are characterized by well-understood convex moment problems or the corresponding support function. The same applies to the heterogeneous bounds considered in this paper. We derive the asymptotic distribution of the nonparametric heterogeneous DML based estimator for the identified set. This nests the univariate generalized Lee bounds by semenova2023generalized as a special case. In addition, we also provide the bounds and inference theory for separate subgroups of always-takers defined by their monotonicity type, i.e. their effect sign of treatment on selection. In contrast to semenova2023debiased, the moment functions are nonlinear in the outcome. The use of machine learning (random forests) for Lee-type bounds has also been heuristically discussed by cornelisz2020addressing. They do, however, not provide any formal theory for estimation and inference. olma2021nonparametric also discusses nonparametric estimation of generic truncated conditional expectations based on similar conditional moments as the ones used in this paper. He suggests kernel estimation for both stages and provides pointwise linearization and distribution results. In contrast, our approach can handle generic first-stage learners that also work in setups with high-dimensional confounding. Moreover, olma2021nonparametric does neither consider local power improvements, misspecification, nor uniform inference.

Inference for partially identified parameters in sample selection models is a non-trivial task. Relying on quantiles of large sample distributions of the bounds can be overly conservative for the actual effect of interest. The relevant uncertainty for the latter depends on the actual width of the identified set. If small, then deviations from a null value are likely to occur in both positive and negative direction, i.e. the problem is effectively two-sided. If large, then uncertainty in one direction dominates, rendering the testing problem close to one-sided. imbens2004confidence consider related confidence intervals for partially identified parameters. Their method is not uniformly valid with regards to width of the underlying identified set due to an implicit superefficiency assumption. stoye2009more suggests to artificially impose superefficiency via shrinkage. andrews2010inference provide a more general approach within a moment inequality framework. andrews2013inference and andrews2017inference study inference based on (many) conditional moment inequalities in the presence of nuisance parameters. While the latter can be infinite dimensional as in this paper, their inference target is essentially parametric. Our approach is closest in spirit to andrews2014nonparametric who consider nonparametric estimation based on conditional moment inequalities. Their method can be adjusted to conduct pointwise inference on conditional effect bounds. However, as they study a relatively general setup, they do neither address local power improvements for effects, estimated nuisances, misspecification, nor uniform inference. chernozhukov2019inference also consider inference using generic moment inequalities in possibly high dimensions allowing for nonparametric sieve-type parameters. In contrast, we exploit the specific shape of the identified set as well as the restriction to low-dimensional heterogeneity variables to obtain stronger asymptotic results that can be used for power improvements.\footnote{chernozhukov2019inference, Appendix B also heuristically discusses approximate moment inequalities/misspecification arising from estimation using for parametric nuisance models.}

When analyzing heterogeneity, we are concerned with effect bounds at potentially many points and thus there is an additional risk of local misspecification compared to unconditional effects. For instance, even when the conditional moments are correctly specified and ordered, their linear predictors might intersect. In addition, there can be global misspecification, e.g. due to misspecified selection probabilities, that can lead to a reversed ordering of the orthogonalized local bounds in both the population and in finite samples. Under such misspecification, aforementioned bounds and moment inequality methods can produce empty or very narrow confidence regions suggesting spuriously precise inference andrews2019inference. andrews2019inference propose inference methods that extend the notion of coverage to pseudo-true parameter sets in general moment inequality frameworks, see also stoye2020simple for a power improvement for generic regular parametric bound estimators. This paper adapts the inference method of stoye2020simple to nonparametric heterogeneous bounds with machine learning in the first stage.

Methodology

Model and Identification Assumptions

In this section, we introduce the sample selection model and the main identification assumptions followed by a brief review of the construction of effect bounds. We then show how to exploit the latter to construct, estimate, and conduct inference on the heterogeneous partially identified effect parameters.

Assume for $i=1,\dots,n$ we observe iid data $W_i = (X_i',D_i,S_i,Y_iS_i)$ where $X_i$ is a vector of predetermined covariates supported on $\mathcal{X} \subseteq \mathbb{R}^d$, $D_i$ is a treatment indicator, $S_i$ indicates whether the outcome is observed, i.e. we only observe $Y_iS_i$ and not $Y_i$. We would like to evaluate the average causal effect of $D_i$ on outcome $Y_i$ for the units that are selected under both treatment or control condition. We focus on the case of known treatment propensities $e(x) = P(D_i=1|X_i=x)$ during the exposition, but the methodology can also be extended to observational studies where they are unknown and have to be estimated. The prototypical setup is depicted by the graph in Figure (ref) derived from the nonparametric structural equation model.

figure[figure omitted — 1,517 chars of source]

The model nests the classic sample selection model heckman1979sample where the treatment of interest enters in selection step, see lee2009training for a more restrictive parametric sample selection model representation. It allows for the treatment $D$ to affect the selection $S$. The model does not restrict the relationship between selection and the partially unobserved outcome of interest. In particular, there can be unobservables $U$ related to both selection and potential outcomes. In addition, there are no variables available which only affect selection that could be used to identify causal effects on the partially unobserved outcome via instrumental variable based methods.

Based on this model, we define unit potential outcomes $Y_i(1)$ and $Y_i(0)$ and potential selection indicators $S_i(1)$ and $S_i(0)$ for units $i=1,\dots,n$. They correspond to a unit's value in outcome or selection when the treatment is exogeneously set to $D_i = 1$ or $D_i=0$ respectively under the model in Figure (ref). Note that potential selection and potential outcomes can still be dependent. This implies that, even conditional on $X_i$, causal effects can differ for different sub-types defined by their respective $S_i(1)$ and $S_i(0)$ variables. The target is to evaluate the causal effect of $D_i$ for the units that are selected in both control and treatment state, i.e. units for which $S_i(1)=S_i(0)=1$. As $D$ affects $S$, a comparison of treated and control units for selected units does not yield a valid causal comparison for the sub-population $S_i(1) = S_i(0) = 1$. However, the model implies the following conditional independence relationships in Rubin-Neyman potential outcome notation\footnote{These are the marginal versions of the joint independence assumption $Y_i(0), Y_i(1), S_i(0), S_i(1) \rotatebox[origin=c]{90}{$\models$}\ D_i | X_i$ considered in lee2009training that yield the same identification results.}

align[align omitted — 118 chars of source]

Condition (ref) cannot be exploited in a standard selection on observables strategy for point identification as, even conditional on $X_i$, we only observe the selected subset of potential outcomes for which $S_i = 1$. We now impose a monotonicity assumption on the selection equation. In particular, treatment can affect selection positively or negatively. The sign, however, must be uniquely determined by the vector of covariates $X_i$. This is a weak or conditional monotonicity assumption based on ex-ante unknown partitions:

ass[Weak/Conditional Monotonicity] There exist partitions of the covariate space $\mathcal{X} = \mathcal{X}_+ \cup \mathcal{X}_-$ such that $P(S_i(1)\geq S_i(0)|X_i\in\mathcal{X}_+) = P(S_i(1)\leq S_i(0)|X_i\in\mathcal{X}_-) = 1$.

This assumptions is consistent with additive separability of observables and unobservables in an otherwise unrestricted selection equation. We also assume that, with positive probability, there are comparable units between selected treated and selected non-treated as well as between treated and non-treated overall:

ass[Multiple Overlap] For all $x\in\mathcal{X}$ and $d \in \{0,1\}$ we have that $0 < P(S_i=1|X_i=x,D_i=d) < 1$ and $0 < P(D_i=d|X_i=x) < 1$.

We also rule out the case where treatment does not affect selection as in this case effects are point-identified:

ass[Margin] There is a constant $c > 0$ and a set $\mathcal{\bar{X}} \subset \mathcal{X}$ with $P(\mathcal{X}\backslash\mathcal{\bar{X}}) = 0$ such that $\inf_{x\in\mathcal{\bar{X}}}\ |P(S_i=1|X_i=x,D_i=1)-P(S_i=1|X_i=x,D_i=0)| > c$.

This assumption is most plausible when there are discrete variables and/or continuous variables with bounded support. For example, a standard binary response model for selection where $D_i$ enters with a fixed, non-zero parameter then obeys this condition. Assumption 3.3 could be weakened to a margin condition that restricts the behavior of the treatment on selection effect distribution around zero with some technical modifications semenova2023generalized. We abstract from these issues for simplicity.

At last, we need a standard continuity condition in each treatment arm:

ass[Continuity] The conditional outcome distributions $P(Y_i\leq y|D_i=1,S_i=1,X_i=x)$ and $P(Y_i\leq y|D_i=0,S_i=1,X_i=x)$ are continuous almost everywhere.

For identification of the sharp effect bounds, we first outline the case of strong or unconditional monotonicity as considered in zhang2003estimation and lee2009training, i.e. when $\mathcal{X}^- = \{ \emptyset \}$ or $S_i(1)\geq S_i(0)$ with probability one. In this case, selected units within the control group must be always-takers or inframarginal in the sense of $S_i(1) = S_i(0) = 1$ while within the treated group there is a mixture of always-takers and compliers or marginal units who are induced to be selected by the treatment. Let the conditional causal effect of the always-takers be defined as

align[align omitted — 70 chars of source]

This is the expected causal effect for a unit with covariates $x$ whose outcome would be observed independently of its treatment status. Conditional on $X_i=x$, let $q_1(u,x)$ be the $u$-th quantile of the treated selected units, $s(d,x) = P(S_i=1|X_i=x,D_i=d)$ the selection probability, and $p_0(x) = {s(0,x)}/{s(1,x)}$ the share of always-takers relative to always-takers and compliers. zhang2003estimation show that the sharp bounds are given by

align[align omitted — 82 chars of source]

where

align[align omitted — 252 chars of source]

While these $x$-specific or “personalized bounds” can provide some insights, they are not very useful in applications with continuous variables and/or many discrete cells to evaluate. Moreover, in such cases it is often no longer possible to consistently estimate the conditional bounds and to construct asymptotically valid confidence intervals without further restrictions. This is equivalent to estimation of personalized treatment effects under unconfoundedness in high dimensions chernozhukov2018generic. Therefore, we propose to instead consider heterogeneous effect bounds conditional on a smaller, pre-specified (policy relevant) variables $Z_i$ supported on $\mathcal{Z}$. These can be (functions of) any observable variable possibly affecting treatment, selection, and outcome. Thus, without loss of generality, we can treat $Z_i = f(X_i)$ where $f(\cdot)$ is a measurable, possibly multi-valued function. As $\mathcal{Z}$ is of low dimension, this allows for conducting asymptotically valid inference for heterogeneous causal effects even if the original confounding dimension of $\mathcal{X}$ is large. The main parameter of interest is then given by

align[align omitted — 148 chars of source]

This can be interpreted as expected causal effect for the always-takers at “summary group” $f(X_i) = Z_i = z$.

Heterogeneous Effect Bounds

We now demonstrate how to obtain sharp bounds and moment estimators for (ref) using (strong) monotonicity. Using the sharp upper bounds $\theta_U(x)$ from zhang2003estimation or lee2009training yields

align[align omitted — 210 chars of source]

by monotonicity and the law of iterated expectations. The same steps apply analogously to the lower bound. Thus, for identification, it is sufficient to have two (non-centered) moment functions: (i) $\psi_B^+(W_i,\eta)$ that identifies the conditional always-taker bound scaled by the density of always-takers and (ii) $\psi_{S_0}(W_i,\eta)$ that identifies the conditional always-taker share itself. $\eta$ here is a vector of nuisance quantities, see Section (ref), Table (ref) for a complete definition. In particular, we need that these are valid conditional on $X_i$ at the true nuisance parameters $\eta = \eta_0$, i.e. that

align[align omitted — 121 chars of source]

In principle, many moment functions have these properties. For estimation and inference, however, we focus on the moment functions $\psi_B^+(W_i,\eta)$ as developed by semenova2023generalized in combination with the well-known doubly robust/augmented inverse probability weighting moment for the expected potential selection robins1994estimation. This will be important to control the large sample properties with unknown nuisance parameters, see Section (ref). We can then express the sharp heterogeneous bounds as

align[align omitted — 212 chars of source]

The corresponding estimators can be obtained by replacing all the conditional expectations with their (nonparametric) projections, see Section (ref).

Weak Monotonicity

We now extend the identification to the weak/conditional monotonicity case. We consider the separate effects for both positive and negative monotonicity partitions as well as the conventional aggregate always-taker effect in the presence of both. Denote $\mathbbm{1}^+(x) = \mathbbm{1}(p_0(x) < 1)$ and $\mathbbm{1}^-(x) = \mathbbm{1}(p_0(x) > 1)$. Given $X_i=x$ the partitions are identified. We can write the conditional always-takers effect as a piece-wise combination

align*[align* omitted — 85 chars of source]

We can correspondingly define partition-specific, truncated estimands as

align[align omitted — 245 chars of source]

and equivalently for $\theta_B^+(x)$ and $\theta_B^-(x)$. As in the previous section, we can now apply monotonicity within the different partitions to bound $\theta_{AT}(z)$. In particular

align[align omitted — 462 chars of source]

and analogously for the lower bound. Thus, for identification, it is sufficient to have four moment functions that again identify the within partition always-taker shares and partition specific effect bounds scaled with the corresponding density, i.e.

align[align omitted — 349 chars of source]

Overall, we obtain that

align[align omitted — 89 chars of source]

with

align[align omitted — 436 chars of source]

For the different monotonicity types, separate bounds are given by

align[align omitted — 167 chars of source]

with

align[align omitted — 520 chars of source]

The aggregate and partition-specific estimands are ratios. They differ by the denominator as, for the aggregate bounds, each partition-specific effect is weighted with its respective conditional probability. These bounds are sharp. In the case of $Z_i = 1$, the aggregate bounds (ref) are identical to the one-dimensional generalized Lee bounds by semenova2023generalized, i.e. our theory covers this as a special case.

Estimation

We now focus on the particular moment functions with desirable properties. For any required $B \in \{L^+,U^+,L^-,U^-,S_0,S_1,S_0^+,S_1^-\}$, let the (uncentered) moment function $\psi_B(W_i,\eta_0)$ be defined according to Table (ref).

table[table omitted — 1,971 chars of source]

Note that all bounds in (ref) and (ref) are functions of various ratios of different conditional expectations $E[\psi_B(W_i,\eta_0)|Z=z]$. Thus, they can be obtained by separate regressions of the respective $\psi_B(W_i,\eta_0)$ onto the spaces spanned by their $k_B$-dimensional transformations of $Z_i$, $b_B(Z_i)$

align[align omitted — 109 chars of source]

with conditional mean errors $E[\varepsilon_{i,B}|Z_i] = 0$ and approximation errors

align[align omitted — 76 chars of source]

Based on this, we construct population estimands for the bounds $B=L,U$ using the linear predictors at each $Z=z$ as

align[align omitted — 177 chars of source]

and analogously for $\theta_{B}^{+,LP}$ and $\theta_{B}^{-,LP}$. Replacing the unobserved true scores in (ref) by their sample counterparts yields estimators

align[align omitted — 137 chars of source]

where the scores of the effect bounds with estimated nuisance quantities $\psi_B(W_i,\hat{\eta})$ serve as pseudo-outcomes in separate least squares regression on their respective basis functions.\footnote{In principle, the equations based on (ref) are only seemingly unrelated and could thus also be estimated using a system based approach that takes into account the correlation structure of the conditional mean errors to increase efficiency. While this is straightforward in the case of a finite-dimensional parametric mean functions, it introduces additional dependencies in the two-step estimation that might offset potential gains in efficiency in the nonparametric case. Thus, we leave an extension along this line for future work.} The estimated aggregate heterogeneous effect bounds can then be calculated by combining the point predictions of four models

align[align omitted — 194 chars of source]

and analogously for $\hat{\theta}_{B}^+(z)$ and $\hat{\theta}_{B}^-(z)$. In principle the two components in both numerator and denominator in (ref) can also be obtained joint regressions that use $(\psi_{B^+}(W_i,\hat{\eta}) + \psi_{B^-}(W_i,\hat{\eta}))$ and $(\psi_{S_0^+}(W_i,\hat{\eta}) + \psi_{S_0^-}(W_i,\hat{\eta}))$ respectively as outcomes in (ref). This requires overall fewer parameters and only two choices of basis functions which could be beneficial in finite samples but is overall less adaptive to different complexity and trade-offs in estimating the relevant conditional expectations.

For the remainder, let $\theta(z) \in \{\theta_{AT}(z),\theta_{AT}^+(z),\theta_{AT}^-(z)\}$ denote a generic parameter of interest. The total sum of basis functions over all required regressions $k^*$ here will depend on this target parameter, see Table (ref) for an overview using a separate basis for each component.

table[table omitted — 667 chars of source]

When all $b_B(Z_i)$ for the relevant $B$ consist only of a constant, the estimator is similar to one for the unconditional generalized Lee bounds in semenova2023generalized with an augmented inverse probability weighting (AIPW) denominator.\footnote{This is possible due to the stronger margin assumption 3.3 compared to semenova2023generalized. Under this strong margin condition, $\rho_N$ in semenova2023generalized, Equation (5.7) can be set to $0$.}

Estimation of nuisance parameters $\hat{\eta}$ can be done via modern machine learning such as random forests, deep neural networks, high-dimensional sparse likelihood and regression models, or other non- and semiparametric estimation methods with good approximation qualities for the nuisance functions at hand. The influence of their learning bias/approximation on the functional bound estimator is limited by the Neyman-orthogonality of the chosen moment functions chernozhukov2018double,semenova2023generalized. In particular, under typical basis choices and suitable regularity conditions, the estimation does not affect the large sample distribution if the $L_2$-approximation rates for all nuisance quantities $\hat{\eta} - \eta_0$ are of order $o((nk_B)^{-1/4})$ and the selection probabilities are $L_1$-consistent at rate $o({k_B}^{-1/2})$. The first condition is identical to standard DML estimation of conditional average treatment effects using Neyman-orthogonal moment functions for nonparametric projection semenova2021debiased. When $z$ is one-dimensional, popular $k_B$ choices for many bases under weak smoothness assumptions are of rate $O(n^{1/5})$ leading to an overall RMSE convergence requirement for the nuisances of $o(n^{-3/10})$. This is a rate achievable by many nonparametric and machine learning estimators such as forests, deep neural networks, or high-dimensional sparse models under moderate complexity and/or dimensionality restrictions, see semenova2023generalized, semenova2021debiased, and heiler2021effect for examples. The additional weak $L_1$ condition is required to control the variance of the trimming indicators. More details regarding the technical assumptions can be found in Section (ref).

We also require that all components in $\hat{\eta}$ are obtained via $K$-fold cross-fitting, see Definition 3.1 in chernozhukov2018double. The use of cross-fitting controls potential bias arising from over-fitting using flexible machine learning methods without the need to evaluate complexity/entropy conditions for the function class that contains true and estimated nuisance quantities with high probability. If finite dimensional parametric models such as linear or logistic are assumed and estimated for the nuisance quantities, the proposed methodology can be applied without the need for cross-fitting.

Under suitable assumptions, the heterogeneous bound estimators are jointly asymptotically normal at each $z \in \mathcal{Z}$

align[align omitted — 508 chars of source]

For the complete definitions of the variance terms consider Appendix (ref). The variance can depend on the sample size and is generally increasing in norm if we allow the basis functions $b_B(z)$ to grow with the sample size. Thus, convergence is slower than the parametric rate equivalently to conventional nonparametric series regression belloni2015some. If the approximation errors $r_{B}(z)$ are rather large, then this distributional result is still valid if centered around the linear predictors of the true bounds (ref) similar to standard regression estimation under misspecification.

Inference

Inverting the quantiles of the distribution in (ref) could in principle be used to construct confidence intervals for the upper and lower bounds. However, this would be overly conservative for the actual effect of interest $\theta(z)$. Inference on this partially identified parameter should be adaptive to the underlying true width of the interval. In addition, if we allow for (local) misspecification, corresponding confidence regions could be empty or very narrow suggesting overly precise inference andrews2019inference. Robustness to misspecification is important in our setup as, in contrast to unconditional bounds, the heterogeneous bound functions are estimated at potentially many points with different variances and varying strength of identification in the sense of different widths of the underlying true identified set. A flexible estimator that is chosen e.g. by a global goodness-of-fit criterion for the effect bound curves could well be locally misspecified at some points. A limiting case of this type of misspecification would be to “overfit” the effect of a treatment that is fully independent of selection for some units instead of imposing local point identification. Corresponding confidence regions should be adaptive to such cases. To do so, we introduce the notion of a pseudo-true parameter $\theta^*(z)$ and its corresponding standard deviation $\sigma^*(z)$ that can be estimated as

align[align omitted — 269 chars of source]

This is a variance weighted version of the upper and lower bound. The pointwise $(1-\alpha)\%$-confidence intervals for the true always-taker effect $\theta(z)$ can then be obtained by the union of two intervals

align[align omitted — 289 chars of source]

where the critical value $\hat{c}(z)$ uniquely solves

align[align omitted — 189 chars of source]

with $(u_1,u_2)$ being jointly normal with unit variances and covariance $\hat{\rho}(z)$. This is an adaptation of the method by stoye2020simple who considers simple parametric estimators for generic bounds without two-stage estimation, cross-fitting, or additional nuisance functions. Interval (ref) is robust against misspecification that yields reverse ordering of the bounds as it is never empty and guarantees at least nominal coverage over an extended parameters space. In particular, it has at least $(1-\alpha)\%$ asymptotic coverage uniformly for all widths $\theta_U(z) - \theta_L(z)$ pointwise at each $z\in\mathcal{Z}$. For uniform confidence bands consider Section (ref).

In principle, population moments (ref) should not cross or be arbitrarily close to point identification under Assumptions 3.1 to 3.4. The same applies to the nonparametric projections at any $z$ under correct specification. However, when estimating the unknown conditional probabilities and quantile functions as well as the final heterogeneity projections, (local) misspecification is increasingly likely. Thus, this problem could also be analyzed from a model selection perspective. Our inference approach is agnostic whether narrow or empty identified sets are a result of a violation of the identification assumptions or of such model selection based misspecification. However, the pseudo-true parameter should be interpreted from that perspective. Note that relaxing the bounds pointwise could potentially be crude. Optimally, an inference procedure should incorporate common information across all bounds. Our misspecification robust intervals do this indirectly: When the density of $z$ and conditional moments are smooth, bounds and confidence intervals will be smooth as well. Thus, there is an implied smoothness on the relaxation parameter of the underlying conditional moment inequalities in the sense of andrews2019inference. For more details regarding the misspecification framework, consider Appendix (ref).

Large Sample Theory

Assumptions and Pointwise Limiting Distribution

In this section we provide the assumptions for asymptotic normality, validity of the confidence intervals (ref), and some more technical discussion. We again present the case with known propensity scores.\footnote{With propensity scores estimated, one has to augment the bias-correction in the moment functions and the nuisance parameter space as in semenova2023generalized as well as Assumption A.6 by additional terms that equivalently depend on the (squared) $L_p$ error rate of the propensity score estimator.} Denote $||\cdot||_p$ as the $L_p$ norm. Let $E[\psi_B(W_i,\eta_0)|Z_i=z] \in \mathcal{G}_B$ where $\mathcal{G}_B$ is a space of functions (possibly depending on $n$) that map from $\mathcal{Z}$ to the real line. Note that $E[\psi_B(W_i,\eta_0)|Z_i=z] = b_B(z)'\beta_{B,0} + r_{B}(z)$ with basis transformations $b_B(z) \in \mathcal{S}^{k_B} := \{b \in \mathbb{R}^{k_B}: ||b|| = 1\}$ and $\beta_{B,0}$ being the parameter of the linear predictor defined as root of equation $E[b_B(Z_i)(\psi_B(W_i,\eta_0) - b_B(Z_i)'\beta_{B,0})] = 0$. Define basis bound $\xi_{k,B} = \sup_{z\in\mathcal{Z}}||b_{B}(z)||$. Let $\eta \in T$ where $T$ is a convex subset of some normed vector space. Denote the realization set $\mathcal{T}_n = (\mathcal{S}_{0,n}\times\mathcal{S}_{1,n}\times \mathcal{Q}_n) \subset T$ as the set that with high probability contains estimators $\hat{\eta} = \{\hat{s}(0,X_i), \hat{s}(1,X_i),\hat{q}_0(u,X_i),\hat{q}_1(u,X_i)\}$ for nuisance quantities $\eta_0 = \{s(0,X_i),s(1,X_i),q_0(u,X_i),q_1(u,X_i)\}$. Let their corresponding $L_p$ error rates be

align*[align* omitted — 254 chars of source]

where $\tilde{U}$ is a compact subset of $(0,1)$ containing the relevant quantile trimming threshold support unions $([supp(p_0(X_i)) \cup supp(1-p_0(X_i))] \cap \mathcal{X}_{+}) \cup ([supp(1/p_0(X_i)) \cup supp(1-1/p_0(X_i))] \cap \mathcal{X}_{-})$. All of the following assumptions are uniformly over $n$ if not stated differently for the required elements $B \in \{L^+,U^+,L^-,U^-,S_0,S_1,S_0^+,S_1^-\}$:

enumerate[itemsep=0pt] • (Identification) $Q_B = E[b_B(Z_i)b_B(Z_i)']$ has eigenvalues bounded above and away from zero. • (Regular outcome) The outcome has bounded conditional moments $E[Y_i^m|X_i=x,D_i=d,S_i=1]$ for some $m > 2$ and a continuous density $f(y|X=x,D=d,S_i=1)$ that is bounded from above and away from zero with bounded first derivative for any $x\in \mathcal{X}$ and $d\in \{0,1\}$. • (Strong multiple overlap) There exist constants $\underline{e},\underline{s} \in (0,1/2)$ such that \begin{align*} e < \inf_{x\in\mathcal{X}} e(x) \leq \sup_{x\in\mathcal{X}} e(x) < 1- e, \end{align*}\begin{align*} s < \inf_{x\in\mathcal{X},d\in{\{0,1\}}} s(d,x) \leq \sup_{x\in\mathcal{X},d\in{\{0,1\}}} s(d,x) < 1- s. \end{align*} • (Approximation) For any $n$ and $k_B$, there are finite constants $c_{k,B}$ and $l_{k,B}$ such that for each $E[\psi_B(W_i,\eta_0)|Z_i=z] \in \mathcal{G}_B$ \begin{align*} ||r_{B}||_{P,2} &:= \sqrt{\int_{z\in\mathcal{Z}}r_{B}^2(z)dP(z)} \leq c_{k,B}, \\ ||r_{B}||_{P,\infty} &:= \sup_{z\in\mathcal{Z}}|r_{B}(z)| \leq l_{k,B}c_{k,B}. \end{align*} • (Basis growth) Let $\sqrt{n}/\xi_{k,B} - l_{k,B}c_{k,B} \rightarrow \infty$ such that \begin{align*} \sqrt{\frac{\xi_{k,B}^2\log k_B}{n}}\bigg(1 + \sqrt{k_B}l_{k,B}c_{k,B} \bigg)&= o(1). \end{align*} • (Machine learning bias) Let $e_n = o(1)$. For all folds, the nuisance parameters obtained via cross-fitting belong to a shrinking neighborhood $\mathcal{T}_n$ around $\eta_0$ with probability of at least $1-e_n$, such that \begin{align*} \xi_{k,B}(\lambda_{q,n,1} + \lambda_{s,n,1} +\lambda_{q,n,2} + \lambda_{s,n,2}) &= o(1) \end{align*} and (i) either \begin{align*} \sqrt{nk_B}(\lambda_{q,n,4}^2 + \lambda_{s,n,4}^2) &= o(1) \end{align*} or (ii) the basis is bounded $\sup_{z\in\mathcal{Z}}||b(z)||_{\infty} < C$ and \begin{align*} \sqrt{nk_B}(\lambda_{q,n,2}^2 + \lambda_{s,n,2}^2) &= o(1). \end{align*} For $B\neq S_0,S_1$ under conditional monotonicity, we also require that \begin{align*} \sqrt{n\xi_{k_B}^2}\lambda_{s,n,2}^2 &= o(1). \end{align*}

}

Assumption A.1 excludes collinearity of the basis transformations of the heterogeneity variables. A.2 puts restrictions on the tails and the smoothness of the distribution for the observed outcome distribution in different selection and treatment states. A.3 assures that there are comparable units between units of different selection and/or treatment status. A.2 and A.3 together imply a continuously differentiable conditional quantile function for the selected observed units that is almost surely bounded. This, together with the strong overlap for the treatment and selection probabilities, assures that the effect bounds are regularly identified khan2010irregular,HEILER2021valid.

A.4 defines $L_2$ and uniform approximation error bounds for function class $\mathcal{G}_B$. This is a typical characterization in the literature on nonparametric series regression without nuisance functions belloni2015some. We say the model is correctly specified if the basis is sufficiently rich to span $\mathcal{G}_B$, i.e. $c_{k,B} \rightarrow 0$ as $k_B \rightarrow \infty$. However, the distributional theory also allows for the case of misspecification, i.e. $c_{k,B} \not\rightarrow 0$. A.5 controls the approximation error from linearization of the estimator with unknown design matrix $Q_B$. This is equivalent to the condition required for localization in general least squares series regression belloni2015some.\footnote{For more specific series methods such as splines huang2003local or local partitioning estimators cattaneo2020large, this rate can be improved to $\sqrt{\xi_{k,B}^2 \log k_B/n}(1+ \sqrt{\log k_B}l_{k,B}c_{k,B})$, see also belloni2015some, Section 4 and cattaneo2020large, Supplemental Appendix Remark SA-4 for a related discussion.}

A.6 is crucial: It says that the machine learning estimators for the conditional quantiles and selection probabilities have sufficiently good approximation qualities in an $L_p$ sense. In the case of a bounded finite basis, the conditions reduce to $L_1$ and $L_2$ consistency as well as the well-known requirement in the semiparametric/DML literature that the nuisance functions have root have mean squared error rates of order $o(n^{-1/4})$ chernozhukov2018double. In the more general case, the higher $L_p$ rates are identical to the ones in semenova2021debiased required for nonparametric estimation of conditional average treatment effects. The additional condition for conditional monotonicity is due to potential classification error. In the typical case of well-behaved basis functions such splines, wavelets, and local partitioning, we have that $\xi_{k,B} \lesssim \sqrt{k_B}$ belloni2015some and thus this requirement is subsumed by the prior $\sqrt{nk_B}\lambda_{s,n,p}^2$ condition. When $B \in \{S_0,S_1\}$, the $L_1$ requirement is omitted as the AIPW moment functions do not contain any trimming indicator.

For the estimation of the asymptotic variance, we also assume that A.V holds:

enumerate[itemsep=0pt] • (Asymptotic variance) The conditions in Appendix (ref), Assumption (ref) hold, i.e.\\ $||\hat{\Omega}_n(z) - \Omega(z)|| = o_p(1)$ pointwise at each $z \in \mathcal{Z}$.

} A.V can require somewhat stronger outcome moment/tail and basis growth conditions. The corresponding primitive assumptions and discussion can be found in Appendix (ref) and are omitted for brevity. We obtain the following Theorem:

thm[Asymptotic Normality] Suppose Assumptions 3.1 - 3.4 and A.1 - A.6 hold and $\hat{\theta}_B(z_0)$ for $B = L, U$ and $\hat{\Omega}_n(z_0)$ are estimators according to (ref) and (ref) respectively. Let $\theta_B^{LP}(z)$ be the the population predictor according to (ref). Then, for any sequence $z_0 = z_{0,n}$, \begin{align*} \sqrt{n}&\hat{\Omega}_n^{-\frac{1}{2}}(z_0)\begin{pmatrix} \hat{\theta}_L(z_0) - \theta_L^{LP}(z_0)\ \\ \hat{\theta}_U(z_0) - \theta_U^{LP}(z_0) \end{pmatrix} \overset{d}{\rightarrow} \mathcal{N}\left(\begin{pmatrix} 0\\0 \end{pmatrix}, \begin{pmatrix} 1 &0 \\ 0 &1 \end{pmatrix}\right). \end{align*} Moreover if $\sup_{B}n^{1/2}k_B^{-1/2}l_{k,B}c_{k,B} = o(1)$, then \begin{align*} \sqrt{n}&\hat{\Omega}_n^{-\frac{1}{2}}(z_0)\begin{pmatrix} \hat{\theta}_L(z_0) - \theta_L(z_0) \\ \hat{\theta}_U(z_0) - \theta_U(z_0) \end{pmatrix} \overset{d}{\rightarrow} \mathcal{N}\left(\begin{pmatrix} 0\\0 \end{pmatrix}, \begin{pmatrix} 1 &0 \\ 0 &1 \end{pmatrix}\right). \end{align*}

Theorem (ref) shows that nonparametric heterogeneous bounds using DML are jointly asymptotically normal. It allows for the case of misspecification when centered around the linear predictor. It is most useful under the additional undersmoothing condition that makes any misspecification bias vanish sufficiently fast.\footnote{For example, when $\mathcal{G}_B$ is in a $s$-dimensional ball on $\mathcal{Z}$ of finite diameter, then the condition simplifies to $n^{1/2}k_B^{-(\frac{1}{2} + \frac{s}{d})}\log(k_B) \rightarrow 0$. See belloni2015some, Comment 4.3 for additional details. Note that undersmoothing does in general not admit mean-squared error optimal $k_B$ choices, see e.g. cattaneo2020large for multiple bias-correction alternatives methods for local partitioning estimators.}

Pseudo Parameterization and Pointwise Inference

Theorem (ref) in principle is sufficient to construct confidence intervals for the heterogeneous effect bounds. They are, however, too wide for the actual effect parameter of interest $\theta(z)$ depending on the width of $\theta_U(z) - \theta_L(z)$ imbens2004confidence. If the difference between upper and lower bound is large, the testing problem is essentially one-sided compared to the case of only having a small difference. This raises the question of how to conduct inference that is uniformly valid with respect to the underlying difference between upper and lower effect bound and has more power compared to using simple two-sided critical values. Moreover, there are multiple sources of misspecification that could invalidate or suggest spuriously precise inference when upper and lower bound estimates are close to each other or reverted. In particular, even when the population bounds do not intersect, their linear predictors still might. Handling such misspecification is particularly relevant for heterogeneous bounds as here we estimate two functions at potentially many points which makes (local) misspecification and/or reordering a much larger concern compared to a single partially identified parameter as in lee2009training or semenova2023generalized.\footnote{Note also that we allow for the use of different basis functions across when estimating upper and lower effect bounds. This seems reasonable as it could well be that lower and upper effect bounds are curves of e.g. different smoothness and not equally difficult to approximate. Alternatively one could impose the largest complexity of any bound for both models, i.e. paying the price of a potentially higher estimation variance in the effect bound with a higher degree of smoothness. While this will generally yield a better approximation for each separate bound, it could also contribute to a reverse ordering in the sense that $\hat{\theta}_L(z) \geq \hat{\theta}_U(z)$ at some $z$. } We relax the notion of coverage following andrews2019inference to include an artificial pseudo-true target parameter that guarantees non-empty confidence intervals. For more details on misspecification and corresponding formalization consider Appendix (ref).

Our method is based on stoye2020simple who demonstrates that, in the simple case of two regular, jointly asymptotically normally distributed parameters, uniformly valid intervals with good power properties can be obtained by concentrating out the unobserved true difference in effect bounds. We adapt his framework to the nonparametric heterogeneous effect bounds with estimated nuisances. For any $z \in \mathcal{Z}$ let the identified set be $\Theta_{z} = [\theta_L(z),\theta_U(z)]$ Now denote the pseudo-true parameter $\theta^*(z)$ and its variance as

align[align omitted — 245 chars of source]

The pseudo-true identified set is then given by $\Theta_{z}^* = \Theta_z \cup \{\theta^*(z)\}$. This set is implicitly defined by an estimand corresponding to setting a test statistic for the heterogeneous effect at $z$ to (i) zero (meaning no rejection of the null) if the interval is nonempty and the null value is inside the estimated intervals and to (ii) the larger of the two $t$-statistics for upper and lower bound for the null of the bound being equal to the hypothesized value. In particular, it chooses the test statistic under the null as $\max\{(\theta(z)-\hat{\theta}_U(z))/\hat{\sigma}_U(z),(\hat{\theta}_L(z)-\theta(z))/\hat{\sigma}_L(z),0\}$. Case (ii) collapses to a standard test for the parameter lying in $\Theta_z$ in the case of a well-defined interval with $\theta_U(z) > \theta_L(z)$ and to a test on the pseudo-true parameter under misspecification $\theta_U(z) < \theta_L(z)$. This inference procedure can also be interpreted as resulting from a moment inequality problem that allows for misspecification by adding slacks/deviations to equations defining upper and lower effect bounds, see Appendix (ref). If such slacks are large, $\theta_U(z) < \theta_L(z)$ can be admitted as solution. In principle, alternative definitions for the pseudo-true set $\Theta^*_z$ that use different weighting compared to (ref) could be considered as well. This would change the definition of the pseudo-true parameter. Not all choices of such pseudo-true parameter are of equal use. This is similar to e.g. GMM models under misspecification where the pseudo-true parameter is defined implicitly as the maximizer of a GMM population criterion using a particular weighting matrix. Therefore, the choice of the pseudo-true parameter and its usefulness should be balanced versus the robustness against spurious precision under misspecification andrews2019inference. The particular choice in (ref) leads to a convenient solution in terms of critical value adjustment due to the specific asymptotic bivariate normal distribution stoye2020simple. We obtain the following Theorem:

thm[Misspecification Robust Inference] Suppose Assumptions 3.1 - 3.4 and A.1 - A.6 hold and $\hat{\theta}_B(z_0)$ for $B = L, U$ and $\hat{\Omega}_n(z_0)$ are estimated according to (ref) and (ref) respectively. Let the confidence interval be constructed according to (ref). Then, for any sequence $z_0 = z_{0,n}$, \begin{align*} \liminf_{n\rightarrow\infty}\inf_{{\theta}^{LP}(z_0)\in{\Theta}^*_{z_0}}\ P({\theta}^{LP}(z_0) \in CI_{{\theta(z_0)},1-\alpha}) \geq 1 - \alpha, \end{align*} where ${\Theta}^*_{z_0} = [{\theta}_L^{LP}(z_0),{\theta}_U^{LP}(z_0)] \cup {\theta}^{LP,*}(z_0)$ where ${\theta}^{LP,*}(z)$ is defined analogously to (ref) with linear predictors instead of bounds. Furthermore if $\sup_{B}n^{1/2}k_B^{-1/2}l_{k,B}c_{k,B} = o(1)$, then \begin{align*} \liminf_{n\rightarrow\infty}\inf_{\theta(z_0)\in{\Theta}^*_{z_0}}\ P(\theta(z_0) \in CI_{\theta(z_0),1-\alpha}) \geq 1 - \alpha, \end{align*} where ${\Theta}^*_{z_0} = [\theta_L(z_0),\theta_U(z_0)] \cup {\theta^*}(z_0)$ as defined in (ref).

Theorem (ref) demonstrates the asymptotic validity of the confidence intervals proposed in (ref) for both the linear predictor as well as the correctly specified case. In particular, we achieve at least nominal coverage independently of the actual width of the true identified region for the heterogeneous effect bound. Coverage is uniform with respect to the width of the true identified region $\theta_U(z) - \theta_L(z)$ pointwise at each $z \in \mathcal{Z}$ or along sequences therein. The coverage notion over the augmented parameter set $\Theta_{z}^*$ assures non-emptiness of the interval and avoids spuriously precise inference in regions where the estimated lower bound might be too large relative to the estimated upper bound. Coverage will be closer to the nominal level when correlations between conditional mean errors are small. In particular, near one-sided critical values apply for the usual levels of confidence when $\rho(z) = 0$. Note that $\rho(z)$ stems from the correlation between two linear combinations of residuals from up to four separate nonparametric regressions. Hence, $\rho(z)$ can generally vary over different values of $z$. Thus, power and size properties of the corresponding tests will depend on the location of the local effect bound. In particular, they are driven by the share of missing outcomes at given $z$. In the case of point identification (no missing outcomes), upper and lower bounds and pseudo true parameter are equivalent and $\rho(z) = 1$. Thus, the confidence interval collapses to one using standard two-sided critical values.\footnote{Alternatively, one could employ an additional data splitting step for estimating the two bounds on different subsamples to assure that $\rho(z)$ is equal to zero. We refrain from this approach to avoid inaccuracies in finite samples as the first estimation step already requires cross-fitting and potential data splitting for tuning of the machine learning methods within folds.}

Uniform Inference

We now consider inference for the whole process $\{\theta(z)\}_{z\in \mathcal{Z}}$ without misspecification. Under stronger approximation requirements, the series approach allows for a uniform linearization that can be used to obtain a process analogue of Theorem (ref). We exploit this to construct a simple multiplier bootstrap procedure that guarantees uniform coverage and avoids retraining the first-stage estimators. Inference is then based on a joint process results for the bounds. This does not exploit locally varying widths of the identified sets and convergence in distribution. As a result, inference will be generally be more conservative. The bootstrap works as follows: For each $b=1,2,\dots$, generate $n$ independent standard exponential random variables $h_1,\dots,h_n$. Then, for the relevant $B \in \{L^+,U^+,L^-,U^-,S_0,S_1,S_0^+,S_1^-\}$, use the same $h_i$ in multiple weighted regressions

align[align omitted — 155 chars of source]

with estimated nuisance parameters to obtain bootstrap estimates

align[align omitted — 149 chars of source]

Then, for all $z$, calculate the centered bootstrap $t$-statistic for bounds $B = L,U$ as

align[align omitted — 93 chars of source]

Denote $\bar{t}_B^b = \sup_{z\in\mathcal{Z}} t_B^b(z)$ and $\underline{t}_B^b = \inf_{z\in\mathcal{Z}} t_B^b(z)$ and, for any $\alpha$, let $c_{n,\alpha}(t^b)$ be the $\alpha$-quantile of $t^b$ over the bootstrap samples $b=1,2,\dots$. The $(1-\alpha)$ confidence band for $\theta(z)$ for all $z \in \mathcal{Z}$ is then given by

align[align omitted — 243 chars of source]

Theorem (ref) contains the joint strong approximation result for both series and bootstrap processes. Theorem (ref) contains the uniform inference result. The assumptions are generally stronger with regards to the basis function growth and existence of moments compared to the pointwise case belloni2015some.

thm[Strong Approximation of Series and Bootstrap Process] Let $k^*$ be defined according to Table (ref). Under the assumptions of Theorem (ref) with (i) $m > 3$ and (ii) $\sup_B \sqrt{{\xi_{k,B}^2\log^3n}/{n}}(n^{1/m_B}\log^{1/2}n + \sqrt{k_B}l_{k,B}c_{k,B}) = o(a_n^{-1})$, (iii) $a_n^3n^{-1/2}\sup_Bk_B^2$ $\xi_{k,B}(\log n + \log^2 n /k_B) = o(1)$, and (iv) $\sup_{B}n^{1/2}$ $k_{B}l_{k,B}c_{k,B}\log^2n = o(a_n^{-1})$, there exists a random standard normal vector $\mathcal{N}_{k^*}$ of length $k^*$ such that \begin{align*} \sqrt{n}{\hat{\Omega}_n(z)^{-1/2}}\begin{pmatrix} {\hat{\theta}_L(z) - \theta_L(z)} \\ {\hat{\theta}_U(z) - \theta_U(z)} \end{pmatrix} =_d {\Omega(z)^{1/2}}\mathcal{N}_{k^*} + o_P(a_n^{-1}) in \ell^{\infty}(\mathcal{Z}) \end{align*} and \begin{align*} \sqrt{n}{\hat{\Omega}^b_n(z)^{-1/2}}\begin{pmatrix} {\hat{\theta}^b_L(z) - \hat{\theta}_L(z)} \\ {\hat{\theta}^b_U(z) - \hat{\theta}_U(z)} \end{pmatrix} =_d {\Omega(z)^{1/2}}\mathcal{N}_{k^*} + o_{P^*}(a_n^{-1}) in \ell^{\infty}(\mathcal{Z}), \end{align*} where $P^*$ is the conditional probability computed given the data $\{W_i\}_{i=1}^n$.
thm[Uniform Inference] Denote $\mathcal{P}$ as the set of probability measures obeying the assumptions of Theorem (ref). The confidence bands (ref) have uniform coverage \begin{align*} \underset{n\rightarrow\infty}{\lim \inf} \inf_{P\in\mathcal{P}} P(\theta(z) \in CB_{\theta(z),1-\alpha} for all z\in\mathcal{Z}) \geq 1 - \alpha. \end{align*}

Monte Carlo Study

In this section, we analyze size and power properties of the proposed misspecification robust confidence intervals in finite samples. In particular we look at the size along a grid of $z$-values that vary with respect to the share of always-takers and thus the width of the identified set in the population. We consider a generalized Roy model with a random binary treatment and missing responses:

table[table omitted — 1,170 chars of source]

Similar designs have been considered for simulation for point-identified effects under exclusion restrictions, see e.g. heckman2005structural or heiler2022efficient. Note that the true model here assumes strong monotonicity of the effect of the treatment on response. This knowledge, however, is not imposed onto the estimation procedures. We chose $Z_i = X_{i,1}$. Thus, the true response function for the always-takers is given by

align[align omitted — 195 chars of source]

The parameters are chosen such that non-response rates vary from 84.13% to 99.00% and are monotonically increasing in $z$. This will allow us to discover potential differences in coverage rates for varying setups including scenarios close to point identification. We consider both the continuous case, i.e. the size of confidence intervals for the continuous $\theta(z)$ as well as power curves for a simple discretized version where $z$ is integrated from 0 to 0.5 and 0.5 to 1 respectively. The nuisance quantities are estimated using honest probability random forests and honest quantile regression forests from the \verb|grf| package athey2019generalized with default tuning parameters and two-fold cross-fitting. The design is sufficiently sparse for the forests to achieve the required convergence rates for estimating (ref) according to Assumption A.6. For the heterogeneity analysis, we use basis splines with nodes and order selected via leave-one-out cross-validation for the continuous case and indicator functions for the discrete case. We analyze the size and power of the misspecification robust confidence intervals (ref) at a 95% confidence level.

figure[figure omitted — 703 chars of source]

Figure (ref) depicts the simulated coverage rates in the case of $p=10$ and $p=100$ regressors for total sample sizes of $n=400$ and $n=2000$. The rates for $n=400$ can sometimes drop to around 90% but overall coverage is still close to nominal. For $n=2000$, the confidence intervals have at least nominal coverage ranging from 95% to 100% depending on $z$. Thus, the theoretical large sample guarantee in Theorem (ref) seems to approximate the finite sample behavior reasonably well in these designs.

figure[figure omitted — 1,373 chars of source]

Figure (ref) depicts the power curves for two different heterogeneous effect parameters $\theta(z_0)$ and $\theta(z_1)$ for $p=10$ and $p=100$ at sample sizes $n=400$ and $n=2000$. $\theta(z_0)$ corresponds to an area with larger uncertainty regarding the partially identified parameter, i.e. it is integrated over the range of the heterogeneity with the largest share of unobserved outcomes while $\theta(z_1)$ integrates over the range with the largest share of observed outcomes. The difference in uncertainty can also be seen in Figure (ref). The results show that the power curves are close to zero around the null for all sample size and designs reflecting the conservativeness of the inferential method. However, power converges quickly to 100% when moving away from the null. Power is lower for smaller sample sizes and a larger amount of possible confounding variables as expected. Moreover, in the $z=z_0$ case for which the share of missing outcome is larger, confidence intervals have lower power compared to the $z=z_1$ case that is closer to the case of point identification. Overall, intervals seem to perform reasonably well for different degrees of identification in the overall heterogeneous effect but being closer to point identification tends to yield more power in finite samples.

Empirical Study

Literature on Social Media and Political Polarization

There is a large developing literature regarding the effects of social media on political polarization. However, the evidence on the effects of social media and news exposure via social media on polarization is very mixed, see haidtONGOINGsocialmediaPolitical (ongoing) for a comprehensive collaborative review. Effects and channels are often context-dependent and can vary by time, country, platform, algorithm, research design, subgroup, outcome measure, and more. Here we focus on research regarding affective polarization, Facebook, online news, and questions of selection and heterogeneity.

There are significant changes in measures of affective polarization in many countries since the 1980's with steepest increase in the US boxell2020cross. The precise role of online news consumption and social media is still under debate: suhay2018polarizing show that online news that contains partisan criticism that derogates political opponents increases affective polarization. munger2020null use different Amazon Turk and Facebook Ad samples with clickbait and conventional news headline treatments. They find that older non-democrat users read clickbait articles more often but no effects on polarization. cho2020search demonstrate that political videos selected via the YouTube algorithm increase affective polarization. nordbrandt2021affective uses Dutch Survey data to argue that higher levels of affective polarization lead to increased social media consumption but finds no evidence for the reverse channel. However, there seems to be significant user-level heterogeneity. beam2018facebook show that, over the course of the US 2016 presidential election campaign, attitudes towards political opponents remained relatively stable. Moreover, they find that Facebook leads to modest depolarization due to an increase in counter-attitudinal news exposure in this period. However, they look at partisan measures based on identity formation and not at standard affective polarization metrics. bail2018exposure study the effects of following counter-attitudinal Twitter bots. Their findings suggest a heterogeneous increase in polarization from counter-attitudinal exposure on Twitter. In particular, Republicans adopted more conservative attitudes after being exposed to a liberal bot while for Democrats the increase in liberal attitudes after following a conservative bot is insignificant. allcott2020welfare show that deactivation of Facebook for one month in Fall 2018 decreased exposure to polarizing news and polarization of political views. Their point estimate on affective polarization is negative but insignificant. However, the study is underpowered to detect small effects and they use a consumption based metric of polarization instead of attitudinal measures. feezell2021exploring do not find evidence that algorithmic or non-algorithmic news sources contribute to higher levels of partisan polarization using explorative survey data. di2021does study varying social media status treatments for inside and outside echo chamber units. They find that Twitter, in particular when allowing for interactions, increases polarization for groups which are already classified as being inside an echo chamber. Subjects in the outside echo chamber group show no significant increases. yarchi2021political argue that there are important cross-platform differences. They provide evidence that Facebook is the least homophilic social media network in terms of interactions, positions, and expressed emotions.

levy2021social employs a large field experiment regarding news consumption, polarization, and algorithmic news selection on Facebook. He shows that a counter-attitudinal nudge towards subscribing to an outlet with an opposing political ideology decreases affective polarization but does not affect political opinions. He argues that the small nudge setup reflects a realistic user experience on Facebook. A channel seems to be that a shock to the selection of news consumption has lasting effects as units do not re-optimize their feed much afterwards. However, levy2021social also provides evidence that the Facebook algorithm limits exposure to precisely these counter-attitudinal news outlets and thus can increase polarization overall. This study has a clear stratified randomization design but suffers from large differential attrition rates in the endline survey, in particular when considering pro- and counter-attitudinal treatments separately. Conventional unconditional Lee bounds for the pro- and counter-attitudinal treatment effects are relatively wide and include zero as well as moderate polarization effects. levy2021social does not find significant heterogeneity in treatment effects when using simple interacted linear regression models that ignore attrition. We re-analyze the specific question of the effect of a counter-attitudinal nudge on Facebook on affective polarization by levy2021social using the refined DML based heterogeneous bounding method developed in this paper.

Experiment, Data, and Attrition

In levy2021social, users were recruited via Facebook ads and filled out a baseline survey between February--March 2018. Units were then stratified by self-reported political ideology and randomly allocated into one of three different treatment arms 1) Liberal, 2) Conservative, or 3) Control. The treatments in 1) and 2) consisted of a nudge to subscribe to (“like”) a selection of four potential (liberal or conservative) outlets. It was explained that a subscription could provide new perspectives, but there were no other incentives or rewards offered. Likes on Facebook make posts from the corresponding outlet more likely to appear on the user's feed, thus exposes them to potentially new information and opinions. The liberal outlets were HuffPost, MSNBC, The New York Times, and Slate. The conservative outlets were Fox News, The National Review, \textit{The Wall Street Journal}, and \textit{The Washington Times}. Based on this, the \textit{counter-attitudinal treatment} is defined as nudge towards outlets with ideological leanings contrary to the leaning of the user.\footnote{Individual leanings are based on party affiliation. If units do not identify as Democrats or Republicans, it is according to self-reported ideology. If they neither identify as liberal nor conservative, support of the candidate in the 2016 elections is used. This excludes about 3% of the total sample that provide no information on leaning.}

Around two months after the baseline survey, participants were asked to fill out an endline survey where political opinions and measures of affective polarization were recorded. levy2021social does not find effects from any treatment on political opinions, Without controlling for attrition, there is a significant decrease in affective polarization for the counter-attitudinal treatment group.

Thus, we consider the index for affective polarization constructed by levy2021social as outcome in what follows. In particular, we analyze the causal effects of nudging users towards counter-attitudinal subscriptions on affective polarization. The parameters can also be interpreted as intent-to-treat effects of subscription. The outcome measure is standardized such that all coefficients are measured in terms of standard deviations in what follows.

The covariates collected via the baseline survey and Facebook contain information on political ideology, party affiliation, voting behavior, approval of President Trump, baseline polarization, news consumption, and socio-demographic variables such as age and gender. For more information and descriptive statistics consider levy2021social, Section II. The final sample (including missing endline survey units) consists of 24230 units of which 12126 are in the treatment and 12104 in the control group. levy2021social estimates the effects of the intervention on affective polarization using only the units which replied to the endline survey, i.e. for which the outcome variable is observed. He argues that the main estimates are likely to generalize beyond the selected population.

table[table omitted — 1,091 chars of source]

The experiment, however, suffers from large differential attrition rates: Table (ref) contains the differential attrition rates between treatment and control group stratified by political ideology. Attrition in the endline survey is large with rates between 46.23% to 65.35%. Moreover, there is significant heterogeneity when looking at the difference in attrition rates between treatment and control condition. In addition, the association between treatment and attrition is not homogeneous across all subgroups. This indicates that there are potential interactions between treatment and baseline characteristics that can lead to heterogeneous response rates. In this case, looking only at unconditional attrition rates between treatment and control group provides a distorted view on the potential bias introduced by ignoring the selection into response. We provide a more thorough analysis of the effect bounds for the counter-attitudinal treatment fully accounting for heterogeneity in treatment effects and attrition rates and potential effects of the treatment on selection into the endline survey. In particular, we re-analyze the unconditional effect bounds using the method suggested in this paper as well as heterogeneous bounds in terms of relevant pre-treatment characteristics.

Parameter, Estimation, and Inference Methods

Replication of unconditional point estimates and unconditional Lee bounds provided by levy2021social, Table A.12(b) use the same methods. We bound the (un)conditional average treatment effect(s) for the always-takers using the methods developed in this paper. Always-takers in this experiment are units that would be taking part in the endline survey regardless of whether they have been nudged or not. They are estimated to make up around 46.18% of the total study population.

We employ two versions of the DML based generalized Lee bounds that differ in terms of nuisance parameter models: The first (DML parametric) uses logistic regression with all confounding variables interacted with the treatment for the response selection probabilities and linear quantile regression with all confounding variables for the conditional quantile of the selected treated and selected controls. These models are more likely to be misspecified. The second (DML forest) uses probability forests and quantile forests with 1000 trees and honest splitting athey2019generalized. For both specifications, all categorical variables are coded as flexible dummy variables leaving us with 36 confounders in total. Conditional quantile trimming levels are rounded towards the closest value on a grid from $(0.01, 0.02,\dots,0.99)$. Cross-fitting is based on 10 folds. For the heterogeneity analysis we use the estimated signals $\psi_B(W_i,\hat{\eta})$ provided by the forest-based DML specification and basis splines for the continuous variables with node and order selection via leave-one-out cross-validation. Confidence intervals are based on the misspecification robust method (ref), confidence bands on the multiplier bootstrap (ref). We report both at 90% due to the convservativeness of the methods.

Results

Unconditional Effects

In this subsection, we provide the unconditional effect analysis of the counter-attitudinal nudge on affective polarization. Table (ref) contains the estimates by levy2021social, the naive Lee bounds assuming (strong) monotonicity as well as the two DML-based methods.\footnote{Table A.12(b) by levy2021social contains an error as his bounds $[-0.172, 0.060]$ (no CI provided) are calculated from $\theta_L(x) = E[Y_i|D_i=1,S_i=1, Y_i \leq q(p_0(X_i)),X_i=x] - E[Y_i|D_i=0,S_i=1, Y_i \leq q(p_0(X_i)),X_i=x]$ (and equivalently for $\theta_U(x)$) and not based on (ref). Moreover, in the specifications with controls, he assumes constant effect bounds and a linear form of the truncated mean which is heavily restrictive and generally does not identify the true bounds. Columns (4) and (5) produce the correctly calculated bounds under the more credible assumptions allowing for non-linearity and weak monotonicity.}

table[table omitted — 1,135 chars of source]

The conventional Lee bounds (3) as well as the DML bounds (4) and (5) contain a null effect in the estimated identified set. All contain the point estimates provided by levy2021social. Both DML bounds are much shorter than the Lee bounds, ruling out even small to moderate polarization effects. The applicability of the conventional Lee bounds (3) under strong monotonicity is questionable. In particular, our estimates for the conditional selection probabilities based on the specification for Column (5) suggest that the effect of the treatment on attrition is negative ($\hat{p}_0(x) > 1$) for 65.42% and positive ($\hat{p}_0(x) < 1)$ for 34.58% of the sample indicating violation of strong monotonicity. This is in line with the heterogeneous attrition rates observed in Table (ref).

Heterogeneous Effects

In this subsection, we analyze the effect bounds of the counter-attitudinal nudge on affective polarization as functions of heterogeneity variables. In particular, we look at political ideology (categorical) and age (continuous). We select political ideology from the baseline survey as it was used for block-randomization of the original experiment and could provide insight regarding potential asymmetries in the effect of counter-attitudinal nudges in terms of partisanship. Age can easily be used for targeting of such an intervention based on social media information only (no survey required) and has been shown to be an important determinant of aggregate affective polarization levels in the US phillips2022affective. These variables were also suggested by levy2021social for heterogeneity analysis.

figure[figure omitted — 605 chars of source]

Figure (ref) contains the heterogeneous bounds sorted by ideology plus $90\%$ confidence intervals. The identified sets for Liberal, Moderate, Conservative, Extremely Conservative and Haven't thought much do not include zero. However, even for these groups, the statistical uncertainty dominates and we cannot rule out a null effect with precision. This reflects the fact that the experiment is not powered enough to detect small effects within smaller subgroups with high probability.

figure[figure omitted — 695 chars of source]

Figure (ref) contains the estimated effect bounds and confidence intervals as functions of age.\footnote{Note that here we have excluded the units for which age information was missing ($n = 791$). Thus, the estimate for the lower bound differs slightly from the unconditional estimate in Table (ref), Column (5). Estimating the bounds separately for the omitted category yields interval $[0.096,~0.167]$ with $90\%$-CI $(-0.779,~1.007)$.} We can see that the identified set is widening monotonically in age reflecting a larger sample selection problem with older users. The identified set suggest a depolarization effect for 18-43 year old users ranging from $[-0.091,-0.041]$ to $[-0.087,-0.001]$. The estimates rule out moderate to large polarization effects for these younger users. The negative bounds, however, cannot statistically reject any positive effects due to relevant standard errors. As we know that inference is conservative, we interpret this as weak evidence in favor of a depolarization effect on young users with magnitudes close to the unconditional point estimate suggested by levy2021social that does not correct for attrition.

Discussion

Overall, the findings do not contradict the conclusion by levy2021social that a counter-attitudinal nudge can decrease affective polarization. There is weak evidence in favor of effect heterogeneity. Tight depolarization bounds can be obtained for users from various political ideologies but they are not statistically significant. Regarding age differences, the identified set excludes non-negative values for young users. This could be due to measurement or, if accurate, dose-related as very young users spend significantly more time online. Their overall activity is also likely contributing to lower attrition rates which leads to the tighter bounds reported. Alternatively or additionally, affective polarization is usually understood through the lens of social identity theory: People internalize partisan affiliation as part of their sense of self. The latter tends to be more malleable for younger people, in particular in their formative years phillips2022affective, which could explain larger effects.

To put the size of identified sets into perspective, we compare the estimates to experimental estimate by allcott2020welfare. The bounds for the youngest age levels suggest that, for these groups, the impact is around 0.4 to 0.9 times as large as the effect of deactivating Facebook for a whole month. As a limitation, note that we are comparing to an unconditional baseline from allcott2020welfare. The effect of deactivation could potentially be much larger (or smaller) for these groups as well. For further research, targeting interventions at particular young age and particular ideological groups could provide more definitive evidence.

Concluding Remarks

This paper provides a method for estimation and inference for bounds for heterogeneous treatment effects under sample selection. We make the general point that heterogeneity in partially identified problems requires special attention as both effect parameters as well as identified sets can be subject to heterogeneity. Exploiting the latter can yield more precise inference in empirical applications compared to crude, unconditional approaches. There are also multiple extensions possible: In many applications where the method could be useful, the i.i.d. assumption is overly restrictive. In particular, in social experiments units are often clustered within groups such as schools or regions. It would also be interesting to see under which conditions the methodology in this paper can be extended to more general moment inequality problems.

\addcontentsline{toc}{section}{References} {\setstretch{1} }