Corrections to: Dyadic Clustering in International Relations

The problem

The passage on p. 187 states a direction: when non-incident dyads are dependent, DCRSEs underestimate uncertainty. From this direction the article infers that its reanalysis is a lower bound on the consequences of inadequate clustering, and the Conclusion repeats the inference. The direction does not hold in general. The DCR variance estimator treats as zero the covariances between observations whose dyads share no member. In the variance of an estimated coefficient, each omitted covariance is weighted by the product of the weights of its two observations, and the sum of these weighted terms can be positive, zero or negative. When the sum is negative, the true variance of the estimated coefficient is smaller than the population dyadic sandwich value, and DCRSEs can overstate uncertainty. The article states no condition on the sign of the sum. The condition in Section 2.1, that $a$ and $c$ are of the same sign, belongs to the toy example. A positive covariance between the errors of non-incident dyads is also not enough, since the weights of two observations in a slope can have opposite signs. Without a sign condition, the comparison of DCRSEs with the original standard errors measures the effect of dyadic clustering and bounds the effect of further dependence in neither direction. The expectation of a larger variance on p. 191 requires a positive omitted weighted covariance sum; the aside there that DCRSEs still improve over the naive approach is not guaranteed.

Evidence an exact counterexample

The evidence is a derivation and six exact examples. Let $X$ be a fixed $n\times p$ matrix of full column rank, and let $Y=X\beta+u$ with $\mathbb{E}[u]=0$ and $\mathbb{E}[uu^T]=\Omega$. Fix $c\in\mathbb{R}^p$ and set $h=X(X^TX)^{-1}c$. Since $\widehat\beta-\beta=(X^TX)^{-1}X^Tu$, I have $c^T\widehat\beta-c^T\beta=h^Tu$, and the true variance is $V=\operatorname{Var}(c^T\widehat\beta)=h^T\Omega h$. Let $d_r$ denote the dyad of observation $r$, and let $\Omega_D$ impose the restriction of Table 4, with entries $(\Omega_D)_{rs}=\Omega_{rs}$ if $d_r\cap d_s\neq\varnothing$ and $(\Omega_D)_{rs}=0$ otherwise. The population dyadic sandwich value is $V_D=h^T\Omega_Dh$. Subtracting, I obtain

where the sum is over ordered pairs of observations. I denote this sum by $\Delta$. This identity implies that $V_D\le V$ if and only if $\Delta\ge0$. Note that positive definiteness of $\Omega$ does not restrict the sign of $\Delta$.

Examples 1 to 3 have two observations, on two dyads with no common member. In Examples 1 and 3 the coefficient is the intercept, with $h=(1/2,1/2)^T$. In Example 2 it is the slope on the single regressor $x=(-1,1)^T$, with $h=(-1/2,1/2)^T$. Example 4 has one observation on each of the six dyads among four members, and Examples 5 and 6 have one observation on each of the dyads $(A,B)$, $(A,C)$, $(D,E)$ and $(D,F)$. In Examples 4 to 6 the coefficient is the intercept. Every error has variance 1. In Examples 1 to 5 the matrix $\Omega$ is positive definite, and in Example 6 it is the covariance matrix of an explicit distribution.

Examples 1, 2, 4 and 5 have $V<V_D$. Example 2 has a positive covariance between the two errors, and Examples 4 and 5 have positive dependence between incident dyads. Examples 3 and 6 show that the other direction also occurs.

The fitted estimator also overstates variance in the limit. In Example 4 the errors can be built from independent signs and are bounded. For $m$ independent copies of this block on disjoint members, I have $n=6m$ and $V=1/n$, and the DCR variance estimator, computed from residuals without a small-sample correction, has expectation $3/(2n)-5/n^2$. The ratio of this expectation to $V$ is $3/2-5/n$, which equals $2/3$ for $m=1$, $13/12$ for $m=2$ and $17/12$ for $m=10$. By the law of large numbers, the ratio of the estimator to $V$ converges to $3/2$. For $m$ independent copies of Example 5, I have $V=3/(16m)$, and the estimator has expectation $(10m-3)/(32m^2)$, which exceeds $V$ for every $m\ge1$. The ratio of the estimator to $V$ converges in probability to $5/3$. In the limit, DCRSEs exceed the standard deviation of the estimate by the factors $\sqrt{3/2}$ and $\sqrt{5/3}$. In Example 4 the variance with every covariance between distinct observations set to zero, as under the naive approach, equals $V$. In Example 5 it is $1/4$, which lies between $V=3/16$ and $V_D=5/16$.

What survives

The toy example of Section 2.1 with its condition on the signs of $a$ and $c$, the sandwich identity, the restriction of Table 4 and the decomposition on p. 191 are unaffected. The ratios of DCRSEs to the original standard errors, the comparisons of $p$-values and the software do not use the direction of the omitted dependence, and this finding changes none of them. I have not independently recalculated those reported values. The sentence in the third paragraph of the Conclusion, that the problem “may be even more severe”, is correct as a statement of possibility (Examples 3 and 6). What changes is the interpretation. The reanalysis is not shown to be a lower bound. A population variance lower bound requires $\Delta\ge0$ for each coefficient, which has not been examined for the 22 studies. The result of Aronow, Crawford and Zubizarreta (2018) cited on p. 191 is not in question.

The problem

Beside the DCRSE, the command exports a comparison standard error from a second fit of the model, and the ratio of the two is the SER of the reanalysis. Three branches choose the second fit. The first runs when crse is given and clusters on the named variable. The second is meant to run when rse is given and to use the robust variance of the estimation command. Its condition repeats that of the first. The branch is reached only when crse is empty and entered only when crse is not empty, which cannot happen. No other line tests rse. A call with rse and without crse then falls to the third branch, which fits the model without adding a variance option. For the reanalysis calls, the comparison standard error, its $p$-value and the SER are then computed from the default variance of the estimation command. The reanalysis has nine such calls. Six are stcrreg models, whose default variance is already robust, and for them the defect need not change anything. Three are models whose default variance is not robust and whose original fit asked for a robust one: the logit model of Kinne (2018), Figure 8, Model 1, and the stcox models of Weisiger (2016), Table 3, Model 1 and Table 4, Model 1. Their seven KEVs, two of the six of Kinne and five of the 28 of Weisiger, are compared with a standard error that does not come from the original variance estimation strategy.

Evidence an exact counterexample

I checked the code and computed two exact examples. In dcr.ado the word rse occurs only on line 16, where the option is declared, and on line 36, where it is removed from the options passed to the estimation command. Lines 310 and 321 carry the same test.

In the example, $y=(0,0,1,3)$ is regressed on an intercept and $x=(-1,0,0,1)$, the four dyads share no member, and the command is called with rse and dofundo(reglike). Then $\widehat\beta=(1,\,3/2)$, the residuals are $(1/2,-1,0,1/2)$, and the three variance matrices are diagonal.

A second exact example is a logit model on the fifteen dyads among six members. Its slope has default variance $47/84$, robust variance $2045/5488$ and DCR variance $195/392$. The squared SER is $585/658<1$ as exported and $546/409>1$ with the robust comparison. In this example the defect moves the ratio from one side of one to the other.

I computed both examples from the formulas of the command and of the estimation commands; I did not run Stata. The seven KEVs are listed below with the values that the command exported. They are not corrected values. I did not fit the three models again.

What survives

The definition of the SER on p. 193 is correct. The DCR variance is computed before this block and is unchanged by this repair, as are the coefficients, the DCRSEs and the $p$-values based on them. The R package, the calls with crse and the calls with neither option are not affected. The comparison standard errors, their $p$-values and the SERs of the seven KEVs need to be recomputed; I have not measured their numerical changes. They enter the rows Kinne (2018) and Weisiger (2016) of Supplementary Table A.2, the 2016, 2018, Security and All rows of Table 5, Figures 1–3 and, through the small-sample scripts, Supplementary Table A.5. Under the inverse study frequency weighting the seven KEVs carry $(2/6+5/28)/22$, about 2%, of the total weight; without an upper bound on their corrected ratios, this gives no upper bound on the corrected mean SER. No rounded figure of the article, the 1.74, the 0.22 and the 67% among them, is shown to change. The DCRSEs of the Weisiger models are the subject of P032-N004.

The problem

The R package is the released implementation of the estimator on p. 191,

Both routines convert the row names of the fitted model frame to integers and use them as row positions in the data. Row names label rows; subsetting or reordering a data frame can preserve names that no longer equal row positions. This can make dcr() build its list of dyad members from the wrong rows, so that $N$ and the sum over members can be wrong; dcr_sandwich() can also attach cluster labels to the wrong observations. Separately, create_dyadid() labels a dyad by sorting its two member labels and joining them with an underscore. The pairs (A_B, C) and (A, B_C) both receive the label A_B_C. Two dyads that share no member are then treated as one repeated dyad in $\widehat V_D$. In each case the routine returns a covariance matrix and gives no warning. The inputs are ordinary. The replication script itself passes a subset of the data for models 1 and 2 of Pratt (2018), where the result happens to be unaffected, and nothing in the documentation of the package excludes underscores in member labels.

Evidence an exact counterexample

The evidence is three exact examples, computed in rational arithmetic and reproduced by running the unmodified package. Each has four observations, responses $y=(1,1,-1,-1)$ and a least-squares fit with an intercept only. The fitted intercept is $0$, the residuals equal $y$, and every variance below is a sum of products of residuals divided by $16$.

Since the dyads are distinct in the first two examples, $\widehat V_D=\widehat V_0=4/16=1/4$. The clusters of members A and D each give $(2^2+1+1)/16=3/8$, and the other four members give $1/4$ each. By the decomposition, I have

With row names 3 to 6, dcr() reads rows 3 and 4 of the four-row frame and two rows that do not exist. It finds the members D, E and F only, and returns

With row names 1, 3, 2, 4, dcr_sandwich() assigns member A to observations 1 and 3 and member D to observations 2 and 4. The residuals in each of these clusters sum to zero, so the clusters of A and D give $2/16=1/8$ each, and the routine returns

In the third example the eight members are distinct. Then no two observations share a member, and I find $\widehat V_r=\widehat V_0=1/4$. Since the first two dyads are both labelled A_B_C, the routines compute $\widehat V_D=\{(1+1)^2+1+1\}/16=3/8$ in place of $1/4$ and return

Row names $101$ to $104$, or names that are not numbers, make both routines stop with an error. With the two changes given below, both routines return $1/2$, $1/2$ and $1/4$ on these inputs.

What survives

This finding does not concern the estimator, its decomposition on p. 191, the Stata command or the simulations. No number reported in the paper is shown to change. For fits 1 to 5 of Pratt (2018), the covariance matrix without the optional small-sample and positive-semidefinite adjustments agrees with an independent calculation to within $10^{-8}$, all $16$ members (fits 1 and 2) and all $25$ members (fits 3 to 5) are counted, and no distinct dyads share a label. The three Arel-Bundock (2017) fits and the sixth Pratt fit could not be checked.

The problem

Twenty-four models in three studies used neither robust nor clustered standard errors. The Stata script runs the DCR command on each of them twice. The first run uses the comparison standard error of the original study, which is the nonrobust one, and the second uses a standard error clustered on repeated dyads. The script writes the two sets of results to two folders. The block of the merge script that builds Table A.3 reads the first folder for all 24 models. The column Avg. SER then repeats ratios of the primary analysis. The entries for Copelovitch and Putnam and for Fang, Johnson, and Leeds are rounded versions of those in their rows in Supplementary Table A.2, 0.898 and 1.691. The table does not report the comparison that its note and footnote 9 describe. The difference changes how the table is read. For Kinne (2018) the average ratio is below one and not above it. For Fang, Johnson, and Leeds (2014) the DCRSE is about 22% larger than the repeated-dyad clustered standard error, in place of the roughly 69% implied by the printed ratio.

Evidence a recalculation from the source

The evidence is a recalculation from the saved result files of the replication package. No model was fitted again. I applied the rules of the merge script to the 50 KEVs that the block selects, once in the 24 files that it reads and once in the 24 files that the Stata script writes for the clustered comparison. The rules round each ratio to three decimals, average the rounded ratios, round the mean to three decimals and display it to two decimals. The coefficients and DCRSEs agree up to numerical rounding; the comparison standard errors and their p-values change.

The files that the block reads reproduce the printed column exactly. For each of the 50 KEVs the significance category at the 5% level is the same in both sets of files, so the four proportions of each row are unchanged. As an instance, the one KEV of Fang, Johnson, and Leeds has DCRSE $0.2308$, nonrobust standard error $0.1365$ and repeated-dyad clustered standard error $0.1889$, to four decimals; the two ratios round to $1.691$ and $1.222$.

What survives

The numbers of KEVs and the four significance proportions of Table A.3 are unchanged, and the table note and footnote 9 stand once the column is corrected. Table 5, with the average SER of 1.74 and the significance-loss proportion of 0.22, Figures 1–3 and Supplementary Tables A.1, A.2 and A.4 are built from other objects of the script and do not change under this correction. The merge script for the small-sample correction makes the same choice of folder in a block whose table is not printed.

The problem

The decomposition is an identity when all of its pieces are sums over the contributions of the same observations, where an observation is one row of the data. In $\widehat V_{c,i}$ all observations whose dyad contains member $i$ form one cluster and each other observation is its own cluster, and in $\widehat V_0$ every observation is its own cluster. The Stata command builds the member-specific pieces in this way. It obtains $\widehat V_0$ by fitting the model again with vce(robust), and it subtracts $(N-2)$ times the result. For stcox and stcrreg, on data declared with stset and an id() variable, the variance that Stata returns for vce(robust) is clustered on the id. The contributions of all observations with one id are summed before they are squared. The reanalysis declares such an id for the three stcox models of Efrat and Newman (2016) and for the ten models of Weisiger (2016), two stcox and eight stcrreg. These 13 models hold 32 KEVs. For them the command returns, before the eigenvalue step of footnote 3 and when every piece carries the same scaling,

where $\widehat V_{\mathrm{id}}$ is the variance matrix clustered on the id. The Weisiger calls set dofundo(asylike), which multiplies the piece by $(n-1)/n$, with $n$ the number of observations, while the factor that Stata documents for a clustered variance is $g/(g-1)$, with $g$ the number of ids. The last term is then $(N-2)\bigl(\widehat V_0-\alpha\widehat V_{\mathrm{id}}\bigr)$ with $\alpha=((n-1)/n)(g/(g-1))$. The difference from $\widehat V_r$ has no fixed sign, and any difference between the two matrices is multiplied by $N-2$. I did not fit the 13 models again, and I do not know the size or the sign of the difference in the two studies.

Evidence an exact counterexample

The evidence has three parts: a check of the code, the documentation of the two Stata commands, and an exact example.

In dcr.ado, lines 108–119 give all observations that contain a member the cluster value $-99$ and give every other observation its own row number. Lines 126 and 202 fit the repeated-dyad piece and the member-specific pieces with explicit cluster variables. Line 164 fits the piece called vhet with vce(robust), and line 190 subtracts $(N-2)$ times it.

The Stata Survival Analysis Reference Manual, Release 15 (StataCorp 2017), in the entry [ST] stcox, technical note under “Robust estimate of variance” (p. 17), says of vce(robust) on data with an id that “stcox knew to specify vce(cluster id)”, and it lets an explicit cluster variable override the id. The entry [ST] stcrreg, last paragraph of “Methods and formulas” (p. 24), states that the clustered version of the variance is used when a cluster variable is given “or an ID variable has been stset.” The description of noadjust in [ST] stcox (p. 4) gives the default multiplier as $n/(n-1)$ or $g/(g-1)$. I read these entries, and the current ones, which agree with them. I did not read Release 17, which the replication files name, and Stata was not run.

The example has six subjects, one for each dyad among four members, with ten observations, one covariate $x$ and no weights. Status 1 is the event of interest, and status 2 is a competing event, treated as censoring.

With Breslow ties the Cox partial likelihood is maximized at $\widehat\beta=0$, the information is $12$, and the scores of the ten observations are those of the last column. Multiplying the sums of squared scores by $1/144$, I obtain the pieces below.

The printed decomposition returns $5/144$, which is also the sum of the products of all pairs of observations that share a member. With the variance matrix clustered on the id, it returns a variance $23/15$ times as large, and with the factor of the Weisiger calls a variance $103/75$ times as large. All three values are positive, so the eigenvalue step changes none of them.

The identity in the corrected statement holds by counting. The product $z_rz_s^\top$ of two observations enters $\sum_i\widehat V_{c,i}$ once for each member $i$ for which the two observations lie in one cluster. If $r=s$, this happens $N$ times, and subtracting $\widehat V_D$ and $(N-2)\widehat V_0$ leaves $N-1-(N-2)=1$. If $r\ne s$ and the two observations have the same dyad, it happens twice, and subtracting $\widehat V_D$ leaves $1$. If the dyads differ and share one member, it happens once and nothing is subtracted. If the dyads share no member, the product does not appear.

What survives

The printed decomposition is correct as an identity, and the 20 studies that are not fitted with these survival commands are not touched by this point. What changes is the standing of the two rows of Supplementary Table A.2. They are not shown to be DCRSEs as p. 191 defines them. The two rows enter the 2016, IOs/International Law, Security and All rows of Table 5, the upper end of the range “0.90 to 4.16” on p. 196, Figures 1–3 and, through the small-sample scripts, Supplementary Table A.5. No corrected value has been computed, and none of these numbers is shown to change. The following averages show how much rests on the two rows and are not corrections. The mean of the Avg. SER column of Table A.2 over the other 20 studies is 1.546, against 1.743 over all 22. The mean for IOs/International Law without Efrat and Newman is 1.701, and the mean for Security without Weisiger is 1.338.

The problem

Supplementary Table A.2 lists four studies in which the share of KEVs that go from significant to insignificant is one half or more: Colgan (2014) with $8/9$, Efrat and Newman (2016) with $2/4$, Fang, Johnson, and Leeds (2014) with $1/1$ and Weisiger (2016) with $17/28$. Three of them are strictly above one half.

Evidence a historical Lean attribution (Generated status: clause support not certified in this release)

The evidence is a recount from the column Sig.$\to$Insig. of Table A.2, where every other study is at $0.375$ or below. The four counts agree with the saved model outputs in the replication files. The spelling is a check of the page.

What survives

The numerical entries of Table A.2 and every aggregate are unchanged. The neighbouring counts in the paragraph, $18$ of $22$ studies with an average SER above one and $7$ of $22$ that lose no significance, agree with the table.

The problem

A symmetric matrix whose diagonal entries are all positive can have a negative eigenvalue. Both packages run the eigenvalue step of footnote 3 only if some diagonal entry of the DCR variance matrix is negative. If the matrix is indefinite and its diagonal is positive, the step is skipped. The matrix that the option returns then assigns a negative variance to some linear combination of the coefficients. The help texts of both packages state without this restriction that the option forces the matrix to be positive semidefinite. Footnote 3 speaks of negative parameter variances, and the code agrees with it. The defect is in the guarantee of the software option.

Evidence an exact counterexample

The evidence is an exact counterexample, computed in rational arithmetic and reproduced by running both R functions.

I regress $y$ on an intercept and $x$ by least squares on the ten dyads among five members, taken in the order $(1,2)$, $(1,3)$, $(2,3)$, $(1,4)$, $(2,4)$, $(3,4)$, $(1,5)$, $(2,5)$, $(3,5)$, $(4,5)$, with

Then

With $N=5$, the decomposition of p. 191 gives, in the order (intercept, slope),

Since both diagonal entries are positive, neither package applies the eigenvalue step. The determinant of $\widehat V_r$ is $-70529/308915776$, which implies that one eigenvalue is negative. For $c=(1,6)$, I have

and this is the estimated variance of the intercept plus six times the slope. Both functions of the R package, installed without modification, return this matrix when posdef = TRUE, with eigenvalues $0.12975$ and $-0.00176$. For the Stata command, the same matrix follows from the source for the call dcr regress y x, dm1() dm2() dofundo(reglike) psd, up to the order of the coefficients. Stata was not run.

The corrected option returns a positive semidefinite matrix. In the notation of the corrected statement below, with $q_j$ the columns of $Q$, I have $c^\top V_+c=\sum_j\max\{\lambda_j,0\}\,(q_j^\top c)^2\ge 0$ for every vector $c$.

What survives

Footnote 3 needs no change. The decomposition and the coefficient estimates are unaffected. The corrected option can change individual DCRSEs when it now leaves an indefinite matrix unchanged. A matrix that is already positive semidefinite is unchanged. Where the option already applies the eigenvalue step, the spectral formula is unchanged. The reanalysis scripts set this option. I did not examine whether any of the 174 models has an indefinite matrix with a positive diagonal, and no published number is shown to change.

The problem

The Stata script names three KEVs for Table 4, Model 7 of Bermeo and Leblang (2015). Its line 100 reads * Table 4, Model 7; KEV: listock_lag MigrationPolicy_lag policy_stock_lag. The saved results of that model have one row for each. The selector compares whole strings, and the string "listock_lag MigrationPolicy_lag" equals no parameter name. Only policy_stock_lag is kept, and two designated KEVs are left out of the sample. The merge script for the small-sample correction has the same line. The sentence on p. 192 counts correctly the sample that the script produced. The omission is in the sample.

Evidence a recalculation from the source

The evidence is a recalculation from the saved estimates of the replication files. I ran the authors' extraction function with the printed selector and with the three names separated. I find 29 rows for the study with the first and 31 with the second, in 11 models both times. The two restored rows have SER 2.218 (listock_lag) and 2.331 (MigrationPolicy_lag), and both are significant at the 5% level before and after. The counts follow by addition: $667+6+18=691$ becomes $669+6+18=693$, where 6 and 18 are the KEVs of the two studies analyzed in R, and the counts for 2015 and for IPE become 124 and 272. For the study, the script averages ratios rounded to three decimals, and $63.811/29\approx 2.200$ becomes $68.360/31\approx 2.205$. The overall average SER lies in $[1.74304,\,1.74314]$ before and in $[1.74325,\,1.74335]$ after, and it rounds to 1.743 both times. These are intervals because the means of the two R studies are known to three decimals.

What survives

The 174 models and the 22 studies are unchanged, as are all fitted coefficients and variances. Every printed SER and proportion of Table 5 keeps its two-decimal value, including 1.74, 0.22 and 0.05, and the overall average SER remains 1.743. What changes is the number of KEVs in Table 5 and in Supplementary Table A.5, the study row and three entries of the last row in Supplementary Table A.2, and Figures 1–3, which incorporate the two additional KEV estimates.

The problem

The plotting code maps the weight to the size of the point without setting a scale that makes point area proportional to weight. In the published figure a point with a larger weight is drawn larger, but the areas are not in proportion to the weights.

Evidence a recalculation from the source

The evidence is a count of pixels in the image embedded in the published PDF. The point for Horowitz and Stam (2014), variable combat, has weight $1/3$ and covers $512$ fully black and $605$ nonwhite pixels. The point for Fang, Johnson, and Leeds (2014) has weight $1$ and covers $1131$ and $1263$. However the gray edge pixels are counted, the ratio of the two areas is at most $1263/512\approx 2.47$, which is less than the ratio $3$ of the weights. In each legend the circles labelled $1.00$ and $0.25$ have an area ratio below $4$: at most $635/214$ in the top panel and $1271/430$ in the bottom panel.

What survives

The plotted $p$-values, the weighted averages and proportions in Table 5 and Supplementary Table A.2, and Figures 1 and 2 do not use the point-size scale and are unchanged.

The problem

The note describes two steps for the original standard errors, and the R script for this reanalysis departs from each. For the three Arel-Bundock (2017) models the original standard error is the HC0 standard error, and no factor is applied. Under HC0 each observation is its own cluster. This implies that $N_c$ is the number $n$ of observations in the fit, and the note calls for the factor $\sqrt{n/(n-1)}$. For the Pratt (2018) models the degrees of freedom are one less than the number of clusters in the data frame passed to the fit. The fit drops rows with missing values, and the clusters that remain are fewer in models 1, 2, 3 and 5.

Two other points qualify the comparison and are not counted as errors. The original standard errors for Pratt (2018) come from sandwich::vcovCL, whose defaults already supply $G/(G-1)$ for $G$ clusters, and $(n-1)/(n-k)$ in the least-squares fits. For the Stata studies the command multiplies the refitted original covariance by $N_c/(N_c-1)$, as the note says. For Shelef (2016) that covariance already carries $\{G/(G-1)\}\{(n-1)/(n-k)\}$, and the cluster factor enters twice.

Evidence an exact counterexample

The evidence is an exact example run through the script’s own function and a refit of the Pratt models. Let $x=(0,1,2,3)$ and $y=(0,1,1,3)$, and consider the least-squares fit with an intercept. The coefficients are $(-1/10,\,9/10)$ and the residuals are $(1/10,\,1/5,\,-7/10,\,2/5)$. A direct computation gives

The function reg_table of the script returns the standard errors $0.1503$ and $0.1435$ from the first matrix, while the note with $N_c=n=4$ requires $0.1736$ and $0.1657$ from the second. In general the squared ratio of the two is $(n-1)/n$. For the first model of Arel-Bundock (2017), with $n=3971$ (Supplementary Table A.1), the missing factor is $\sqrt{3971/3970}\approx 1.000126$.

For Pratt (2018), refitting models 1 to 5 on the data in the replication files and counting the clusters among the fitted rows, I obtain the following.

The $p$-value of the KEV in model 1 is $0.07748$ with $237$ degrees of freedom and $0.07766$ with $208$. No KEV in these five models crosses the 5% level when the degrees of freedom are corrected. Model 6 was not refitted.

What survives

This finding does not contradict the DCRSEs, the primary analysis in Table 5, or the statement of footnote 12 that the small-sample reanalysis shows no substantive differences. The two corrections bear on Table A.5 alone. There the original standard errors of the first model of Arel-Bundock (2017) grow by a factor of about $1.0001$, and the original $p$-values of the KEVs in models 1 to 5 of Pratt (2018) move by less than $0.0002$. A corrected Table A.5 has not been computed.

The problem

The passage points to a section instead of the covariance table used in its comparison. The sentence stands in Section 3 and introduces the covariance structure assumed under clustering by repeated dyad, which Table 3 displays.

Evidence a check of the printed page

I checked the page. The next paragraph begins “In Table 3”, and the parallel sentence on p. 189 reads “Table 2 demonstrates”.

What survives

The covariance comparison is unaffected by the wrong cross-reference.

The problem

Volume 22, issue 4 of Political Analysis is the Autumn 2014 issue. The publisher’s page for the article also carries the line “Published online by Cambridge University Press: 04 January 2017”, and that date is not the year of the journal issue.

Evidence a check of the printed page

The evidence is a check of the publisher’s page for the article, https://doi.org/10.1093/pan/mpt051, which gives Volume 22, Issue 4, Autumn 2014, pp. 457–463. In the same reference list, Political Analysis 23 (4) is dated 2015 and 24 (4) is dated 2016.

What survives

The article is identified correctly by its authors, title, volume, issue and pages; only the year changes, in the reference and in its three citations.

The problem

Each of the nine Shelef regressions behind this row has one set of indicators, i.regionfe_region_decade, and no indicators for the state that loses territory (loseccode) or for the state that gains it (gainccode). The first of them reads

reg midpost homelandclaim priormid coethnicity oil_gas capital_distance_log dyaddem allydummy cincrat i.regionfe_region_decade if colony2==0, cluster(dyadid4)

The options dm1(loseccode) and dm2(gainccode) of the DCR command name the dyad members for the variance and add no regressors. Region-by-decade indicators do not contain the member effects. These models are not models with fixed effects for both members of the dyad, and the subgroup that p. 192 and footnote 10 describe has two studies.

Evidence an exact counterexample

The evidence is an exact check on the data of the replication files. Observations 3 and 7 of the Shelef data file are in the estimation sample of all nine models and have the same level, 7, of regionfe_region_decade. Their members are $(255,315)$ and $(315,260)$, with the losing state first. Let $w$ be the vector with $+1$ on observation 3, $-1$ on observation 7 and $0$ elsewhere. Then $w$ is orthogonal to every indicator of regionfe_region_decade. Its inner product with the indicator of losing state 255 is $1$, and the same holds for the indicator of gaining state 315 and for the indicator that state 255 is either member. This implies that none of these three indicators is a linear combination of the region-by-decade indicators. The other regressors cannot supply the member effects either. In the three groups of models (Tables 3, 4 and 5 of the study) the design matrix has rank at most 34, 31 and 33. The indicators of the losing and of the gaining states have ranks 114 and 114, 105 and 104, and 83 and 83.

What survives

The numbers of the Shelef row are those of the nine models as fitted. I did not recompute them, and no error in them is shown. Shelef (2016) stays in the primary analysis of 22 studies, and Table 5 and Supplementary Tables A.1 and A.2 do not change. What changes is the subgroup with fixed effects for both members. It has two studies and 9 KEVs ($8+1$), not three studies and 18. The article draws no conclusion from this subgroup.

The problem

The program nclogit, supplied with dcr for conditional-logit models, begins with the statement version 16.1, which is executed before any argument is read. Stata 15 and 15.1 cannot run a program that requests a later version. Every call of dcr nclogit stops at that statement, although dcr.ado itself declares version 15. The replication uses this path for three models, one of Weisiger and Yarhi-Milo (2015) and two of Kinne (2018).

Evidence a derivation

The evidence is a reading of the source of dcr.ado and of nclogit.ado, whose two released copies are identical, together with the rule for Stata’s version command. Stata was not run.

What survives

Results produced with Stata/MP 17.0, the version named in the replication README, are unaffected, and no estimate, table or figure changes. The other estimator paths of dcr are not shown to fail on Stata 15.

The problem

The note to the table says that the statistics come from a representative reanalyzed model. The three reanalyzed models of this study are fitted to 952 donor–recipient–year observations. The printed 4760 equals $5\times952$, and the code in the replication files requests five imputations. The entry has the size of five completed copies of the sample and does not count the observations once. The entry 756 is the number of directed donor–recipient pairs among the 952 observations. With direction ignored, as in the clustering of the reanalysis (p. 192), there are 741 pairs, and the table does not say which convention it uses.

Evidence a recalculation from the source

The evidence is a recalculation from the replication files. The data file has 4624 rows, of which 952 satisfy the condition g7==1 of the three models. Each of these is a distinct donor–recipient–year, and the years run from 1976 to 2009. The 96 saved $p$-values that use a $t$ distribution (two columns for each of 16 coefficients in three models) are reproduced with $n=952$ and 936 degrees of freedom to within $2.8\times10^{-16}$, and not with $n=4760$. The 952 observations contain 756 distinct ordered pairs and 741 unordered pairs: 15 pairs of countries occur in both orders, and six of the 741 pair a country with itself. The code that filled the summary workbook is not in the replication files, and how the entry 4760 was produced cannot be traced.

What survives

The correction changes this row of Table A.1 and the corresponding cells of the summary workbook in the replication files. The years 1976 and 2009 are correct, and 756 is correct as a count of directed pairs. The estimates, standard errors, standard error ratios and significance comparisons are computed from the model outputs and do not use these two columns.

The problem

Each of the three calls runs every component of the decomposition on p. 191 under the prefix mi estimate, post:. Under the standard behaviour of this prefix, the matrix e(V) returned for a component is the covariance pooled over the five imputations by Rubin’s rules: the average of the within-imputation covariances plus $(1+1/M)B$, where $B$ is the between-imputation covariance of the coefficients and $M=5$. The command then multiplies the returned matrix by the factor that removes the finite-sample correction of that component. This factor belongs to the within-imputation covariance only. Applied to the pooled matrix, it also multiplies $(1+1/M)B$. The components carry different factors, and the between-imputation term enters the result with a coefficient $c$ that differs from 1 in general. The DCR covariance computed within each imputation and then pooled has coefficient 1. The decomposition on p. 191 is correct for one complete data set if no same-country dyad appears in more than one observation.

Evidence an exact counterexample

The evidence is a derivation from the code and an exact counterexample. No Stata run was made, and the three models were not refitted. Let $K$ be the number of dyad members (the $N$ of p. 191), $n$ the number of observations and $M$ the number of imputations. The decomposition has one component for each member $i$, the repeated-dyad component $D$ and the heteroskedasticity-consistent component $0$, with signed weights $s_i=1$, $s_D=-1$ and $s_0=-(K-2)$, so that $\sum_js_j=1$. For component $j$, let $\bar U_j$ be the average over the imputations of the within-imputation covariance without finite-sample correction, and let $a_j$ be the factor by which the command multiplies. The corrected within-imputation covariance is the uncorrected one divided by $a_j$. Suppose that no same-country dyad appears in more than one observation, that the sample, the clusters and the factors are the same in every imputation and that no eigenvalue adjustment is made. Since the fitted coefficients do not depend on the clustering, every component has the same $Q=(1+1/M)B$. The matrix returned for component $j$ is $T_j=\bar U_j/a_j+Q$, and the command computes

The DCR covariance computed within each imputation and then pooled is $T_*=\sum_js_j\bar U_j+Q$. Subtracting, I find $T_{\mathrm{code}}-T_*=(c-1)(1+1/M)B$. With the option dofundo(asylike), the factors are $a_i=1-1/G_i$, $a_D=1-1/G_D$ and $a_0=1-1/n$, where $G_i$ and $G_D$ are the numbers of clusters in the component for member $i$ and in the repeated-dyad component. Substituting these factors into the definition of $c$, I obtain

which is not zero in general.

The decomposition includes each observation paired with itself; it omits only pairs of different observations on a common same-country dyad. In the source data, I count each of the six same-country dyads just once among the 952 observations with g7==1, so none contributes such an omission.

For the counterexample, I consider 24 observations, four on each of the six dyads among four members, an intercept-only logit and five imputations. Three imputations complete the data to 12 successes and two to 18. Then $K=4$, $n=24$, $G_D=6$ and $G_i=13$ for every member, and $c=4(12/13)-5/6-2(23/24)=49/52$. The fitted intercepts are $0$ and $\log3$, and I have $B=3(\log3)^2/10$ and $Q=9(\log3)^2/25$. The DCR variances within the two kinds of imputation are $2/9$ and $8/81$, with average $14/81$. A direct computation gives $T_*=14/81+(9/25)(\log3)^2$ and $T_{\mathrm{code}}=14/81+(441/1300)(\log3)^2$, and the difference is $-(27/1300)(\log3)^2<0$. All the variances are positive, and the eigenvalue adjustment is not active.

What survives

The decomposition on p. 191 is unaffected. The reanalyses of the other 21 studies do not use this multiple-imputation prefix. The logit coefficients of the three models are unaffected, since the covariance estimator does not change them. The DCR covariance of the three models depends on this code, and with it the DCRSEs, the standard error ratios and the DCR $p$-values of the 18 KEVs in the row of Supplementary Table A.2. These enter the 2020, IPE and overall rows of Table 5 and Figures 1–3, and the same calls are repeated in the small-sample reanalysis behind Supplementary Table A.5. No value in these tables is shown to be wrong. The size and direction of the difference for the three models have not been computed, because the replication files do not keep the coefficient and covariance matrices of each imputation.