Corrections to: Estimating Causal Effects in Experiments with Network Spillovers

The problem

Recall that a variance estimator is conservative if its expectation over the randomization is at least the variance of the effect estimator. The Horvitz–Thompson variance estimator has this property when it is computed with the exact exposure probabilities of the design and every unit has positive probability of both contrasted exposures. The chapter computes it with exposure probabilities simulated from $R=10{,}000$ assignments, and the property does not carry over. The estimator divides by the simulated probability that two units are jointly in given exposures, and for rare joint exposures this frequency rests on few assignments. In the regenerated network, among pairs of units that are both in isolated direct exposure, the smallest positive count in the 10,000 simulated assignments is 1.

Finite simulation introduces numerical error here, just as a finite number of bootstrap resamples approximates the ideal bootstrap calculation. This familiar issue can be addressed by controlling simulation error; a draw count alone does not show that the resulting Horvitz–Thompson variance estimates remain conservative.

The worked example uses an assignment from the probability simulation. Both calls to make_tr_vec_permutation set the seed 4224, which implies that the realized assignment of p. 217 is the first of the 10,000 simulated assignments of p. 214. The text does not say so; it says “We also need to assign treatment”. In an experiment analyzed by this recipe, the realized assignment is ordinarily drawn independently of the simulated assignments. The two cases differ in the small-world network regenerated from the printed code. When the realized assignment is uniform among the simulated assignments, the expectation of the Horvitz–Thompson variance estimator is 1.673, 1.687 and 0.996 times the variance of the Horvitz–Thompson estimator for the contrasts $(d_{11},d_{00})$, $(d_{10},d_{00})$ and $(d_{01},d_{00})$. When the realized assignment is independent of them, the expectation is 0.617, 0.035 and 0.995 times that variance, and the observed shares of negative variance estimates are 9%, 40% and 0%. Inclusion raises the first two ratios above 1 but leaves the third below 1. Applied to an experiment, the recipe can return a Horvitz–Thompson variance estimate that is a small fraction of the variance of the estimator, or negative. Exact counterexamples show that neither a large $R$ nor a realized assignment among the simulated ones guarantees that the estimator is conservative.

Evidence an exact counterexample

I checked the claim with a recalculation in the chapter’s small-world design, exact counterexamples and a derivation, supplemented by a second recorded 10,000-assignment table and a finite-bootstrap comparison.

The small-world design. The small-world network was regenerated with the printed code ($N=1000$, nei = 5, rewiring probability 0.1, seed 42), the 10,000 simulated assignments (100 of 1000 units treated, seed 4224) and the potential outcomes (seed 547). The package then computed the exposure probabilities. The regenerated network need not be the chapter’s realization: its degree variance is 1.87 and its skewness is 0.07, against 2 and 0.12 on p. 214. Conditional on the 10,000 simulated assignments, the expectation of each variance estimator and the variance of each effect estimator were computed under two distributions of the realized assignment. In the first, the realized assignment is uniform on the 10,000 simulated assignments. In the second, it is drawn by complete randomization, independently of them. For the Horvitz–Thompson estimator, both moments are exact sums over the joint exposure probabilities of the distribution in question, with an empty exposure sum defined as zero; the package instead returns an undefined value when an exposure group is absent. The shares of negative variance estimates, the coverage and the Hajek ratios use all 10,000 simulated assignments in the first case and are estimated from 30,000 independent assignments in the second; all evaluated assignments have all four exposure groups present.

Averaging also over the simulated assignments by numerical integration, with a negligible bound on the approximation of sampling without replacement, the Horvitz–Thompson ratios are 1.664, 1.667 and 0.996 in the first case and 0.632, 0.072 and 0.980 in the second. With exact exposure probabilities the ratios are 1.745, 1.696 and 1.008. With exact marginal probabilities and simulated joint probabilities they are 0.648, 0.020 and 1.006 for an independent realized assignment. This shows that the shortfall persists with exact marginal probabilities. The ratios for the degree-corrected block model were not computed.

Exact counterexamples. In each, the exposure mapping is the chapter’s, treatment is by complete randomization, and every unit has positive probability of each of the four exposures. The contrast is $(d_{01},d_{00})$. The simulated assignments are a uniform sample of distinct assignments, as make_tr_vec_permutation draws them, and the moments are exact rational numbers.

1. Twelve units form three disjoint cliques of four, two units are treated, all potential outcomes equal 1, and the probabilities are simulated from 65 of the 66 assignments. Over all such samples, the variance of the Horvitz–Thompson estimator is 0.48766 and the expected variance estimate is 0.44710 when the realized assignment is independent of the sample. They are 0.48734 and 0.44874 when it is one of the 65. With exact probabilities both equal $6050/11907$.
2. Ten units form a cycle, two units are treated, $y_i(d_{00})=1$ and the other potential outcomes are 0, and the probabilities are simulated from 20 of the 45 assignments. Over all such samples, the expected variance estimate is $-0.0711$ against a variance of 0.0370 when the realized assignment is independent of the sample, and 0.0173 against 0.0334 when it is one of the 20. With exact probabilities the expected variance estimate is $32/245$ and the variance is $3/98$.
3. Forty units form eight disjoint cliques of five, four units are treated, all potential outcomes equal 1, and the probabilities are simulated from one recorded sample of $R=10{,}000$ of the 91,390 assignments. Given this sample, the expected variance estimate is 0.12980 against a variance of 0.13262 when the realized assignment is independent of the sample, and 0.13013 against 0.13036 when it is one of the 10,000.

For a second probability table for the same forty-unit design, using Python seed 20261019 for its 10,000 distinct assignments, I obtain an expected variance estimate of $0.130527822999$ against variance $0.132664519808$ for an independent realized assignment, and $0.130907208515$ against $0.131130209715$ when the realized assignment is included. Every positive relevant joint has a positive empirical count, with minimum 845; exact probabilities give equality at $7430620902683355/55988922226386688\approx0.132715912491$. These moments condition on this probability table rather than averaging over all 10,000-assignment samples. Both tables concern this forty-unit design rather than the chapter’s published networks, and neither calculation establishes a numerical error in a published table entry.

Finite simulation and accuracy. The literal bootstrap analogy concerns a computational approximation. For fixed observed data $(0,1)$, the ideal variance of a two-observation bootstrap mean is $1/8$. With $B\ge2$ independent bootstrap means, their sample variance with denominator $B-1$ has expectation $1/8$, but is zero with probability $2(1/4)^B+(1/2)^B>0$; denominator $B$ gives expectation $(B-1)/(8B)$. This illustrates finite computational error without equating its bias mechanism to the nonlinear Horvitz–Thompson plug-in or establishing bootstrap validity for an original sampling distribution.

The package samples distinct assignments by default. If $M$ assignments are possible and an exposure event has exact probability $q$, its count $C$ in $R$ uniform draws without replacement is hypergeometric, with

Thus its relative standard deviation depends on $Rq$, not on $R$ alone. At $R=10{,}000$ and negligible sampling fraction, the illustrative probabilities $.01$, $.001$ and $.0001$ give relative standard deviations about $.0995$, $.3161$ and $.9999$; these are not measurements of the chapter’s graph. The routine also uses $(C_i+1)/(R+1)$ for marginals, $C_{ij}/(R+1)$ for same-exposure off-diagonal joints and $C_{ij}/R$ for cross-exposure joints. Smoothing and these distinct normalizers must be included in an accuracy argument, along with support and zero-count handling. A bound holding with high probability does not alone prove conservative expectation without control of the failure event and any needed margin or inflation.

For independent with-replacement simulation from a fixed finite design, fixed finite outcomes and mapping, positive contrast marginals and finitely many relevant positive joints, the strong law and vanishing smoothing give almost-sure convergence of the coefficient formulas and their fixed-design conditional moments to their exact-probability counterparts once support is correctly distinguished. This is no finite-$R$ inequality, no uniform guarantee over changing designs or shrinking positive joints, no growing-$N$ result, and no coverage guarantee. For the default without-replacement scheme, $R\le M$ at fixed design; even a full census retains the routine’s $R+1$ smoothing. Analytic probabilities or unsmoothed census frequencies are exact alternatives. Sufficiently accurate simulation with justified error accounting can also be adequate.

Derivation for exact probabilities. Let $\pi_i(d)=\Pr(D_i=d)$ and $\pi_{ij}(d,d')=\Pr(D_i=d,\,D_j=d')$, with $\pi_{ii}(d,d)=\pi_i(d)$, and write $I_i(d)=\mathbf{1}\{D_i=d\}$ and $Y_i=y_i(D_i)$. The Horvitz–Thompson estimator of the total under $d$ is $\widehat T_d=\sum_i I_i(d)\,Y_i/\pi_i(d)$, and $\widehat\tau_{HT}(d,d')=(\widehat T_d-\widehat T_{d'})/N$. The variance estimator is $\widehat V_{HT}(d,d')=N^{-2}\bigl(\widehat W_d+\widehat W_{d'}-2\,\widehat C_{dd'}\bigr)$, where

Each term of the two sums over positive joint probabilities has expectation equal to the corresponding term of $\operatorname{Var}(\widehat T_d)$ or of $\operatorname{Cov}(\widehat T_d,\widehat T_{d'})$. For a pair with $\pi_{ij}(d,d)=0$, the term of the variance is $-y_i(d)\,y_j(d)$, and the second sum of $\widehat W_d$ has expectation $\bigl(y_i(d)^2+y_j(d)^2\bigr)/2$ for that pair. For a pair with $\pi_{ij}(d,d')=0$, the term of the covariance is $-y_i(d)\,y_j(d')$, and the second sum of $\widehat C_{dd'}$ has expectation $-\bigl(y_i(d)^2+y_j(d')^2\bigr)/2$. Collecting terms gives

The first step of this computation uses that the probabilities inside $\widehat V_{HT}$ are those of the distribution of the realized assignment. With simulated probabilities in their place the identity fails, as the counterexamples show.

What survives

The exposure mapping, the estimands and the Horvitz–Thompson variance estimator computed with exact exposure probabilities under positive marginal support are unaffected. No entry of Table 16.1 is shown to be wrong. Its Horvitz–Thompson standard errors and $t$-statistics were computed with the realized assignment among the simulated assignments, and they describe that case. Figure 16.3 reports the bias and the RMSE of the point estimators and does not depend on the variance estimators. In the regenerated small-world design, the Hajek variance ratios are approximately between 1.01 and 1.82 in the first case and are estimated to lie approximately between 1.02 and 1.82 in the independent case. These results concern that design and give no general guarantee. What is lost is the assurance that the Horvitz–Thompson variance estimator computed with simulated exposure probabilities is conservative.

The problem

The chapter defines the Bernoulli design and gives a recipe for exposure probabilities: simulate assignments from the design, and count how often each unit falls in each exposure condition. Independent draws suffice, but dependent rows with the right marginal law or suitable sampling weights can also yield unbiased estimates. The function make_tr_vec_bernoulli of the package does not generally preserve that marginal law. A proposed profile that equals an earlier row is discarded and drawn again, both in the unit-level sampler (lines 181–185) and in the cluster sampler (lines 117–152).

Let $q(\mathbf{z})=p^{\sum_i z_i}(1-p)^{N-\sum_i z_i}$ be the probability of the profile $\mathbf{z}$ under the unit-level design, and let $S$ be the set of rows already accepted, with $q(S)=\sum_{\mathbf{s}\in S}q(\mathbf{s})$. While $q(S)<1$, the next row to be accepted equals a profile $\mathbf{z}\notin S$ with probability $q(\mathbf{z})/(1-q(S))$, and it never equals a member of $S$. The first row has the law $q$; later rows are conditioned to exclude earlier profiles. The routine make_exposure_prob weights the returned rows equally. This can overrepresent low-probability profiles and distort exposure probabilities. Inverse probability weights and variance estimates that use the distorted probabilities can be affected, and the package gives no warning. The documentation says correctly that sampling is without replacement; it does not explain that the returned rows need not retain the design’s marginal law with parameter p.

The distortion is not confined to small examples. Among distinct profiles, at most one has no treated unit and at most $N$ have exactly one treated unit. With $N=100$, $p=0.01$ and $R=1000$, at least $899$ of the returned rows have two or more treated units, and the rows contain at least $100+2\cdot 899=1898$ treated entries among $100{,}000$. The treatment frequency is at least $0.01898$, where the design has $0.01$. When repeated profiles are unlikely over the $R$ draws, the loop usually has no effect.

Evidence an exact counterexample

I checked the counterexample in exact rational arithmetic and by running the unmodified package functions. Since each proposal has the law $q$, for $\mathbf{z}\notin S$ and $q(S)<1$ the probability that $\mathbf{z}$ is the next row to be accepted is

Let the network consist of two units joined by an edge, so that neither is an isolate, and let $p=1/5$. The profiles $(0,0)$, $(0,1)$, $(1,0)$ and $(1,1)$ have probabilities $16/25$, $4/25$, $4/25$ and $1/25$, and they place unit 1 in the conditions no, ind1, isol_dir and dir_ind1, respectively. With $R=4$, the largest value that the function accepts when $N=2$, each profile is returned exactly once. The routine make_exposure_prob reports the marginal probability of a condition as $(C+1)/(R+1)$, where $C$ is the number of rows that place the unit in the condition.

For two different conditions, the routine reports the plain frequency. The reported probability that unit 1 is in isol_dir and unit 2 is in ind1 equals $1/4$, where the design has $4/25$. Since no adjustment enters this entry, the difference comes from the discarded repeats alone.

With $R=2$ on the same network, there are fewer rows than profiles. In the first row, unit 1 is treated with probability $1/5$. Conditioning on the first row, which is $(0,0)$, $(0,1)$, $(1,0)$ or $(1,1)$, I obtain for the second row

For the expected frequency of treatment of unit 1 over the two rows, I find $(1/5+257/630)/2=383/1260$, where the design has $1/5$. The cluster sampler puts both units of this network in one cluster. With $R=2$, it must return $(0,0)$ and $(1,1)$ once each. The condition dir_ind1 then has frequency $1/2$, where the design has $1/5$. With a single unit, $p=1/5$ and $R=2$, the function returns one $0$ and one $1$, and the unit is treated in the second row with probability $4/5$.

What survives

This defect does not alter the chapter’s definition of the Bernoulli design, its exposure mappings or its design-based identities when the exposure probabilities are those of the design. The worked example summarized in Table 16.1 uses the function make_tr_vec_permutation, which carries out complete random assignment, and the simulation in Figure 16.3 continues that design. Under a design that is uniform on its support, each row remains marginally uniform when repeats are discarded, so this point does not implicate those results. A single row ($R=1$) from either Bernoulli sampler is unaffected by duplicate rejection. What can change is the calculation of exposure probabilities from several rows of make_tr_vec_bernoulli under a nonuniform design, with the weights and variance estimates that use them.

The problem

The chapter defines the effect $\tau(d,d')=N^{-1}\sum_{i=1}^N(y_i(d)-y_i(d'))$ as an average over the whole population, and it says on p. 212 that the exposure probabilities give unbiased Horvitz–Thompson estimators of it. This holds only under a condition that the chapter does not state. Let $\pi_i(d)=\Pr(f(\mathbf{Z},\theta_i)=d)$ be the probability that unit $i$ receives exposure $d$ under the design. The condition is that $\pi_i(d)>0$ and $\pi_i(d')>0$ for every unit $i$. If $\pi_j(d)=0$ for some unit $j$, then $y_j(d)$ is never observed. It can be changed without changing the distribution of anything that is observed, while $\tau(d,d')$ changes with it. In that case no estimator is unbiased for $\tau(d,d')$ over all potential outcomes, whatever weights it uses.

Random assignment does not guarantee the condition: there exist networks and designs such that some unit has probability zero of an exposure. Under complete random assignment, for example, a unit with enough peers cannot have all of them untreated. The only related guidance in the chapter is Note 1, which recommends that isolates be removed. The note does not claim that this is sufficient, and it is not sufficient: the example below has no isolates. The sentence on p. 216 needs the same condition when it is read as a general rule. The condition is the usual proviso of the Horvitz–Thompson result. The error is that the chapter states the result in general without it.

Evidence an exact counterexample

The evidence is an exact counterexample, obtained by enumerating all assignments in rational arithmetic, and a derivation of the corrected statement. I consider four units in a star: unit 1 is joined to units 2, 3 and 4, and there are no other edges. The degrees are $(3,1,1,1)$, so no unit is an isolate. Exactly two units are treated, chosen uniformly from the six pairs. This is complete random assignment with $N=4$ and $N_1=2$. Under the chapter’s mapping with four conditions, the exposure probabilities are as follows.

Unit 1 never receives $d_{10}$ or $d_{00}$. Indeed, if it is treated, then one of its peers is treated, and if it is untreated, then two of its peers are treated. I compare two schedules of potential outcomes. In the first, $y_i(d)=1$ for every unit and every exposure. The second differs from the first only in that $y_1(d_{00})=5$. Under each of the six assignments, the observed outcomes are $(1,1,1,1)$ in both schedules. The joint distribution of assignments, exposures and outcomes is then the same in both, and any estimator has the same expectation under both. However, $\tau(d_{11},d_{00})=0$ in the first schedule and $\tau(d_{11},d_{00})=(1-5)/4=-1$ in the second. No estimator is unbiased under both. If each inverse-weighted sum is taken over the units observed in the exposure and divided by $N=4$, then the expected contrast is $1-3/4=1/4$ under both schedules. The bias is $1/4$ in the first schedule and $5/4$ in the second.

Under the condition, the estimator is unbiased. Let $D_i=f(\mathbf{Z},\theta_i)$, and let $\widehat\tau_{HT}(d,d')$ be the difference between the averages over $i$ of $\mathbf{1}(D_i=d)\,Y_i/\pi_i(d)$ and of $\mathbf{1}(D_i=d')\,Y_i/\pi_i(d')$. Since $Y_i=y_i(d)$ on the event $D_i=d$ and $\mathbb{E}[\mathbf{1}(D_i=d)]=\pi_i(d)$, I have

Averaging over $i$ and subtracting the same identity for $d'$, I obtain $\mathbb{E}[\widehat\tau_{HT}(d,d')]=\tau(d,d')$. No condition on joint exposure probabilities, and no independence across units, is used.

For complete random assignment of $M$ of $N$ units and the mapping with four conditions, a unit with $k$ peers has

where a binomial coefficient is zero when its lower index is negative or exceeds its upper index. For $2\le M\le N-1$, all four exposures have positive probability exactly when $1\le k\le N-M-1$. If isolates are removed, the lower end of this range is secured, and the upper end is not. In the star, $N=4$, $M=2$ and $k=3$ for the centre. In the chapter’s design, $N=1000$ and $M=100$, and the range is $1\le k\le 899$.

What survives

The exposure mapping and the estimand are well defined, and the Horvitz–Thompson estimator is unbiased wherever the condition holds. This counterexample does not establish an error in the estimates, the simulation or the comparison of precision for the two example networks. The degrees displayed in Figure 16.2 are compatible with the range $1\le k\le 899$, and no failure of the condition is established for either network. Note 1 is not itself incorrect. What changes is the general statement on p. 212, and the sentence on p. 216 when it is read as a general rule. An application in which some unit has probability zero of one of the two exposures must change the design or the target.

The problem

The sentence gives the wrong outcome model for the only simulation of the chapter. The printed call passes make_corr_out to make_dilated_out with hop = 1, and the baseline outcome of unit $i$ is then generated as $B_i=(k_i+q_i+k_iq_i)\,|G_i|+E_i$. Here $k_i$ is the degree of unit $i$, and $q_i$ is the number of other units that a walk of two edges reaches from $i$. The variable $G_i$ is standard normal, and $E_i$ is normal with mean 1 and standard deviation $1/4$. The baseline outcome is not the absolute value of a standard normal variable, and it can be negative. Its mean given the network is $(k_i+q_i+k_iq_i)\sqrt{2/\pi}+1$, while the absolute value of a standard normal variable has mean $\sqrt{2/\pi}\approx0.80$. The first exceeds the second by at least 1 for every unit of a network without isolates. With the multipliers $(2,1.5,1.25)$, the three average effects are the mean of the $B_i$, one half of it and one quarter of it. Under the printed description they would be about 0.80, 0.40 and 0.20 in expectation. Table 16.1 reports estimates between 115.79 and 891.08. The sentence also does not say how strongly the outcomes depend on degree. The simulation cannot be rebuilt from the sentence alone, although the printed call and the package determine it.

Evidence an exact counterexample

I checked the function that the printed call invokes, using an exact counterexample and a short derivation. The function computes

(degree + degree2 + degree * degree2) * abs(rnorm(length(degree))) +

rnorm(length(degree), 1, 0.25)

where degree2 counts the other units $j$ for which the $(i,j)$ entry of the squared adjacency matrix is positive. For the counterexample, let three units form a triangle. Then $k_i=q_i=2$ and $B_i=8\,|G_i|+E_i$ for each unit. With the draws $G=(1/64,1,-1/2)$ and $E=(-3/4,1,5/4)$, I obtain $B=(-5/8,\,9,\,21/4)$. The functions of the package, run with these draws in place of the random ones, return these values and the corresponding multiples for the other exposures. On the event that $|G_1|<1/32$ and $-1<E_1<-1/2$, which has positive probability, I have $B_1<-1/4$. The absolute value of a standard normal variable is never negative. For the moments, recall that $\mathbb{E}|G_i|=\sqrt{2/\pi}$ and $\operatorname{Var}(|G_i|)=1-2/\pi$. Since $G_i$ and $E_i$ are independent given the network, I find $\mathbb{E}[B_i]=a_i\sqrt{2/\pi}+1$ and $\operatorname{Var}(B_i)=a_i^2(1-2/\pi)+1/16$ with $a_i=k_i+q_i+k_iq_i$.

Two separate recalculations illustrate the scale in the chapter’s designs. On the small-world network regenerated from the printed code and seeds under R 4.4.2 and igraph 2.2.2, which need not be the chapter’s realization (P044-N004), I obtain the three true effects 467.75, 233.88 and 116.94 from the generated potential outcomes. On the degree-corrected block model they are 623.80, 311.90 and 155.95 (P044-N009).

What survives

The multipliers $(2,1.5,1.25)$ agree with the code. The exposure mapping, the estimators and the design-based argument do not depend on the description of the baseline. The correction changes the sentence and the matching sentence of the package’s help page; it changes no code and no reported number. I did not recompute Table 16.1 or Figure 16.3 for this correction, and neither is shown to be wrong by it.

The problem

The text reads the argument rho = 0.7 of the generator as the proportion of edges within communities. In randnet::BlockModel.Gen the argument has a different meaning. Between-community connection probabilities are controlled by beta, relative to within-community connection probabilities. The printed call does not pass beta, and its default is $0$. The block probability matrix is then diagonal, and no edge joins two communities. Under the default simple = TRUE, the argument rho is the probability that a unit is given the degree weight $0.2$ in place of $1$, and alpha is not used. The second example network is therefore misdescribed, and it cannot be recovered from the sentence.

Evidence a check of the printed page

I checked the source code of the generator and a recalculation. In randnet 0.5, the arguments of BlockModel.Gen have the defaults beta=0, rho=0, simple=TRUE and alpha=5. The function sets the block matrix to diag(w) and fills its off-diagonal entries only if(beta > 0), and in the simple branch it sets node.degree[runif(n)<rho] <- 0.2. The manual describes beta as the “out-in ratio: the ratio of between-block edges over within-block edges”. The code of the generator is the same in versions 0.4, 0.5 and 0.7. The chapter does not say which version it used.

When the printed function is run with seed 42, the loop accepts the 193rd draw. The accepted network has all of its edges within communities and ten connected components. A share $0.723$ of its units have weight $0.2$, its mean degree is $9.858$, its degree variance is $75.87$ and its skewness is $1.075$. It reproduces the six effect estimates and the six variance estimates that are printed for the DCBM on p. 218, to their printed precision.

What survives

The degree variance of 75 and the skewness of 1.07 quoted in the same paragraph approximate those of the network that the printed code generates. The printed DCBM estimates follow from that network. The definitions, the estimand and the design-based arguments do not depend on the proportion of edges within communities, and the small-world network is unaffected. If a network with 70% of its edges within communities was intended, the settings of the generator would have to be changed and the DCBM example recomputed.

The problem

The statement is correct for the true exposure probabilities $\pi_i(d)=\Pr(D_i=d)$, when they are positive. The routine that computes “these exposure probabilities” does not return them. Let $c_i(d)$ denote the number of the $R$ supplied assignments under which unit $i$ has exposure $d$, and let $c_{ij}(d,d')$ denote the number under which unit $i$ has exposure $d$ and unit $j$ has exposure $d'$. The routine returns $(c_i(d)+1)/(R+1)$ for the probability that unit $i$ has exposure $d$, $c_{ij}(d,d)/(R+1)$ for the joint probability of two different units in the same exposure, and $c_{ij}(d,d')/R$ for two different exposures. It does so even when the $R$ assignments are all the assignments of the design, a case that p. 214 allows (“the set of all possible (or a large random sample of all possible) treatment assignments”). The four returned probabilities of a unit sum to $(R+4)/(R+1)$, which exceeds 1. Weighting by their inverses does not give an exactly unbiased estimator. The chapter does not mention the smoothing.

Evidence an exact counterexample

The evidence is an exact counterexample, which I also ran through the unchanged functions of the package. Let eight units form a cycle, and let the design be uniform on the eight rotations of the assignment $(1,1,0,1,0,0,0,0)$. Every assignment places two units in $d_{11}$, one in $d_{10}$, three in $d_{01}$ and two in $d_{00}$, and every unit has the true exposure probabilities $(1/4,1/8,3/8,1/4)$. Supplied with all eight assignments, the routine returns $(1/3,2/9,4/9,1/3)$, which sum to $4/3$. Let every unit have potential outcomes $(2,3/2,5/4,1)$ under $(d_{11},d_{10},d_{01},d_{00})$. Then $\tau(d_{11},d_{00})=1$, $\tau(d_{10},d_{00})=1/2$ and $\tau(d_{01},d_{00})=1/4$. With the returned probabilities, the Horvitz–Thompson contrasts equal $3/4$, $3/32$ and $39/128$ under every assignment, and these values are their expectations. With the true probabilities they equal $1$, $1/2$ and $1/4$.

The general computation is short. Let $q_i(d)$ be positive weights used in place of $\pi_i(d)$. Since $\mathbb{E}[\mathbf{1}\{D_i=d\}]=\pi_i(d)$, I have

For a design of $R$ equally likely assignments, all supplied, $\pi_i(d)=c_i(d)/R$ and $q_i(d)=(c_i(d)+1)/(R+1)$. The ratio $\pi_i(d)/q_i(d)$ then equals $c_i(d)(R+1)/\bigl(R(c_i(d)+1)\bigr)$, which differs from 1 whenever $0<c_i(d)<R$.

What survives

The exposure mapping and the estimands are unaffected, and the Horvitz–Thompson estimator remains unbiased when the exposure probabilities are the true ones and are positive. The chapter already says that the Hajek estimator is biased. The estimator computed from the package’s probabilities has no guarantee of exact unbiasedness, even when all assignments are enumerated. The counterexample does not measure the bias at the chapter’s $R=10{,}000$. The recalculation made for P044-N009 measures it on the degree-corrected block model generated by the printed code. There, conditional on the chapter’s simulated probabilities and with an independent assignment, the expectation of the Horvitz–Thompson estimator exceeds the true effect by 0.5%, 1.5% and 0.9% for $(d_{11},d_{00})$, $(d_{10},d_{00})$ and $(d_{01},d_{00})$. No entry of Table 16.1 and no comparison of Figure 16.3 is shown to change.

The problem

The package documents yT_ht as the Horvitz–Thompson estimator of the total of the potential outcomes under an exposure, $\widehat T_d=\sum_i\mathbf{1}\{D_i=d\}\,Y_i/\pi_i(d)$. When no unit receives $d$, every term of the sum is zero and $\widehat T_d=0$. Line 219 replaces this zero by NA. The contrast tau_ht is then missing, although the estimator that the chapter calls unbiased (pp. 212 and 216) equals a finite number. An average over the assignments for which the function returns a value conditions on a random event and can have a different expectation. A missing value is appropriate for the Hajek mean, which is $0/0$ for an empty exposure.

Evidence an exact counterexample

The evidence is an exact counterexample, which I also ran on the unchanged function. Let six units form a cycle, and let two of them be treated, uniformly over the 15 assignments. Every unit has exposure probabilities $(2/15,1/5,7/15,1/5)$ for $(d_{11},d_{10},d_{01},d_{00})$. Let every unit have potential outcomes $(2,3/2,5/4,1)$, so that $\tau(d_{11},d_{00})=1$. Under the assignment $(1,0,1,0,0,0)$ no unit receives $d_{11}$ and one unit receives $d_{00}$. The totals are $\widehat T_{d_{11}}=0$ and $\widehat T_{d_{00}}=5$, and the Horvitz–Thompson contrast is $(0-5)/6=-5/6$. The function returns NA. Over the 15 assignments, the contrast equals $10/3$ six times (the two treated units are adjacent), $-5/6$ six times and $0$ three times, and its mean is 1. The function returns a value only in the first six cases, and the mean of the returned values is $10/3$. With the overwrite removed, the function returns the exact totals and contrasts under all 15 assignments.

The general statement takes one line. Since $\mathbf{1}\{D_i=d\}\,Y_i=\mathbf{1}\{D_i=d\}\,y_i(d)$ and $\pi_i(d)>0$, I have $\mathbb{E}\bigl[\mathbf{1}\{D_i=d\}\,Y_i/\pi_i(d)\bigr]=y_i(d)$ for each unit, whether or not the exposure occurs under a given assignment. Summing over $i$ gives $\mathbb{E}[\widehat T_d]=\sum_i y_i(d)$, which implies that $\mathbb{E}\bigl[(\widehat T_d-\widehat T_{d'})/N\bigr]=\tau(d,d')$.

What survives

The Horvitz–Thompson estimator remains unbiased when the exposure probabilities are the true ones and are positive. The estimands are unaffected, and so are the point estimates under any assignment in which every exposure occurs. Table 16.1 prints finite estimates, and no entry of it and no result of Figure 16.3 is shown to change. I did not estimate how often an exposure is empty on the chapter’s networks of 1000 units. The correction concerns the point totals only and implies nothing about the variance estimators.