Corrections to: Randomization-based confidence sets for the local average treatment effect

The problem

Step 4 follows the simulated AR function that currently realizes the quantile, with index $m_i$, to its next intersection $e_i$, and then chooses a new index for the interval to the right of $e_i$. The printed step does this in both cases: when later intersections exist it evaluates $\mathcal I$ between $e_i$ and the next one, and when $e_i$ is the largest of all intersections it evaluates $\mathcal I$ at $e_i+1$. The code carries out the first case. In the second case it appends the current index s, which is $m_i$, for the final interval and leaves the loop. If the crossing changes the function that realizes the quantile, the code compares the observed AR function with the wrong simulated one on $(e_i,\infty)$. The returned set then differs from the intended inversion where the observed function is strictly above the smaller of the two critical functions and at or below the larger. The difference between the two sets consists either only of added points or only of deleted points, and it can be unbounded.

The defect is confined to the interval to the right of the largest intersection of simulated AR functions. Its effect on the reported numbers is not known. When the last intersection is a unique simple crossing of two singly occurring functions, the kept function is one rank away from the correct one, but this bounds neither how often the branch is reached nor how far an endpoint moves: in the example below the kept function is exactly one rank away at $M=1000$, $5000$ and $10^4$, and the right endpoint moves from $\infty$ to $12$.

Evidence an exact counterexample

The evidence is an exact counterexample, run through the unmodified code of both sources, and a derivation. I take eight units, four of them treated in every assignment, with the observed assignment $Z$ and four other assignments $A$, $B$, $H$ and $L$:

For an assignment $j$ let $t_j(\beta)$ be the difference of the arm means of $Y-\beta D$, and let $T(\beta)=2(\beta-1)^2+38$. The squared AR function of $j$ is $F_j=16t_j^2/(T-2t_j^2)$, with positive denominator, and

Since $F_j$ is increasing in $t_j^2$, every comparison reduces to one between squares of these linear functions. The simulated AR functions intersect at $-8,-3,0,2/3,1,2,5$ and $8$, and the largest intersection, at $\beta=8$, is a crossing of $F_A$ and $F_B$. Further,

Let $r$ be the rank of the quantile, and let the $M$ simulated assignments be $r-2$ copies of $L$, one copy each of $A$ and $B$, and $M-r$ copies of $H$. The squared critical value at rank $r$ is then $\min\{F_H,\max(F_A,F_B)\}$, and $\max(F_A,F_B)$ equals $F_A$ for $|\beta|\le8$ and $F_B$ otherwise. For (13), $r=\lceil0.95M\rceil$, and I obtain

Just left of $\beta=8$ the quantile is realized by $A$, and just right of it by $B$. The code keeps $A$ on $(8,\infty)$ and returns $[-4,-2]\cup[2,12]$. At $\beta=16$ the kept function equals $32/57$, the observed function $18/13$ and the squared critical value $128/45$.

The rows with a covariate use $32$ units, obtained by repeating each unit four times, with a covariate taking the values $-1,-1,1,1$ within each group of four and with $Y$ replaced by $Y-16D$. The endpoints in these rows agree with the stated values to $10^{-12}$. In these rows the returned set excludes the value $0$, which the intended set contains. With the same construction, the stated sets also hold at $M=5000$ and $M=10^4$, with ranks $4750$ and $9500$ for the research code and $4751$ and $9501$ for the package; the code was run at $M=1000$ only. A second, independent example with eight simulated assignments has set $(-\infty,-3]\cup[-1,\infty)$, and both routines return only $(-\infty,-3]$.

With correctly computed intersections and quantile indices, the repair is justified as follows. Under Assumption 6, or Assumption 7 with covariates, the squared AR functions are continuous, and two distinct ones differ by a rational function with at most four real zeros. To the right of the largest intersection no difference changes sign, so the order of the simulated values is constant there. Any index in $\mathcal I(e_i+1)$ realizes the quantile on the whole final interval, including its left endpoint by continuity.

The effect of the defect on the reported numbers is not established. The simulated assignments of the published runs are not in the replication files, and neither are the data of the application to Angrist et al. or of four of the six cities. An argument that the branch is reached with probability of order $2/M$ does not hold in general: for four units with $Y=(0,1,2,4)$, $D=(1,1,0,0)$ and two treated, under complete randomization, the code keeps a wrong index with probability at least $1-2\times10^{-106}$ at $M=1000$. Coverage over $5000$ data sets moves in steps of $0.0002$, and one changed decision can change a printed third decimal.

What survives

Step 4 as printed, Lemma 1, the definitions (13)–(17), Theorem 3 and the theory for the confidence sets are unaffected by this code defect. This defect does not change membership at or to the left of the largest intersection of simulated AR functions. A value $\beta$ to the right of it changes status only if the observed AR function is strictly above the smaller of the kept function and the correct one and at or below the larger; in particular a changed coverage decision at the true value $0$ in Section 6 requires both a negative largest intersection and this bracket at $0$. The p-value of the rilate package at $\beta=0$ is computed separately and is unaffected. The procedures that do not use randomization critical values are unaffected. The simulation and application scripts of the research code call the affected routine. I have not established whether any reported entry changes: this applies to the coverage rates and lengths of the randomization sets in Table 1 and Supplementary Tables 3–5 and 12, to their endpoints in Table 2 and Supplementary Tables 6 and 7, and to the relative lengths in Supplementary Table 10.

The problem

Let $\Delta_m^{*2}=N_m/Q_m$ be the squared AR function of a simulated assignment $m$, where $N_m$ and $Q_m$ are the quadratics of (16). Under Assumption 6, two such functions meet exactly at the real zeros of $P_{ms}=N_mQ_s-N_sQ_m=\sum_{j=0}^4v_j\beta^j$, and they are the same function exactly when all five coefficients $v_j$ are zero. This is Remark 7. The code replaces “$v_j=0$” by “$|v_j|<10^{-8}$”, both to decide whether two functions are the same and to choose the degree of the polynomial that it solves. The coefficients $v_j$ have no fixed scale. Multiplying the coefficient vectors of both functions by $t>0$ multiplies every $v_j$ by $t^2$ and leaves the functions unchanged, and in the binary examples below the $v_j$ decrease as $n$ grows. The code discards nonzero terms that are below the threshold, which can change the set of real zeros. When all five coefficients of a pair are below the threshold, it reports no intersection for that pair. When this happens for every pair, the grid of intersections is empty and the code for Algorithm 2 returns the whole real line (P002-E015).

The threshold is reached on unrescaled data in the settings of the paper. In recalculations with the same simulated assignments, nonzero terms were discarded in one saved data set for each of the fourteen combinations of design and allocation of Section 6 ($M=1000$) and in the Bridgeport and Columbus data, with and without the covariate. In one recalculation for Columbus without the covariate at $M=1000$, $51\,705$ pairs of distinct functions are classified as identical and the code reports $621\,548$ real roots, compared with $977\,137$ for the full polynomials, counted separately for each pair. Every endpoint returned in these runs still agrees with the exact set (14), or (F.4) with the covariate, to three decimals. The largest error is approximately $0.0000336$, at the lower endpoint for Columbus with the covariate at $M=10^4$, where the code returns approximately $-0.041607$ and the exact value is approximately $-0.041574$. The failure is of a different kind when outcome and treatment are binary and the sample is large: the code returns the whole real line where (14) is a short interval. This occurs from $n=1176$ in a constructed family and at $n=19\,392$ for the composition of the Columbus data with independent simulated assignments.

Evidence an exact counterexample

I use exact counterexamples, run through the unmodified R and C++ code, and recalculations using exact zeros and roots of the calculated AR functions.

First, consider the coefficient vectors $(1,0,0,1,0,1)$ and $(1,-2,1,1,0,1)$, which represent $\beta^2/(\beta^2+1)$ and $(\beta-1)^2/(\beta^2+1)$. Their only intersection is $\beta=1/2$, and $P=(2\beta-1)(\beta^2+1)$. After both vectors are multiplied by $10^{-5}$ each function is unchanged, the coefficients of $P$ are $(-10^{-10},2\times10^{-10},-10^{-10},2\times10^{-10},0)$, and the root finders of the rilate package and of the research code (R and C++) return no intersection. The interval routine of the package then returns the whole real line where the first function is at most the second exactly on $(-\infty,1/2]$.

Second, I take eight units with $Y=(0,1,2,3,4,5,6,8)$ and $D=Z=(0,0,0,0,1,1,1,1)$, and $M=40$ simulated assignments consisting of $20$ copies each of $(1,0,0,0,0,1,1,1)$ and $(0,1,0,0,0,1,1,1)$. All denominators are positive, the two simulated AR functions intersect only at $\beta=5$, and the set (14) at level $0.95$ is $[3,14/3]$. The package returns this interval to within $10^{-12}$. When the outcomes are multiplied by $10^{-7}$, the set (14) is $[3\times10^{-7},\,(14/3)\times10^{-7}]$, the code loses the intersection at $5\times10^{-7}$, and AR_algo2 returns the whole real line.

Third, no rescaling is needed. Let $n$ be a multiple of $8$ such that $n>8$, and consider $n/4$ units of each type $(Y,D)\in\{(0,0),(1,0),(0,1),(1,1)\}$ with $n/2$ treated. Let the observed treatment group contain $0$, $n/8$, $n/8$ and $n/4$ units of these types, the first simulated one $n/8-1$, $n/8+1$, $n/8$ and $n/8$, and the second $n/8-1$, $n/8$, $n/8+1$ and $n/8$, with $20$ copies of each simulated assignment. In the representation used by the code, the numerators are $N_o=(\beta-1)^2/4$, $N_1=16/n^2$ and $N_2=16\beta^2/n^2$, the denominators are $Q_j=(\beta^2+1-N_j)/n>0$, and

The two simulated AR functions intersect at $\pm1$. The code discards $P_{12}$ when $16/n^3<10^{-8}$, that is, for $n\ge1170$. At $n=1168$ it finds both intersections and returns a bounded set. At $n=1176$ the set (14) is $[146/147,\,147/146]$, and the package returns $(-\infty,\infty)$.

Fourth, independent simulated assignments lead to the same state. For the pooled composition of the Columbus data ($2075$, $190$, $130$ and $29$ units of the four types) replicated $8$ times, $n=19\,392$, with balanced assignment and $M=1000$ or $M=10^4$ independent simulated assignments, the largest $|v_j|$ over all pairs is below $10^{-8}$, while (14) equals $[-14/323,\,50/629]$ and $[-14/323,\,50/621]$ respectively. At $n=9696$ some pair is still above the threshold. There is no threshold in $n$ alone: for equal numbers of the four types the first tested size is $51\,200$ at $M=1000$ and $102\,400$ at $M=10^4$.

The first and second rows were run through the unmodified code. The fourth used the same root and inversion decisions with optimized storage and grid searches. In the third row the returned set follows from the classification of all pairs, which was computed exactly. In the fourth row the set is the regression-adjusted analogue (F.4) of (14).

What survives

Remark 7, Lemma 1, Algorithms 1 to 4 as printed and Theorem 3 use the actual intersections and are not contradicted by this code defect. The definitions of the confidence sets and the finite-sample and asymptotic theory are unaffected. No reported value is shown to change: every endpoint checked in the Section 6 and city recalculations agrees with exact inversion at the three decimals printed in Supplementary Tables 6 and 7. The checks do not cover the application to the data of Angrist et al., the cities Raleigh, Minneapolis, Detroit and St. Paul, or the Bridgeport subsamples, whose inputs are not in the replication files, and they do not use the simulated assignments of the published runs. The research code also contains a C++ routine in exact integer arithmetic. It tests for exact zero and is not affected. A user who applies the replication code to a large experiment with binary outcome and treatment can obtain the whole real line where (14) is bounded.

The problem

Lemma 20 concerns a sequence of models with covariates $\theta_n\in\Xi_n$, and its claims hold almost surely with respect to the randomness in $Z$. The symbols $\mathbb P^*$ and $\operatorname{var}^*$ refer to reassignment. The feasible statistic $\widehat\tau^*_{\mathrm{adj}}$ estimates the regression coefficients anew in each arm for every reassignment (equation (9) of the article). The infeasible statistic $\tilde\tau^*_{\mathrm{adj}}$ of (D.4) holds them fixed at the coefficients of the regression on the full observed population, and its variance is stated in (D.6). Items 1–3 of the lemma, (E.58)–(E.60), are normalized by $\sqrt{\operatorname{var}^*(\tilde\tau^*_{\mathrm{adj}})}$, and the proof of Theorem 2 needs item 4 on the same scale. As printed, item 4 asserts instead that $\widehat\sigma^2_{\mathrm{adj},*}$ is consistent for the variance of the feasible statistic. This is a different claim, and it is false under the stated assumptions. Assumption 5 controls the feasible statistic in probability, but it does not control what rare reassignments with ill-conditioned arm regressions contribute to its second moment.

Evidence an exact counterexample

The evidence is an exact counterexample to the printed statement and a derivation of the corrected one. Recall that $\widehat\sigma^2_{\mathrm{adj},*}=RSS_1/n_1^2+RSS_0/n_0^2$, where $RSS_a$ is the residual sum of squares of the regression in arm $a$ after reassignment.

Counterexample. Let $k=1$, $\delta=0$, $\tilde\varepsilon=1/4$, $A=16$, $B=5$ and $r=1/3$. For even $n=2m\ge8$, let $t_j=(2j-(m-1))/(m-1)$ for $j=0,\ldots,m-1$, $\eta_n=2^{-n^2}$ and $v_m=m^{-1}\sum_j t_j^2$. The population has two units for each $j$, indexed by $s\in\{1,-1\}$, with

and $n_1=n_0=m$ under complete randomization. Every unit is a complier, the local average treatment effect is $\beta=0$, Assumptions 1, 2, 4 and 5 hold with the constants above, and every arm regression has full rank. On the reassignment that treats the whole block $s=1$, both fitted slopes equal $1/\eta_n$ and $\widehat\tau^*_{\mathrm{adj}}=-2/\eta_n$; on its complement, $\widehat\tau^*_{\mathrm{adj}}=2/\eta_n$. Exchanging the two arms changes the sign of $\widehat\tau^*_{\mathrm{adj}}$, so its mean over reassignments is zero, and

Since the fit with zero intercept and zero slope is available in each arm and $|y_i|\le11/10$, every reassignment has $\widehat\sigma^2_{\mathrm{adj},*}\le121/(25n)<5/n$. Combining the two bounds, I find that the ratio in (E.61) is at most $(5/(8n))\,2^{n-2n^2}<1/2$ on every reassignment. This implies that

for every even $n\ge8$ and every original assignment, while (E.61) asserts that this probability tends to $0$.

Derivation of the corrected statement, taking Lemma 19 as given. I write $e_i$ for the residuals of the regression on the full observed population, and I set $H=n^{-1}\sum_{i=1}^ne_i^2$ and $V=\operatorname{var}^*(\tilde\tau^*_{\mathrm{adj}})$. By (D.6), I have $V=(1/n_1+1/n_0)\,nH/(n-1)$. Let $Q_a$ and $h_a$ be the averages of $\tilde x_i\tilde x_i'$ and of $\tilde x_ie_i$ over the units reassigned to arm $a$. When $Q_a$ is invertible, I obtain from the normal equations of least squares that

Let $Q=n^{-1}\sum_{i=1}^n\tilde x_i\tilde x_i'$. Its smallest eigenvalue is at least $\min\{1,\tilde\varepsilon\}$ by Assumption 5(iii) and (v). Using Assumption 5(iv), the variance formula for sampling without replacement and the identity $\sum_{i=1}^n\tilde x_ie_i=0$, I find that $Q_a-Q\to0$ and $\|h_a\|_2^2/H\to0$ in $\mathbb P^*$-probability. By Lemma 19 and Chebyshev’s inequality, I have $n_a^{-1}\sum_{i\in\mathcal T^*(a)}e_i^2/H\to1$ in probability. Then $RSS_a/(n_aH)\to1$ in probability for each $a\in\{0,1\}$, and I conclude that

in probability, which is the corrected (E.61).

What survives

With the tilde in the denominator, (E.60) and (E.61) are on the same scale, and Slutsky’s theorem, applied conditionally, gives the normal critical values that this step of the proof of Theorem 2 uses. The correction adds no assumption and changes neither the test statistic, nor the confidence set, nor the coverage statement. Equations (9)–(11) of the article, (F.1)–(F.4) and Algorithms 3 and 4 use residual sums of squares and need no change. Theorem 1 and the unadjusted procedure do not use Lemma 20. The counterexample does not contradict the coverage conclusion of Theorem 2, since all its treatment effects are zero.

The problem

The quantity that is centered is $w_i=y_i(d_i(1))-\beta_nd_i(1)$, whose population mean is $\mu_{y(d(1))}-\beta_n\mu_{d(1)}$. By (B.3), I have $v_i(1)=w_i-\mu_{y(d(1))}+\beta_n\mu_{d(1)}$. The term subtracted in the bracket is then $\mu_{y(d(1))}-\beta_n\mu_{d(1)}$. With the printed sign the bracket is the error of the treated sample mean minus $2\beta_n\mu_{d(1)}$. The identity is then false, and the expectation of the numerator in (C.17) exceeds the printed value by $4\beta_n^2\mu_{d(1)}^2$, which need not vanish. Both displays also print $y_i(d(1))$ for $y_i(d_i(1))$.

Evidence a derivation

The evidence is a derivation and an exact example. Set $\mu_w=\mu_{y(d(1))}-\beta_n\mu_{d(1)}$ and $\bar w_1=n_1^{-1}\sum_iZ_iw_i$. For a treated unit, $\widehat v_i(1)=w_i-\bar w_1$ by equation (3) of the article and $v_i(1)=w_i-\mu_w$ by (B.3). Expanding the square, I obtain

Under complete randomization, I have $E[\bar w_1]=\mu_w$ and

which is the right side of (C.17) after division by $S_1^2(\beta_n)$. For the example, set $n=4$, $n_1=2$ and $\beta_n=2$, let all units be compliers with $y(0)=(-1,-1,1,1)$ and $y(1)=y(0)+2$, and treat the first two units. The left side of the identity is $0$. The printed right side is $1-(-5)^2=-24$, and the corrected right side is $1-(-1)^2=0$.

What survives

With the corrected sign the identity holds and (C.17) holds with its printed right side. Lemma 7 and the rest of its proof are unchanged, and no assumption is added.

The problem

The ratio divides an estimated standard deviation by a variance. In the balanced, bounded-outcome example below, the denominator is of order $1/n$ and $\widehat\sigma^*$ is of order $n^{-1/2}$, so the printed ratio is of order $n^{1/2}$ and does not approach $1$. The sentence that introduces the lemma, the parallel statement in Lemma 7 and the use of the lemma to studentize in Lemma 13 all concern the ratio of the estimated variance to the variance.

Evidence a derivation

The evidence is an exact example in which the printed statement fails, and a derivation of the corrected one. Let $n$ be even and $n_1=n_0=n/2$. Suppose that every unit is a complier with $y_i(1)=y_i(0)$, that half of these outcomes equal $-1$ and half equal $1$, and that $\beta_n=0$. Let $q$ be the mean outcome in a reassigned treatment group, so that the control mean is $-q$. Computing the three quantities directly, I find

The printed ratio equals $2(n-1)n^{-1/2}(1-q^2)^{1/2}$, which diverges in $\mathbb P^*$-probability for every original assignment. The ratio with the squared numerator is $R=(n-1)(1-q^2)/n$. Since $0\le R\le1$ and $E^*(1-R)=2/n$, Markov’s inequality implies that $\mathbb P^*(|R-1|\ge\kappa)\le2/(n\kappa)$. With $\kappa=1/2$ and $n=8,32,128,512$, the exact probability for the printed ratio equals $1$ up to rounding error in each case, and for the corrected ratio it is $0.029$, $4.9\times10^{-5}$, $6.5\times10^{-17}$ and $1.6\times10^{-63}$.

In general, let $a_i=Y_i-\beta_nD_i$ minus its mean, $S^2=(n-1)^{-1}\sum_ia_i^2$ and $M^2=\max_ia_i^2$. Lemma 11 implies that $M^2/(nS^2)\to0$ for almost every sequence of original assignments. In a reassigned arm of size $m$, the mean of the $a_i^2$ has expectation $(n-1)S^2/n$ and variance at most a constant times $M^2S^2/n$, and the squared arm mean has expectation $O(S^2/n)$. By Chebyshev’s and Markov’s inequalities, the centred second moment in each arm divided by $S^2$ tends to $1$ in $\mathbb P^*$-probability. Substituting into the definition of $\widehat\sigma^{*2}$ on p. 6 of the article and comparing with (C.25), I obtain the corrected statement.