Corrections to: Estimating Average Causal Effects Under General Interference, with Application to a Social Network Experiment

The problem

Proposition 6.2 states that, under Conditions 3, 5 and 6, Wald intervals built from the variance estimator of equation (11) are asymptotically conservative. Condition 3 bounds the potential outcomes and the inverse marginal exposure probabilities $1/\pi_i(d_k)$. Through equations (7) and (10), the variance estimator also divides by the positive joint exposure probabilities $\pi_{ij}(d_k)$, $\pi_{ij}(d_l)$ and $\pi_{ij}(d_k,d_l)$, and the three conditions do not keep these away from zero. The step of Appendix A.7 that takes every $a_{ij}(D_i,D_j)$ to be bounded “by Condition 3” is therefore not justified. The bound $O(N^{-1})$ on the variance of the scaled variance estimator rests on that step, and it can fail. When it fails, the variance estimator can be unbiased and still converge in probability to a value below the true scaled variance, because the estimate is large on rare assignments. The interval then covers less often than its nominal level in the limit. This is so even when Condition 5 is read as a dependency graph of bounded degree, that is, after the correction of P013-E002.

Evidence an exact counterexample

The evidence is an exact counterexample to the printed statement and a derivation of the corrected one.

For $N=2n$, let the units form $n$ pairs that are randomized independently. Within a pair, the assignments $(0,0)$ and $(1,1)$ each have probability $\epsilon_n=1/(4n^3)$, and $(1,0)$ and $(0,1)$ each have probability $q_n=1/2-\epsilon_n$. There is no interference: $D_i=Z_i$, with $d_k=1$ and $d_l=0$. Both potential outcomes equal $1$ for the first unit of each pair and $-1$ for the second, and the average effect is $0$. The outcomes are bounded and every marginal exposure probability is $1/2$, and Condition 3 holds with $c_1=1$ and $c_2=2$. Since the pairs are independent, Condition 5 holds with $m=2$, for a genuine dependency graph. Further, I have $N\operatorname{Var}[\widehat{\tau}_{\mathrm{HT}}]=8(1-2\epsilon_n)\to8$, and Condition 6 holds. All joint exposure probabilities of distinct units are positive, and those of equal exposures within a pair equal $\epsilon_n\to0$.

Let $\widehat{V}_N$ denote the estimator of equation (11). Substituting into (7), (10) and (11), I obtain $N^2\widehat{V}_N=\sum_{b=1}^{n}B_b$, where $B_b=2/\epsilon_n$ if the two units of pair $b$ receive the same assignment and $B_b=16-2/q_n$ otherwise. Then $\mathrm{E}[B_b]=16(1-2\epsilon_n)$, and $\widehat{V}_N$ is exactly unbiased. The probability that some pair receives equal assignments is at most $2n\epsilon_n=1/(2n^2)\to0$, and on the complementary event I have $N\widehat{V}_N=8-1/q_n\to6$. I conclude that $N\widehat{V}_N\to6$ in probability, while $N\operatorname{Var}[\widehat{\tau}_{\mathrm{HT}}]\to8$. Further, I find that

contradicting the concentration bound in Appendix A.7. Since $\sqrt{N}\,\widehat{\tau}_{\mathrm{HT}}$ is asymptotically normal with mean $0$ and variance $8$, the coverage of the nominal $100(1-\alpha)\%$ interval tends to

which is $0.9104$ for $\alpha=0.05$. A variant of the example with $\epsilon=N^{-2}$ has the same limits.

For the corrected statement, suppose in addition that every positive joint exposure probability that enters (7) and (10) is at least $b>0$, and let Condition 5 hold for a symmetric dependency graph with row sums at most $m$. I add every self-edge, writing $\widetilde h_{ij}=1$ when $i=j$ or $h_{ij}=1$. This preserves the dependency condition and bounds every row sum by $m+1$. I write $N^2\widehat{V}_N=\sum_{i,j}B_{ij}(D_i,D_j)$, where $B_{ij}$ collects the terms that involve units $i$ and $j$. A term with positive joint probability is bounded in absolute value by $(1+1/b)c_1^2c_2^2$, and the corrections for zero joint probabilities divide only by marginal probabilities. Then there exists a constant $C$, which does not depend on the population size, such that $|B_{ij}|\le C$ for every pair of units. If $\widetilde h_{ij}=0$, then the exposures of $i$ and $j$ are independent, the joint probabilities are products of marginal probabilities, and $B_{ij}=0$. Then at most $N(m+1)$ of the $B_{ij}$ are nonzero. Two of them can be correlated only if a unit of one is adjacent under $\widetilde h$ to a unit of the other, which leaves at most $4N(m+1)^3$ nonzero covariances, each at most $4C^2$ in absolute value. I obtain

Let $c_N=N\operatorname{Var}[\widehat{\tau}_{\mathrm{HT}}(d_k,d_l)]\to c>0$. Proposition 5.6 implies that $\mathrm{E}[N\widehat{V}_N]\ge c_N$. By Chebyshev’s inequality, I have $P\{N\widehat{V}_N<c_N-\delta\}\to0$ for every $\delta>0$. The statistic $T_N=\sqrt{N}(\widehat{\tau}_{\mathrm{HT}}(d_k,d_l)-\tau(d_k,d_l))/\sqrt{c_N}$ is asymptotically standard normal by the central limit theorem that Appendix A.7 invokes, which involves no joint probability. For $0<\delta<c$, on the event $N\widehat{V}_N\ge c_N-\delta$ the interval covers whenever $|T_N|\le z_{1-\alpha/2}\sqrt{1-\delta/c_N}$. Then the lower limit of the coverage is at least $2\Phi(z_{1-\alpha/2}\sqrt{1-\delta/c})-1$. Letting $\delta\downarrow0$, I obtain the lower bound $1-\alpha$. The argument does not require that $\mathrm{E}[N\widehat{V}_N]$ converge.

The added condition holds automatically for the following designs described in the paper. Consider first Bernoulli assignment at a fixed rate $p$ with first-order interference on a graph of degree at most $\Delta$, a fixed-rate version of the illustration given after Condition 5. A joint exposure event of two units is determined by at most $2(\Delta+1)$ independent assignments, and every positive joint probability is at least $\min(p,1-p)^{2(\Delta+1)}$. In the design of Section 9, a fixed fraction of each school is treated uniformly and schools are randomized independently, and the bound follows from Condition 3. In the block-randomized design of Section 10, the blocks are randomized uniformly at fixed quotas and a school has at most 64 eligible students. Condition 3 then implies that every positive joint probability is at least $c_2^{-2}2^{-129}$. For arbitrary known designs the condition is a further restriction. The counterexample above, which has independent groups of two units, does not satisfy it.

What survives

The estimand, the Horvitz–Thompson point estimator and the variance estimator of equation (11) with its corrections for zero joint probabilities are unchanged. The results on unbiasedness in finite samples, the exact variance and covariance formulas and Proposition 5.6 are unaffected; in the counterexample Proposition 5.6 holds with equality. The point estimator is consistent and asymptotically normal under a dependency graph of bounded degree, and these properties do not involve joint probabilities. With the added condition the coverage conclusion holds, and the condition is automatic for the Bernoulli illustration with a fixed rate and for the designs of Sections 9 and 10 wherever Conditions 3 and 5 hold for them. No entry of Table 1 or Table 2 changes, and no reported interval is shown to undercover. What is lost is the claim for arbitrary known designs under Conditions 3, 5 and 6 alone; the refinements of Section 7, which use equation (11), need the same control in addition to their own assumptions. The example has limiting coverage $0.9104$, and it is not established that coverage can tend to zero.

The problem

Proposition 4.1 bounds the relative bias of the Horvitz–Thompson total computed with simulated exposure probabilities by $(1-\pi_0(d_k))^{R+1}$, and it places no restriction on the signs of the potential outcomes. With independent replication, the expectation identities of Appendix A.3 are correct: the bias equals $-\sum_i y_i(d_k)(1-\pi_i(d_k))^{R+1}$. The last step of the proof bounds the absolute value of this weighted sum by $|y^T(d_k)|(1-\pi_0(d_k))^{R+1}$. That step can fail when the outcomes have mixed signs and the exposure probabilities differ across units, because positive and negative outcomes can cancel in the total $y^T(d_k)$ without cancelling in the weighted sum. When $y^T(d_k)=0$ the relative bias is not defined. The context does not exclude this case, since the paper considers outcomes centred on zero. For mixed signs, the rule $R\ge\log(b)/\log(1-\pi_0(d_k))-1$ does not in general guarantee a relative bias of at most $b$, and no number of replicates chosen from $\pi_0(d_k)$ and $b$ alone can guarantee it for all signed outcomes with nonzero total.

Evidence an exact counterexample

The evidence is an exact counterexample and a derivation of the corrected bound. Let $N=2$, with no interference, and let the two units be exposed to $d_k$ independently with probabilities $\pi_1(d_k)=1/3$ and $\pi_2(d_k)=2/3$. Let $y_1(d_k)=2$ and $y_2(d_k)=-1$, and let $R=1$ replicate be drawn independently of the experiment. Then $y^T(d_k)=1$. Summing over the 16 pairs of an actual and a simulated assignment, I obtain $\mathrm{E}[\widehat{y^T_{\mathrm{HT},R}}(d_k)]=2/9$. The relative bias is $7/9$, and the printed bound is $(1-1/3)^2=4/9$. For the target $b=1/2$ the printed rule selects $R=1$, and the relative bias $7/9$ exceeds the target. With outcomes $(1,-1)$ in the same design, the total is $0$ and the expected estimate is $-1/3$.

For the corrected bound, I write $q_i=(1-\pi_i(d_k))^{R+1}$ and $q_0=(1-\pi_0(d_k))^{R+1}$, so that $0\le q_i\le q_0$. By the identity of Appendix A.3 and the triangle inequality, I have

If all $y_i(d_k)$ have the same sign and $y^T(d_k)\ne0$, then $\sum_i|y_i(d_k)|=|y^T(d_k)|$ and the printed bound follows. In the example, the corrected rule with $b=1/2$ gives $R=4$, for which the corrected bound is $32/81$ and the relative bias is $7/27$.

What survives

The expectation identities of Appendix A.3 hold under independent replication, and the bias of the approximate total still tends to zero as $R$ grows in a fixed finite population. The unbiasedness of the Horvitz–Thompson estimator with exact probabilities (Lemma 4.1), Proposition 4.2 and the later variance and asymptotic results do not use Proposition 4.1. The arithmetic of the printed example is right: the threshold is $10{,}592.99$, and 10,593 is the smallest admissible integer. The outcomes of the simulation (counts multiplied by positive factors) and of the application (a binary outcome) are nonnegative, and the printed bound applies when the corresponding potential-outcome total is nonzero. It need not apply to totals of signed residuals. What changes is the scope of the proposition and of the rule for $R$.

The problem

I could not establish that this implementation was used for computations in the 2017 article. The loop on lines 181–185 redraws any assignment vector that repeats an earlier one. The $R$ vectors kept are then distinct, and they are not independent draws from the Bernoulli design. The function make_exposure_prob estimates the exposure probabilities by the unweighted frequencies over these vectors, as in the estimators $(\widehat{\mathbf{I}}_k\widehat{\mathbf{I}}_k'+\iota_N)/(R+1)$ and $\widehat{\mathbf{I}}_k\widehat{\mathbf{I}}_l'/R$ of Section 3. Unweighted frequencies estimate the probabilities under the design when the draws are independent, or when all assignment vectors are equally likely. Under a Bernoulli design with $p\ne1/2$ the assignment vectors are not equally likely, and the sampler keeps each of them at most once. A vector whose probability exceeds $1/R$ is therefore under-represented in every sample, and less likely vectors take its place. This happens in populations of ordinary size. For $N=100$ and $p=0.01$ there are $2^{100}$ assignment vectors, but the vector with no treated unit has probability $(99/100)^{100}\approx0.366$. Two independent draws coincide with probability about $0.135$. In this design the less likely vectors have more treated units, and the estimated probabilities of treatment are too large for the numbers of replicates examined in the evidence. The estimated probabilities enter the Horvitz–Thompson estimates and their variance estimates as weights. The check on lines 169–177 follows from the same choice. The distinct vectors are exhausted at $R=2^N$, and the limit $R\to\infty$ of Proposition 3.1 is not available to this sampler. The article names Bernoulli designs on a graph of bounded degree as a setting for its asymptotic results (p. 1926). Conditions 1–6 do not exclude the examples below, since a sequence of Bernoulli designs with $p$ fixed at $0.01$ and bounded outcomes can satisfy them, with or without first-order interference on disjoint edges.

Evidence an exact counterexample

The evidence is an exact calculation for a design of ordinary size, a run of the unchanged code on that design, and an exact counterexample with two units.

Let $N=100$ and $p=1/100$, and let the exposure of a unit be its own treatment. Among $R$ distinct assignment vectors, at most one has no treated unit and at most $100$ have exactly one treated unit. Every further vector has at least two. For $R=1000$ the vectors kept then contain at least $100+2\cdot899=1898$ treated entries. I call the fraction of the vectors kept in which a unit is treated its raw frequency of treatment. Averaged over the units, the raw frequency of treatment is at least $1898/(100\cdot1000)=0.01898$. This bound holds for every sample that the sampler can return, whatever the seed. The same count gives the other rows of the table.

Since the sampler treats the units symmetrically, the expected raw frequency of each unit obeys the same bounds. The diagonal entries that the package returns add one to each count and to $R$, and for $R=1000$ their average is at least $1998/100100\approx0.01996$. The raw joint frequency of treatment of a pair of units is the fraction of the vectors kept in which both units are treated. Counting the ordered pairs of treated units in the same way, I find that its average over the ordered pairs of units is at least $1798/(1000\cdot100\cdot99)\approx0.0001816$ for $R=1000$, against $p^2=0.0001$. The value $R=1500$ is the number of replicates in the replication script of the package, and $R=10{,}593$ is the number in the example after Proposition 4.1. The rows apply these numbers to the Bernoulli design. They do not recalculate any result of the article. For small $R$ the bound is not informative: for $R=100$ it is $0.0099$. The distortion is present there as well. For $R=2$ the expected raw frequency of treatment is $0.011057$ to six decimals, which is $10.57\%$ above $0.01$.

I ran the unchanged sampler with $N=100$, $p=0.01$, $R=1000$ and seed $101$. The $1000$ vectors kept are distinct and contain $2240$ treated entries. The average raw frequency of treatment is $0.0224$, and every unit is treated in at least $12$ of the vectors. On a graph of $50$ disjoint edges, the first-order exposure mapping of the package, which has the four exposures of Section 9, gives the following frequencies.

Let every outcome under treatment equal $1$, so that the total under treatment is $100$. Suppose that this total is estimated by the Horvitz–Thompson estimator with the probabilities that the package returns for this run. Given the vectors kept, the expectation of the estimator over the Bernoulli design is $44.31$, with an empty sum counted as zero. With independent draws, the identity of Appendix A.3 implies that the expectation over the replicates and the design is $100[1-(99/100)^{1001}]$.

For an exact counterexample with two units, let $N=2$ and $R=4$, and let each unit be treated with probability $p=1/5$. The sampler then keeps each of the four assignment vectors once. Each unit is treated in half of them, although its probability of treatment under the design is $1/5$. When the exposure of a unit is its own treatment, the package returns $3/5$ and $1/5$ for the marginal and joint probabilities of treatment. Their values under the design are $1/5$ and $1/25$. The distortion does not require that the sampler exhaust the assignment vectors: for $R=2$ the expected raw frequency of treatment is $383/1260$. At the observed assignment $(1,0)$ with both outcomes equal to $1$, estimates() returns a Horvitz–Thompson total of $5/3$ for the treated exposure with variance estimate $10/9$, and the probabilities under the design give $5$ and $20$. I confirmed these values by running the unchanged code.

The size of the distortion depends on the design. For $p=1/2$ all assignment vectors are equally likely, and the raw frequencies are unbiased. If the first $R$ independent draws are distinct, the sampler returns exactly these draws. For $N=1000$, $p=0.01$ and $R=1500$, the probability that two of them coincide is at most $0.0024$, and the expectations of the raw frequencies differ from those under independent draws by at most the same amount.

What survives

Section 3 and Proposition 3.1 are unaffected when the replicates are drawn independently from the randomization plan, and the two estimators of Section 3 are implemented as printed. The designs, the exposure mappings, the estimands and the results that take the exposure probabilities as known are unchanged. The simulation of Section 9 treats a fixed number of units chosen uniformly; the later package’s Table 1 replication script uses the permutation sampler, not this function, and does not establish which program generated the published table. That sampler also rejects repeats by default, but this distortion does not arise when equally likely vectors are sampled without replacement. No entry of Table 1 or Table 2 is shown to be affected by this error. I could not establish when the loop was introduced, who introduced it, or a causal connection between the code and Appendix A.1. I could not identify the code used to compute the exposure probabilities for the application of Section 10 or determine how its replicates were drawn. What is lost is the assurance that the probabilities computed from this sampler are those of the Bernoulli design. In sparse designs such as the one above they are not, and the Horvitz–Thompson estimates and variance estimates that use them inherit the error.

The problem

The article describes the outcome as whether the student wore a program wristband, measured by self-report (p. 1934), and it gives no rule for a student whose report is absent. On the public data release of the field experiment (ICPSR study 37070, version 2), the variable that reproduces Table 2 is an indicator for a recorded orange-wristband report. It is $0$ or $1$ for every student. The general wristband item of the second-wave survey, “I wore a Roots wristband this past week”, has no recorded answer for 342 of the 2,050 students of the reconstructed sample. Of these 342 students, 151 have no recorded second-wave survey, and 191 have a survey and no answer to the item. All 342 have the value $0$ on the indicator. Among the three policies compared here, the twelve point estimates of Table 2 are reproduced to three decimals only when these students are kept in the sample with outcome $0$. Without the 342 students the sample has 1,708 students, and without the 151 alone it has 1,899. In both cases none of the twelve estimates is reproduced. The article does not distinguish an observed negative answer from a zero recorded without an answer. For one student in six in this reconstruction, the zero stands where no answer to the general wristband item is recorded. Under the exposure assumptions for this recorded outcome, the reconstructed point estimates concern effects on a recorded report of wearing, with an absent report counted as zero. In the outcome model below, effects on wearing require a further condition that the article does not state: the students whose reports would be absent must contribute equally to the mean of wearing under the two exposures compared. The sizes that the text quotes, about 30 percentage points for direct exposure and 15 to 16 for indirect exposure, depend on the outcome definition. ICPSR notes that some released variables were changed. The authors' own outcome file is not public, and I could not check whether it codes these students as the public file does.

Evidence a recalculation from the source

The evidence is a recalculation on the public data, a derivation and an exact example. The recalculation uses the released code’s unblocked school-assignment generator, seed 123 and 10,000 draws, with the sampler that R used before version 3.6; the article describes school assignment as blocked. It keeps the eligible students whose five unsmoothed estimated exposure probabilities exceed $0.01$, and it fits the WLS rows with one slope on network degree. This reconstruction has $N=2{,}050$ students and, with the orange-wristband indicator as released, reproduces all twelve point estimates of Table 2 to three decimals. None of the eleven scripts released with the data computes Table 2, and the reconstruction is not the authors' calculation. Recall that the article limits the analysis to students with nonzero probabilities and does not state the cut-off of $0.01$. Whether the published analysis used it is not determined.

Every student removed in the last two columns has outcome $0$. With the exposure probabilities held fixed, the numerators of the Horvitz–Thompson means are then unchanged, and only the divisor $N$ changes. Each HT contrast with these students kept equals $1708/2050=854/1025$ of its value without the 342 students, and $1899/2050$ of its value without the 151. These ratios describe a change of divisor. They do not measure a bias.

To isolate the effect of absent reports, let $R_i(d)$ indicate that usable wearing information for student $i$ is recorded under exposure $d$, let $y_i(d)$ be the wearing outcome, and suppose the analyzed outcome is $g_i(d)=R_i(d)\,y_i(d)$. This model assumes that recorded reports equal wearing. Suppose also that $y_i$ and $R_i$ depend on the assignment only through the exposure and that the exposure probabilities used in the estimator are exact and positive. By Lemma 4.1 applied to $g$, I have

where the estimator is computed from the recorded outcomes and $\mu(d)$ is the mean of the wearing outcomes $y_i(d)$. Subtracting the same identity for $d_{000}$, I obtain

This identity implies that the estimator with absent reports counted as zero is unbiased for the effect on wearing if and only if $L(d_k)=L(d_{000})$. It suffices that under both exposures no student whose report would be absent would wear the wristband. Randomization, bounded outcomes and positive exposure probabilities do not imply this condition.

The exact example has two units with independent fair assignments and no interference. The first unit always reports, and its wearing outcome is $0$ without treatment and $1$ with it. The second unit never reports. Its wearing outcomes without and with treatment may be $(0,0)$, $(0,1)$ or $(1,0)$, and the recorded data are the same in the three cases under all four assignments. The effect on wearing is $1/2$, $1$ or $0$. The estimator with the absent report counted as zero has expectation $1/2$ in each case, and its bias is $0$, $-1/2$ or $1/2$. The estimator on the reporting unit alone has expectation $1$. The recorded data therefore do not determine the sign of the bias.

Two limits apply. In the whole public file, 4,108 of the 24,471 records have no recorded answer to the general item, and four of them have the value $1$ on the orange-wristband indicator. An absent answer to the general item does not prove that no information on colour was recorded. The released scripts do not show how the zeros arose: the cleaning script sets codes above 4 to missing and otherwise leaves the indicator as it is.

What survives

All twelve point estimates of Table 2 are reproduced on the public data under this coding, and so is the mean of zero in the no-exposure condition (0 of 1,034 students). The broad order of the effects is the same under each of the three treatments and for each estimator: school-only exposure, then indirect exposure, then the two direct exposures. This outcome-coding issue does not affect the estimators or the results of Sections 3 to 8, and the simulation of Section 9 has no absent outcomes. What changes is the definition of the outcome, and with it the estimand of Table 2 and the meaning of the sizes quoted on pp. 1936–1937. No bias in the effects on wearing is identified. This finding concerns point estimates only; a different treatment of these students would need its own standard errors.

The problem

As printed, the two sets of exposures must be independent only when no edge of $h$ lies inside either set or leaves either set for the outside; edges between the two sets are allowed. This is not the condition that no edge joins the sets, and the requirement can be vacuous. If $h$ is the identity matrix, or the adjacency matrix of disjoint triangles, the premise holds only when both sets are empty. Every joint distribution of the exposures then satisfies the definition, with row sums $1$ or $2$. The printed definition therefore supports neither the statement on p. 1926 that Condition 5 subsumes Condition 4 nor the central limit theorem used for Proposition 6.2. The sentence before the definition, which announces “joint independence of observations rather than pairwise independence”, and the description on p. 1926 of a dependency graph of limited degree show what is meant.

Evidence an exact counterexample

The evidence is two exact counterexamples to the printed definition. First, let every unit receive the same fair Bernoulli exposure, $D_i=X$. With $h$ the identity matrix, the printed definition is satisfied and Condition 5 holds with $m=1$. But every pair of units is dependent, which implies that $\sum_{i=1}^{N}\sum_{j=1}^{N}g_{ij}=N^2$, and Condition 4 fails. Second, for $N=36r^2$ let the number $S$ of treated units take the values $18r^2-9r$, $18r^2$ and $18r^2+9r$ with probabilities $1/18$, $8/9$ and $1/18$, let the treated units be chosen uniformly given $S$, and let $D_i=Z_i$, $y_i(1)=1$ and $y_i(0)=0$. Distinct units are pairwise independent, and Conditions 1–4 hold. Condition 5 holds as printed, with $h$ the adjacency matrix of disjoint triangles and $m=2$. Condition 6 holds, since $N\operatorname{Var}[\widehat{\tau}_{\mathrm{HT}}]=1$ for every $r$. The estimator standardized by its true standard deviation takes the values $-3$, $0$ and $3$ with probabilities $1/18$, $8/9$ and $1/18$ for every $r$, which is not a normal limit. The interval of Proposition 6.2 with critical value $2$ has coverage exactly $8/9$ for every $r\ge5$, which is below its nominal level. Neither example satisfies the corrected definition.

Applying the corrected definition to two singletons, I find that any two dependent units are adjacent. Since the row sums are at most $m$, I obtain $\sum_{i}\sum_{j}g_{ij}\le N(m+1)=O(N)$, which is the statement on p. 1926 that Condition 5 subsumes Condition 4.

What survives

The unbiasedness results, the exact variance formulas, Proposition 5.6 and Proposition 6.1, which uses Conditions 3 and 4, do not involve $h_{ij}$. With the four signs corrected the definition is the dependency-graph condition that p. 1926 describes and that Appendix A.7 uses, and no reported number depends on it. The correction does not by itself settle Proposition 6.2, which is the subject of P013-E001.