Difference-in-Differences with Interference

Difference-in-Differences with Interference

Introduction

According to the stable unit treatment value assumption (SUTVA), potential outcomes only depend on one's own treatment assignment. In many cases, SUTVA fails due to an unknown interference structure among neighbors. In the fields of environmental economics, urban economics, labor economics, criminal justice, and many other fields of social sciences, place-based policies often generate spillover effects. One example is minimum wage increase in Seattle studied by jardim2022boundary. Through the channels of competition in the regional labor market for workers and the possibility of relocation of businesses, they find that significant spillover effects on wages and hours are seen up to a 40-minute drive from Seattle city limits.

When spillover effects are of interest, one often needs to observe the entire population. For example, we can typically collect information about all counties in the United States. In the example above, jardim2022boundary use administrative employment records in the state of Washington. In this way, we do not need to deal with partial observation of other units' treatment status and can allow for diverse interference patterns based on access to the entire treatment assignment vector. If we take sampling from the superpopulation/infinite population approach literally, what we are estimating turns out to be the spillover effect in a researcher's sample with missing neighbors unless interactions are restricted within clusters of friends or household members and one randomly samples clusters.\footnote{See further explanation below equation (7) on page 537 in manski1993identification. }

In this paper, I study the “difference-in-differences" (DID) type estimators that allow interference from a finite population perspective, where inference is conditional on covariates and the whole population is observed. This approach is closest to the conditional inference discussed by abadie2014inference and jin2023tailored. Conditional treatment effect parameters have also been mentioned in abadie2002simple, imbens2004nonparametric, and balzer2015targeted. Recently, viviano2024policy adopts the same inference framework when studying optimal treatment allocation under network interference.

The conditional inference approach adopted here allows arbitrary spatial correlation and nonstationarity of covariates. Meanwhile, stochastic potential outcomes allow the possibility of modeling the conditional means of the outcome variables and incorporate more uncertainty into the inference framework. Consequently, the proposed estimators below are more robust to model specification with a straightforward causal interpretation. One could legitimately argue that researchers should not stick to a single inference framework. That said, many attribute variables containing locational information and neighborhood characteristics, such as landlocked status, are deemed non-stochastic for spatial data. I therefore consider the current approach a natural starting point for studying population interference/spillover effects.\footnote{The conditional inference framework only affects the definition of the parameter of interest and the inference later on. Estimation is not affected by whether the current conditional inference, design-based inference, or superpopulation framework is adopted.}

One challenge of incorporating interference is the modeling of spillover functions. Following the interference literature I introduce exposure mapping, which is a function that maps an assignment vector to an exposure value aronow2017estimating, manski2013identification. One leading example is specifying exposure as the average value of treatment statuses of neighbors within a distance $\bar{\rho}$ of a unit $i$. I start with the case where exposure mapping is assumed to be known and correctly specified. Later on, I generalize the analysis to the case where exposure mapping could be misspecified. Regardless of the correct specification of the exposure mapping, the same estimator and inference procedure is proposed for the direct average treatment effect on the treated (DATT) and the spillover effect defined in the current paper. As a result, practitioners can still use the spillover function/exposure mapping they choose based on their domain/institutional knowledge. Inspired by savje2023causal, the estimand remains well-defined under misspecification of the exposure. Only its interpretation needs to be modified. In addition, the assignment variables are allowed to be spatially correlated as is often the case in practice with spatial data.

Putting all the pieces together, I propose doubly robust estimators for the direct treatment effect and spillover effect. The proposed doubly robust estimator is a modified version of the augmented inverse probability weighting (AIPW) estimator, which only requires correct specification of either the propensity scores of treatment/exposure or the conditional mean of the outcomes. The conditional inference approach in the current paper leads to a different variance-covariance matrix which may require a new variance estimator when necessary.

Besides the main contribution above, I study the identification of canonical DID estimators available in the literature in Section (ref) below. I provide conditions under which canonical estimators can still identify meaningful causal parameters. This discussion alone would be of interest to practitioners. The proposed doubly robust estimators are applied in Section (ref) to study the policy effect of special economic zones (SEZ) in China. Appendix (ref) summarizes the detailed steps for direct effect estimation.

Related Literature: delgado2015difference and butts2021difference allow interference in DID estimation in a two-way fixed effects (TWFE) estimating equation (often without covariates) from a superpopulation perspective. huber2021framework also propose a DID approach to estimate spillover effect and total effect within a superpopulation framework.\footnote{Their potential outcomes are defined as functions of individual and regional treatments, where individual treatment status is a function of the regional treatment. Therefore, huber2021framework is more applicable to studies of local equilibrium effects.} I instead propose doubly robust estimators with flexible choice of exposure mappings and adopt a finite population framework. In addition, all three papers mentioned above share some or all the restrictions of the general interference literature. More specifically, most methodological literature studies spillover effects in a single cross section of experimental data and assumes partial interference or limits interference to immediate neighbors. Additionally, they assume that the function of dependence on neighbors' treatments is known and correctly specified. See, for instance, hudgens2008toward and aronow2017estimating. I relax these assumptions in a DID context in Section (ref) below discussing misspecification of the exposure mapping, adapting tools from savje2023causal and leung2022causal. Design-based DID estimation has been studied by athey2022design and rambachan2020design, but they keep the SUTVA. sant2020doubly have proposed AIPW estimators in the DID context, maintaining SUTVA and the superpopulation perspective.

Setup

Environment

I start with the relatively simple setting of panel data with two time periods; $t=1,2$ stands for the time period before and after treatment respectively. Consider a sequence of lattices of (possibly) unevenly placed locations in $\mathbb{R}^d$, $\{D_M\subseteq \mathbb{R}^d, d\geq 1\}$, where $M$ indexes the sequence of finite populations. In a finite population setting, the parameters of interest are defined within the finite population. Because I consider the case where the sample coincides with the population for spatial data, I let the population size $|D_M|$ diverge to infinity in deriving the asymptotic properties of the proposed estimators and conducting inference, where $|V|$ denotes the cardinality of a finite subset $V\subseteq D_M$.

I briefly summarize the notation used throughout the paper. I adopt the metric $\rho(i,j)=\max_{1\leq l\leq d}|j_l-i_l|$ in space $\mathbb{R}^d$, where $i_l$ is the $l$-th component of $i$. This metric could capture geographic distance or some economic distance between units $i$ and $j$. The distance between any subsets $K,V\subseteq D_M$ is defined as $\rho(K,V)=\inf\{\rho(i,j): i\in K\mbox{ and } j\in V\}$. For any random vector $X$, $\left\lVertX\right\rVert_p=\big(\mathbbm{E}\left\lVertX\right\rVert^p\big)^{1/p}$, $p\geq 1$, denotes its $L_p$-norm. Lastly, $C$ denotes a generic positive constant that may vary under different circumstances.

For each unit $i$ in the population, there is a stochastic assignment variable $W_i\in \{0, 1\}$, a vector of fixed attributes $z_i=(z_i^{ind},z_i^{neigh})$ that possibly includes attributes of $i$'s neighborhood $z_i^{neigh}$ in addition to individual characteristics $z_i^{ind}$, and a vector of stochastic unobservables $U_{it}$. The potential outcome function for any $i\in D_M$ is defined as $h_{it}(\cdot): \{0,1\}^{|D_M|}\times \mathbb{R}^{\dim(z_i)}\times \mathbb{R}^{\dim(U_{it})} \to \mathbb{R}$. I emphasize the treatment vector of the entire population by denoting the potential outcomes as $y_{it}(w_i,\bm{w}_{-i})=h(w_i,\bm{w}_{-i}, z_{i}, U_{it})$, where $\bm{w}_{-i}=\{w_j, j\in D_M, j\neq i\}$.\footnote{As in manski2013identification, the potential outcome function defined here can be considered as the response function, namely the reduced form of structural equations where the structural potential outcome may depend on other units' treatments as well as outcomes.} The dependence of the potential outcomes on the fixed attributes and stochastic unobservables is indicated by its $i,t$ subscript. The realized potential outcomes are denoted by $Y_{it}=y_{it}(\bm{W})$. Notice that $(\bm{W}, \bm{z}, \bm{Y},\bm{U})=\{(W_i, z_i, Y_{it}(\cdot),U_{it}), i\in D_M, M\geq 1\}$ are triangular arrays of random fields defined on a probability space $(\Omega, \mathcal{F}, P)$. Exposure mapping is defined by the function $G_i=G(i, \bm{W}_{-i})\in \mathcal{G}$, where $\mathcal{G}$ is a discrete set.\footnote{Given our finite population setting and a binary individual treatment variable $W_i$, $\mathcal{G}$ is a finite set. This is also the approach taken by aronow2017estimating and leung2022causal.} Therefore, $G(\cdot)$ maps the treatment status of all units except $i$ to an exposure value.

In line with common practices, empirical researchers can construct the exposure mapping $G(\cdot)$ in the following manner, which is assumed to capture the true interference structure. Given a fixed $K$, define the $K$-neighborhood of the unit $i$ as \[ \mathcal{N}(i,K)=\{j\in D_M: \rho(i,j)\leq K, j\neq i \} \] Let $\bm{w}_{\mathcal{N}(i,K)}=(w_j: j\in \mathcal{N}(i,K))$ be the treatment vector of units within $i$'s $K$-neighborhood. There exists $K<\infty$ such that for all $\bm{w}_{-i}$ and $\bm{w}'_{-i}$ such that $\bm{w}_{\mathcal{N}(i,K)}=\bm{w}'_{\mathcal{N}(i,K)}$, $G(i,\bm{w}_{-i})=G(i,\bm{w}'_{-i})$. As a result, the specified exposure mapping function restricts spillover effects within the immediate $K$-neighborhood of each unit.\footnote{With correctly specified exposure mapping, asymptotic distribution of the proposed estimators below can be derived for general exposure mapping not restricted to $K$-neighborhood interference under a different set of local dependence assumptions than the ones listed below. Later on, I relax the immediate $K$-neighborhood interference to allow for potential misspecification of the exposure mapping. With the current approach of constructing the exposure mapping, the asymptotic results can be derived in a unified framework accommodating both correctly or incorrectly specified exposure mapping.} A leading example is $G_i=\sum_{j\in D_M, j\neq i} A_{ij}W_j/\sum_{j\in D_M,j\neq i} A_{ij}$, where $A_{ij}=1$ if the distance between units $i$ and $j$ is within a cutoff $K$. Given a correctly specified exposure mapping with $G(i,\bm{w}_{-i})=g$, $\tilde{y}_{i2}(w_i,g)=y_{i2}(w_i,\bm{w}_{-i})$.

Estimands of Interest

This paper is interested in the expected finite population average, i.e., the average of the expected potential outcome across all units in the finite population. In other words, I focus on conditional inference given fixed attributes $z_i$; see abadie2014inference and jin2023tailored for detailed discussion of conditional parameters and conditional inference.

There are two types of estimands of interest. For the main part of the paper, I will focus on the first type, the direct treatment effect. There is more than one way to define the parameter of interest. As an analogy to the expected average treatment effect in savje2021average, the overall direct effect can be defined as \[ \tau=\frac{1}{|D_M|}\sum_{i\in D_M}\mathbbm{E}\big[y_{i2}(1,\bm{W}_{-i})-y_{i2}(0,\bm{W}_{-i})|W_i=1, z_i\big], \] which marginalizes over the treatment assignment vector. The overall direct effect is a natural extension of the average treatment effect on the treated (ATT) as it coincides with the ATT when units do not interfere ($y_{i2}(W_i,\bm{W}_{-i})$ reduces to $y_{i2}(W_i)$).

Often the time, in addition to a summary of the direct effects, researchers can also be interested in direct effect at different exposure levels. The overall direct effect is highly related to the direct average treatment effect on the treated (DATT) at exposure levels $g\in \mathcal{G}$ defined in equation ((ref)) below.\footnote{Their relationship is explained in equation ((ref)) in Appendix (ref).}

Without interference, $\tau(g)$ becomes a conditional ATT with $G_i=g$ serving as another characteristic of unit $i$. I focus on the identification and estimation of DATT as the estimation of the overall direct effect follows using weighted averages when the weights are correctly specified.\footnote{As shown in Section 3.1, the canonical DID estimand, which ignores interference, fails to identify $\tau$ for spatially correlated assignments. Notice that I do not propose an estimator for the overall direct effect in this paper. savje2021average show when common estimators, such as the Horvitz-Thompson and H\'ajek estimators, are consistent for the expected average treatment effect without proposing a new estimator.} Also, by contrasting different exposure levels, the definition of DATT facilitates the discussion of the second estimand, the spillover effect. In the interest of space, this is delegated to Appendix B.

The all-or-nothing effect, $\frac{1}{|D_M|}\sum_{i\in D_M}\mathbbm{E}\big[y_{i2}(\bm{1})-y_{i2}(\bm{0})| z_i\big]$, where $\bm{1}$ and $\bm{0}$ are unit and zero vectors, cannot be consistently estimated (basse2018limitations). Instead, direct effects and spillover effects summarize different aspects of policy effects. For the DATT I consider in the main text, it captures the direct treatment effect at different exposure levels. In a vaccination example, if the direct effect of vaccinating an additional individual is almost zero given that a certain fraction of the population are already vaccinated, this can serve as an indicator of herd immunity being achieved.

I use the empirical application in Section (ref) below to illustrate the relevant variables and estimands. The data I use come from lu2019place, who study how China's SEZ policy impacts various outcomes $Y_{it}$, such as the logarithm of firm output at the village level. If village $i$ is located within the boundaries of a SEZ, direct treatment variable $W_i$ is equal to one and otherwise, it is equal to zero. Exposure mapping $G_i$ is a binary variable equal to one if the leave-one-out ratio of SEZ villages to the total number of villages in the county $c$ in which village $i$ is located is greater than the mean ratio among all counties. In this case, the DATT captures the direct effect of establishing a SEZ in village $i$ given the fraction of SEZ villages within a county. It is possible that the direct effect is lower when there is a higher proportion of neighboring SEZ villages. This can help determine whether establishing an additional SEZ is cost-effective.

Identification

The first question when relaxing SUTVA is what the canonical DID estimator identifies if spillover effects are incorrectly ignored. Namely, will the canonical DID estimator still consistently estimate ATT in the presence of interference? forastiere2021identification discuss bias of the difference-in-means estimator when SUTVA is wrongly assumed in observational studies on networks. To my knowledge, the literature has not yet investigated DID type estimators. To facilitate the discussion of identification, I impose the following assumptions.

To simplify notation, I assume that the overlap assumption applies to every unit in the population. With certain exposure mapping specifications, this might not be plausible. An easy fix is to change the estimand by averaging over the subpopulation where $G_i$ can take on the value $g$. Failure to satisfy the overlap condition for $p(z_i)$ is trickier. If one is willing to move the goalpost by redefining the population, one can drop units that always or never take treatment. The good news is that for the redefined population, we can still observe the treatment assignment vector of the original population since the treatment status of the dropped units is fixed and known. This way, dropping the always or never takers will not affect the exposure mapping. On the other hand, to deal with weak overlap conditions in practice without changing the population or estimand, one can consider approaches proposed by ma2020robust and man2023doubly to trim propensity scores and correct the resulting bias simultaneously.

Assumption (ref) requires that the potential outcome in the first time period prior to treatment is always equal to the potential outcome without treatment nor spillover ($\bm{w}_{-i}=\underline{0}$). The no-anticipation assumption is quite standard in the literature, sometimes implicitly assumed.

With correctly specified exposure mapping, I impose the following parallel trends assumption:

Notice that in the parallel trends assumption, as no one is treated at $t=1$ there is no spillover in the potential outcome function in the first time period. Namely, moving from the first to the second time period, in the absence of direct treatment the conditional mean of the potential outcomes for the treated and the untreated with the same level of exposure in the second time period follows the same trend. Equation ((ref)) serves as our starting point for identification.

There is a growing literature on justification and falsification of the parallel trends assumption under SUTVA; see, for instance, ghanem2022selection and roth2023parallel. When parallel trends might be violated, rambachan2023more present confidence sets for the identified set of treatment effects. The extension of these analyses to parallel trends with interference is outside the scope of the current paper and left as future research.

Canonical DID

The usual ATT under the SUTVA is \[ \tilde{\tau}=\frac{1}{|D_M|}\sum_{i\in D_M}\mathbbm{E}\big[y_{i2}(1)-y_{i2}(0)|W_i=1,z_i\big]. \] Here, the potential outcomes are determined solely by unit $i$'s own treatment. Suppose the canonical DID estimator consistently estimates

Examples include the TWFE linear estimating equation in Remark 1 in sant2020doubly under the additional restrictions of the data generating process therein, as well as the inverse probability weighting (IPW) estimator in abadie2005semiparametric. If the usual (conditional) parallel trends assumption holds without interference, $\tau_{canonic}$ would be equivalent to $\tilde{\tau}$.

If SUTVA is violated, $\tilde{\tau}$ is not well defined as the potential outcome should depend on the entire assignment vector. Also, DATT is generally determined by the specified exposure level. As a result, I use the overall direct effect as a benchmark for comparison. Using the law of iterated expectations, $\tau$ and $\tau_{canonic}$ can be decomposed in the following way:

The equality of the conditional probabilities holds if $G_i\perp \!\!\! \perp W_i\ |\ z_i$. However, conditional independence can be easily violated if either of the following is true: (\romannumeral 1) $G_i$ and $W_i$ are linked through covariates not included in $z_i$; (\romannumeral 2) neighbors' behavior affects unit $i$'s treatment uptake; (\romannumeral 3) similar neighborhood characteristics drive the assignment mechanism; see forastiere2021identification for a parallel discussion allowing interference on networks under unconfoundedness. As a consequence, the overall direct effect can be either underestimated or overestimated by the canonical DID estimators.

Following Remark (ref), Table (ref) below provides numerical results from a simple simulation study. The detailed data generating process is summarized in Appendix (ref). The exposure mapping $G_i$ is correlated with $W_i$ in the design and there is no first period outcome heterogeneity associated with $G_i$. In the first design, when there is no spillover effect on the non-directly treated units, i.e., $\mathbbm{E}(Y_{i2}|W_i=0,G_i=1,z_i)-\mathbbm{E}(Y_{i2}|W_i=0,G_i=0,z_i)=0$, $\forall i$, $\tau_{canonic}$ is almost identical to $\tau$. By contrast, in the second design, there is spillover effect on the non-directly treated units with $\mathbbm{E}(Y_{i2}|W_i=0,G_i=1,z_i)-\mathbbm{E}(Y_{i2}|W_i=0,G_i=0,z_i)=1$, $\forall i$. With the direct effects $\tau(1)=-1$ and $\tau(0)=0$, the spillover effect is comparable to the direct effect in magnitude. In this case, $\tau_{canonic}$ can be substantially different from $\tau$.

Modified Two-Way Fixed Effects

One way to estimate the spillover effect suggested in the existing literature is to augment the TWFE DID estimating equation with another binary indicator $S_{i}$ equal to one if a unit is close to a treated unit; see, for instance, attempts in di2004police and butts2021difference. Using the notation in the current paper, I modify the estimating equation to be

where $W_{it}=W_i*\mathbbm{1}\{t=2\}$ and $S_{it}=S_i*\mathbbm{1}\{t=2\}$. $\hat{\beta}_1$ and $\hat{\beta}_1+\hat{\beta}_3-\hat{\beta}_2$ estimated from equation ((ref)) would be consistent for the DATT defined by

and \[ \bar{\tau}(1)=\frac{1}{|D_M|}\sum_{i\in D_M}\Big[\mathbbm{E}\big(y_{i2}(1,\bm{W}_{-i})-y_{i2}(0,\bm{W}_{-i})|W_i=1,S_i=1\big)\Big] \] respectively, under some parallel trends assumptions.

We can see that given the estimating equation of the augmented TWFE, the specified exposure mapping is fixed as $\mathbbm{1}\{A_s\bm{W}>0\}=S_i$, where $A_s$ is the spatial weighting matrix with units being neighbors if their distance is less than or equal to $\bar{\rho}$. Only when the interference structure coincides with the indicator function $\mathbbm{1}\{A_s\bm{W}>0\}$ along with the correct distance cutoff, can the augmented TWFE identify the exact direct ATT. In contrast, when the true interference structure is not $\mathbbm{1}\{A_s\bm{W}>0\}$, the proposed estimands in Section (ref) below can still identify the exact direct ATT by choosing correct specification of the exposure mapping with potentially multiple levels of neighborhood exposure $g$. Meanwhile, covariates can be flexibly accounted for in the proposed estimands by assuming conditional parallel trends.

Doubly Robust Estimand

Since ignoring the spillover effect is only harmless under special scenarios, new estimators need to be proposed for the DATT. Under parallel trends and overlap assumptions, the DATT can be identified by inverse weighting using propensity scores.

To simplify notation, I use $\mathbbm{E}_D$ to denote the finite population average conditional on the attributes $\bm{z}$ from now on. Without the $G$ indicator and the additional propensity scores for spillover, the IPW-DID estimand is the same as the estimand proposed in abadie2005semiparametric.

Alternatively, the DATT can also be identified through regression adjustment. Define the conditional means of the potential outcome as

The regression adjustment estimand is

To allow for more robustness against misspecification of the propensity scores or the conditional means of the outcomes, the IPW-DID estimand can be extended to an AIPW estimand. Let $m_{t,wg}(z_i)$ denote the model for equation ((ref)). Denote $\Delta m_{wg}(z_i)=m_{2,wg}(z_i)-m_{1,wg}(z_i)$. Furthermore, let $\eta(z_i)$, $\eta_{1g}(z_i)$, and $\eta_{0g}(z_i)$ be the models for the propensity scores in equations ((ref))-((ref)), respectively. The doubly robust estimand is

As a result, we can consistently estimate $\tau(g)$ as long as either the models for the propensity scores or the models for the conditional means of the outcome are correctly specified in the doubly robust estimand.

Misspecified Exposure Mapping

In this section, I show how to proceed with DID estimation with a chosen $G(\cdot)$ function that is potentially misspecified. For more discussion on some common choices of $G(\cdot)$ and how they can be accommodated in the current framework, I refer readers to Appendix (ref).

In an ideal scenario, empirical researchers would like to come up with a functional form of $G(\cdot)$ that captures actual interactions among neighbors as well as conveys clear causal explanations. Because of the unknown interference structure and the high dimensional treatment assignment vector of the entire population, choosing a $G(\cdot)$ that achieves both goals is challenging.\footnote{manski2013identification and basse2018limitations formally point out that there exist no consistent treatment effect estimators under arbitrary interference. It is therefore necessary to make dimension reduction assumptions about the interference structure in order to identify meaningful treatment effect parameters.} Recall the example in Section (ref), $G_i=\sum_{j\in D_M, j\neq i} A_{ij}W_j/\sum_{j\in D_M,j\neq i} A_{ij}$, where $A_{ij}=1$ if the distance between units $i$ and $j$ is within a certain cutoff $\bar{\rho}$. Besides the somewhat arbitrary cutoff $\bar{\rho}$, the impact of $i$'s neighbors may not be exchangeable in reality, e.g., unit $l$ may have greater influence on $i$ compared to unit $m$. That said, the specification above might still capture meaningful policy effects. Consequently, using domain knowledge in the context of each specific empirical question, coming up with a $G(\cdot)$ that summarizes interesting and relevant policy implications for both direct treatment and spillover effects might be a good starting point.

To overcome potential misspecification of the spillover pattern, I consider the expected direct treatment effect at certain neighborhood exposure levels as the parameter of interest inspired by savje2023causal. The causal estimands I define coincide with the exact direct treatment effect when the exposure mapping is correctly specified and remain well-defined even under misspecification.

The parameter of interest is now the expected direct average treatment effect on the treated (EDATT) at exposure levels $g\in \mathcal{G}$, which identifies the direct ATT that would realize in expectation at the specified exposure level.\footnote{One can similarly aggregate EDATT to the overall direct effect.}

The key ingredient of the definition is the expected potential outcome at exposure level $g$,

where, in addition to the stochastic potential outcomes, the expectation is taken over all possible realizations of $\bm{W}_{-i}$ given the specified exposure mapping $G(i,\bm{W}_{-i})=g$. The definition of the expected potential outcome is different from what is initially proposed in savje2023causal, in which the potential outcome is fixed in an experimental setting and the expectation is with respect to the assignment variables only. Not only that I split the entire treatment vector into $w_i$ and $\bm{w_{-i}}$, but also the stochastic nature of the potential outcomes needs to be taken into account. The randomness of the potential outcome function brings up challenge to causal interpretation of the spillover effect estimand. Appendix (ref) provides more detailed reasoning.

In order to provide a general framework for identifying EDATT, I impose the following assumption instead. It reduces to Assumption (ref) under correctly specified $G(\cdot)$.\footnote{I consider Assumption (ref) as an attempt to accommondate misspecification of exposure mapping in practice. Its plausibility depends on the true interference pattern, the chosen exposure mapping, and the design of the individual assignments and hence needs to be examined case by case. If Assumption (ref) does not hold exactly with misspecified $G(\cdot)$, one might also consider sensitivity analysis along the lines in rambachan2023more.}

Corollary (ref) implies that the same identification results for DATT also hold for EDATT. Consequently, one can use the same estimator proposed in Section (ref) below. Notice that the identification results here and those in Section (ref) hold for any arbitrary $G(\cdot)$, either correctly or incorrectly specified. To show the asymptotic properties of the estimator in a unified framework, I apply the device of approximate neighborhood interference (ANI) in leung2022causal to spatial data. The ANI device implies that treatments of units outside of $i$'s $K$-neighborhood can legitimately influence $i$'s potential outcome as long as treatments assigned to units further from $i$ have a smaller, but possibly nonzero, effect on $i$’s response. This device is used to make the asymptotic derivation more tractable but is not required for the identification of the parameter.

In summary, practitioners can still use the spillover function they choose according to the construction of the exposure mapping $G(\cdot)$ in Section (ref) based on their domain/institutional knowledge. Having said that, the actual potential outcome function does not restrict the interference structure exactly as $G(\cdot)$. Based on the same causal effect estimator and a unified asymptotic distribution, practitioners do not need to change their estimation and inference procedure based on their stance on the specification of the spillover pattern. If misspecification of the exposure mapping is a concern, one only needs to modify the interpretation of DATT to EDATT.

Pre-trends

In empirical research, tests for pre-trends remain common despite the caution described in roth2022pretest. A placebo DID is typically applied to multiple periods observed before treatment by imposing a hypothetical period of adoption of treatment. Something similar can be done in the context of interference. Since Assumption (ref) nests Assumption (ref) when the exposure is correctly specified, I focus on the pre-trends of Assumption (ref). In the simplest case, suppose there are two time periods $t=0, 1$ prior to treatment, by imposing the placebo treatment between time periods 0 and 1, we would like to test

The potential outcome $y_{i1}(0, \bm{W}_{-i})$ is not observable because no unit is treated in time period 1. Nevertheless, under the no anticipation assumption, equation ((ref)) reduces to

Namely, one can test whether subgroups defined by the combination of direct treatment status and exposure level have different trends before actual treatment. In reality, the testing equation ((ref)) is a test of both no anticipation and parallel pre-trends.

Asymptotic Properties of the Parametric Estimator

I am primarily concerned with estimating the DATT (EDATT) in this section. Spillover effects are defined in Appendix (ref). Their estimation is similar to that of the DATT. I propose a GMM estimator combining equation ((ref)) with moment conditions for the propensity scores and conditional means of outcomes chosen by the empirical researcher. Since the inference is only conditional on the covariates $\bm{z}$, all unobservables are identically distributed conditional on $z_i$.\footnote{It suffices that the first conditional moment is identical.} The inference framework implies that the individual propensity score function and the individual conditional mean function of the outcome remain the same across units.

To make the estimators more robust to misspecification of these functions, one can use various moment conditions to identify the propensity scores. One option is the covariate balancing propensity scores (CBPS) in imai2014covariate or similarly the inverse probability tilting estimator in graham2012inverse, which can be locally more robust than the propensity scores based on maximum likelihood estimation (MLE).\footnote{The alternative would be estimating all functions semiparametrically or nonparametrically, which is left as future work.}

I denote a generic moment condition for propensity scores as

and

where $z_i$ can contain neighbors' attributes. For instance, the moment conditions for CBPS are

and for $g=1,2,\dots,G-1$,

where $P(W_i=1|z_i)$ is some probability for a binary response, such as $\frac{exp(z_i\gamma^*_1)}{1+exp(z_i\gamma^*_1)}$, and $P(G_i=g|W_i,z_i)$ is some probability for discrete choices. Similarly, generic conditional moment conditions are denoted by

and

Alternatively, one can model the conditional mean for $\Delta Y_i=Y_{i2}-Y_{i1}$ and formulate the moment condition as

Leading cases for outcome regression are moment conditions from (nonlinear) least squares. Lastly, the moment condition for $\tau(g)$ is a restatement of equation ((ref))\footnote{In practice, it is recommended to normalize the weights for IPW type estimators. Changing the moment condition with normalized propensity scores -- where the weights sum to unity -- does not affect asymptotic normality of the GMM estimator. In fact, estimators with normalized weights consistently show better finite sample performance in the simulations below.}. Denote $\theta^*_M=({\gamma^*_1}',{\gamma^*_2}',{\gamma^*_3}',{\gamma^*_4}',\tau(g))'$.

Let $X_i=\{Y_{it},W_i,G_i, z_i\}$, $q(X_i,\theta)=(q_1'(\gamma_1),q_2'(\gamma_2),q_3'(\gamma_3),q_4'(\gamma_4),q_5(\theta))'$, and $\widehat{\Psi}$ be the weighting matrix with dimensions larger or equal to that of $\theta$.

The GMM estimator is the solution to the finite population minimization problem in equation ((ref)). And the estimator of $\tau(g)$ is the last element of $\hat{\theta}$. Notice that the same estimation procedure applies to $\tau^*(g)$. I impose the following assumptions to study the asymptotic distribution of the GMM estimator.

I adopt the increasing domain asymptotics with spatial data as implied by Assumption (ref). The assumption of the minimum distance ensures the expansion of the finite population region in our asymptotic framework when the population size keeps growing. This rules out the case where population size grows within a bounded population region in the sense that units become arbitrarily dense in a given region.

Assumption (ref) is a modified version of Assumption 4 in leung2022causal. leung2022causal varies $\bm{W}'_{D_M\backslash \mathcal{N}(i,s)}$ in an arbitrary way but these treatments outside of the $s$-neighborhood are fixed at zero here. Assumption (ref) essentially implies that treatments of units from $s$ distance away from $i$ should have a minimal impact as the distance $s$ gets larger. This way, we can allow interference from outside the immediate $K$-neighborhood while still being able to derive the asymptotic properties of the proposed estimators. With that said, interference is restricted in a way that a unit's exposure is primarily, but not entirely, determined by the assignments of neighbors closer to it.

Needless to say, Assumption (ref) is automatically satisfied under correctly specified exposure mapping within the $K$-neighborhood. We can also accommodate misspecified exposure mapping that satisfies Assumption (ref) in our asymptotic framework. Appendix (ref) gives an overview of different approaches to modeling interference taken in the literature and compares them to ANI.

I adopt $\psi$-dependence in kojevnikov2021limit as the notion of weak dependence throughout the paper. Notice that $\alpha$-mixing is a special case of $\psi$-dependence. Let $\mathcal{L}_{\nu, h}$ denote the collection of bounded Lipschitz real functions $f(\cdot)$ on $\mathbb{R}^{\nu\times h}$ with the Lipschitz constant $\text{Lip}(f)<\infty$ and $\left\lVertf\right\rVert_\infty<\infty$, where $\left\lVertf\right\rVert_\infty=\sup_x|f(x)|$. Denote the collection of subset pairs as \[ \mathcal{P}_M(h,h';s)=\{(H,H'):H,H'\subseteq D_M, |H|=h,|H'|=h', \rho(H,H')\geq s\}. \]

I require $\tilde{\kappa}_{M,s}$ to approach zero as $s$ grows. $\psi$-dependence bounds the covariances of any two subsets of observations distant from each other.

On top of possible interference, Assumption (ref) allows assignment variables to be spatially correlated as well.

Lemma (ref) shows that Assumption (ref) together with Assumption (ref) ensure that the data is weakly dependent even with spatially correlated assignment and spillover effect, regardless of the correct specification of the exposure mapping.

To adapt the limit theorems in kojevnikov2021limit to spatial data, I replace the network denseness with the cardinality of the spatial sets implied by Lemma A.1 in jenish2009central. As a result, Assumption 3.2 in kojevnikov2021limit is modified as

Assumption (ref) is in line with Assumption 3(b) in jenish2009central for $\alpha$-mixing random fields. It imposes a trade-off between the sizes of $s$-neighborhood boundaries and the rate of decay of spatial dependence. Let $\sigma_M^2=Var\big[\sum_{i\in D_M}\lambda'q(X_i,\theta_M^*)|\bm{z}\big]$ for a nonzero vector $\lambda$. Similarly, Assumption 3.4 in kojevnikov2021limit is modified as

The rate of $\tilde{\kappa}_{M,s}$ is implicitly implied by Assumption (ref), which restricts the spatial sets and the cross-sectional dependence of random variables. A sufficient condition for the first part of the assumption is \[ \frac{|D_M|}{\sigma_M^{2+k}}r_M^{kd}\sum^{\infty}_{s=0}s^{d-1}\tilde{\kappa}_{M,s}^{1-\frac{2+k}{p}}\to 0. \] Analogous conditions can be found in jenish2009central as equations (B.18) and (B.19) therein.

The notation used in the asymptotic distribution of the GMM estimator is introduced as follows. Define

where

and

$\Delta_{ehw,M}$ and $\Delta_{spatial,M}$ account for heteroskedasticity and spatial correlation respectively, whereas $\Delta_{E,M}$ and $\Delta_{ES,M}$ are their finite population counterparts. Denote \[ R_M^*=\mathbbm{E}_D\big[\nabla_\theta q(X_i,\theta^*_M)\big] \] and

where $\widehat{\Psi}-\Psi_M\overset{p}\to \textbf{0}$.

With the consideration of interference and potentially spatially correlated assignments, we need to make inference robust to spatial correlation. As a common approach to adjust the variance estimator for spatial correlation, the usual spatial heteroskedasticity and autocorrelation consistent (SHAC) variance estimator is defined as \[ \hat{V}=\big(\hat{R}'\widehat{\Psi}\hat{R}\big)^{-1}\hat{R}'\widehat{\Psi}\tilde{\Omega}(\hat{\theta})\widehat{\Psi}\hat{R} \big(\hat{R}'\widehat{\Psi}\hat{R}\big)^{-1}, \] where \[ \hat{R}=\frac{1}{|D_M|}\sum_{i\in D_M}\nabla_\theta q(X_i,\hat{\theta}) \] and \[ \tilde{\Omega}(\theta)=\frac{1}{|D_M|}\sum_{s=0}^{\infty}\omega\Big(\frac{s}{b_M}\Big)\sum_{i\in D_M}\sum_{j\in D_M, s\leq \rho(i,j)< s+1}q(X_i,\theta)q(X_j,\theta)' \] with $b_M$ being the bandwidth parameter. I impose the following assumption for the estimation of the variance-covariance matrix.

Assumption (ref)(\romannumeral 1) is satisfied by common choices of kernels, including the Bartlett and Parzen kernels. Assumption (ref)(\romannumeral 2) requires that the kernel weights $\omega\big(\frac{s}{b_M}\big)$ converge to one sufficiently fast as $M\to\infty$. Part (\romannumeral 3) of Assumption (ref) regulates the growth rate of the bandwidth $\{b_M\}$.

The conservativeness of the usual variance estimator for conditional variance has also been investigated in abadie2014inference under the independence assumption for the heteroskedasticity-robust variance matrix. I extend it to the case with spatial correlation here when $\Omega_E$ is PSD based on a PSD kernel weighting matrix. An exception to Remark (ref) is when $\mathbbm{E}\big[q(X_i,\theta^*_M)|\bm{z}\big]=\bf{0}$ for all $i\in D_M$. In this case, the usual variance-covariance matrix estimator is no longer conservative as $V_E=\bf{0}$. With heterogeneous direct treatment effect or misspecification of either the propensity scores or conditional means, $\mathbbm{E}\big[q(X_i,\theta^*_M)|\bm{z}\big]\neq\bf{0}$.

That said, I would like to highlight a few points. First, because $\tilde{\Omega}(\hat{\theta})$ is a conservative estimator for $\Omega_M$, even if we choose $\Psi_M$ as the optimal weighting matrix $\Omega_M^{-1}$, using $\widehat{\Psi}=\tilde{\Omega}(\hat{\theta})^{-1}$ in estimation is not going to achieve the most efficient GMM estimator. The usual variance estimator is therefore conservative not only because of the neglect of the additional terms in the variance-covariance matrix but also because the optimal weighting matrix is not consistently estimated. Of course, when the model is just identified, the weighting matrix choice is irrelevant.

Second, unlike the finite population variance-covariance matrix in xu2022design, the conditional spatial-correlation robust variance matrix is consistently estimable because it is no longer conditional on the unobserved potential outcomes. There are different approaches one can take. However, since the usual SHAC variance estimator is known to suffer from downward bias especially when the spatial correlation is high, it is not always necessary to estimate the smaller conditional variance matrix.

Simulations

In the simulation, I show the finite sample performance of the proposed estimators for DATT (EDATT). I consider an irregularly spaced lattice with $|D_M|=400$ units. The locations ($s_{1,iM}$,$s_{2,iM}$) are drawn once and kept fixed across replications. Each of $s_{1,iM}$ and $s_{2,iM}$ is independently drawn from $\mathcal{U}(0,20)$. The distance between units $i$ and $j$ is measured by $\rho(i,j)=\max\{|s_{1,iM}-s_{1,jM}|,|s_{2,iM}-s_{2,jM}|\}$. Units are considered neighbors if $\rho(i,j)\leq 0.3$ with the neighborhood structure summarized by the normalized contiguity matrix, $A$. After ruling out units without neighbors, the effective size of the subpopulation eligible for spillover reduces to 350.

I consider two time period panel data. The potential outcome function in the first time period remains the same across different designs. \[ y_1(0,\underline{0})=1+z+e_1, \] where $z$ is the individual covariate independently drawn from the standard normal distribution and kept fixed, while $e_1$ is the first time period unobservable. There is a single binary treatment variable $W=\mathbbm{1}\{p(z^*)>u\}$ with $u_i\overset{i.i.d}\sim \mathcal{U}(0,1)$. I vary the second time period potential outcome function and the assignment probability $p(z^*)$ in different designs summarized in Table (ref) below. $z^*=(z,z_u)$, where the vector of $z_u$ in the assignment probability is drawn from a multivariate normal distribution with mean zero and a variance-covariance matrix equal to 0.5 raised to the power of the distance between units. Thus, $z_u$ is a spatially correlated locational covariate that stands for neighborhood similarity, which might be neglected in naive estimation assuming away spillover effect. Along with the individual second time period unobservable $e_2$, $e_{i1}|W_i,W_{-i},z_i\sim \mathcal{N}(W_i*z_i,1)$, $e_{i2}|W_i,W_{-i},z_i\sim \mathcal{N}(W_i*z_i,1)$ and $e_{i1}\perp \!\!\! \perp e_{i2}|W_i,W_{-i},z_i$, $\forall\ i$. The specified exposure mapping is denoted by $G=\mathbbm{1}\{A\bm{W}>0\}$, which may or may not coincide with the true interference structure.

I compare the following estimators: the canonical TWFE, Abadie's IPW estimator, the augmented TWFE, regression adjustment, IPW estimator with either MLE or CBPS moment condition for the propensity scores, and the proposed AIPW estimator with either MLE or CBPS moment condition for the propensity scores. Appendix (ref) contains the standard deviation of the proposed estimators and the coverage rate of the 95% confidence intervals based on the usual SHAC standard errors for the doubly robust estimators.

For the canonical Abadie's IPW estimator, I only include $z$ in the logit model of $W$ as interference is assumed away when employing the canonical DID estimators. As an illustration of Proposition (ref), I also report Abadie's IPW estimator with $z$, $A\bm{z}$, and $z_u$ included in the logit model, which leads to conditional independence of $W$ and $G$. The estimation of the augmented TWFE follows equation ((ref)) with $S_i=G_i$. For the estimation of the proposed IPW, regression, and AIPW estimator accounting for spillover effect, the propensity scores for $W$ and $G$ are estimated based on a logit model on $z$, $A\bm{z}$, and $z_u$ and a logit model on $W$, $z$, $A\bm{z}$, and $z_u$, respectively. For the first time period data, I regress $Y_1$ on $W$, $z$, and $W*z$. As for the second period data, I regress $Y_2$ on $W$, $z$, $W*z$, and $G$. All estimators involving weighting are weighted by the normalized propensity scores. The results are summarized across 10,000 replications.

According to the population generating process, the direct effects are $\tau(1)=\tau(0)=1$ in designs 1-5 and $\tau(1)=1$, $\tau(0)=0$ in design 6. In the last design, the overall direct effect is approximately 0.607. The point estimates for the direct effect are summarized in Table (ref). In designs 1 and 2, neighborhood similarity does not drive treatment assignment. As a result, the canonical Abadie's IPW estimator with covariate $z$ closely estimates the overall direct effect. The canonical TWFE only performs well in design 1 as the estimating equation of TWFE rules out $z$-specific time trends, which is violated in all other designs. The augmented TWFE estimators suffer from the same linearity restriction in their estimating equation as the regular TWFE. With the inclusion of both $z$ and $z_u$, Abadie's IPW estimates are very close to the overall direct effect.

The proposed estimators accounting for the spillover effects all perform relatively well. Due to the specific exposure mapping functional form, the overlap condition holds better for exposure level one than zero. Consequently, the point estimates for the direct effect estimator at exposure level one are slightly more accurate than the results for the estimator at exposure level zero. It is worth mentioning that the propensity score model for $G$ is always misspecified. The outcome regressions are also misspecified in designs 4-6. Nevertheless, the estimates from the proposed IPW and doubly robust estimators are all quite close to the truth and much more accurate than the TWFE type estimators.

The doubly robust estimators improve upon regression adjustment and IPW alone, especially at exposure level zero. The only exception is design 5. Since the outcome regression is more severely misspecified than in other designs, we do not see improvement moving from IPW to AIPW. Nevertheless, the AIPW estimates are still better than the regression adjustment estimates. Estimators with CBPS moment condition slightly improve upon estimators with MLE moment condition. When the overlap condition holds weaker in other population generating processes, for instance, changing the assignment probability to $p(z^*)=\frac{exp(z+2z_u)}{1+exp(z+2z_u)}$, we can see more noticeable improvement from using the CBPS moment condition instead of the MLE moment condition. Moreover, the doubly robust estimator can perform substantially better than the proposed IPW estimator at exposure level zero.

Empirical Illustration

I evaluate the effects of China's special economic zones (SEZ) policy using the proposed doubly robust estimators. SEZs are a prominent development strategy that aims to foster agglomeration economies. The benefits of SEZs include corporate tax concessions, customs duty exemptions, discounts on land use fees, and special bank loan programs. SEZs are likely to affect neighboring non-SEZ areas through, for instance, firm relocation or knowledge spillover.

The data for the empirical illustration come from lu2019place. There are five waves of SEZ establishment in China. Each wave is different in nature and targets different regions with earlier waves creating more national-level economic zones.\footnote{As a result, SEZ establishment in China cannot simply be considered as staggered adoption.} Since detailed village level data is only available starting from 2004, lu2019place focus on the latest wave of SEZs established between 2005 and 2008. China established 663 SEZs at the provincial level in 2006, accounting for 42 percent of the country's SEZs. These SEZs cover the coastal, central, and western regions and are considered small-scale regional SEZs. As a result, the policy effect is interpreted as the treatment effect on villages that had not yet been treated prior to this wave. This means that areas covered by zones from earlier waves are not included in the finite population.

lu2019place collect comprehensive data on China's economic zones based on the economic censuses conducted by China's National Bureau of Statistics in 2004 and 2008 covering all manufacturing firms. Consequently, the entire finite population of village-level data in 2004 and 2008 is observed, where 2004 is the period prior to the treatment and 2008 post the treatment. The units of observation are villages, which are the most disaggregated geographical units and smaller than an SEZ. Treated villages are referred to as SEZ villages. Unfortunately, the publicly available data from lu2019place do not contain an identifier of villages nor distances among villages. Nevertheless, I can match counties in which villages are located from separate datasets published by lu2019place. In the organized dataset, there are 3,963 SEZ villages and 99,259 non-SEZ villages. It would be ideal to set certain neighborhood boundaries for each village based on the distance between villages. The exposure mapping is then a function of neighboring villages' SEZ assignment status. In the absence of detailed geographical data, I consider each village's neighborhood to be its corresponding county, with its neighbors being the other villages within the same county. Without distance measures, standard errors are clustered at the county level, which can be considered a special case of spatial-correlation robust inference.

According to equation ((ref)), the outcome variables $Y_{it}$ include the logarithm of capital, employment, and output of firms in a village. The direct treatment variable $W_i$ is equal to one if village $i$ is located within the boundaries of SEZs and zero otherwise. As an illustration, exposure mapping is defined in the following way. I define SEZ ratio as the fraction of other SEZ villages over the total number of other villages in a county. $G_i$ is a binary variable equal to one if this ratio in county $c$ in which village $i$ is located is above the mean of the ratio among all counties.\footnote{Counties on average contain 148.3 villages. On average, there are 5.8 SEZ villages in a county.} In this way, exposure is characterized when there are sufficient number of SEZ villages in a county. It is in line with the initial goal of SEZ establishment, namely promoting agglomeration economies. Villages are considered as intensively exposed to neighbors' economic zones if the corresponding county contains a relatively high fraction of SEZ villages.

There are four baseline village characteristics including logs of a village’s distance from an airport and port, log of the capital-to-labor ratio, and log of the number of firms in the village in 2004. The distance from an airport and port can be considered as fixed attributes, while the capital-to-labor ratio and the number of firms can be treated as stochastic.\footnote{In the conditional inference framework adopted here, we can consider a subset of attributes fixed with the rest stochastic. The identification and estimation procedure stays the same. Only the inference will change in a way that the additional uncertainty resulting from a subset of stochastic attributes needs to be accounted for.} These baseline characteristics and their interactions with the direct treatment variable are included as regressors in moment conditions ((ref)) and ((ref)). The four baseline characteristics and their leave-one-out means at the county level are covariates in moment conditions ((ref)) and ((ref)) for the propensity scores. The spillover effects are estimated analogously using equations ((ref)) and ((ref)) in Appendix (ref).

The first row of Table (ref) below reports DID estimates using the IPW approach in abadie2005semiparametric with the four baseline village characteristics as covariates. Because of potential spillover effects, these canonical estimates are difficult to interpret causally. When interference has been considered, the direct effects are mostly smaller than the canonical DID estimates, especially for SEZ villages with exposure level one. To summarize, SEZ establishment has positive and statistically significant direct effects at the 1% level. This implies that the SEZ villages benefit from the program by gaining investment, employing more labor, and producing more output.

In terms of spillover effects, having relatively high ratio of neighboring SEZ villages does not significantly affect economic activities in SEZ villages. The only exception is for log output, where the spillover effect under more intensive exposure is quite negative and marginally statistically significant. By contrast, non-SEZ villages benefit from SEZ neighbors when there are sufficient number of SEZ villages in the same county. These patterns of direct and spillover effects are not found in lu2019place.

As suggested in Section (ref), I also examined pre-trends with the exposure mapping as well as classical pre-trends for canonical DID estimation. Unfortunately, there is only one economic census period prior to treatment. As a result, the pre-trends are tested with China’s Annual Surveys of Industrial Firms (ASIF) used by lu2019place, which contain more data in the pre-treatment period but only cover firms with relatively large sizes. I use ASIF data from 2004 and 2005. Given this is not the same dataset used for the main analysis, the results presented in Table (ref) are only demonstrative. For canonical DID, only the differential pre-trend for log of employment is statistically significant at the 10% level with small magnitude. None of the doubly robust estimates for placebo direct effects are significant.

Conclusion

I propose doubly robust estimators for the direct treatment effect and spillover effect in a DID context. I later generalize to the case where exposure mapping could be misspecified and interference is not restricted within a fixed boundary of neighborhoods. Given the general spillover effect, one needs to account for spatial correlation when conducting inference. With the entire population observed, the usual spatial-correlation robust variance estimator could be conservative.

I provide identification results of the direct and spillover effect for the IPW estimand, outcome regression estimand, and the doubly robust estimand. From here, researchers can approach these estimands using various parametric, semiparametric, or nonparametric estimation methods. In the current paper, I proved the asymptotic properties of GMM-type parametric estimators as an illustration of estimation. Given the inclusion of neighbors' treatments and attributes in the propensity score and the conditional mean functions, other nonparametric estimation methods are left as future work.