Corrections to: Design-Based Inference for Spatial Experiments under Unknown Interference

The problem

Throughout, $d$ is fixed, and I write $R_i=\mu_i(\mathbf Y;d)$, $p_i=p(O_i\mid\theta_0)$, $m_{ia}=E[\mu_i(a;d)]$ for $a\in\{0,1\}$, $\tau_i=m_{i1}-m_{i0}$ and $\delta_i=O_i'(\beta_{1,N}-\beta_{0,N})$. Let

so that $\beta_{a,N}=H_N^{-1}N^{-1}\sum_ip_i(1-p_i)m_{ia}O_i$. The upper limits $n$ in the two residual sums of the proposition are read as $N$, and “C1” in its hypotheses is the condition C 1 of A.6 (see P004-N016).

The first-order term of the expansion (42) has the coefficient $N^{-1}\sum_iZ_iR_i(1-p_i)O_i'/p_i$, whose expectation is $N^{-1}\sum_i(1-p_i)m_{i1}O_i'$. Equation (43) concerns a different sum, in which $Z_iR_i/p_i$ is multiplied by $p_i(1-p_i)O_i$ and not by $(1-p_i)O_i$, and it is this sum that converges to $H_N\beta_{1,N}$. The display after “Hence (42) becomes” puts $\beta_{1,N}$ in place of the coefficient of (42) and moves a factor $1/p_i$ into the score sum. Lemma A.6 expresses $\widehat\theta_{MLE}-\theta_0$ through $H_N^{-1}S_N$, in which every node has the same weight, so no factor $1/p_i$ can arise there. The two expressions agree when all $p_i$ are equal and differ in general. Expanding the definition, I find that the printed linear form is

whereas the IPW estimator satisfies $\widehat\tau_{IPW}=\widehat\tau_{HT}-S_N'\gamma_N+O_p(N^{-1})$ with $\gamma_N$ as in the headline. In general the two differ at order $N^{-1/2}$, and then (36) fails. The variance calculation (37)–(40) is based on the printed linear form rather than the correct first-order form of the IPW estimator; the HAC and SAH estimators built on the residuals $R_i-O_i'\widehat\beta_{a,N}$ can estimate the wrong quantity, and the normal limit and the two coverage statements for the IPW estimator do not follow. As printed, Appendix A.6 does not supply the inferential theory that Section 6.4 announces.

Evidence an exact counterexample

The evidence is an exact counterexample to (36), computed in rational arithmetic, and a derivation of the corrected statement.

The expansion of the IPW estimator. For fixed outcomes and assignments,

Since $Z_i$ is independent of the other assignments, I have $E[Z_iR_i]=p_im_{i1}$ and $E[(1-Z_i)R_i]=(1-p_i)m_{i0}$. The expected gradient of $\widehat\tau_{IPW}$ at $\theta_0$ is therefore $-D_N$, where $D_N=N^{-1}\sum_i[(1-p_i)m_{i1}+p_im_{i0}]O_i$. Under C2, C3, C4, C6 and C 1, the summands are bounded and each shares assignments with a bounded number of others, which implies that the gradient differs from $-D_N$ by $O_p(N^{-1/2})$. The second derivatives are bounded near $\theta_0$. Inserting the expansion $\widehat\theta_{MLE}-\theta_0=H_N^{-1}S_N+O_p(N^{-1})$ from the logistic score equation, with $\theta_0$ an interior point (P004-E003), I obtain $\widehat\tau_{IPW}=\widehat\tau_{HT}-S_N'\gamma_N+O_p(N^{-1})$ with $\gamma_N=H_N^{-1}D_N$.

The counterexample. Let $N=3K$. The nodes are the points $(10i,0)$, $i=1,\dots,N$, and the outcome points are the same points, so that at $d=0$ the circle average of a node is the outcome at its own point. The nodes have the types $t=-1,0,1$, repeated in this order. A node of type $t$ has $O_t=(1,t)'$ and the outcome $q_t$ at its own point under every assignment, with $q_{-1}=1$, $q_0=21/20$ and $q_1=1$. The assignments are independent with $\theta_0=(0,\log2)'$ and $\Theta_0=[-2,2]^2$. Then $p_t=2^t/(1+2^t)$, that is, $1/3$, $1/2$ and $2/3$. Every hypothesis of the proposition holds. The outcomes are bounded. There is no interference, and $\mathcal B(i;0)=\{i\}$. The nodes are 10 apart, the probabilities lie in $(1/10,9/10)$, the Gram matrix of the confounders is $\operatorname{diag}(1,2/3)$, and $\theta_0$ is an interior point of $\Theta_0$. Further, $\tau_i=0$ and $\beta_{1,N}=\beta_{0,N}$, which implies that the sum in C 2 is 0, in the printed form and in the corrected one. The variance condition of the proposition also holds, since $V_P>0$ below. A direct computation gives

Since the two arms have the same outcome, both linear forms are averages over the nodes of multiples of $Z_i-p_i$. At a node of type $t$, the first-order form of the IPW estimator has the coefficient $a_t$ and the printed linear form has the coefficient $b_t$, where

Note that the two sets of coefficients have opposite signs. Writing $w_t=p_t(1-p_t)$, I find that $N$ times the variances of the two linear forms and their difference are

for the first-order form of the IPW estimator, for the printed linear form, and for their difference, respectively. By the expansion above, $\sqrt N(\widehat\tau_{IPW}-\widehat\tau^{L}_{IPW})=U_K+o_p(1)$, where $U_K$ is a normalized sum of $K$ independent and identically distributed terms, one for each block of three nodes. Its second moment is $V_D$ and its fourth moment is at most $2205448981/2250000000000$ for every $K$. By the Cauchy–Schwarz inequality, I have

Hence the difference in (36) is not $o_p(N^{-1/2})$. Further, $\sqrt N\,\widehat\tau_{IPW}$ converges in distribution to a normal variable with variance $3/625$, and the AME is 0. Then the ratio that Proposition 4 standardizes by the printed variance has limiting variance $V_A/V_P=3600/1753$, not 1. In this example the printed HAC estimator, multiplied by $N$, converges to $V_P$. A normal interval with critical value $z$ then has limiting coverage $2\Phi(z\sqrt{1753/3600})-1$, which is about $0.83$ at the nominal level $0.95$ and about $0.75$ at the nominal level $0.90$.

The corrected linear form. With $L_i=O_i'\gamma_N$, the fits $f_{i1}$ and $f_{i0}$ of the corrected statement satisfy $f_{i1}-f_{i0}=\delta_i$ and $f_{i1}/p_i+f_{i0}/(1-p_i)=L_i$. For every assignment,

Summing over $i$, I obtain $N^{-1}\sum_i\delta_i+N^{-1}\sum_iT_i=\widehat\tau_{HT}-S_N'\gamma_N$ as an identity, with $T_i$ as in the corrected statement, and $E[T_i]=\tau_i-\delta_i$. The corrected form of (36) then follows from the expansion of the IPW estimator.

The variance estimators. Let $V_N$ be the variance of $N^{-1}\sum_iT_i$, and let $W$ be the sum of the products $T_iT_j$ over $i$ and $j\in\mathcal B(i;d)$, divided by $N^2$. By C3, the covariance of $T_i$ and $T_j$ vanishes when $j\notin\mathcal B(i;d)$, which implies that

The last sum is nonnegative under the corrected C 2, which is the homophily condition of the paper with $O_j$ in its second factor. For the SAH estimator, the inequality $\operatorname{Cov}(T_i,T_j)\le(\operatorname{Var}T_i+\operatorname{Var}T_j)/2$ and the symmetry of the neighbourhoods give $\operatorname{Var}(\sum_iT_i)\le\sum_ic_i(d)E[T_i^2]$, with no homophily condition, and $T_i^2$ is the bracket in the corrected SAH estimator. The $T_i$ are bounded, and each depends on the assignments of a bounded number of nodes. By Chebyshev’s inequality, $N$ times each of the two sums concentrates around its expectation. The estimated fits $\widehat f_{ia}$ are consistent uniformly in $i$, since $\widehat\theta_{MLE}$, $\widehat\gamma$ and $\widehat\beta_{a,N}$ have errors of order $N^{-1/2}$. Then replacing $T_i$ by $\widehat T_i$ changes $N$ times each sum by $o_p(1)$. Finally, a sum of bounded variables whose dependency graph has bounded degree is asymptotically normal when its variance is of order $N$, which is the condition that $NV_N$ is bounded below.

Three remarks on the correction. First, the index in C 2 matters. In a second exact example, with three types on a path, the printed sum of C 2 equals $7/339>0$ and the corrected sum equals $-775/38307<0$. The corrected HAC estimator, multiplied by $N$, then converges to $6053/153228$, below the value $27/452$ of $NV_N$. Second, a positive printed variance does not imply that $V_N$ is positive. With $q=(1,9/8,1)$ in the counterexample, the first-order form of the IPW estimator vanishes identically, while $N$ times the printed variance is $1753/120000$. Third, the fits are not unique. The fits $g_{i1}=p_iO_i'\gamma_1$ and $g_{i0}=(1-p_i)O_i'\gamma_0$, with $\gamma_1=H_N^{-1}N^{-1}\sum_i(1-p_i)m_{i1}O_i$ and $\gamma_0=H_N^{-1}N^{-1}\sum_ip_im_{i0}O_i$, yield the same linear form. With them, however, the bound for the HAC estimator needs homophily of $\tau_i-g_{i1}+g_{i0}$, which C 2 does not imply: in a third example the sum for $g$ equals $-9/2000$ while the sum of C 2 is 0. The fits $f$ avoid this additional condition.

What is not established. Equation (36) does not fail in every design with unequal probabilities. With two types of node, a saturated logistic model and outcomes that depend only on the type and the arm, (36) holds although the probabilities differ. No published figure, estimate or interval of the observational application has been shown to change, and the corrected conditions have not been checked for that application.

What survives

The definition of the AME, its identification under independent assignment, and the Horvitz–Thompson estimator with known probabilities (with the weight corrected in P004-E004) are not affected. The theory of Sections 4 and 5 for experiments with a constant probability does not use Proposition 4 of A.6. The IPW estimator remains consistent for the AME. The defect concerns its first-order approximation and the variance derived from it. The printed expansion is correct when the logistic model has an intercept only, since all probabilities are then equal. What changes is the inference for the IPW estimator with a fitted logistic model: the linear form, its variance, the two variance estimators and the conditions are replaced as below. The correction keeps the estimand, the estimator and the homophily condition in its intended form. It adds two regularity conditions, an interior $\theta_0$ and a positive lower bound on $N$ times the corrected variance.

The problem

The residual $\widehat\epsilon_i$ defined above the display already contains the inverse propensity weight. I write $\widehat p_i=p(O_i|\widehat\theta_{\mathrm{MLE}})$ and $\widehat r_{ia}=\mu_i(\mathbf Y;d)-O_i'\widehat\beta_{a,N}$, and I read the unbound $O_j$ of the residual as $O_i$ (see P004-E015). For a treated node, $\widehat\epsilon_i=(N_1/N)\,\widehat r_{i1}/\widehat p_i$. The SAH display divides $\widehat\epsilon_i^2$ by $\widehat p_i^2$ again. Substituting the residual into the display, I obtain

whereas the analogue of (20) has $Z_i\widehat r_{i1}^2/\widehat p_i^{\,2}$ and $(1-Z_i)\widehat r_{i0}^2/(1-\widehat p_i)^2$ in the brackets. The two differ by the factors $(N_1/(N\widehat p_i))^2$ and $(N_0/(N(1-\widehat p_i)))^2$. These factors equal $1$ for an intercept-only logistic model with an interior maximum-likelihood estimate, where $\widehat p_i=N_1/N$, and they can be far below $1$ when the propensities vary across nodes. The example below shows that the printed estimator need not be a conservative bound for the variance, and the intended IPW interval of Proposition 4(i) of A.6 can cover less often than its nominal level. A.6 defines no unweighted residual for the SAH display, and the page offers no other reading of $\widehat\epsilon_i$ there.

Evidence an exact counterexample

The evidence is an exact counterexample. Let $N$ be a multiple of $8$. Under independent assignment, a quarter of the nodes have propensity $1/2$ and treated circle averages $a_i=+1$ or $-1$, in equal numbers; the others have propensity $1/10$ and treated circle average $0$. All control circle averages are $0$, there is no interference at $d=0$, and $c_i(0)=1$. The covariates are the two group indicators, for which the logistic model holds, and the parameter space is the compact set on which the fitted propensities lie in $[1/20,19/20]$. The AME equals $0$ and the coefficients $\beta_{1,N}$ and $\beta_{0,N}$ equal $0$. Let $S$ be the sum of $a_iZ_i$ over the first group and $\widehat p_H$ the fitted propensity of that group. Then $\widehat\tau_{\mathrm{IPW}}=S/(N\widehat p_H)$. A direct computation gives $\widehat\tau_{\mathrm{IPW}}-2S/N=O_p(N^{-1})$ and $N\operatorname{Var}(\widehat\tau_{\mathrm{IPW}})\to1/4$. The printed estimator satisfies $N\widehat V^{\mathrm{obs}}_{\mathrm{SAH}}\to2/25$ in probability, a ratio of $8/25$ to the variance. For the corrected estimator the limit is $1/2$. By an exact calculation over the assignment distribution at $N=32$, the variance of $\widehat\tau_{\mathrm{IPW}}$ is $0.01212$, the expectation of the printed estimator is $0.00256$, and that of the corrected estimator is $0.01298$. At $N=64$ the three values are $0.00483$, $0.00145$ and $0.00764$. The nominal $95\%$ interval built on the printed estimator has limiting coverage $2\Phi(1.96\sqrt{8/25})-1\approx0.73$.

What survives

The experimental SAH estimator (20) of the article, which is (24) of Supplement A, uses unweighted residuals with a common $p$; it is a different formula and is unaffected. The experimental simulations, the AME, the point estimators and the randomization tests are also unaffected by this error, and the observational HAC display is not touched by this point. The forest application of Section 8.2 reports HAC intervals, not observational SAH intervals; this error does not establish that any reported interval changes. What changes is the SAH display of A.6 and, with it, the intended IPW interval of Proposition 4(i) there. The corrected display removes the double weighting and adds no assumption. The general interval also requires the corrected linearization and fitted values in P004-E001.