Using Multiple Outcomes to Adjust Standard Errors for Spatial Correlation

titlepage\begin{spacing}{1} \begin{abstract} \smalltonormalsize{ Empirical research in economics often examines the behavior of agents located in a geographic space. In such cases, statistical inference is complicated by the interdependence of economic outcomes across locations. A common approach to account for this dependence is to cluster standard errors based on a predefined geographic partition. A second strategy is to model dependence in terms of the distance between units. Dependence, however, does not necessarily stop at borders and is typically not determined by distance alone. This paper introduces a method that leverages observations of multiple outcomes to adjust standard errors for cross-sectional dependence. Specifically, a researcher, while interested in a particular outcome variable, often observes dozens of other variables for the same units. We show that these outcomes can be used to estimate dependence under the assumption that the cross-sectional correlation structure is shared across outcomes. We develop a procedure, which we call Thresholding Multiple Outcomes (TMO), that uses this estimate to adjust standard errors in a given regression setting. We show that adjustments of this form can lead to sizable reductions in the bias of standard errors in calibrated U.S.\ county-level regressions. Re-analyzing nine recent papers, we find that the proposed correction can make a substantial difference in practice. } \\ \\ Keywords: Spatial Correlation, Clustering, Large-Scale Inference \\ JEL: C10, C21, C31 \end{abstract} \end{spacing}

spacing{1.4} \section{Introduction} Empirical research in economics often considers data indexed by locations in space. In 2023 alone, nearly half---61 of 128---of all empirical, observational papers published in five leading, general interest journals primarily consider data with this characteristic.\footnote{We hand-labeled the 370 papers published in the American Economic Review, Econometrica, Journal of Political Economy, Review of Economic Studies, and Quarterly Journal of Economics in 2023 according to whether they are, primarily, empirical, observational, and consider data indexed by locations in space. Further details concerning the criteria used to define these categories are given in (ref). Examples of empirical analyses without a spatial aspect include those where units are products, or students in a single educational institution, or patents.} The interpretation of inferences produced using these data is complicated by a ubiquitous form of dependence: economic outcomes are typically linked across locations in potentially complex ways, through direct causal relationships or a common underlying influence. Productivity shocks to one region may propagate to other areas through trade linkages or labor market flows greenstone2010identifying,bustos2016agricultural,giroud2024propagation, fiscal policies in one municipality may affect both the design of symmetric policies and economic outcomes in other locations case1993budget,dellavigna2022policy, and political events in one jurisdiction may determine voter behavior in others besley1995incumbent,madestam2013political. Spatial dependence has the potential to induce substantial biases in statistical inferences that are made under the assumption that outcomes are independent across space. Existing methods for correcting standard errors fall into two broad categories. The first set of methods cluster standard errors based on a partition of the space. A prototypical instance of this approach allows for dependence between pairs of U.S.\ counties in the same state (for discussion, see moulton1986random, moulton1990illustration). The second set of methods parameterize dependence as a monotonic, decreasing function of geographic distance conley1999gmm, muller2022spatial, muller2023spatial. Both approaches impose strong, ex ante restrictions on how spatial dependence varies with geography, in one case with sharp boundaries, and in the other with smooth functions of distance. However, the correlation between economic outcomes at two locations need not be driven, in large or substantial part, by their geographic distance. For instance, similarities in economic or political behavior could instead arise from factors like urbanization, income, or education levels, some of which may be unobserved. \begin{figure}[p] \caption{Correlation with Two California Counties} \begin{tabular}{l | l} \textbf{Correlation w/ San Francisco County} & \textbf{Correlation w/ Modoc County} \\ \hline \textit{California} &\\ & \\ \textit{North Dakota} & \\ & \\ \textit{New York} & \\ & \\ \multicolumn{2}{c} \\ \end{tabular} \justifying { Notes: (ref) displays heat-maps giving the correlation of a collection of economic outcomes between pairs of U.S.\ counties. These outcomes are listed in (ref) and have been normalized to have mean zero and unit variance. The first column gives the correlation between outcomes in San Francisco County, California, with each county in the states of California, North Dakota, and New York. The second column is analogous, but gives correlations relative to Modoc County, California. } \end{figure} To illustrate this point, we collect a set of 91 U.S. county-level economic outcomes (including the unemployment rate, the average income, the poverty rate, etc.,) for each of the 3,144 U.S. counties.\footnote{We list these outcomes, and give further details concerning their construction, primarily from the U.S.\ Census, in (ref).} We normalize each outcome to have mean zero and unit variance across the set of counties. The left column of (ref) displays the correlation between the economic outcomes of San Francisco County, California and the outcomes of every other county in the states of California, North Dakota, and New York. As one might expect, San Francisco is more correlated with other densely populated areas, such as Los Angeles County and New York County (which contains the borough of Manhattan), than with the less populated areas in the rural state of North Dakota. The same, intuitive, patterns hold among less populated areas. Take Modoc County, California, for instance. Modoc County is the third least populous county in the state of California, with an economy substantially driven by the agricultural sector modoc2023. The right column of (ref) displays the correlation between the economic outcomes of Modoc County, and every other county in California, North Dakota, and New York. Modoc County is highly correlated with other rural, agricultural areas. In fact, Modoc County is more correlated with very distant less populated areas in North Dakota and New York than with closer, but more populated, areas in California.\footnote{These results generalize beyond Modoc county. (ref) displays the correlation in outcomes between each pair of counties in California, North Dakota, and New York. Counties are sorted by state and population. It is clear that both population and geography are important determinants of the correlation between two counties. (ref) gives an analogous visualization for each pair of counties in the United States. The same pattern holds.} Methods based on the assumption that cross-sectional dependence is determined by geographic proximity alone cannot capture the many of the main features of the estimates displayed in (ref), and so may not perform well in many regression problems considered in empirical economics. Indeed, in (ref), we calibrate a simulation to the estimates of the correlation between pairs of U.S.\ counties across the outcomes introduced above. We find that, in this setting, many standard adjustments for spatial dependence yield standard errors that are badly biased. In this paper, we introduce a method that uses observations of multiple outcomes to adjust standard errors for spatial dependence. We detail the proposal in (ref). Often, in a single study, researchers observe many outcome variables for the same set of units. We give a condition under which collections of outcomes, of the type considered above, can be used to estimate dependence across units. In particular, rather than directly restricting cross-sectional dependence to be determined by, say, geographic proximity, we assume that the structure of cross-sectional correlation is shared across outcomes. This is a highly stylized restriction---so too are alternative assumptions. The main appeal is that, in settings where this assumption is reasonable, i.e., where the main features of cross-sectional correlation are captured by a collection of available outcomes, spatial dependence can be estimated directly from the data at hand. We propose a procedure, which we call Thresholding Multiple Outcomes (TMO), for using estimates of spatial dependence to adjust standard errors in a given regression setting. The method consists of three steps. First, we estimate the correlation across locations using a collection of auxiliary outcomes. Second, we use the empirical distribution of these estimates to determine a threshold; pairs of locations whose correlation exceeds this threshold are treated as correlated, while pairs below the threshold are treated as uncorrelated. Third, we construct standard errors by allowing for dependence among location pairs treated as correlated. This approach is analogous to the construction of clustered standard errors. There, pairs of locations in the same cluster are assumed to correlate \emph{a priori}. Here, we use the auxiliary outcomes to \emph{determine} the pairs of locations that correlate. Revisiting the calibrated simulation from (ref), we find that this procedure exhibits substantially smaller bias, and more accurate rejection rates, than existing methods. Our procedure is premised on the hypothesis that, in many settings of practical interest in applied economics, pairs of locations can be effectively categorized as either “null” or “non-null.” Outcomes at non-null location pairs exhibit meaningful correlation, whereas outcomes at null pairs do not. The main idea is that the empirical distribution of the estimates of the correlation between each pair of locations can be used to distinguish null pairs from non-null pairs. We appeal to the intuition that the center of this distribution should predominately consist of the correlations associated with null pairs. Thus, we estimate the null distribution using the center, and then extrapolate this estimate into the tails to infer the analogous distribution for the non-nulls. We rely on a Gaussian approximation to the null distribution. These inferred distributions can then be used to select a threshold that appropriately balances the contributions of null and non-null pairs to the error when constructing the resultant standard errors. The main non-standard feature of our setting is that the data under consideration are correlated across both locations and outcomes. In (ref), we consider how both sources of correlation impact the performance of the TMO estimator. We present two results. First, we give a quantitative description of the extent of the correlation across locations and outcomes under which spatial dependence is consistently estimable. Second, we characterize the quality of a Gaussian approximation to the empirical distribution of the null correlation estimates in terms of the extent of the correlation across outcomes. In both cases, the main takeaway is that our ability to recover spatial dependence is adversely impacted by the use of highly correlated outcomes. Thus, when applying the TMO estimator, researchers face a trade-off: pairs of outcomes with the most similar cross-sectional correlation are likely to be, themselves, more correlated. We provide recommendations on how to best balance these considerations in the remaining two sections. In (ref), we illustrate the application of TMO standard errors to nine recent papers. We begin by documenting that the potential for spatial correlation is widespread in empirical economics. In a set of 370 economics papers published in 2023 in five leading journals, we identify 126 observational papers, 48% of which have the potential for spatial correlation to affect a main estimate in the paper. The most common geographic unit in such papers is the U.S.\ county, employed in 15 papers. We re-analyze seven of these papers, which differ in their economic context and underlying methodology. For each paper, we identify a set of county-level outcomes, taken from each paper's replication packages as well as from external sources. We use these outcomes to construct the TMO estimator. When applied to these papers, we find that the TMO estimator can substantively change the estimated standard errors, with a median increase of 37%. We then compare TMO to several leading alternatives based on modeling dependence in terms of geographic proximity: clustered standard errors, conley1999gmm, and muller2022spatial,muller2023spatial. In some cases, the TMO standard error leads to a similar adjustment. In others, the resultant standard errors are qualitatively different. To show that the method is applicable to other units of observation, we apply the TMO estimator to two additional papers---one that studies commuting zones and another that studies countries. We conclude in (ref) by detailing our recommendations for practice. We emphasize rules-of-thumb concerning what types of data our method is best suited for, diagnostics for determining that the method is performing as intended, and further illustrations of principals for choosing a relevant collection of outcomes. Auxiliary Figures and Tables are displayed in (ref). Further details and extensions are given in (ref). Proofs for all results stated in the main text are given in (ref); Proofs for supporting Lemmas are given in (ref). Details concerning our empirical applications are given in (ref). Code implementing the method developed in this paper is available at the link \url{https://github.com/wjnkim/tmo}. \subsection{Related Literature} We contribute to an extensive econometric literature concerning quantification of uncertainty in data that exhibit cross-sectional dependence. Much of this research focuses on methods that partition units into clusters according to an observable characteristic, allowing for correlation among observations within the same cluster. liang1986longitudinal introduce the first heteroskedasticity-consistent clustered standard error (see also mackinnon1985some, bell2002bias, and imbens2016robust for alternatives that exhibit better finite-sample performance). bertrand2004much provide a compelling demonstration of the practical importance of these methods in panel data. cameron2015practitioner give a comprehensive review. Like this approach, we assume that only a small proportion of pairs of units correlate.\footnote{Available methods for constructing clustered standard errors with large clusters are either conservative ibragimov2010t,ibragimov2016inference or rely an assumption that error covariance matrices are identical across clusters canay2017randomization,canay2021wild.} By contrast, we do not assume that these pairs are know \textit{ex ante}, but rather, can be learned from observations of auxiliary outcomes. A smaller, but widely applied, subset of this literature develops methods premised on the assumption that dependence can be modeled as a function of the distance between units. This line of work was initiated by conley1999gmm, who emphasizes that dependence should be modeled by general measures of “economic distance” (see also kelejian2007hac, kim2011spatial, bester2016fixed, pollmann2020causal, colella2023acreg and muller2022spatial,muller2023spatial for alternative estimators based on this program). Despite this, most applications of this framework parameterize dependence in terms of geographic proximity alone. Separately, barrios2012clustering, muller2024spatial, and conley2025standard demonstrate that applying standard variance estimators to spatially correlated data often yields spurious findings. The simulation experiments conducted in (ref) generate results that are consistent with this view. Despite this, our outlook is pragmatic. The core of our proposal is that, due to the widespread availability of additional outcome data, at the very least, some of the cross-sectional correlation that is common across measured outcomes can be controlled. Our approach for identifying pairs of locations that correlate draws on an extensive literature in multiple hypothesis testing, particularly as developed for DNA micro-array analyses. See efron2012large for a textbook treatment. As our primary aim is to construct a single standard error, the criteria that we use to distinguish nulls from non-nulls differs from the methods developed in that literature (see also bailey2016two and bailey2019multiple for connections between multiple hypothesis testing and the construction of spatial standard errors). Moreover, we emphasize simple methods that require minimal nonparametric modeling. More sophisticated methods for estimating null distributions, as well as optimality theory, are discussed in langaas2005estimating, efron2007size, jin2007estimating, and cai2010optimal, among other places. Finally, we contribute to a growing literature that considers how standard statistical procedures can be improved by incorporating measurements of multiple outcomes. For example, ludwig2019machine, sun2023using, and abadie2024doubly each consider how collections of outcomes can be used to improve estimates of causal effects. Our primary focus, by contrast, is on obtaining more accurate estimates of uncertainty in regression problems. \section{Standard Adjustments for Spatial Dependence are Biased} In this section, we use a simulation experiment to measure the performance of several methods for adjusting standard errors for spatial dependence. We consider settings with the following structure. A researcher observes the data $(Y_i, W_i,X_i)$ for each unit $i$ in $1,\ldots,n$, where $Y_i$ is a scalar outcome, $W_i$ is a scalar treatment, and $X_i$ is a vector of covariates. For the sake of concreteness, $Y_i$ could measure the change in the income per capita in U.S.\ county $i$ over a given time period, $W_i$ could denote the change in the percent of the population in county $i$ that has completed a college degree, and $X_i$ could collect a set of state fixed effects. The researcher assumes that the data arise from the linear model \begin{equation} Y_i = \alpha + \tau W_i + \theta^\top X_i + \varepsilon_i,\quad\mathbb{E}[\varepsilon_i|W,X]=0,\quad i = 1,\ldots,n , \end{equation} where $W=(W_i)_{i=1}^n$ and $X=(X_i)_{i=1}^n$ collect the treatment and covariate measurements over the full sample of units. Let $\hat{\tau}$ denote the ordinary least squares estimate of the parameter $\tau$. Our interest is in estimating an appropriate standard error for $\hat{\tau}$. Often, it is unreasonable to assume that the components of the residual vector $\varepsilon = (\varepsilon_i)^n_{i=1}$ are uncorrelated. For most settings of interest in applied economics, factors that determine outcomes in San Francisco, say, plausibly also determine outcomes in Berkeley. There are two common ways of adjusting standard errors for spatial dependence. One set of methods is based on clustering standard errors based on a partition of the space, thereby imposing geographic borders on the scope of spatial dependence moulton1986random, moulton1990illustration. A second set of methods parameterize spatial dependence smoothly in terms of geographic distance conley1999gmm, muller2022spatial, muller2023spatial. How well do these methods address the types of spatial dependence exhibited by the data considered in applied economics? An answer to this question requires the specification of a “typical” correlation structure. In (ref), we measure the main features of the spatial dependence exhibited by the types of U.S. county-level outcomes considered in applied economics.\footnote{We focus our discussion on U.S.\ counties because this is the most common level unit of aggregation considered in leading general interest economics journals in 2023. In particular, 12% of all observational papers in this sample take U.S.\ counties as their unit of observation. This proportion increases to 25% if the sample is restricted to papers whose unit of observation has the potential for spatial correlation. Further details are given in (ref).} In (ref), we use these measurements to calibrate a simulation and document the performance of several widely applied approaches to adjusting standard errors for spatial dependence. \subsection{Measuring The Structure of Spatial Dependence} What types of spatial dependence should we expect in the economic outcomes of U.S.\ counties? The basic premise of this paper is that non-trivial correlation structures can be inferred from measurements of multiple outcomes. To operationalize this idea, we collect 91 U.S. county-level outcomes. This is the same set of outcomes considered in (ref). These outcomes, which include measurements of income, poverty, employment, crime, and health, and are intended to capture the main features of the types of outcomes considered in applied economics. We list these outcomes and give further details concerning their construction in (ref). Given a “representative” outcome $Y_i$ and treatment $W_i$, and a pair of counties $i$ and $i^\prime$, should we expect the residuals $\varepsilon_i$ and $\varepsilon_{i^\prime}$ to correlate? We argue that an answer to this question can be recovered by averaging the product of the empirical versions of these residuals across the outcomes in our sample. Formally, let $Y^{(1)}_i, \ldots, Y^{(d)}_i$ denote the $d$ outcomes measured for county $i$. As a working example, we take the treatment of interest $W_i$ to be the change in the percentage of the population in county $i$ that has completed a college degree from 1980 to 2009. For each outcome $Y_i^{(j)}$, let $\hat{\varepsilon}^{(j)}_i$ be the empirical residuals from the least squares regression of $Y^{(j)}_i$ on $W_i$ and a vector of state fixed effects. Let $\tilde{\varepsilon}^{(j)}_i$ denote a version of $\hat{\varepsilon}^{(j)}_i$ that has been standardized to have variance one across units, given by \begin{equation} \tilde\varepsilon_{i}^{(j)} = \hat{\varepsilon}^{(j)}_i \Biggl/\sqrt{\frac{1}{n} \sum_{i=1}^n (\hat{\varepsilon}^{(j)}_i)^2} . \end{equation} The empirical covariance and correlation between the empirical residuals for counties $i$ and $i^\prime$ \emph{across outcomes} are given by \begin{equation} \hat{\lambda}_{i,i^{\prime}} = \frac{1}{d} \sum_{j=1}^d (\tilde{\varepsilon}^{(j)}_i - \bar{\varepsilon}_i) (\tilde{\varepsilon}^{(j)}_{i^{\prime}} - \bar{\varepsilon}_{i^{\prime}}) \quad \text{and} \quad \hat{\rho}_{i,i^\prime} = \frac{\hat{\lambda}_{i,i^{\prime}}}{\sqrt{\hat{\lambda}_{i,i} \hat{\lambda}_{i^{\prime},i^{\prime}}}} , \end{equation} respectively, where $\bar{\varepsilon}_i=\frac{1}{d} \sum_{j=1}^d \tilde{\varepsilon}^{(j)}_i$ is the average of the normalized residuals $\tilde{\varepsilon}^{(j)}_i$ across outcomes for county $i$. Intuitively, the statistic $\hat{\rho}_{i,i^\prime} $ captures the tendency for residuals associated with counties $i$ and $i^\prime$ to be correlated. (ref) displays a histogram of the correlation estimates $\hat{\rho}_{i,i^\prime}$ over each pair of distinct U.S.\ counties. A Gaussian density function, centered at zero, is overlaid in blue. The scaling of this density function has been chosen to best fit the quartiles of the distribution of the measurements $\hat{\rho}_{i,i^\prime}$. The left and right tails of the distribution are displayed as enlarged insets. Two features stand out: The Gaussian distribution does an exceptional job at approximating the center, but under-estimates the right-tail, of the empirical distribution of the correlation estimates. \begin{figure}[t] \begin{centering} \caption{The Distribution of Correlations Between Pairs of U.S.\ Counties} \begin{tabular}{c} \end{tabular} \end{centering} \justifying {Notes: (ref) displays a histogram of the estimates $\hat{\rho}_{i,i^\prime}$ over each pair of distinct U.S.\ counties. A mean-zero, Gaussian density function is overlaid in blue. The scaling of this density function has been chosen to best fit the corresponding quartiles of the distribution of the estimates $\hat{\rho}_{i,i^\prime}$. The left and right tails of the distribution are displayed as enlarged insets.} \hrulefill \end{figure} These observations can be rationalized as follows. Suppose that the correlation estimates $\hat{\rho}_{i,i^\prime}$ are approximately Gaussian. This will hold if the number of outcomes $d$ is sufficiently large and there is not too much correlation across residuals. (We will have more to say about this in (ref).) The excellent approximation of a mean-zero, Gaussian distribution to the central quantiles of the correlation estimates suggests that most pairs of counties are not systematically correlated. By contrast, the poor approximation of the Gaussian distribution to the right tail of the correlation estimates suggests that there is a relatively small proportion of pairs of counties that are systematically positively correlated. It is worth scrutinizing, however, whether the right-tail of the distribution of correlation estimates captures pairs of counties that are genuinely correlated. Consider the set of county pairs $(i,i^\prime)$ in which the correlation estimate $\hat{\rho}_{i,i^\prime}$ exceeds the positive threshold $\delta$. (ref) reports, for several values of $\delta$, the proportion of county pairs $(i,i^\prime)$ in this set in which county $i$ is among the closest 10% of counties to $i'$ on a collection of observable characteristics. Proximity on observable characteristics is very predictive of having highly correlated residuals. Geographic proximity does not tell the whole story, however. Demographic and socioeconomic proximity, e.g., income, population, and population density, are similarly predictive of correlations in residuals. \begin{table}[t] \begin{centering} \caption{Predictors of Highly-Correlated County Pairs} \begin{tabular}{p{10cm} c c c} \toprule \multirow[b]{2}{10.5cm}{\textit{Among county pairs $(i,i')$ with $\hat{\rho}_{i,i'} > \delta$,\\ percentage in which $i$ ranks among the 10% closest to $i'$ in:}} & \multicolumn{3}{c}{$\delta$} \\ \cmidrule(lr){2-4} & 0 & 0.2 & 0.4 \\ \midrule Population & 13 & 21 & 41 \\ Urban % & 14 & 20 & 41 \\ Geographic Distance & 12 & 17 & 32 \\ Median income & 13 & 17 & 31 \\ Non-white % & 12 & 15 & 24 \\ Vote-share & 12 & 14 & 22 \\ Any of the above & 51 & 62 & 83 \\ \bottomrule \end{tabular} \end{centering} \justifying {Notes: (ref) shows the percent of county-pairs $(i,i^\prime)$ in which county $i$ is within the top 10% of counties closest to $i'$ on a collection of observable characteristics, among the county-pairs that have a correlation greater than a threshold $\delta$. The threshold $\delta$ increases across columns. The percentages of county-pairs with estimated correlations above 0, 0.2, and 0.4 are 49.3%\unskip, 6.4%\unskip, and 0.3%\unskip, respectively.} \hrulefill \end{table} \subsection{Calibrated Simulation} The estimates displayed in (ref) suggest that the outcomes of at least some pairs of U.S.\ counties are systematically correlated. The next question is whether this matters. We assess, through a calibrated simulation, whether ignoring these correlations leads to biased standard error estimates. Specifically, we use the estimates $\hat{\rho}_{i,i^\prime}$ to parameterize a covariance matrix $\Sigma$ for the residual vector $\varepsilon$ in the specification (ref). We assume that the residual vector $\varepsilon$ has a Gaussian distribution with mean zero and variance $\Sigma$. Let $\sigma_{i,i^\prime}$ denote the $i,i^\prime$th component of $\Sigma$. Our goal is to assess the performance of standard methods that ignore spatial dependencies, as well as the performance of conventional methods for adjusting standard errors that are based on modeling spatial correlation in terms of geographic distance. To do this, we must meet two objectives. First, the non-zero off-diagonal elements of $\Sigma$ should primarily consist of pairs of counties whose correlation estimate $\hat{\rho}_{i,i^\prime}$ is likely to reflect real dependence. Second, the covariance matrix $\Sigma$ must be a proper, positive-definite, covariance matrix, so that it can be used to draw simulation replicates. We develop a simple approach that achieves these goals. We set $\sigma_{i,i^\prime} = 0$ for all pairs of counties with $\vert \hat{\rho}_{i,i^\prime} \vert < 0.45$. We then implement an algorithm, described in (ref), to produce a positive-definite matrix by forming “clusters” among the remaining pairs of counties that have $\vert \hat{\rho}_{i,i^\prime} \vert \geq 0.45$. Within these clusters, which contain 39% of these highly correlated pairs, we set $\sigma_{i,i^\prime} = \hat{\rho}_{i,i^\prime}$. The simulation is structured as follows. For each of 1,000 simulation replicates, we draw a random vector $\varepsilon$ from $\mathsf{N}(0,\Sigma)$, set $Y = \beta W + \varepsilon$, where $\beta = 0$ and $W$ is a specified treatment vector, and apply various methods to (i) estimate the standard error associated with the least squares regression of $Y$ on $W$ and (ii) test the null hypothesis that $\beta = 0$ at level $\alpha = 0.05$. (ref) displays the results of this exercise. We consider two choices of the treatment vector $W$: The change in the percentage of the population in county $i$ that has completed a college degree from 1980 to 2009 and the change in the per-acre value of farm-land. The first row displays the bias for heteroskedasticity-consistent standard errors (HC1). Across the 1,000 simulation draws, the median estimated standard error is only 0.44 of the true standard error, leading to a rejection rate of 0.39. Conventional methods to address spatial correlation, for example clustering by state (row 2) or using conley1999gmm standard errors with a 150-mile radius (row 3), shows little improvement, with the median standard error at still less than half of the true standard error.\footnote{conley1999gmm standard errors are consistent only if the number of location pairs within the bandwidth is small relative to the sample size. In our setting, a bandwidth of 150 miles includes 3.1% of county pairs. Setting the bandwidth to 300 miles, say, includes 11.3% of county pairs, moving beyond the regime where we should expect consistency to hold. The methods proposed in bester2016fixed and muller2022spatial, muller2023spatial, by contrast, control size under various restrictions in strongly correlated data.} These adjustments have no impact on the bias in the standard error estimates (This can be attributed to the fact that $\Sigma$ is based on the covariance of residuals that have already been been demeaned by state). The fourth row displays the performance of the Spatial Correlation Principal Components (SCPC) method of muller2022spatial, muller2023spatial, which is based on parameterizing spatial correlation in terms of geographic distance. The results here exhibit a sizable improvement over the previous methods. The bias of the median estimated standard error is reduced to around 40% of the the true standard error and the rejection rate falls to 0.24. \begin{table}[t] \begin{centering} \caption{Comparison of Methods for Adjusting for Spatial Dependence} \begin{tabular} {l c c c c} \toprule & \multicolumn{2}{c}{$ \Delta$ % College Grad} & \multicolumn{2}{c}{$ \Delta$ Per-acre Farm Value} \\ \cmidrule(lr){2-3} \cmidrule(lr){4-5} & $ \text{Mean}\left(\frac{\text{Est. SE}}{\text{True SE}}\right)$ & Rej. rate & $ \text{Mean}\left(\frac{\text{Est. SE}}{\text{True SE}}\right)$ & Rej. rate \\ Method & (1) & (2) & (3) & (4) \\ \hline (1) HC1 & 0.44 & 0.39 & 0.54 & 0.29 \\ (2) Cluster by state & 0.47 & 0.36 & 0.57 & 0.25 \\ (3) Conley (150mi) & 0.47 & 0.37 & 0.57 & 0.26 \\ (4) SCPC & 0.64 & 0.24 & 0.71 & 0.17 \\ \hline (5) TMO & 0.77 & 0.14 & 0.76 & 0.12 \\ \bottomrule \end{tabular} \end{centering} \justifying {Notes: (ref) reports the performance of various methods in the calibrated simulation experiment outlined in (ref). Columns (1) and (3) report the mean standard error estimate, over simulation replications, as a proportion of the true standard error. Columns (2) and (4) reports the rejection rate at level $ 0.05$ of the associated test that $ \beta=0$. Results are from 1,000 simulation draws. The treatments of interest are the change in the percent of college graduates from 1980-2009 (Columns 1-2) and the change in the per-acre value of farmland (Columns 3-4). \say{Robust} refers to HC1 heteroskedasticity consistent standard errors hinkley1977jackknifing. \say{Conley (150mi)} refers to the standard error estimator proposed by conley1999gmm, implemented with a bandwidth of 150 miles. \say{SCPC} refers to the Spatial Correlation Principal Components method proposed by muller2022spatial, muller2023spatial. \say{TMO} refers to the Thresholding Multiple Outcomes estimator proposed in (ref).} \hrulefill \end{table} \section{Using Multiple Outcomes to Adjust Standard Errors} The results reported in (ref) are concerning. In a setting calibrated to approximate the spatial dependence exhibited by outcomes frequently studied in applied economics, we find that standard methods can be severely biased. In this section, we propose an alternative method. Our method takes as input a collection of outcomes, of the sort constructed in (ref). This collection is used to recover an estimate of spatial dependence. This estimate is used to adjust standard errors. As before, we consider the linear model \begin{equation} Y^{(0)}_i = \alpha + \tau^{(0)} W_i + \varepsilon^{(0)}_i,\quad i = 1,\ldots,n , \end{equation} where $Y^{(0)}_i$ measures some outcome and $W_i$ denotes some treatment of interest. Extensions to settings with covariates, instruments, and panel data are treated in (ref). Again, the statistic $\hat{\tau}^{(0)}$ denotes the ordinary least squares estimate of $\tau^{(0)}$ in the specification (ref) and we are interested in estimating an appropriate standard error for $\hat{\tau}^{(0)}$. For the purposes of this paper, we interpret this problem as desiring an estimate of the conditional variance \begin{equation} V(\Sigma^{(0)}) = \Var(\hat{\tau}^{(0)} \mid W) = S^{-2}_n\left( W^\top \Sigma^{(0)} W\right) , \end{equation} where $W = (W_i)^n_{i=1}$ collects the observations $W_i$, $\Sigma^{(0)} = \Var(\varepsilon^{(0)} \mid W)$, and $S_n = W^\top W$. As a consequence, estimation of the conditional variance (ref) reduces to estimation of the error covariance matrix \begin{equation} \Sigma^{(0)} = (\sigma^{(0)}_{i,i^\prime})_{i,i^\prime=1}^n,\quad\text{where}\quad\sigma^{(0)}_{i,i^\prime} = \Cov(\varepsilon^{(0)}_{i}, \varepsilon^{(0)}_{i^\prime} \mid W) . \end{equation} This conditional perspective is shared by, e.g., muller2022spatial, muller2023spatial. Design based settings, where interest is in estimating the randomness in the estimator generated by the treatment $W$, are worth further consideration abadie2020sampling,abadie2023should. \subsection{Multiple Outcomes and Proportionality} In absence of further assumptions on its structure, the residual covariance matrix $\Sigma^{(0)}$ is not recoverable. As a consequence, in practice, researchers impose a variety of restrictions on the structure of cross-sectional correlation. For example, in many settings, it is common to assume that all of the off-diagonal elements of the matrix $\Sigma^{(0)}$ are equal to zero white1980heteroskedasticity. In cases where this is less credible, researchers often assume that units can be partitioned into clusters on the basis of some shared characteristic moulton1986random,moulton1990illustration,liang1986longitudinal, or that cross-sectional correlation can be parameterized smoothly in terms of geographic proximity conley1999gmm,muller2022spatial,muller2023spatial. We take an alternative route. Our approach is premised on the observation that, often, in a single study, researchers observe many outcome variables for the same set of units. For example, a researcher, interested in the correlation across U.S.\ counties between average income and the proportion of the population with a college education, could also feasibly observe all of the outcomes considered in (ref). In (ref), we give several concrete examples from empirical economics with this characteristic. Formally, for each unit $i$, in addition to the outcome of interest $Y_i^{(0)}$, we assume that we also observe a collection $Y^{(1)}_i,\ldots,Y^{(d)}_i$ of $d$ auxiliary, post-treatment outcomes. Let $\varepsilon_{i}^{(j)}$ and $\hat{\varepsilon}_{i}^{(j)}$ denote the population and empirical residuals associated with unit $i$ when outcome $Y^{(j)}_i$ replaces $Y^{(0)}_i$ in the linear model (ref), respectively. The basic idea is to use the auxiliary outcomes to estimate the residual correlation matrix associated with the outcome of interest. To rationalize this approach, we must impose some conditions. First, we impose a restriction that implies that the residual \emph{correlation} matrix is shared across outcomes. We call a collection of outcomes with this property \emph{proportional}. \begin{assumption}[Proportionality] There exist matrices $\Lambda_{n}=(\lambda_{i,i^{\prime}})_{i,i^{\prime}=1}^{n}$ and $\Gamma_{d}=(\gamma_{j,j^{\prime}})_{j,j^{\prime}=0}^{d}$ such that \begin{equation} \sigma^{(0)}_{j,j^\prime} = \Cov(\varepsilon_{i}^{(j)}, \varepsilon_{i^\prime}^{(j^\prime)} \mid W) = \gamma_{i,i^\prime} \lambda_{j,j^\prime} , \end{equation} for each pair of units $i,i^\prime$ and outcomes $j,j^\prime$, respectively. By convention, we set $\gamma_{0,0}=1$. \end{assumption} The restriction (ref) means that the covariance between the residuals associated with two units and a pair of, potentially distinct, outcomes depends only on the outcomes through a constant factor.\footnote{(ref) is testable. guggenberger2023test propose a computationally efficient test and illustrate its application to a setting related to identification robust inference for linear instrumental variables regression.} As a consequence, the empirical residuals for the auxiliary outcomes can be used to estimate the residual covariance matrix $\Sigma^{(0)}$. In particular, let $\tilde\varepsilon_{i}^{(j)}$ denote a version of the empirical residual $\hat{\varepsilon}^{(j)}_i$ that has been normalized to have variance one across units. That is, the normalized residual $\tilde\varepsilon_{i}^{(j)}$ is given by \begin{equation} \tilde\varepsilon_{i}^{(j)} = \hat{\varepsilon}^{(j)}_i \hat{\gamma}_{j}^{-1/2} , \quad\text{where}\quad \hat{\gamma}_{j} = \frac{1}{n} \sum_{i=1}^n (\hat{\varepsilon}^{(j)}_i)^2 \end{equation} denotes the sample variance of the residuals for the $j$th outcome. The covariance between these, normalized, residuals for units $i$ and $i^\prime$, across outcomes, is given by \begin{equation} \hat{\lambda}_{i,i^{\prime}} = \frac{1}{d} \sum_{j=1}^d (\tilde{\varepsilon}^{(j)}_i - \bar{\varepsilon}_i) (\tilde{\varepsilon}^{(j)}_{i^{\prime}} - \bar{\varepsilon}_{i^{\prime}}) , \quad\text{where}\quad \bar{\varepsilon}_i = \frac{1}{d} \sum_{j=1}^d \tilde\varepsilon_{i}^{(j)} \end{equation} denotes the average normalized residual for $i$th unit. The associated correlation is given by \begin{equation} \hat{\rho}_{i,i^\prime} = \frac{\hat{\lambda}_{i,i^{\prime}}}{\sqrt{\hat{\lambda}_{i,i} \hat{\lambda}_{i^{\prime},i^{\prime}}}} . \end{equation} The statistics $\hat{\rho}_{i,i^\prime}$ are exactly the correlation estimates considered in (ref). Under (ref), the correlation estimates $\hat{\rho}_{i,i^\prime}$ recover the residual covariances $\sigma^{(0)}_{j,j^\prime}$ up to a scale factor. Proportionality is a highly stylized condition and should be viewed as only ever holding approximately for any given collection of outcomes. The same holds, of course, for assumptions that decompose spatial correlation into clusters or parameterize spatial correlation in terms of geography. From a practical perspective, by imposing (ref), we appeal to the intuition that the cross-sectional correlation exhibited by the auxiliary outcomes captures the main features of the cross-sectional correlation of the outcome of interest. On the other hand, we are ignoring aspects of spatial correlation that might be particular to any given outcome. Moreover, we are assuming that these particularities are idiosyncratic and, when considered in aggregate, negligible. Assumptions equivalent to proportionality have been considered in other fields of applied statistics. In particular, a literature in the analysis of micro-array experiments, initiated by efron2010correlated, develops procedures for multiple hypothesis testing in proportional data. See, e.g., muralidharan2010detecting, allen2012inference, and hoff2016limitations for further discussion.\footnote{More generally, dawid1981some gives a classical discussion of applications of models satisfying (ref) to Bayesian inference. gupta2018matrix gives a textbook treatment of related distributions. See also zhou2014gemini and hornstein2018joint for more recent analyses of related settings.} Proportionality, however, is not enough. In fact, in the absence of further restrictions, any estimator of the conditional variance (ref) can still be made arbitrarily inaccurate. We formalize this statement as follows. \begin{theorem} Let $\widecheck{V}(Y,W)$ be any measurable, locally bounded function of the outcomes $Y = (Y^{(j)})_{j=0}^d$ and treatments $W$. For any constants $0<\eta<1$ and $0<K$ and $W$ with $W^\top W \neq 0$, there exists a conditional distribution for $Y$, that satisfies (ref), such that either \begin{align} P\left\{ \widecheck{V}(Y, W) - V(\Sigma^{(0)}) \geq K \mid W \right\} \geq 1-\eta \quad\text{or}\quad P\left\{ V(\Sigma^{(0)}) - \widecheck{V}(Y, W) \geq K \mid W \right\} \geq 1-\eta \end{align} holds. \end{theorem} This result can be understood with a simple example. If we observe a single realization of a random variable $X$, then there is no way to estimate its variance. Translated to our context, without further conditions, there is no way to recover changes in residual covariance matrix $\Sigma^{(0)}$ that are proportional to $WW^\top$. Thus, some further restriction must be made for the standard error of the least squares estimator $\hat{\tau}^{(0)}$ to be estimable. \subsection{Sparsity and Thresholding} Our approach is based on the hypothesis that, in many settings of practical interest in applied economics, the residuals for \emph{most} pairs of units are \emph{not} meaningfully correlated. That is, pairs of units can be effectively categorized as “nulls” and “non-nulls.” The components of the residual covariance matrix $\Sigma^{(0)}$ that correspond to the null pairs are equal to zero, the components that correspond to non-nulls can be large, and the nulls outnumber the non-nulls. In other words, the covariance matrix $\Sigma^{(0)}$ is sparse. This is the same idea that underlies clustering. There, pairs of units in the same cluster are known to be non-null \textit{a priori}. Here, we use the auxiliary outcomes to \emph{determine} the pairs of units that are non-nulls. There are many ways that one might make such a determination. We propose to choose the pairs of units whose empirical residual correlation $\hat{\rho}_{i,i^\prime}$ is large in absolute value.\footnote{An alternative approach, not pursued in this paper, is to use the estimates $\hat{\rho}_{i,i^\prime}$ to learn clusters of units. This is worth further consideration, although it has the disadvantage that, necessarily, some geographically proximate units are determined to not correlate. cao2024inference pursue a related idea using a single outcome.} There are two reasons for this. First, correlation coefficients are self-normalized. That is, the covariances $\hat{\lambda}_{i,i^\prime}$ and $\hat{\lambda}_{i,i^{\prime\prime}}$ could be very different simply because the residuals for unit $i$ are are much noisier than the residuals for unit $i^{\prime\prime}$. The correlations $\hat{\rho}_{i,i^\prime}$ and $\hat{\rho}_{i,i^{\prime\prime}}$, by contrast, are comparable. Second, if the pair of units $i,i^\prime$ is a null, \textit{i.e.,} the associated residuals are not correlated, then the distribution of the statistic $\hat{\rho}_{i,i^\prime}$ can be characterized. Knowledge of the null distribution, then, allows us to discriminate between nulls and non-nulls. We treat this point in detail in (ref). Before turning to this, however, we spell out how we can use the identified non-nulls to adjust standard errors. Suppose that we were willing to determine that all pairs of units whose absolute residual correlation is greater than some threshold $\delta$ are non-nulls. That is, we say that the pair of units $i,i^\prime$ is a non-null if $\vert \hat{\rho}_{i,i^\prime}\vert \geq \delta$. How should this inform how we construct standard errors for the regression problem (ref)? One sensible estimator of the covariance is given by \begin{equation} \hat{\sigma}^{(0)}_{i,i^\prime}(\delta) = \begin{cases} \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_i ,& i = i^\prime, \\ \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}\text \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \geq \delta \} , & i\neq i^\prime . \end{cases} \end{equation} Recall that the outcomes $Y_i^{(0)}$ are not used to construct the correlations $\hat{\rho}_{i,i^\prime}$. Collect these estimates into the matrix $\widehat{\Sigma}^{(0)}(\delta)$ and plug this matrix into the expression (ref), giving the estimator \begin{equation} \widehat{V}(\delta) = S^{-2}_n \left( W^\top \widehat{\Sigma}^{(0)}(\delta) W\right) , \end{equation} where we recall that $S_n = W^\top W$. Written differently, the components of the covariance matrix $\widehat{\Sigma}^{(0)}$ for the null pairs are set to zero; the components for the non-null pairs are determined entirely by the empirical residuals associated with the outcome of interest. This is analogous to standard heteroskedasticity or cluster robust standard errors which use term like $ \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime} $ to estimate non-zero elements of residual covariance matrices white1980heteroskedasticity,liang1986longitudinal. The only difference is that, here, the auxiliary outcomes are used to determine which elements of the residual covariance matrix will be non-zero. As a consequence, the estimator (ref) is minimally sensitive to (ref)---its value is affected by the auxiliary outcomes only through the choice of which covariances to zero out. We refer to the estimator (ref) as the “TMO” variance estimator, for “Thresholding Multiple Outcomes.” In the formulation (ref), the TMO estimator (ref) does not threshold diagonal elements of the residual covariance matrix. That is, in effect, we augment the standard white1980heteroskedasticity heteroskedasticity consistent variance estimate. The same idea can be applied to wrap our approach around other existing variance estimators, particularly those that explicitly model cross-section dependence in terms of distance. For instance, an analogous estimator could be constructed, where elements that correspond to pairs of units in the same cluster, e.g., state, are never thresholded to zero. The spatial standard errors proposed by conley1999gmm and muller2022spatial, muller2023spatial can be augmented analogously. See (ref) for further discussion. We give a detailed comparison of the practical performance of these alternatives in (ref). \subsection{Threshold Choice} How should we choose the threshold $\delta$? Many estimators of sparse covariance matrices considered in the literature involve thresholding each element of an initial estimate at critical values analogous to standard corrections for multiple hypothesis testing, such as the Bonferroni procedure cai2010optimalconv,cai2011adaptive,bailey2019multiple. That is, in these formulations, $\delta$ is chosen so that the probability of mistakenly retaining even one null correlation is kept below a fixed level. If the objective is to consistently estimate a growing covariance matrix, then this is reasonable. If too many nulls are never thresholded, the aggregate error in the estimate will never converge. In this section, we argue that, in fixed samples, thresholds constructed in this way can be too large to be practically useful. That is, rather than targeting consistent estimation of $\Sigma^{(0)}$, per se, we seek to choose a threshold $\delta$ that improves the downstream estimate of the regression variance $V(\Sigma^{(0)})$ for a fixed collection of units and outcomes. This task entails balancing the contributions of the null and non-null pairs of units to the error of the resultant estimator. To see this, consider the decomposition \begin{flalign} \widehat{V}(\delta) - V(\Sigma^{(0)})& = S^{-2}_n \left( \sum_{i=1}^n W_i^2 (\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_i - \sigma_{i,i}^{(0)}) + \mathsf{T}(\delta) + \mathsf{F}(\delta) \right) ,\quad\text{where} \\ \mathsf{T}(\delta) & = \sum_{i,i^\prime \in \mathcal{T}}W_i W_{i^\prime} \left( (\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime} - \sigma_{i,i'}^{(0)}) \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \geq \delta \} - \sigma_{i,i'}^{(0)} \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \leq \delta \}\right)\quad\text{and}\nonumber \\ \mathsf{F}(\delta) &= \sum_{i,i^\prime \in \mathcal{F}}W_i W_{i^\prime} \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime} \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \geq \delta \} .\nonumber \end{flalign} Here, the sets $\mathcal{T}$ and $\mathcal{F}$ collect non-null and null pairs of units, respectively. Optimization of the expectation, or any higher order moment, of the error (ref) requires knowledge of several unknown quantities, e.g., the distribution of the statistics $\hat{\rho}_{i,i^\prime}$ and the values of the true residual covariances $ \sigma_{i,i}^{(0)}$. To make this problem tractable, we apply several, heuristic, approximations. Taken in aggregate, these approximations will bias our analysis toward smaller values of $\delta$, giving, more often than not, more conservative downstream standard errors. Let $p_0$ denote the proportion of pairs of units that are nulls. Let $f_0(\cdot)$ denote the density the absolute empirical correlation $\vert \hat{\rho}_{i,i^\prime}\vert$ for a randomly selected null pair $i,i'$. Let $f_1(\cdot)$ denote the analogous density for the non-null pairs. Consider the loss function \begin{align} L(\delta) &= \widetilde{\mathsf{T}}(\delta) +\widetilde{\mathsf{F}}(\delta),\quad\text{where}\\ \widetilde{\mathsf{T}}(\delta) & = (1 - p_0) \int_0^\delta t\text{ d} f_1(t) \quad\text{and}\quad \widetilde{\mathsf{F}}(\delta) = p_0 \int_\delta^\infty t \text{ d} f_0(t) .\nonumber \end{align} The function (ref) abstracts away, or ignores, several features of the error (ref). First, any heterogeneity across pairs $i,i'$ generated by differences in the term $W_i W_{i'}$ is dropped. Moreover, for the non-nulls, when $\vert \hat{\rho}_{i,i^\prime}\vert \geq \delta$, the difference $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime} - \sigma_{i,i'}^{(0)}$ is taken, in effect, to be zero.\footnote{As the statistics $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$ and $\hat{\rho}_{i,i^\prime}$ are likely to be positively correlated, on the event $\vert \hat{\rho}_{i,i^\prime}\vert \geq \delta$, the statistic $\vert \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}\vert $ is likely to be be positively biased for $\vert \sigma_{i,i'}^{(0)}\vert$. We are ignoring this bias, appealing to the intuition that when $\hat{\rho}_{i,i^\prime}$ is reasonably precisely estimated, its correlation with $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$ should not be very big, and so the resultant bias will not overwhelm the other contributions to the error.} In turn, when $\vert \hat{\rho}_{i,i^\prime}\vert \leq \delta$, we ignore the fact that the distributions of $\hat{\rho}_{i,i^\prime}$ and $\sigma_{i,i'}^{(0)}$ are not the same. For the nulls, likewise, we ignore the fact that the distributions of $\hat{\rho}_{i,i^\prime}$ and $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$ are not the same.\footnote{The correlation estimates $\hat{\rho}_{i,i^\prime}$ are more disperse than the true covariances $\sigma_{i,i'}^{(0)}$, but less disperse than the statistics $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$. Thus, in our approximation, we overestimate the contribution from the non-nulls and underestimate the contribution from the nulls. Both factors make the threshold smaller than would be optimal, and thereby yield, more often than not, a more conservative standard error.} Finally, we drop the signs of the terms $\sigma_{i,i'}^{(0)}$ and $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$, letting each contribute to the loss only in absolute value.\footnote{By letting the terms $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$ contribute only in absolute value, we rule out situations where positive and negative contributions from the nulls cancel each other out, as it would be undesirable to rely on this sort of cancellation in practice.} Nevertheless, (ref) captures, in a parsimonious and tractable way, many features of the structure encoded in (ref). Including a non-null reduces the error in proportion to the population covariance $\sigma_{i,i'}^{(0)}$. Including a null increases the error in proportion to the empirical covariance $\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$. Under (ref), these objects can be proxied by the correlation $\hat{\rho}_{i,i^\prime}$. The following Theorem characterizes the threshold $\delta$ that minimizes the loss (ref), under a collection of simple conditions. \begin{theorem} Assume that the densities $f_0(\cdot)$ and $f_1(\cdot)$ are continuously differentiable. Moreover, assume that the distribution associated with $f_1(\cdot)$ is stochastically larger than the distribution associated with $f_0(\cdot)$. Define the local false discovery rate \begin{equation} \mathsf{fdr}(\delta) = P\left\{ (i,i^\prime)\in \mathcal{F} \mid \vert \hat{\rho}_{i,i^\prime}\vert = \delta \right\} , \end{equation} where $\mathcal{F}$ denotes the set of null pairs of units. The loss $L(\delta)$, defined in (ref), has a unique minimum at the threshold $\delta^\star$ that solves the equality $\mathsf{fdr}(\delta) = 0.5$. \end{theorem} We refer to the point that $\delta^\star$ solves $\text{fdr}(\delta) = 0.5$ as “the point of equalized classification.” We aim to set the threshold $\delta$ at $\delta^\star$. The intuition underlying (ref) is simple. Suppose that we choose a very large value of $\delta$, e.g., 0.9. This far out in the tail, almost all of the pairs with $\vert\hat{\rho}_{i,i^\prime}\vert \geq \delta$ will be non-nulls. Most of the non-null pairs, however, won't have correlations this large and so won't have been included in the variance estimate. Consider the decision to incrementally reduce the threshold. We are willing to do this as long as we are letting in more non-nulls than nulls. Iterating this process, we stop at the point where the proportion of new non-nulls and nulls is equalized. It remains, then, to specify a procedure for estimating $\delta^\star$. For the time being, assume that we know of the distribution function $F_0(\cdot)$ of the empirical correlations $\hat{\rho}_{i,i^\prime}$ for the null pairs of units. We treat the construction of an estimate of this quantity in the following subsection. Let $F^+_0(\cdot) = 2(1-F_0(\cdot))$ denote the associated right distribution function for the absolute correlations $\vert\hat{\rho}_{i,i^\prime}\vert$. Analogously, let $F^+(\delta)= P\{\vert \hat{\rho}_{i,i^\prime} \vert \geq \delta \} $ denote the unconditional right distribution function of the absolute correlations. Observe that \begin{flalign} P\left\{ (i,i^\prime) \in \mathcal{T} \mid \vert \hat{\rho}_{i,i^\prime} \vert \geq \delta \right\} & = (F^+(\delta) - p_0 F^+_0(\delta)) / F^+(\delta) \quad\text{and} \nonumber \\ P\left\{ (i,i^\prime) \in \mathcal{F} \mid \vert \hat{\rho}_{i,i^\prime} \vert \geq \delta \right\} &= p_0 F^+_0(\delta) / F^+(\delta) \end{flalign} give the probabilities of being non-null and null, conditional on the absolute correlation being larger than $\delta$, respectively. The difference between these probabilities is proportional to the function \begin{equation} Q(\delta) = (F^+(\delta) - p_0 F^+_0(\delta)) - p_0 F^+_0(\delta) = F^{+}(\delta) - 2 p_0 F^+_0(\delta) . \end{equation} The function $Q(\delta)$ is increasing as $\delta$ decreases if and only if more non-nulls are being let in than nulls. Thus, the point of equalized classification can be recovered by maximizing the function $Q(\delta)$. An an unbiased, but infeasible, estimate of $Q(\delta)$ is given by \begin{equation} \tilde{Q}_{n,d}(\delta) = \widehat{F}_{n,d}(\delta) - 2 p_0 F^+_0(\delta) , \end{equation} where \begin{equation} \widehat{F}_{n,d}(\delta) = \frac{2}{n (n-1)} \sum_{i = 1}^n \sum_{i^\prime < i} \text \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \geq \delta \} \end{equation} denotes the empirical right distribution function of the absolute correlations. Constructing a feasible version of the estimator (ref) requires an estimate of the proportion of nulls $p_0$. There is a large literature that studies estimators of this quantity (see e.g., langaas2005estimating, efron2007size, jin2007estimating, and cai2010optimal). For the most part, these estimators are quite complicated. Again, we take a simpler approach. In many applications of interest, the proportion of null pairs $p_0$ is close to one. In this case, the feasible estimator \begin{equation} \widehat{Q}_{n,d}(\delta) = \widehat{F}_{n,d}(\delta) - 2 F^+_0(\delta) , \end{equation} will be accurate. Our preferred approach chooses the threshold $\hat{\delta}^\star$ as the maximizer of (ref).\footnote{Setting $p_0$ to one when approximating the local false discovery rate is adopted, implicitly, by benjamini1995controlling, proposed explicitly by efron2004large, and has subsequently become commonplace efron2021computer. An infeasible threshold chosen by maximizing the function (ref) would be smaller than $\hat{\delta}^\star$, as $p_0 \leq 1$.} \subsection{Estimating the Null Distribution} To close the procedure, it remains to specify an estimator of the null distribution $F_0(\cdot)$. If the outcomes $Y_i^{(1)}, ..., Y_i^{(d)}$ were independent, this would be straightforward. In particular, in this case, if the pair $i,i^\prime$ is a null, then the correlation coefficient $\hat{\rho}_{i,i^\prime}$ will be approximately Gaussian with mean zero and variance $1/d$ lehmann2022Testing. In practice, of course, outcomes are unlikely to be independently distributed. Correlation among the outcomes reduces the effective sample size associated with the estimator $\hat{\rho}_{i,i^\prime}$, increasing its variance. To see this, we return to the collection of outcomes considered in (ref). Like (ref), Panel A of (ref) displays a histogram of the correlation estimates $\hat{\rho}_{i,i^\prime}$ across pairs of U.S.\ counties. A Gaussian density with mean zero and variance $1/d$ is overlaid in black; the dispersion in the empirical correlations is severely underestimated. \begin{figure}[p] \begin{centering} \caption{Estimating the Null Distribution and Choosing a Threshold} \begin{tabular}{c} \textit{Panel A: Correlation Between Pairs of U.S.\ Counties, Revisited} \\ \\ \textit{Panel B: Approximation to Non-Nulls Minus Nulls} \\ \end{tabular} \end{centering} \justifying { Notes: Panel A of (ref) reproduces (ref). A mean-zero, Gaussian density function with variance $1/d$ is overlaid in black. A mean-zero, Gaussian density function with variance $\hat{v}$ is overlaid in blue. Panel B displays the estimator estimator $\widehat{Q}_{n,d}(\delta)$, defined in (ref), over a grid of values of $\delta$, where we have replaced the the null distribution $F_0(\cdot)$ with its Gaussian approximation $\Phi_{\hat{v}}(\cdot)$. The solid vertical red line displays the threshold $\hat{\delta}^\star$, defined as the maximizer of $\widehat{Q}_{n,d}(\delta)$. The dashed vertical black line displays the Bonferroni threshold $\Phi^{-1}_{\hat{v}}(1-0.05/n^2)$.} \end{figure} Consequently, to estimate a reasonable null distribution, we must account for the correlation across outcomes. One way to do this would involve estimating this correlation directly. We propose a simpler approach, motivated by methods developed in efron2007size, that works by directly inspecting the empirical distribution of the correlations. In particular, we assume the center of the distribution largely consists of correlations associated with null pairs. Accordingly, we fit the variance of a mean-zero Gaussian distribution to the center of the empirical distribution of the correlation estimates. Formally, let $\Phi_v(\cdot)$ denote the cumulative distribution function of a mean-zero Gaussian random variable with variance $v$ and let $\Phi^{-1}_v(\cdot)$ and $\mathsf{IQR}_{v,0} = \Phi^{-1}_v(0.75) - \Phi^{-1}_v(0.25)$ denote the associated quantile function and interquartile range, respectively. In parallel, let \begin{equation} G_{n,d}(\delta) = \frac{2}{n (n-1)} \sum_{i = 1}^n \sum_{i^\prime < i} \text \mathbb{I}\{ \hat{\rho}_{i,i^\prime} < \delta \} \end{equation} denote the empirical distribution of the correlation estimates and let $G_{n,d}^{-1}(\cdot)$ and $\widehat{\mathsf{IQR}} = G_{n,d}^{-1}(0.75) - G_{n,d}^{-1}(0.25)$ denote the associated quantile function and interquartile range. We let $\hat{v}$ denote the value that equates the Gaussian and empirical interquartile ranges, \textit{i.e.,} $\mathsf{IQR}_{\hat{v},0} = \widehat{\mathsf{IQR}}$. We refer to the estimate $\widehat{\mathsf{df}} = 1 / {\hat{v}}$ as the effective sample size or degrees of freedom. There is an extensive statistical literature that develops more sophisticated approaches for estimating null distributions in large-scale multiple testing problems, with an emphasis on applications to analysis of DNA microarray data (see e.g., efron2007size, jin2007estimating, and cai2010optimal). efron2012large gives an influential textbook treatment. We adopt the approach outlined above, because it is computationally simple, \textit{i.e.,} it does not require any nonparametric density estimation, and performs well in our applications. In (ref), we assess the sensitivity of our results to alternative null distribution estimators. The density function associated with a mean-zero Gaussian with variance $\hat{v}$ is overlaid on Panel A of (ref) in blue. The fit at the center of the distribution is remarkably accurate. The right-tail of the Gaussian density underestimates the empirical distribution of the correlation estimates. That is, if our estimate of the null distribution is accurate, then there is evidence of a substantial number of non-null pairs of units whose residuals are significantly positively correlated. Panel B of (ref) displays the values of the estimator $\widehat{Q}_{n,d}(\delta)$, defined in (ref), over a grid of values of $\delta$, where we have replaced the null distribution $F_0(\cdot)$ with its Gaussian approximation $\Phi_{\hat{v}}(\cdot)$. A vertical red line denotes the maximizer of this curve, i.e., $\hat{\delta}^{\star}$, at roughly $0.4$. It is worth highlighting that the threshold $\hat{\delta}^{\star}$ is much smaller than the threshold that would be chosen with a standard correction for multiple hypothesis testing. The threshold that would be used in an application of a Bonferroni correction, i.e., $\Phi^{-1}_{\hat{v}}(1-0.05/n^2)\approx 0.73$, is displayed with a vertical black line. We conclude by noting that, in small samples, correlation coefficients have a highly skewed distribution away from the null. In other words, correlation coefficients are not pivotal---even in Gaussian data their variances depend on the true underlying correlation. fisher1915frequency shows that correlation coefficients can be made approximately Gaussian through the transformation \begin{equation} \tilde{\rho}_{i,i^\prime} = \frac{1}{2} \log\left( \frac{1 + \hat{\rho}_{i,i^\prime}}{1 - \hat{\rho}_{i,i^\prime}} \right) . \end{equation} See hotelling1953new and efron1982transformation for further discussion.\footnote{muralidharan2010detecting shows that the Fisher transformation maintains it normalizing property in data satisfying (ref).} In our applications, we find that estimating a null distribution by using a Gaussian approximation to the Fisher transformed correlations $\tilde{\rho}_{i,i^\prime}$, rather than the empirical correlations $\hat{\rho}_{i,i^\prime}$, is more robust. In practice, we recommend using the Fisher transformed correlation in place of the empirical correlation when implementing the construction proposed in this section. We have centered our discussion around correlation coefficients to ease exposition. \subsection{Performance} We now return to the calibrated simulation considered in (ref). We show that the TMO estimator, implemented with the threshold $\hat{\delta}^\star$, exhibits smaller bias and has more accurate rejection rates than standard methods for constructing spatial standard errors. In particular, the fifth row of (ref) displays estimates of the bias and rejection rate of the TMO estimator when used to construct standard errors for the two regressions considered in (ref). The biases are substantially smaller, and the associated rejection rates are closer to the nominal Type I error rate, \textit{i.e.,} $\alpha = 0.05$. It is worth acknowledging, however, that this simulation environment is particularly suited to the TMO estimator, as the underlying spatial correlation structure is determined by the same set of outcomes used by TMO. Moreover, the spatial correlation structure used in the simulation is quite sparse. More diffuse, geographically structured correlation patterns may be better suited for the alternative methods. Despite this, the good performance of the TMO estimator in this environment can be interpreted as demonstrating that it is able to adjust for forms of spatial correlation that are not captured by existing methods. There may be some concern that these results are driven by the two choices of treatments considered in (ref), \textit{i.e.,} the change in the percentage of the population in county $i$ that has completed a college degree and the change in the per-acre value of farm-land. To address this concern, we replicate the Monte Carlo experiment, setting each of the 91 outcome variables discussed in (ref) as the treatment variable. (ref) displays the distribution functions of the bias and rejection rate associated with each method across these specifications. The TMO standard errors have uniformly better bias, and more accurate rejection rates, than the HC1 standard errors, robust standard errors clustered at the state level, and conley1999gmm standard errors. In Panel A, we display two forms of the SCPC standard errors. These differ according to whether or not the alternative critical value associated with this procedure is incorporated into the standard error.\footnote{Unless otherwise specified, when reporting results associated with SCPC, we incorporate the critical value into the estimate of the standard error.} The SCPC standard errors, once adjusted in this way, more closely approximate the TMO correction, though, unlike the TMO correction, they under-reject in about 20% of the specifications. \begin{figure}[t] \begin{centering} \caption{Distributional Comparison of Performance Across Treatments} \begin{tabular}{cc} \textit{Panel A: Bias in standard errors} & \textit{Panel B: Rejection Rate}\tabularnewline & \tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) displays CDFs of the bias and rejection rates associated with various methods across several specifications of the calibrated simulation experiment outlined in (ref). In particular, the distributions are taken across the specifications where each of the 91 outcome variables discussed in (ref) is set as the treatment variable. For each specification, we take 1,000 simulation draws. In Panels A and B, the vertical red lines display the hypothetical performance of an unbiased estimator and an unbiased test, respectively. In Panel A, we display two forms of the SCPC standard errors. The label “SCPC (Adj.)” denotes that the alternative critical value associated with the SCPC procedure has been incorporated into the standard error.} \hrulefill \end{figure} \section{Consistency and Gaussian Approximation across Correlated Outcomes} The main non-standard feature of our setting is that the data under consideration are correlated across both locations and outcomes. In this section, we consider how both sources of correlation affect the methods developed in (ref). We have two objectives. First, we give sufficient conditions under which the spatial correlation across units is consistently estimable. In particular, we give a quantitative description of the extent of the correlation across outcomes and units allowed by our estimator. Second, we give sufficient conditions for the validity of a Gaussian approximation to the empirical distribution of residual correlations of null pairs of units. Here, we highlight how the extent of the correlation across outcomes impacts the dispersion of the null distribution and the quality of the Gaussian approximation. \subsection{Sparsity and Consistency} We begin by imposing a series of simplifications that will greatly ease the notational burden, while preserving the key features of our setting. First, we focus our attention on inference for the mean $\tau^{(0)} = \mathbb{E} [ Y_i^{(0)} ]$. That is, we set $W_i = 1$ for all units. Second, we impose the normalization $\mathbb{E} [ Y_i^{(j)} ] = 0$ for all units $i$ and outcomes $j$ and consider the simplified estimators \begin{align} \hat{\lambda}_{i,i^\prime} &= \frac{1}{d} \sum_{j=1}^d \tilde{Y}^{(j)}_i \tilde{Y}^{(j)}_{i^\prime} ,\quad\text{where}\quad \tilde{Y}^{(j)}_i = \hat{\gamma}^{-1/2}_j Y^{(j)}_{i} \quad\text{and}\quad \hat{\gamma}_j = \frac{1}{n} \sum_{i=1}^n (Y_{i}^{(j)})^2 . \end{align} That is, relative to the feasible estimators (ref) and (ref), the infeasible estimators (ref) have been centered with the true values of the unknown parameters $\mathbb{E} [ Y_i^{(j)} ]$. The variance estimator proposed in (ref) is premised on the assumption that there is not too much dependence across either locations or outcomes. Here, we describe this assumption formally. That is, we give sufficient conditions on the scope of the dependence across units and outcomes under which we can recover the spatial correlation across units. Our results are expressed in terms of the matrices $\Lambda_n$ and $\Gamma_d$ introduced in (ref). Recall that the matrix $\Lambda_n$ parametrizes the structure of the dependence across units. We assume that the rows of this matrix are sparse, in the following sense. \begin{assumption}[Sparsity] There exists a constant $0\leq q<1$ and a sequence $\kappa_{n}$ such that the components of $\Lambda_{n}$ satisfy \begin{equation} \max_{i\in\{1,\ldots,n\}}\sum_{i^{\prime}=1}^{n}\vert\lambda_{i,i^{\prime}}\vert^{q}\leq\kappa_{n} \end{equation} for each integer $n$, where we adopt the convention that $0^0=0$. \end{assumption} In the extreme case that $q=0$, (ref) means that at most $\kappa_n$ elements of any row of $\Lambda_{n}$ are non-zero. Non-zero values of $q$ accommodate less stringent forms of sparsity.\footnote{A version of (ref) was first considered in the normal-means problem by abramovich2006special. It has been applied to estimation of covariance matrices in independent data by bickel2008regularized and cai2011adaptive, among others.} In turn, the matrix $\Gamma_d$ parametrizes the structure of the dependence across outcomes. We place no \textit{a priori} restrictions on this quantity. In the statement of our results we will reference the matrix \begin{equation} \Omega_d = \Gamma_d^{1/2}\Xi_d \Gamma_d^{1/2} ,\quad\text{where}\quad \Xi_d = \mathsf{diag}(1/\gamma_{1,1},\ldots, 1/\gamma_{d,d}) \end{equation} and $A^{1/2}$ denotes the square root of the matrix $A$. Observe that $\Omega_d$ can be interpreted as the correlation matrix across outcomes. Let $\omega_{j,j^\prime}$ denote the $j,j^\prime$th element of $\Omega_d$. The following Theorem characterizes the rate of convergence of the estimator $\hat{\gamma}_0\hat{\lambda}_{i,i^\prime}$, where we normalize the covariance estimator $\hat{\lambda}_{i,i^\prime}$ by the variance estimate $\hat{\gamma}_0$ to fix the scale. To ease exposition, throughout the main text, we allow the constants $c$ and $C$ to depend on the largest and smallest diagonal terms of the matrices $\Lambda_n$ and $\Gamma_d$ as well as the maximum sub-Gaussian norm of the data $Y^{(j)}_i$. We track the dependence of the rate of convergence on these quantities in the proof. \begin{theorem}[Consistency] Suppose that a strengthened version of (ref), stated in (ref), and (ref) hold. Fix a constant $0<\phi<1$ and define the sequence \begin{equation} \varphi_{n,d} = \left( \frac{1}{d} \sqrt{\sum_{j=1}^d\sum_{j^\prime=1}^d\omega_{j,j^\prime}^2} + \sqrt{\frac{\kappa_n}{n}} \right) \log^{3/2}(nd / \phi) . \end{equation} There exist constants $c < 1$ and $C$ such that if $\varphi_{n,d}<c$, then \begin{equation} \max_{i,i^\prime \in \{1,\ldots,n\}} \vert \hat{\gamma}_0 \hat{\lambda}_{i,i^\prime} - \lambda_{i,i^\prime} \vert \leq C \varphi_{n,d} \end{equation} with probability greater than $1-C\phi$. \end{theorem} (ref) bounds the error in our estimate of cross-sectional dependence in terms of two quantities: \begin{equation} \frac{1}{d}\sqrt{\sum_{j=1}^d\sum_{j^\prime=1}^d\omega_{j,j^\prime}^2} \quad\text{and}\quad \sqrt{\frac{\kappa_n}{n}} . \end{equation} The former parametrizes the extent of the correlation across outcomes. Observe that if the matrix $\Omega_d$ is diagonal, then this term is equal to $d^{-1/2}$, \textit{i.e.,} the parametric rate. As the correlation across outcomes increases, this term becomes larger, and the precision of our estimate decreases. Consequently, when constructing collections of outcomes, researchers face an essential trade-off. Outcomes that more credibly reflect the same underlying correlation structure might tend to be highly correlated. If outcomes are too similar, the resultant estimate of cross-sectional dependence will be imprecise. If outcomes are too different, they might not capture the same cross-sectional correlation, and estimates of cross-sectional dependence will be biased. We test approaches for navigating this trade-off in the applications considered in (ref). We summarize our recommendations in (ref). The second term in (ref) describes the sparsity of the matrix $\Lambda_n$, and arises due to the error in variance estimator $\hat{\gamma}_j$. The interpretation of this quantity is cleanest when (ref) holds with $q=0$. In this case, our estimate of cross-sectional dependence is accurate so long as the maximal number of non-zero elements of any row of $\Gamma_n$ is small relative to the total number of units $n$. As a result, there is less scope to apply our procedure in settings where the number of units is small, e.g., in cases where the units of observation are U.S.\ states. This situation echos results in the literature on estimation of clustered standard errors. For example, hansen2019asymptotic show that standard estimators of clustered standard errors are consistent only if the number of units in each cluster is small relative to the sample size. \subsection{Gaussian Approximation} In (ref), we apply a Gaussian approximation to the empirical distribution of the correlations of the residuals of the null pairs of units. Approximations of this sort are widely applied in independent data (see e.g., efron2007size for discussion). However, in our context, there are two forms of dependence that complicate this program. First, outcomes are likely to be highly correlated with each other. Second, the correlations $\hat{\rho}_{i,i^\prime}$ and $\hat{\rho}_{i,i^{\prime\prime}}$ will be dependent. In this section, we characterize how dependence across outcomes, and across correlations taken with the same unit, impact a Gaussian approximation to the null distribution.\footnote{efron2010correlated and azriel2015empirical give less specialized methods for approximating the empirical distribution of correlated observations.} In particular, continuing with the simplifications introduced in the previous subsection, we give a bound on the difference between the cumulative distribution \begin{equation} G_{n,d}(\delta) = \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \mathbb{I}\{\hat{\gamma}_0 \hat{\lambda}_{i,i^\prime} \leq \delta\} . \end{equation} and an appropriately scaled Gaussian distribution function. For the sake of simplicity, we assume that all pairs of units are nulls. Further, as the distribution function (ref) is taken over covariances, rather than correlations, we assume that the data are homoskedastic. That is, we assume that the matrix $\Lambda_n$ is given by $\lambda I_n$ for some constant $\lambda$, where $I_n$ is the $n\times n$ identity matrix. Moreover, we assume that the underlying data $Y_i^{(j)}$ are Gaussian. Even in this highly stylized case, the calculations involved are non-standard. \begin{theorem}[Gaussian Approximation] Suppose that (ref) holds, $\Lambda_n = \lambda I_n$ for some constant $\lambda$, and that the data $Y_i^{(j)}$ are normally distributed. Let $\eta_1,\ldots,\eta_d$ denote the eigenvalues of $\Omega_d$, defined in (ref). Fix the parameter \begin{equation} v^\star = \frac{\lambda^2}{d^2} \sum_{j=1}^d \sum_{j^\prime=1}^d \omega_{j,j^\prime}^2 \end{equation} and recall that $\Phi_v(\cdot)$ denotes the cumulative distribution function of a mean-zero Gaussian random variable with variance $v$. Fix a constant $0<\phi<1$. For each threshold $\delta$ in $\mathbb{R}$, it holds that \begin{equation} \vert G_{n,d}(\delta) - \Phi_{v^\star}(\delta) \vert \leq C\left(\frac{\sum_{j=1}^d \eta^3_j}{(\sum_{j=1}^d \eta_j^2)^{3/2}} + \sqrt{\frac{\log^{3}(dn/\phi)}{n}}\right) \end{equation} with probability greater than $1-\phi$. \end{theorem} (ref) establishes that the empirical distribution of the null covariances fluctuates around a Gaussian distribution whose variance $v^\star$ is given by (ref). Observe that the square-root of $v^\star$ appears in the bound (ref). Written differently, given a collection of outcomes with correlation matrix $\Omega_d$, the true “effective sample-size” or “degrees of freedom” is given by \begin{equation} \mathsf{df}^\star = \frac{1}{v^{\star}} = \frac{1}{\lambda^2} \frac{d^2}{\sum_{j=1}^d \sum_{j^\prime=1}^d \omega_{j,j^\prime}^2} . \end{equation} If the outcomes are mutually independent, then $\mathsf{df}^\star = d/\lambda^2$. As the correlation across outcomes increases, the effective sample size decreases, and the dispersion of the null distribution increases. Correlation across outcomes also reduces the quality of the Gaussian approximation. Consider the term \begin{equation} \frac{\sum_{j=1}^d \eta^3_j}{(\sum_{j=1}^d \eta_j^2)^{3/2}} \end{equation} appearing in the bound (ref). If the outcomes are mutually independent, then all of the eigenvalues of the correlation matrix $\Omega_d$ are equal to one, and the error (ref) is equal to $d^{-1/2}$. Indeed, this is the order of the error that we should expect in independent data. As the correlation across outcomes increases, this term becomes larger. In the extreme, when $\Omega_d$ is rank one, having a single eigenvalue equal to $d$, the error (ref) is equal to one. In sum, correlation across outcomes impacts both our ability to discriminate nulls from non-nulls, through the dispersion of the null distribution, and our ability to accurately recover the null distribution through a Gaussian approximation. \section{Applications} We now examine the impact of applying the TMO variance estimator, as well as alternative spatial standard error estimators, to a sample of recently published papers in applied economics. We find that the proposed estimator can make a substantial difference in practice. \subsection{Spatial Correlation in Applied Economics} We start off by documenting the prevalence of the potential for spatial correlation in a set of recently published papers. We review all 370 papers published in 2023 in the \textit{American Economic Review}, \textit{Econometrica}, the \textit{Journal of Political Economy}, the \textit{Quarterly Journal of Economics}, and the \textit{Review of Economic Studies}. We hand-code each paper as belonging to four (overlapping) categories: Econometrics (26 papers), Economic Theory (108 papers), Macroeconomics (66 papers), and Empirical work (197 papers). We take an expansive view of the latter category, including macroeconomic, econometric, or theoretical papers with empirical applications. Within the 197 papers with empirical content, there are 17 laboratory experiments, 28 field experiments (RCTs), and 23 descriptive papers (e.g., measures of inequality or inter-generational mobility). Empirical studies in these three categories do not typically give rise to the issue of spatial correlation in the errors, due to the experimental design or to the nature of the issue being studied. \begin{table}[t] \caption{Survey of Papers Published in Leading Economics Journals} \setlength\extrarowheight{-3pt} \begin{tabular}{l c c c c c c} \toprule & All & \textit{AER} & \textit{ECTA} & \textit{JPE} & \textit{QJE} & \textit{RES} \\ \hline All papers ($N$) & 370 & 93 & 66 & 71 & 48 & 92 \\ \multicolumn{7}{l}{\textbf{Fields (not mutually exclusive)}} \\ Economic Theory & 108 & 18 & 23 & 34 & 6 & 27 \\ Econometrics & 26 & 4 & 12 & 2 & 2 & 6 \\ Macroeconomics & 66 & 15 & 12 & 15 & 7 & 17 \\ Empirical & 197 & 67 & 20 & 30 & 33 & 47 \\ \multicolumn{7}{l}{\textbf{Method}} \\ Descriptive & 23 & 5 & 1 & 3 & 4 & 10 \\ Randomized control trial & 28 & 11 & 4 & 5 & 4 & 4 \\ Lab experiment & 17 & 8 & 2 & 2 & 3 & 2 \\ Observational & 128 & 41 & 13 & 19 & 22 & 33 \\ \multicolumn{7}{l}{\textbf{Potential for spatial correlation}} \\ No & 67 & 23 & 10 & 12 & 10 & 12 \\ Yes & 61 & 18 & 3 & 7 & 12 & 21 \\ U.S.\ unit of observation & 31 & 9 & 2 & 2 & 7 & 11 \\ \textit{Most common:} County & 15 & 4 & 1 & 1 & 6 & 3 \\ \textit{Second most common:} State & 6 & 1 & 1 & 1 & 1 & 2 \\ Non-U.S.\ unit of observation & 30 & 9 & 1 & 5 & 5 & 10 \\ \textit{Most common:} Country & 6 & 2 & 1 & 2 & 1 & 0 \\ \textit{Second most common:} Chinese county & 2 & 1 & 0 & 0 & 1 & 0 \\ \bottomrule \end{tabular} \justifying { Notes: (ref) gives the results of a survey of papers published in the \textit{American Economic Review}, \textit{Econometrica}, \textit{Journal of Political Economy}, \textit{Review of Economic Studies}, and \textit{Quarterly Journal of Economics} in 2023. Each paper is categorized by field. Empirical papers are categorized according to their primary methodology. Observational papers are categorized according to whether their primary analyses have the potential to be affected by spatial correlation. For these papers, the geographic level of such correlation is coded.} \hrulefill \end{table} We classify the remaining 128 papers, a third of the original sample, as belonging to an “observational" category, which includes all (non-randomized) studies of how variation in a variable $W_i$ affects an outcome variable $Y_i$. This category includes difference-in-difference designs, structural papers, and shift-share instruments, among other designs. Within these 128 papers, we code whether a main specification of the paper can be affected by spatial correlation of the observations. If so, we record the geographical level of such correlation. In our assessment, spatial correlation of the errors plays a role in a striking number of cases: 61 cases out of 128, or 48 percent, of observational papers. These spatial designs are present in all five journals and especially common in the \textit{Review of Economic Studies} (21 papers), the \textit{American Economic Review} (18 papers), and the \textit{Quarterly Journal of Economics} (12 papers). Our categorization is likely a lower bound of the importance for spatial correlation because in some papers the analysis has a spatial aspect, but it is not emphasized, e.g., the study of establishment-level productivity in cases in which the location of the establishment is not recorded. The most common geographic unit at which the correlation is likely to occur is the U.S.\ county, with 15 papers in this category. For papers examining outcomes in the United States (31 papers overall), the next most common geographic units are states (6 papers), census blocks (2 papers), zip codes (2 papers), and commuting zones (2 papers). For papers with coverage outside the United States (30 papers), the most common levels of geography are the country (6 papers), Chinese counties (2 papers), followed by a variety of other local data units (e.g., French municipalities, Austrian regions, and Ukrainian districts). Given that the most common spatial-based setting in this sample is the county-level spatial design, accounting for 8 percent of all empirical papers, in the next subsection we focus on this category of papers, systematically revisiting the estimates in the 15 county-level papers. \subsection{Results} For the 15 papers that use U.S.\ county-level observations, we query replication packages and, where necessary, authors to access underlying data. For eight out of the 15 papers, we are not able to reproduce a main specification, due to a key variable being confidential. In the remaining seven cases, we identify and reproduce a main empirical specification. Panel A of (ref) displays the point estimate, original standard error, and level of clustering used for each of these papers. For the majority, there is no clustering at a higher level of geographic aggregation. (Online appendices often include additional corrections, typically a conley1999gmm standard error.) The papers differ in a number of ways, ranging from cross-sections to panel data (as indicated by the number of periods in Column 7), from modern to historical outcome data, and span the most common observational designs. For each of the papers, to apply the TMO correction we construct a collection of auxiliary outcome variables, typically from replication packages associated with the papers themselves. In several cases, these outcomes are supplemented with additional variables obtained from external sources. Column 8 reports the number of outcomes obtained for each paper. Further details on how we handle the data associated with each paper, and construct a relevant set of auxiliary outcomes, are given in (ref). Column 1 presents the key finding: the ratio of the TMO standard error (Column 2) to the original standard error (Column 3). For each paper, we use the TMO procedure to augment the original level of clustering (Column 4). That is, a county pair $(i,i^\prime)$ is never thresholded if counties $i$ and $i'$ are elements of the same cluster. Written differently, the TMO procedure augments the original standard error by allowing for dependence between pairs of counties, from different clusters, whose absolute correlation $\vert \hat{\rho}_{i,i^\prime} \vert $ is above the threshold $\hat{\delta}^{\star}$ (Column 10). Column 11 reports the proportion of location pairs that meet this criteria and contribute to the TMO standard error. The proportion varies from 0.18% to 4.06%. TMO leaves the standard error unaltered in one case, and increases the remaining six. The median increase is 37% percent. To illustrate that TMO is applicable to geographic units beyond U.S.\ counties, we consider two additional, prominent, non-county papers: chetty2014land, which consider U.S.\ commuting zone level data, and acemoglu2019democracy, which studies country level data.\footnote{muller2024spatial also consider the data from chetty2014land.} The findings are broadly consistent, although the magnitude of the corrections are both on the lower end of the U.S.\ county-level corrections. As a concrete example, bernini2023race studies the impact of the Voting Rights Act on the racial makeup of local governments in the Southern U.S. We obtain 60 outcome variables from the replication package, which we standardize and use to compute the correlation across counties. Given the correlation across outcomes, find that the null distribution has approximately 26 degrees of freedom (as reported in Column 9), stressing the importance of starting with a large number of outcomes. For this paper, the optimal threshold is set at 0.54, keeping only 0.7% of county-pairs (as reported in Columns 10 and 11). The TMO method increases the standard error by 37% relative to the estimate reported in the paper. \newgeometry{left=0.5in, right=0.5in, top=1in, bottom=1in} \begin{landscape} \begin{table}[t] \caption{Ratio of TMO to Original Standard Errors in Nine Recent Papers} { \begin{tabular} {c c c c c c c c c c c} \toprule SE Ratio (TMO/Orig.) & TMO & Original SE & Clustering & Coefficient & Units $ n$ & Periods & Outcomes $ d$ & $ \widehat{\mathsf{df}}$ & $ \hat{\delta}^\star$ & $ \%\ge\hat{\delta}^\star$ \\ (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) \\ \hline \multicolumn{11}{l}{\textbf{Panel A. County-level papers from 2023 survey}} \\ \multicolumn{11}{l}{\textit{Cook et al. (2023): Effect of White casualties in WWII on Green Book establishments (Table 3 Column 4)}} \\ \textcolor{blue}{\textbf{3.50}} & 0.03 & 0.01 & County & 0.06 & 3104 & 12 & 78 & 36.7 & 0.42 & 3.57 \\ \multicolumn{11}{l}{\textit{Caprettini and Voth (2023): Effect of New Deal grants on purchase of WWII bonds (Table 2 Column 1)}} \\ \textcolor{blue}{\textbf{2.10}} & 0.05 & 0.02 & County & 0.19 & 3022 & 1 & 44 & 26.1 & 0.51 & 1.25 \\ \multicolumn{11}{l}{\textit{Esposito et al. (2023): Effect of The Birth of a Nation screening on patriotic vs. divisive language in newspapers (Table 3 Column 3)}} \\ \textcolor{blue}{\textbf{1.76}} & 0.13 & 0.07 & County & 1.09 & 786 & 132 & 172 & 192.3 & 0.19 & 4.06 \\ \multicolumn{11}{l}{\textit{Bernini et al. (2023): Effect of Black % $\times$ Voting Rights Act on share of black elected officials (Table 2 Column 4)}} \\ \textcolor{blue}{\textbf{1.37}} & 0.06 & 0.04 & Judicial div. & 0.10 & 971 & 1 & 60 & 25.8 & 0.54 & 0.70 \\ \multicolumn{11}{l}{\textit{Bazzi et al. (2023): Effect of Southern Whites % on 2016 Trump vote-share (Table 2 Column 4)}} \\ \textcolor{blue}{\textbf{1.20}} & 0.20 & 0.17 & $60\text{mi}^2$ & 1.03 & 1886 & 1 & 60 & 21.7 & 0.54 & 1.99 \\ \multicolumn{11}{l}{\textit{Calderon et al. (2023): Effect of Black population on Democratic vote-share (Table 2 Column 6)}} \\ \textcolor{blue}{\textbf{1.09}} & 0.49 & 0.45 & County & 1.88 & 1263 & 3 & 78 & 21.4 & 0.55 & 1.59 \\ \multicolumn{11}{l}{\textit{Moscona and Sastry (2023): Effect of extreme temperature exposure $\times$ innovation on price per agricultural land acre (Table 3 Column 1)}} \\ \textcolor{blue}{\textbf{1.00}} & 0.09 & 0.09 & State & 0.25 & 3000 & 2 & 39 & 17.1 & 0.67 & 0.18 \\ \hline \multicolumn{11}{l}{\textbf{Panel B. Prominent non-county examples}} \\ \multicolumn{11}{l}{\textit{Chetty et al. (2014): Correlation between share of African American and upward mobility (Figure 8 Row 1)}} \\ \textcolor{blue}{\textbf{1.11}} & 0.05 & 0.05 & State & -0.36 & 693 & 1 & 158 & 14.8 & 0.68 & 0.35 \\ \multicolumn{11}{l}{\textit{Acemoglu et al. (2019): Effect of democracy on log GDP per capita (Table 2 Column 3)}} \\ \textcolor{blue}{\textbf{1.06}} & 0.24 & 0.23 & Country & 0.79 & 175 & 47 & 79 & 79.1 & 0.27 & 5.97 \\ \bottomrule \end{tabular} } \end{table} \justifying \begin{spacing}{0.9} { Notes: Panel A of (ref) shows the results from applying TMO to seven recent county-level papers surveyed in (ref). Panel B shows the results for two prominent non-county examples, chetty2014land at the commuting-zone level and acemoglu2019democracy at the country level. Column 1 reports the ratio of the TMO standard error estimate (Column 2) relative to the original standard error in each paper (Column 3). The TMO standard error in Column 2 is combined with the original level of clustering in the paper (Column 4). The table also shows the coefficient on the main regressor of interest (Column 5), the number of locations in each sample (Column 6), the number of time periods (Column 7), the number of outcomes used in the TMO procedure (Column 8), the estimated degrees of freedom for the null distribution (Column 9), the optimal threshold $\hat{\delta}^{\star}$ that maximizes (ref) (Column 10), and the percent of location pairs $(i,i^\prime)$ that have a correlation $\hat{\rho}_{i,i^\prime}$ greater than $\hat{\delta}^{\star}$ in absolute value, out of all location pairs $(i,i^\prime)$ where $i$ and $i^\prime$ are in different clusters (Column 11).} \end{spacing} \end{landscape} \restoregeometry (ref) provides further details concerning the application to bernini2023race. Panel A shows the distribution of the correlation estimates $\hat{\rho}_{i,i^\prime}$. The solid curve corresponds to the estimate of the null density, proposed in (ref). Panel B plots the ratio of the TMO standard error and the original standard error, as the threshold $\delta$ varies. At low thresholds ($<0.3$), the standard error ratio hovers around 1, indicating the high prevalence of noise, relative to signal, added to the standard error. In the neighborhood around the optimal threshold $\hat{\delta}^\star = 0.54$, the TMO estimate is relatively stable. As the threshold increases, the standard error reverts to the the value reported in the paper. \begin{figure}[p] \begin{centering} \caption{Details for Application to bernini2023race } \begin{tabular}{c} \textit{Panel A: Correlation Between Pairs of Counties} \\ \\ \textit{Panel B: TMO across Thresholds}\\ \tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) gives details concerning the application of TMO to bernini2023race. In Panel A, we plot the distribution of the pairwise correlation estimates, $\tilde{\rho}_{i,i^\prime}$, constructed using 60 auxiliary outcomes obtained from the paper's replication package. The solid curve denotes the null density estimate, obtained with the procedure outlined in (ref). Panel B shows the ratio of the TMO standard error to the original standard error as the threshold $\delta$ varies. In both panels, the vertical line marks the estimate of the optimal threshold $\hat{\delta}^{\star}$.} \end{figure} \subsection{Comparison} How does TMO compare, in practice, to other popular corrections for spatial correlation? Column 1 of (ref) displays the ratio of the TMO standard error to the original standard error. Here, unlike Column 1 of (ref), TMO does not augment the original level of clustering. That is, TMO is permitted to threshold pairs of counties in the same cluster. Column 2 gives the analogous ratio, where now TMO augments state level clusters, i.e., pairs of counties in the same state are never thresholded. In most cases, the standard error increases moderately, although the increase tends not to be dramatic. Combining TMO with state level clustering, in this way, is well suited for settings where policies are adopted at the state level.\footnote{Clustering by state can be problematic in samples with a small number of states. For example, the data used in bernini2023race cover only 12 states. There, 12% of county pairs are in the same state.} Columns 3 and 5 display analogous results for two leading alternatives: conley1999gmm standard errors, constructed with a 150 mile bandwidth, and muller2022spatial,muller2023spatial SCPC standard errors. For each procedure, we report the ratio of the adjusted standard error to the original standard error. In several cases, the conley1999gmm and SCPC are similar to the TMO. In others, they are qualitatively different. Columns 4 and 6 display analogous ratios, where now TMO has been applied to augment the conley1999gmm and SCPC, respectively. Combining the two forms of correction tends to lead to further increases in the resultant standard error. This is intuitive---distance-based methods are better suited to capture more diffuse, geographically driven, dependencies. The TMO standard errors, by contrast, will better capture strong dependencies that are not driven by geographic proximity. The applications that we have considered thus far exhibit spatial correlation predicted well by both geographic and non-geographic factors. In an exercise similar to (ref), (ref) displays, for each of the U.S.\ county or commuting zone level papers, the proportion of highly correlated pairs of locations that are similar along observable characteristics. Pairs of locations with strong positive correlations tend to be close (within the top 10%) in at least one of these observable characteristics. While geographic proximity is typically the first or second most predictive dimension, there are still a significant proportions of location pairs that are highly correlated, but not necessarily close in geographic distance. \begin{table}[t] \caption{Comparison with Alternative Spatial Standard Errors } \resizebox{\textwidth}{!}{ \begin{tabular} {l c c c c c c} \toprule & \multicolumn{2}{c}{TMO} & \multicolumn{2}{c}{Conley 150mi} & \multicolumn{2}{c}{SCPC} \\ \cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} & & +State Clusters & & +TMO & & +TMO \\ & (1) & (2) & (3) & (4) & (5) & (6) \\ \hline Cook et al. (2023) & 3.50 [3.6%] & 3.57 [6.2%] & 0.97 [3.1%] & 3.40 [7.5%] & 1.37 & 4.09 \\ Caprettini and Voth (2023) & 2.10 [1.3%] & 2.37 [4.1%] & 2.16 [3.2%] & 2.48 [4.2%] & 1.81 & 3.01 \\ Esposito et al. (2023) & 1.76 [4.1%] & 1.88 [7.4%] & 1.79 [3.9%] & 2.16 [7.6%] & 0.86 & 2.29 \\ Bernini et al. (2023) & 1.32 [0.7%] & 1.72 [12.3%] & 1.40 [9.0%] & 1.51 [9.6%] & 2.02 & 2.26 \\ Bazzi et al. (2023) & 1.15 [2.1%] & 1.43 [5.4%] & 1.47 [4.5%] & 1.46 [6.0%] & 1.12 & 1.58 \\ Calderon et al. (2023) & 1.09 [1.6%] & 0.45 [5.5%] & 0.95 [5.2%] & 0.96 [6.2%] & 2.79\textsuperscript{\textdagger} & 2.40\textsuperscript{\textdagger} \\ Moscona and Sastry (2023) & 0.64 [0.3%] & 1.00 [3.3%] & 0.96 [3.2%] & 0.95 [3.4%] & 0.95 & 1.08 \\ Chetty et al. (2014) & 0.77 [0.4%] & 1.11 [3.2%] & 0.97 [2.6%] & 1.07 [3.0%] & 2.70 & 2.71 \\ Acemoglu et al. (2019) & 1.06 [6.0%] & - & 1.03 [2.9%] & 1.07 [8.1%] & 1.01 & 1.10 \\ \bottomrule \end{tabular} } \justifying \begin{spacing}{0.9} { Notes: (ref) compares the behavior of several methods for adjusting standard errors for spatial correlation. Each entry gives the ratio of the standard error produced by the method denoted in each column to the original standard error. Wherever applicable, the proportion of location pairs that are allowed to correlate is displayed in brackets. Column 1 reports the TMO estimates. Here, TMO does not augment the original level of clustering in the paper (which leads to differences from Column 1 in (ref)). Column 2 gives the analogous ratio when TMO augments U.S.\ state-level clusters (which does not apply to the country-level acemoglu2019democracy paper). Column 3 reports the adjustment associated with conley1999gmm standard errors using a 150-mile bandwidth. Column 4 combines TMO with the Conley correction. For the country-level acemoglu2019democracy paper, the bandwidth is set to 650 miles, which covers 3% of all country-pairs. Column 5 displays the results from the muller2022spatial,muller2023spatial SCPC method. Column 6 integrates TMO with the SCPC estimates. Columns 5 and 6 do not report the proportion of location pairs that are allowed to correlate, as SCPC does not directly “zero out” terms as in the (ref) formula. \textsuperscript{\textdagger}calderon2023racial conduct a weighted regression, which is not supported by the current SCPC Stata package; these estimates show the adjustment for the unweighted specification. } \end{spacing} \hrulefill \end{table} \section{Recommendations for Practice} Economic outcomes are often linked across locations. Statistical inferences that do not account for these linkages can exhibit substantial biases. This paper introduces a method for adjusting standard errors for spatial correlation. We call this method “Thresholding Multiple Outcomes” (TMO). The essential input is a collection of auxiliary outcome variables. The main assumption is that the spatial correlation in the residuals for a regression problem of interest is shared by the analogous residuals constructed using these auxiliary outcomes. Under this assumption, the auxiliary outcomes can be used to estimate spatial correlation. We propose to use the empirical distribution of these estimates to determine pairs of locations that are very correlated. Standard errors in the original regression problem can then be adjusted by allowing for correlation among these pairs. The method is summarized in (ref). \begin{algorithm}[t] \caption{Thresholding Multiple Outcomes (TMO)} \setstretch{1.35} \SetAlgoHangIndent{0pt} \KwIn{Outcome $Y_i^{(0)}$ and treatment $W_i$ for units $i$ in $1,...,n$} \KwResult{Estimate of the variance of $\hat{\tau}^{(0)}$, the least-squares regression coefficient of $Y_i^{(0)}$ on $W_i$} Identify a collection of auxiliary outcomes $Y_i^{(1)},\ldots,Y_i^{(d)}$ For each unit $i$ and outcome $j$, compute the normalized residuals $\tilde\varepsilon_{i}^{(j)}$, defined in (ref) For each pair of units $i\neq i^\prime$, measure the correlation $\hat{\rho}_{i,i^\prime}$, defined in (ref) Approximate the null distribution of $\hat{\rho}_{i,i^\prime}$ with $\mathsf{N}(0,\hat{v})$, where $\hat{v}$ solves \begin{equation*} \Phi_v^{-1}(0.75) - \Phi_v^{-1}(0.25) = G_{n,d}^{-1}(0.75) - G_{n,d}^{-1}(0.25) , \end{equation*} the function $\Phi_v(\cdot)$ is the CDF of $\mathsf{N}(0,v)$, and $G_{n,d}(\cdot)$ is defined in (ref) Choose the threshold $\hat{\delta}^\star$ by minimizing \begin{equation*} \hat{Q}_{n,d}(\delta) = (\hat{F}_{n,d}(\delta) - F_0^{+}(\delta)) - F_0^{+}(\delta) \end{equation*} over $\delta>0$, where $F_0^{+}(\cdot) = 2(1- \Phi_{\hat{v}}(\cdot))$ and $\hat{F}_{n,d}(\cdot)$ is defined in (ref) Collect the estimates \begin{equation*} \hat{\sigma}^{(0)}_{i,i^\prime}(\hat{\delta}^\star) = \begin{cases} \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_i ,& i = i^\prime, \\ \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}\text \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \geq \hat{\delta}^\star \} , & i\neq i^\prime , \end{cases} \end{equation*} into the matrix $\widehat{\Sigma}^{(0)}(\hat{\delta}^\star)$ \nonl \textbf{Return} the variance estimate \begin{equation*} \widehat{V}(\hat{\delta}^\star) = \left(W^\top W\right)^{-1} \left( W^\top \widehat{\Sigma}^{(0)}(\hat{\delta}^\star) W\right)\left(W^\top W\right)^{-1} \end{equation*} \nonl \hrulefill \setstretch{1} \nonl {Notes: (ref) details the Thresholding Multiple Outcomes (TMO) estimator proposed in this paper. Extensions of this procedure to settings with covariates, instrumental variables, and panel data are given in (ref). } \end{algorithm} It is worth emphasizing that we view the TMO procedure as as a complement, not a substitute, to standard error corrections based on geographic distance. Indeed, in the examples above, combining the TMO estimator with clustered, conley1999gmm, or SCPC standard errors often impacts the resultant estimate. In practice, combinations of this form are straightforward to implement and are likely to offer the greatest robustness to different forms of spatial correlation. We conclude by detailing some recommendations for practice. We focus our discussion on the construction of a relevant set of auxiliary outcomes. We again use bernini2023race as a running example to contextualize our recommendations. (ref) shows the ratio of the TMO standard error to the original standard error across different choices of the set of auxiliary outcomes. We conduct two experiments. First, we compute the TMO standard error using random subsamples of the 60 outcome variables sourced from the replication package associated with this paper. The circular points denote the average adjustment, where the size of the subsample varies along the $x$-axis. The intervals around each point show the 25-75th percentiles of the adjustment, over 20 random draws. The size of the adjustment increases with the number of outcomes.\footnote{(ref) reproduces Panel A of (ref) using a random subset of 20 outcomes.} With fewer outcomes, the null distribution has a higher variance and fewer degrees of freedom, making the optimal threshold higher, and reducing our ability to detect non-nulls. Ultimately, with a limited number of outcomes, the procedure has low power to distinguish true correlations from noise. This exercise highlights the importance of using a sufficiently large number of outcomes, so that the correlation estimates are sufficiently precise to distinguish pairs of locations with correlated errors. \begin{figure}[t] \begin{centering} \caption{TMO Adjustment for Different Sets of Auxiliary Outcomes} \begin{tabular}{c} \tabularnewline \end{tabular} \end{centering} \justifying {Notes: (ref) shows the ratio of the TMO standard error to the original standard error in bernini2023race, for different choices of auxiliary outcomes. The baseline TMO estimate uses the set of 60 \textit{Internal} outcomes from the bernini2023race replication file. The circular points show the average TMO adjustment across 20 randomly drawn subsets of the \textit{Internal} outcomes, where the number of randomly drawn outcomes in each subset is shown on the $x$-axis. For example, the circular point at $x=30$ indicates the average adjustment across 20 rounds of drawing 30 of the \textit{Internal} outcomes randomly and running TMO with those 30 outcomes. The intervals correspond to the interquartile range of the TMO adjustment across the 20 rounds. \textit{Same area} refers to the subset of \textit{Internal} outcomes that have a racial, demographic, or political characteristic like the main outcome of interest (the change in the share of black elected officials), while \textit{Different area} represents the subset of the remaining \textit{Internal} outcomes. The \textit{Changes} subset includes the \textit{Internal} outcomes that are expressed in differences, whereas the \textit{Levels} set includes those expressed in levels. The set of \textit{External} outcomes are from external sources such as the Census and are based on data from more recent time periods. \textit{Int + Ext} is the union of the \textit{Internal} and \textit{External} outcomes.} \hrulefill \end{figure} Second, we change various characteristics of the auxiliary outcomes. In bernini2023race, the outcome of interest is the change in the share of black elected officials in a long-difference regression. Among the 60 “Internal” outcomes taken from the replication file, around half are “Similar” with either a racial, demographic, or political component (e.g., change in NAACP branches), and the remaining half are “Different” outcomes that cover other areas (e.g., the unemployment rate). Though both sets have roughly 30 outcomes each, using the “Similar” set of outcomes produces a much stronger adjustment. Furthermore, using outcomes that are expressed in changes, like the outcome of interest, also generates a higher adjustment than using outcomes in levels. This exercise demonstrates the importance of identifying outcomes whose underlying spatial correlation is \emph{relevant} for the outcome of interest. These two experiments illustrate a fundamental trade-off between the number of auxiliary outcomes and their relevance to the outcome of interest. Adding more outcomes can lead to less informative estimates of the standard error if the additional outcomes do not reflect the structure of correlation in the outcome of interest. For example, (ref) shows that adding a set of around 60 “External” outcomes (e.g., median age from Census data) reduces the adjustment to the standard error, relative to using only the 60 “Internal” outcomes. In cases where researchers can choose from a large set of outcomes, we recommend that they aim to select the most relevant ones that provide at least 20 degrees of freedom.\footnote{We have found that larger degrees of freedom are necessary for reliable results in panel data settings.} Of course, in the application to bernini2023race, we do not know what the standard error \emph{should} be, i.e., whether larger adjustments are in fact indicative of better performance. (ref) displays the results of a set of experiments analogous to the results displayed in (ref), but in the context of the simulation exercise considered in (ref). Reassuringly, similar patterns emerge: the size of the adjustment increases---toward the true standard error---with both the number of auxiliary outcomes and the relevance of the outcomes to the underlying correlation structure. We also recommend two diagnostic exercises. The first is to plot the distribution of the correlation estimates and to check that the central mass of this distribution is well-approximated by a Gaussian. Heavily skewed distributions suggest that there is too much correlation across outcomes and that a larger sample of outcomes should be used. The second is to ensure that the function (ref) varies smoothly with the threshold. Otherwise, the estimate for the optimal threshold may be unreliable. These issues are more prevalent in cases with fewer locations, such as U.S.\ state or country level analyses. For example, funkeetal2023 study the effect of populist leaders on GDP per capita growth rates for 60 countries. (ref) displays these diagnostics implemented with data from this paper. We find that the Gaussian approximation to the null distribution is poor and that the estimate of the optimal threshold is unstable. Finally, we stress that we see more avenues for further research along the lines outlined in this paper. For example, although our estimator can be applied to panel settings, as we have done in some of the examples discussed above, data sets with panel structures have special features which merit their own analysis. Furthermore, there may be alternative ways to use information on spatial correlations obtained from multiple outcomes. For instance, this information could be used to identify clusters of outcomes which share a similar spatial correlation structure, as opposed to assuming a common correlation among all outcomes considered.

spacing{1.18}

appendix\setcounter{figure}{0} \setcounter{table}{0} \begin{center} {\it Supplemental Appendix to:} \vskip0.2cm \begin{spacing}{1} Using Multiple Outcomes to \\ Adjust Standard Errors for Spatial Correlation\daggerfootnote{Date: \today} \\ \begin{tabular}[t]{c@{\extracolsep{4em}}c} {Stefano DellaVigna} & {Guido Imbens}\\ {UC Berkeley and NBER} & {Stanford University} \\ \\ {Woojin Kim} & {David M. Ritzwoller}\\ {NBER} & {Stanford University} \\ \end{tabular} \end{spacing} \end{center} \begin{spacing}{1.13} \DoToC \end{spacing} \thispagestyle{empty} \setcounter{page}{0} \setcounter{figure}{0} \begin{spacing}{1.4} \section{Auxiliary Figures and Tables} \subsection{Figures} \begin{figure}[h] \caption{Spatial Dependence Across Counties in Selected States } \begin{centering} \begin{tabular}{c} \\ \end{tabular} \end{centering} \justifying { Notes: (ref) displays a correllogram reflecting the measurements $\hat{\rho}_{i,i^\prime}$ for each pair of U.S.\ counties in California, North Dakota, and New York State. Counties are sorted by population within each state.} \end{figure} \begin{figure}[h] \begin{centering} \caption{Spatial Dependence Across U.S. Counties} \begin{tabular}{c} \end{tabular} \end{centering} \justifying { Notes: (ref) displays a correllogram reflecting the measurements $\hat{\rho}_{i,i^\prime}$ for each pair of U.S.\ counties. Counties are sorted by region, state, and population within each state.} \end{figure} \begin{figure}[t] \begin{centering} \caption{Correlation Between Pairs of Counties (20 vs. 60 Auxiliary Outcomes)} \begin{tabular}{c} \tabularnewline \end{tabular} \end{centering} \justifying { Notes: (ref) builds on Panel A of (ref). As in Panel A of (ref), we plot the distribution of the pairwise correlation estimates, $\tilde{\rho}_{i,i^\prime}$, in green constructed using 60 auxiliary outcomes obtained from the replication package associated with bernini2023race. The solid curve denotes the null density estimate, obtained with the procedure outlined in (ref). The analogous histogram, constructed from a subsample of 20 randomly drawn outcomes, is displayed in red. In both cases, the associated optimal thresholds are displayed with vertical lines.} \end{figure} \begin{figure}[t] \begin{centering} \caption{Recommendations for Outcome Selection} \begin{tabular}{c} \tabularnewline \end{tabular} \end{centering} \justifying { Notes: (ref) shows the average ratio of the TMO standard error to the true standard error in the simulation developed in (ref). We vary the set of auxiliary outcomes used to construct the TMO estimator. The treatment variable is the change in the percentage of the population in county $i$ that has completed a college degree from 1980 to 2009. In each of the 1000 simulation rounds, the outcome of interest is a randomly drawn vector from a covariance structure based on recent changes in various economic, health, demographic, agricultural, and environmental variables in the last few decades (see (ref) and (ref) for details). The baseline Recent changes TMO adjustment uses 50 randomly drawn outcomes from the same covariance structure as the outcome of interest. The other Recent changes adjustments use fewer randomly drawn outcomes as shown on the $x$-axis, but still from the same covariance structure. The outcomes for Recent levels are drawn from a covariance structure based on the same outcomes listed in (ref), but the procedure in (ref) is performed on the average levels of these outcomes across years, rather than the changes between the years. The outcomes for Historical changes are drawn from a covariance structure built using the changes in different outcomes from the early to mid 1900s, rather than recent decades. \textit{Recent + Historical changes} is the union of 50 randomly drawn outcomes from the \textit{Recent changes} covariance structure and 50 randomly drawn outcomes from the \textit{Historical changes} covariance structure.} \end{figure} \begin{figure}[t] \begin{centering} \caption{Diagnostic Tests for TMO Assumptions funkeetal2023} \begin{tabular}{c} \textit{Panel A: Non-normal Distribution of Correlations} \\ \\ \textit{Panel B: Approximation to Non-Nulls Minus Nulls} \\ \tabularnewline \end{tabular} \end{centering} \justifying { Notes: (ref) is analogous to (ref), but is constructed using data from funkeetal2023. Panel A plots the distribution of the pairwise correlation estimates, $\tilde{\rho}_{i,i^\prime}$. Panel B shows the shows the value of the estimator (ref) over a range of values of the threshold $\delta$.} \end{figure} \end{spacing} \begin{landscape} \subsection{Tables} \begin{table}[H] \caption{Predicting Significant Positive Correlations Between Locations in U.S. Applications} { \begin{tabular} {l c c c c c c c c c} \toprule & & & \multicolumn{7}{c}{\textit{Among county/CZ pairs $ (i,i')$ with $ \hat{\rho}_{i,i'} \ge \hat{\delta}^\star$,}} \\ & & & \multicolumn{7}{c}{\textit{percentage in which $ i$ ranks among the 10% closest to $ i'$ in:}} \\ \cmidrule(lr){4-10} & $ \hat{\delta}^{*}$ & $ \%\ge\hat{\delta}^\star$ & Distance & Population & Urban % & Median income & Non-white % & Vote-share & Any of $ \leftarrow$ \\ & (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) \\ \hline Cook et al. (2023) & 0.42 & 2.89 & 46 & 23 & 27 & 22 & 21 & 16 & 81 \\ Caprettini and Voth (2023) & 0.51 & 0.88 & 35 & 16 & 16 & 15 & 16 & 15 & 64 \\ Esposito et al. (2023) & 0.19 & 3.22 & 25 & 16 & 15 & 16 & 15 & 14 & 61 \\ Bernini et al. (2023) & 0.54 & 0.67 & 34 & 17 & 18 & 20 & 20 & 16 & 68 \\ Bazzi et al. (2023) & 0.54 & 1.27 & 50 & 23 & 21 & 17 & 16 & 23 & 75 \\ Calderon et al. (2023) & 0.55 & 1.33 & 57 & 22 & 20 & 16 & 22 & 20 & 81 \\ Moscona and Sastry (2023) & 0.67 & 0.25 & 55 & 31 & 29 & 34 & 31 & 23 & 87 \\ Chetty et al. (2014) & 0.68 & 0.58 & 35 & 53 & 47 & 36 & 34 & 28 & 88 \\ \bottomrule \end{tabular} } \end{table} \justifying \begin{spacing}{1} { Notes: (ref) displays the results of exercise analogous to results reported in (ref). Column 1 reports the optimal threshold estimated for each paper. Column 2 shows the percent of county/CZ-pairs that have a positive correlation (i.e., not in absolute value) above the threshold. Columns 3-9 show the percent of county/CZ-pairs $(i,i^\prime)$ in which county/CZ $i$ is within the top 10% of those closest to $i'$ in the observable characteristic listed in the column header, among the pairs that have a correlation greater than a threshold $\hat{\delta}^{\star}$.} \end{spacing} \end{landscape} \begin{spacing}{1.4} \section{Details and Extensions} \subsection{Outcome Selection and Cleaning} We construct a dataset of 91 county variables that reflect typical variables used in empirical research. We take these variables from the following sources: \begin{itemize} • The Social Determinants of Health Database from the Agency for Healthcare Research \scriptsize{\url{https://www.ahrq.gov/sdoh/data-analytics/sdoh-data.html#download}} • Bailey et al. (2016). “U.S. County-Level Natality and Mortality Data, 1915-2007.” Inter-university Consortium for Political and Social Research [distributor]. \scriptsize{\url{https://doi.org/10.3886/E100229V4}} • U.S. Census Bureau. (2012). “Consolidated File: County Data, 1947-1977.” Inter-university Consortium for Political and Social Research [distributor]. \scriptsize{\url{https://doi.org/10.3886/ICPSR07736.v2}} • U.S. Census Bureau. (2011). “USA Counties.” Statistical Compendia. \scriptsize{\url{https://www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}} \end{itemize} The set of variables are listed in Tables (ref)-(ref). They cover a range of policy areas, including economics (e.g., income and labor force participation), health (e.g., life expectancy and health insurance coverage), demographics (e.g., population density and percent black), agriculture (e.g., cropland and percent employed in agriculture), and environment (e.g., water use and air pollution). We perform the following data processing steps. For variables that scale with population (e.g., the number of black residents, public school enrollment, and manufacturing establishments), we normalize by the county population. We take the logarithm for all positive variables that do not have negative values or are not in percentages. Then for each variable, we keep the first and last years for which there are data for at least 3000 counties. We then standardize the variable to have mean zero and unit standard deviation within each year and winsorize at the 0.1 and 99.9 percentiles to prevent results being skewed by outliers. Finally, to construct the long-difference outcomes, we take the difference between the first and last years for which there are data. \begin{landscape} \begin{table}[htbp!] \begin{threeparttable} \newcolumntype{L}{>{\arraybackslash}X} \newcolumntype{C}{>{\arraybackslash}X} \newcolumntype{S}{>{\arraybackslash}X} \caption{County outcomes used in simulation} \begin{tabularx}{\linewidth}{L l c c c} \toprule Description & Source & First year & Last year & No. counties \\ \hline 1. Annual mean of Particulate Matter (PM2.5) concentration (g/m3) & AHRQ & 2009 & 2018 & 3106 \\ 2. Average household size & AHRQ & 2009 & 2018 & 3105 \\ 3. Gini index of income inequality & AHRQ & 2009 & 2018 & 3105 \\ 4. Median distance in miles to the nearest emergency department, calculated using population weighted tract centroids in the county & AHRQ & 2009 & 2018 & 3105 \\ 5. Median distance in miles to the nearest obstetrics department, calculated using population weighted tract centroids in the county & AHRQ & 2009 & 2018 & 3105 \\ 6. Median distance in miles to the nearest pediatric ICU, calculated using population weighted tract centroids in the county & AHRQ & 2009 & 2018 & 3105 \\ 7. Median selected monthly owner costs for houses with a mortgage (dollars) & AHRQ & 2009 & 2018 & 3103 \\ 8. Percentage of families with children that are single-parent families & AHRQ & 2009 & 2018 & 3106 \\ 9. Percentage of households with same-sex unmarried partner & AHRQ & 2009 & 2018 & 3106 \\ 10. Percentage of limited English speaking households & AHRQ & 2009 & 2018 & 3106 \\ 11. Percentage of occupied housing units with utility gas heating & AHRQ & 2009 & 2018 & 3106 \\ 12. Percentage of population reporting American Indian and Alaska Native race alone & AHRQ & 2009 & 2018 & 3106 \\ 13. Percentage of population reporting multiple races & AHRQ & 2009 & 2018 & 3106 \\ 14. Percentage of population that speaks Spanish (ages 5 and over) & AHRQ & 2009 & 2018 & 3106 \\ 15. Percentage of unmarried partner households that received food stamps/SNAP benefits & AHRQ & 2009 & 2018 & 3102 \\ 16. Total cardiovascular disease deaths per 100,000 population (ages 35 and over) & AHRQ & 2009 & 2018 & 3106 \\ 17. Total expenditure (Dollars) per student & AHRQ & 2009 & 2018 & 3100 \\ 18. Total number of days with daily maximum heat index, absolute threshold: 90F & AHRQ & 2009 & 2018 & 3106 \\ 19. Total number of hospital beds per 1,000 population & AHRQ & 2009 & 2018 & 3106 \\ 20. Total number of hospitals per 1,000 population & AHRQ & 2009 & 2018 & 3106 \\ \bottomrule \end{tabularx} \begin{tablenotes} • \textit{Sources:} AHRQ: \href{www.ahrq.gov/sdoh/data-analytics/sdoh-data.html#download}{Agency for Healthcare Research and Quality}, Bailey et al. (2016): \url{www.doi.org/10.3886/E100229V4}, Country Data Book: \url{www.doi.org/10.3886/ICPSR07736.v2}, IHME: \href{www.doi.org/10.6069/3WQ2-TG23}{Institute For Health Metrics and Evaluation}, U.S. Census Bureau: \href{www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}{USA Counties}. USA Counties variables based on complete counts are indicated in the description. \end{tablenotes} \end{threeparttable} \end{table} \addtocounter{table}{-1} \begin{table}[htbp!] \begin{threeparttable} \newcolumntype{L}{>{\arraybackslash}X} \newcolumntype{C}{>{\arraybackslash}X} \newcolumntype{S}{>{\arraybackslash}X} \caption{County outcomes used in simulation} \begin{tabularx}{\linewidth}{L l c c c} \toprule Description & Source & First year & Last year & No. counties \\ \hline 21. Total standardized Medicare costs, fee for service (dollars) & AHRQ & 2009 & 2018 & 3104 \\ 22. Total stroke deaths per 100,000 population (ages 35 and over) & AHRQ & 2009 & 2018 & 3102 \\ 23. Birth Weight by Residence (2500 grams or less) & Bailey et al. (2016) & 1982 & 1988 & 3108 \\ 24. Nonmarital births by place of residence & Bailey et al. (2016) & 1968 & 1988 & 3097 \\ 25. Civilian Labor Force.Pct.Male & County Data Book & 1960 & 1970 & 3094 \\ 26. Divorce Rate Per 1000 Pop. & County Data Book & 1970 & 1975 & 3044 \\ 27. Employed.Pct.In Manfgr & County Data Book & 1950 & 1970 & 3053 \\ 28. Loc.Gov.Dir.Gen.Exp.Pct.Educ & County Data Book & 1967 & 1972 & 3090 \\ 29. Loc.Gov.Dir.Gen.Exp.Pct.Hghway & County Data Book & 1967 & 1972 & 3075 \\ 30. Loc.Gov.Prop.Tax.Per.Capita USD & County Data Book & 1962 & 1972 & 3088 \\ 31. Manfgr.Estab. & County Data Book & 1947 & 1972 & 3088 \\ 32. Marriage Rate Per 1000 Pop. & County Data Book & 1970 & 1975 & 3101 \\ 33. Oasdhi Payments Per Mon. USD1000 & County Data Book & 1971 & 1976 & 3104 \\ 34. Ou.Pct.Owner Occupied & County Data Book & 1940 & 1970 & 3087 \\ 35. Population Pct Foreign Stock & County Data Book & 1960 & 1970 & 3095 \\ 36. Population Pct. Urban & County Data Book & 1960 & 1970 & 3026 \\ 37. Population Rank & County Data Book & 1950 & 1960 & 3092 \\ 38. Pub.Assis.Recipients.Afdc. & County Data Book & 1972 & 1976 & 3079 \\ 39. Retail Trade Estab. & County Data Book & 1954 & 1972 & 3096 \\ 40. Retail Trade Estab.Sales USD1000 & County Data Book & 1948 & 1972 & 3076 \\ \bottomrule \end{tabularx} \begin{tablenotes} • \textit{Sources:} AHRQ: \href{www.ahrq.gov/sdoh/data-analytics/sdoh-data.html#download}{Agency for Healthcare Research and Quality}, Bailey et al. (2016): \url{www.doi.org/10.3886/E100229V4}, Country Data Book: \url{www.doi.org/10.3886/ICPSR07736.v2}, IHME: \href{www.doi.org/10.6069/3WQ2-TG23}{Institute For Health Metrics and Evaluation}, U.S. Census Bureau: \href{www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}{USA Counties}. USA Counties variables based on complete counts are indicated in the description. \end{tablenotes} \end{threeparttable} \end{table} \addtocounter{table}{-1} \begin{table}[htbp!] \begin{threeparttable} \newcolumntype{L}{>{\arraybackslash}X} \newcolumntype{C}{>{\arraybackslash}X} \newcolumntype{S}{>{\arraybackslash}X} \caption{County outcomes used in simulation} \begin{tabularx}{\linewidth}{L l c c c} \toprule Description & Source & First year & Last year & No. counties \\ \hline 41. Workers.Pct.Used Pub.Trans & County Data Book & 1960 & 1970 & 3095 \\ 42. Life expectancy & IHME & 1999 & 2019 & 3113 \\ 43. All persons 18 to 64 years without health insurance, percent & U.S. Census Bureau & 2005 & 2007 & 3108 \\ 44. All persons under 18 years without health insurance, percent & U.S. Census Bureau & 2005 & 2007 & 3108 \\ 45. Average age of farm operators (NAICS) & U.S. Census Bureau & 2002 & 2007 & 3064 \\ 46. Average travel time to work for workers 16 years and over who did not work at home & U.S. Census Bureau & 1990 & 2009 & 3107 \\ 47. Average value of land and buildings per acre (NAICS) & U.S. Census Bureau & 2002 & 2007 & 3058 \\ 48. Average value of land and buildings per farm (NAICS) & U.S. Census Bureau & 2002 & 2007 & 3058 \\ 49. Births per 1,000 population & U.S. Census Bureau & 1970 & 2007 & 3102 \\ 50. Black population (complete count) & U.S. Census Bureau & 1990 & 2010 & 3106 \\ 51. Civilian labor force & U.S. Census Bureau & 1990 & 2010 & 3106 \\ 52. Civilian labor force unemployment rate & U.S. Census Bureau & 1990 & 2010 & 3107 \\ 53. Commercial banks and savings institutions (FDIC-insured) - total deposits & U.S. Census Bureau & 1980 & 2010 & 3103 \\ 54. Cropland - harvested cropland (NAICS) (acres) & U.S. Census Bureau & 2002 & 2007 & 2963 \\ 55. Cropland - total (NAICS) (acres) & U.S. Census Bureau & 2002 & 2007 & 3051 \\ 56. Deaths per 1,000 population & U.S. Census Bureau & 1970 & 2007 & 3102 \\ 57. Earnings in retail trade (NAICS 44-45) & U.S. Census Bureau & 2001 & 2007 & 2989 \\ 58. Educational attainment - persons 25 years and over - percent bachelor's degree or higher & U.S. Census Bureau & 1980 & 2009 & 3105 \\ 59. Educational attainment - persons 25 years and over - percent high school graduate or higher & U.S. Census Bureau & 1980 & 2009 & 3105 \\ 60. Employment in farming (NAICS) & U.S. Census Bureau & 2001 & 2007 & 3077 \\ \bottomrule \end{tabularx} \begin{tablenotes} • \textit{Sources:} AHRQ: \href{www.ahrq.gov/sdoh/data-analytics/sdoh-data.html#download}{Agency for Healthcare Research and Quality}, Bailey et al. (2016): \url{www.doi.org/10.3886/E100229V4}, Country Data Book: \url{www.doi.org/10.3886/ICPSR07736.v2}, IHME: \href{www.doi.org/10.6069/3WQ2-TG23}{Institute For Health Metrics and Evaluation}, U.S. Census Bureau: \href{www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}{USA Counties}. USA Counties variables based on complete counts are indicated in the description. \end{tablenotes} \end{threeparttable} \end{table} \addtocounter{table}{-1} \begin{table}[htbp!] \begin{threeparttable} \newcolumntype{L}{>{\arraybackslash}X} \newcolumntype{C}{>{\arraybackslash}X} \newcolumntype{S}{>{\arraybackslash}X} \caption{County outcomes used in simulation} \begin{tabularx}{\linewidth}{L l c c c} \toprule Description & Source & First year & Last year & No. counties \\ \hline 61. Employment in government - state and local (NAICS) & U.S. Census Bureau & 2001 & 2007 & 3078 \\ 62. Employment in retail trade (NAICS 44-45) & U.S. Census Bureau & 2001 & 2007 & 2988 \\ 63. Federal Government direct loans FY & U.S. Census Bureau & 1983 & 2010 & 3087 \\ 64. Hospital insurance and/or supplemental medical insurance (Medicare) - aged persons enrolled & U.S. Census Bureau & 1998 & 2007 & 3084 \\ 65. Infant deaths per 1,000 live births & U.S. Census Bureau & 1990 & 2007 & 3107 \\ 66. Land in farms (NAICS) (acres) & U.S. Census Bureau & 2002 & 2007 & 3026 \\ 67. Median contract rent of specified renter-occupied housing units paying cash rent (complete count) & U.S. Census Bureau & 1980 & 2009 & 3105 \\ 68. Median household income & U.S. Census Bureau & 1995 & 2009 & 3107 \\ 69. Median selected monthly owner costs of specified owner-occupied housing units with a mortgage & U.S. Census Bureau & 1980 & 2000 & 3106 \\ 70. Median value of specified owner-occupied housing units (complete count) & U.S. Census Bureau & 1980 & 2009 & 3105 \\ 71. New private housing units authorized by building permits - total & U.S. Census Bureau & 2004 & 2010 & 3107 \\ 72. People of all ages in poverty - percent & U.S. Census Bureau & 1995 & 2009 & 3107 \\ 73. People under age 18 in poverty - percent & U.S. Census Bureau & 1995 & 2009 & 3107 \\ 74. Per capita personal income & U.S. Census Bureau & 1969 & 2000 & 3075 \\ 75. Persons 16 to 19 years not enrolled in school and not a high school graduate & U.S. Census Bureau & 1990 & 2000 & 3108 \\ 76. Population per square mile & U.S. Census Bureau & 1980 & 2010 & 3105 \\ 77. Private nonfarm establishments & U.S. Census Bureau & 1990 & 2009 & 3107 \\ 78. Private nonfarm establishments - arts, entertainment and recreation (NAICS 71) & U.S. Census Bureau & 2002 & 2009 & 3108 \\ 79. Public school enrollment Fall & U.S. Census Bureau & 1987 & 2009 & 3088 \\ 80. Related children age 5 to 17 in families in poverty - percent & U.S. Census Bureau & 1995 & 2009 & 3107 \\ \bottomrule \end{tabularx} \begin{tablenotes} • \textit{Sources:} AHRQ: \href{www.ahrq.gov/sdoh/data-analytics/sdoh-data.html#download}{Agency for Healthcare Research and Quality}, Bailey et al. (2016): \url{www.doi.org/10.3886/E100229V4}, Country Data Book: \url{www.doi.org/10.3886/ICPSR07736.v2}, IHME: \href{www.doi.org/10.6069/3WQ2-TG23}{Institute For Health Metrics and Evaluation}, U.S. Census Bureau: \href{www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}{USA Counties}. USA Counties variables based on complete counts are indicated in the description. \end{tablenotes} \end{threeparttable} \end{table} \addtocounter{table}{-1} \begin{table}[htbp!] \begin{threeparttable} \newcolumntype{L}{>{\arraybackslash}X} \newcolumntype{C}{>{\arraybackslash}X} \newcolumntype{S}{>{\arraybackslash}X} \caption{County outcomes used in simulation} \begin{tabularx}{\linewidth}{L l c c c} \toprule Description & Source & First year & Last year & No. counties \\ \hline 81. Renter-occupied housing units (complete count) & U.S. Census Bureau & 1990 & 2010 & 3106 \\ 82. Resident population under 18 years, percent & U.S. Census Bureau & 2000 & 2009 & 3108 \\ 83. Resident population: Black alone, percent & U.S. Census Bureau & 2000 & 2009 & 3059 \\ 84. Resident population: Median age (complete count) & U.S. Census Bureau & 1980 & 2010 & 3105 \\ 85. Resident population: total females, percent & U.S. Census Bureau & 2000 & 2009 & 3108 \\ 86. Total physicians & U.S. Census Bureau & 2004 & 2009 & 3108 \\ 87. Valuation of new private housing units authorized by building permits & U.S. Census Bureau & 2004 & 2010 & 3108 \\ 88. Value of farm products sold - total (NAICS) & U.S. Census Bureau & 2002 & 2007 & 3003 \\ 89. Vehicles available per occupied housing unit & U.S. Census Bureau & 1990 & 2000 & 3108 \\ 90. Vote cast for president- percent Republican & U.S. Census Bureau & 1980 & 2008 & 3105 \\ 91. Water use: per capita use & U.S. Census Bureau & 1990 & 2005 & 3107 \\ \bottomrule \end{tabularx} \begin{tablenotes} • \textit{Sources:} AHRQ: \href{www.ahrq.gov/sdoh/data-analytics/sdoh-data.html#download}{Agency for Healthcare Research and Quality}, Bailey et al. (2016): \url{www.doi.org/10.3886/E100229V4}, Country Data Book: \url{www.doi.org/10.3886/ICPSR07736.v2}, IHME: \href{www.doi.org/10.6069/3WQ2-TG23}{Institute For Health Metrics and Evaluation}, U.S. Census Bureau: \href{www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}{USA Counties}. USA Counties variables based on complete counts are indicated in the description. \end{tablenotes} \end{threeparttable} \end{table} \end{landscape} \subsection{Simulation Design} For the simulation considered in (ref), we construct covariance matrix $\Sigma$ for the residual vector $\varepsilon$. Our aim is to ensure that the off-diagonal elements of $\Sigma$ primarily consist of “non-null” pairs of counties, while maintaining positive definiteness. To meet these objectives, we take as input the correlation matrix of the standardized outcomes. Next, we set a threshold $\delta$, above which, we believe most pairs of counties with larger correlations should be non-nulls. After some experimentation, we set $\delta=0.45$. We then find a “central” county that has the highest number of correlations with other counties that are above 0.45 in absolute value. That central county and all the counties with which it is sufficiently correlated are designated as a cluster. We find the next “central” county, among those that have not been assigned to a cluster, that has the highest number of correlations above 0.45 with other counties that have also yet to be assigned to a cluster. Note that once a county has been assigned to a cluster (as the central one or not), then it cannot be assigned to another cluster. We repeat the process, and once all the correlations above 0.45 have been exhausted, we set all correlations between counties outside the clusters to 0, and keep all the correlations within clusters intact. The resulting correlation matrix, ordered by clusters, has a block-diagonal structure that fulfills the objectives. \subsection{Alternative Estimators and Data Structures} In this section, we detail several extensions to the method proposed in (ref). We keep the same notation. That is, we consider the linear model \begin{equation} Y^{(0)}_i = \alpha + \tau^{(0)} W_i + \theta^\top X_i \varepsilon^{(0)}_i,\quad i = 1,\ldots,n , \end{equation} where $Y^{(0)}_i$ measures some outcome of interest, $W_i$ denotes some treatment of interest, and $X_i$ denotes a $k$-vector of covariates. The variables $Y^{(1)}_i,\ldots,Y^{(d)}_i$ denote $d$ auxiliary, post-treatment outcomes. The variables $\varepsilon_i^{(j)}$ and $\hat{\varepsilon}_i^{(j)}$ denote the population and empirical residuals associated with the regression of the outcome $Y_i^{(j)}$ on the treatment $W_i$, respectively. \subsubsection{Covariates} In the main text, for ease of exposition, we treat the case where there are no covariates. The extension to the case with covariates is straightforward, by the Frisch-Waugh-Lovell theorem. In particular, let $\tilde{W_i}$ and $\tilde{Y}_i^{(j)}$ denote the residuals from regressing $W_i$ and $\tilde{Y}_i^{(j)}$ on $X_i$, respectively. A standard error for $\hat{\tau}$ can then be constructed as before, by using these data. \subsubsection{Weighted regression} Throughout the main text, we study the construction of standard errors for ordinary least squares estimators. The proposed method can be easily adapted to weighted least squares estimators. In particular, suppose that we are interested in the estimator $\hat{\tau}^{(0)}_h$ defined by \begin{equation} \hat{\tau}^{(0)}_h = \arg \min_{\tau} \left\{ \sum^n_{i=1} h_i (Y^{(0)}_i - \tau W_i)^2 \right\} , \end{equation} where $h=(h_i)_{i=1}^n$ is an exogenous weight vector. Observe that, in this case, we have that \begin{flalign} \Var(\hat{\tau}^{(0)}_h \mid W) &= S_n(h)^{-2} (W^\top H \Sigma^{0} H W),\quad\text{where} \nonumber \\ S_n(h) &= W^\top H W\quad\text{and}\quad H = \mathsf{diag}(h) . \end{flalign} In other words, the problem again reduces to the estimation of the residual covariance matrix $\Sigma^{0} = \Var(\varepsilon^{(0)})$ and our method for estimating this matrix applies. \subsubsection{Instrumental Variables Regression} For each unit $i$, suppose that we observe the instrumental variable $Z_i$. If we are interested in constructing a standard error for the two-stage least squares estimator, we can simply apply our method to the second-stage regression. We caution that this approach is not robust to weak identification. The construction of identification robust standard error estimates for settings with spatial correlation is an interesting direction for future research. See andrews2019weak for a recent review of identification robust standard errors for linear instrumental variables regression in non-homoskedastic data. \subsubsection{Panel Regressions} Suppose that we observe the collection of outcomes across $t$ time periods. That is, we observe $(Y_{i,s}^{(0)}, Y^{(1)}_{i,s},\ldots,Y^{(d)}_{i,s})$ for each unit $i$ in $1,\ldots,n$ and period $s$ in $1,\ldots,t$. Here, we are often interested in constructing standard errors for coefficients in regressions of the form \begin{equation} Y^{(0)}_{i,s} = \alpha + \tau^{(0)} W_{i,s} + \theta^\top X_{i,s} \varepsilon^{(0)}_{i,s} , \end{equation} where, in many cases, the covariates $X_{i,s}$ contain various fixed effects. To adjust standard errors for correlation across units, there are at least two options. First, we might account for the correlation between the units $i$ and $i^\prime$ only \textit{within} each time period. That is, the residuals $(\varepsilon^{(j)}_{i,s},\varepsilon^{(j)}_{i^{\prime},s})$ are allowed to be correlate, but the residuals $(\varepsilon^{(j)}_{i,s},\varepsilon^{(j)}_{i^{\prime,s^\prime}})$, for $s{\neq}s^{\prime}$, are not. We note that an analogous assumption is widely applied in treatments of two-way or multi-way clustered data menzel2021bootstrap. An alternative option would be to allow for correlation \textit{across} time periods as well. For the purposes of this paper, we adopt the latter approach. A more thorough consideration of the construction of standard errors for panel data, that accommodates spatial correlation across units is useful direction for further research. In our treatment, we estimate the covariance between residuals associated with units $i$ and $i^\prime$ with the estimator \begin{align} \hat{\lambda}_{i,i^{\prime}} = \frac{1}{dt} \sum_{j=1}^d \sum_{s=1}^t (\tilde{\varepsilon}^{(j)}_{i,t} - \bar{\varepsilon}_i) (\tilde{\varepsilon}^{(j)}_{i^{\prime},t} - \bar{\varepsilon}_{i^{\prime}}) \end{align} That is, in effect, we treat each outcome-period pair as a distinct outcome. These estimates are then be used to determine which pairs of units are “non-nulls” with the procedure developed in the main text. Standard errors are then adjusted, as before, allowing for correlation across units and across time periods. \subsection{Augmenting Distance-Based Estimators} The TMO estimator can be used to augment existing spatial standard errors that are based on modeling dependence in terms of geographic distance. Recall that, for the method detailed in (ref), the covariance $\sigma^{(0)}_{i,i^\prime}$ is estimated by the quantity \begin{equation} \hat{\sigma}^{(0)}_{i,i^\prime}(\hat{\delta}^\star) = \begin{cases} \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_i ,& i = i^\prime, \\ \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}\text \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \geq \hat{\delta}^\star \} , & i\neq i^\prime . \end{cases} \end{equation} That is, the TMO estimator does not threshold diagonal elements of the residual covariance matrix---in effect, augmenting the standard white1980heteroskedasticity heteroskedasticity consistent variance estimate. The same idea can be applied to augment alternative standard error estimates. Let $\mathcal{C} \subseteq [n]^2$ denote a subset of the pairs of units $(i,i^\prime)$ that includes all pairs of the form $(i,i)$. These might be pairs of units in the same cluster, or pairs of units whose geographic distance is smaller than some threshold. Suppose that, moreover, we have at our disposal an alternative estimator of the covariance $\sigma^{(0)}_{i,i^\prime}$ for each pair $(i,i^\prime)$ in $\mathcal{C}$. Denote this estimator by $\tilde{\sigma}^{(0)}_{i,i^\prime}$. This estimator could be simply $\tilde{\sigma}^{(0)}_{i,i^\prime}=\hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}$ or it could be, e.g., the estimator associated with a conley1999gmm standard error with a chosen kernel and bandwidth. We can construct the covariance matrix estimate $\hat{\Sigma}(\delta)$ by collecting the estimates \begin{equation} \hat{\sigma}^{(0)}_{i,i^\prime}(\hat{\delta}^\star) = \begin{cases} \tilde{\sigma}^{(0)}_{i,i^\prime} ,& (i, i^\prime)\in\mathcal{C}, \\ \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime}\text \mathbb{I}\{ \vert \hat{\rho}_{i,i^\prime}\vert \geq \hat{\delta}^\star \} , & (i, i^\prime)\not\in\mathcal{C} . \end{cases} \end{equation} In this case, the optimal threshold $\hat{\delta}^\star$ should be estimated by considering the distribution of only those pairs of units $(i, i^\prime)$ not in $\mathcal{C}$. A similar approach can be taken to use TMO to augment the muller2022spatial,muller2023spatial SCPC standard error estimator. muller2022spatial,muller2023spatial estimate the variance $V(\Sigma^{(0)})$ with an estimator of the form \begin{equation} \hat{V}_{\text{SCPC}}(q) = \frac{S_n^{-2}}{q} \sum_{j=1}^q r_j^\top (W^\top \hat{\varepsilon}^{(0)} (\hat{\varepsilon}^{(0)})^\top W) r_j^\top , \end{equation} where $r_j$ is an approximation to the $j$th principal component of $W^\top \Var(\varepsilon)W$. SCPC can, thus, be augmented with the TMO standard error by taking \begin{equation} \hat{V}_{\text{SCPC} +\text{TMP}}(\delta, q) = \frac{S_n^{-2}}{q} \sum_{j=1}^q r_j^\top (W^\top \hat{\Sigma}_{SCPC}(\delta) W) r_j^\top + S_n^{-2}(W^\top \hat{\Sigma}(\delta) W) \end{equation} where $\hat{\Sigma}_{SCPC}(\delta)$ collects \begin{equation} \hat{\sigma}^{\text{SCPC}}_{i,i^\prime,}(\delta) = \begin{cases} \hat{\varepsilon}^{(0)}_i \hat{\varepsilon}^{(0)}_{i^\prime} ,& \vert \hat{\rho}_{i,i^\prime}\vert \leq \delta, \\ 0 , & \vert\hat{\rho}_{i,i^\prime}\vert \geq \delta , \end{cases} \end{equation} and $\hat{\Sigma}(\delta)$ collects the estimates (ref), as before. \begin{figure}[t] \begin{centering} \caption{Robustness to Alternative Null Distribution Estimators} \begin{tabular}{c} \tabularnewline \end{tabular} \end{centering} \justifying { Notes: (ref) builds on Panel A of (ref). As in Panel A of (ref), we plot the distribution of the pairwise correlation estimates, $\tilde{\rho}_{i,i^\prime}$ using 60 auxiliary outcomes obtained from the replication package associated with bernini2023race. The solid black curve denotes the null density estimate, obtained with the procedure outlined in (ref). The blue and red curves display the fits using the methods outlined in (ref) that give the largest and smallest threshold. A $\chi^2$ approximation is also displayed.} \end{figure} \subsection{Alternative Null Distribution Estimators} In this section, we consider alternative methods for estimating the null distribution. The methods we consider use the entire vector of off-diagonal correlations between the $q$ and $1-q$ percentiles, where we vary $q\in\{0.1,0.2,0.25\}$. We group the correlations into granular bins. For a candidate variance $v$, we calculate the “distance” between the empirical and null distributions in one of two ways. The first computes the difference in the density at the center of each bin versus the null density at the midpoints. The second calculates the difference in the mass within each bin versus the mass within the endpoints of each bin under the null. Differences are calculated using either $\ell^1,\ell^2,$ or $\ell^\infty$ norms. We compute the variance $v$ that minimizes the difference. (ref) is analogous to Panel A of (ref), but additionally displays the fits to the null distribution and maximize and minimize the resultant threshold. We also display a approximation to the null distribution based on a $\chi^2$ distribution. The quality of the approximation, and the resultant threshold are insensitive to how the null distribution is estimated. \section{Proofs for Results Stated in the Main Text} \subsection{Proof of (ref)} Throughout, we let the operator $\sim$ denote equality in distribution. Fix any constant $b>0$ and define the matrix \begin{equation} \Sigma_b = b \cdot R , \quad\text{where}\quad R = W(W^\top W)^{-1} W^\top . \end{equation} Suppose that, conditional on $W$, each pair of population residual vectors $\varepsilon^{(j)}$ and $\varepsilon^{(j')}$ had the joint distribution \begin{equation} { \begin{pmatrix} \varepsilon^{(j)} \\ \varepsilon^{(j')} \end{pmatrix} \sim \mathsf{N}\!\left( \begin{pmatrix} \bm{0}_n \\ \bm{0}_n \end{pmatrix},\, \begin{pmatrix} \Sigma_b & \Sigma_b \\ \Sigma_b & \Sigma_b \end{pmatrix} \right) } . \end{equation} In this case, the data $Y$ are jointly Gaussian conditional on $W$ and (ref) holds. Fix $\tau^{(j)} = 0$ for all $j$ in $0, 1, ..., d$. The parameter of interest is given by \begin{align} V(\Sigma_b) & = (W^\top W)^{-1} W^\top \Sigma_b W (W^\top W)^{-1} = b (W^\top W)^{-1} . \end{align} Observe that the collection \begin{equation} (Y, W) = (Y^{(0)}, Y^{(1)}, \ldots, Y^{(d)}, W) \end{equation} can be reconstructed from \begin{equation} (\hat{Y}, \hat{\varepsilon}, W) = (\hat{Y}^{(0)}, \hat{\varepsilon}^{(0)}, \hat{Y}^{(1)}, \hat{\varepsilon}^{(1)}, \ldots, \hat{Y}^{(d)}, \hat{\varepsilon}^{(d)}, W) \end{equation} where \begin{equation} \hat{Y}^{(j)} = R Y^{(j)} \quad\text{and}\quad \hat{\varepsilon}^{(j)} = (I_n - R)Y^{(j)} . \end{equation} Moreover, conditional on $W$, we have that \begin{equation} { (\hat{Y}^{(j)}, \hat{Y}^{(j')}, \hat{\varepsilon}^{(j”)}, \hat{\varepsilon}^{(j”')})^{\top} \sim \mathsf{N}\!\left( \begin{pmatrix} \bm{0}_{2n} \\ \bm{0}_{2n} \end{pmatrix},\, \begin{pmatrix} b R \otimes \bm{1}_{2\times2} & \bm{0}_{2n\times 2n} \\ \bm{0}_{2n\times 2n} & \bm{0}_{2n\times 2n} \end{pmatrix} \right) , } \end{equation} where we have used the facts that \begin{equation} M \Sigma_b^{(0)} M = bM R M = \bm{0}_{n\times n} \quad\text{and}\quad b R^\top R R = b R . \end{equation} Here, the operator $\otimes$ denotes the Kronecker product. Thus, there exists some locally bounded, measurable function $\widecheck{V}^\prime(\cdot)$ such that \begin{equation} \widecheck{V}^\prime(\hat{Y}, W) = \widecheck{V}(Y, W) \end{equation} almost surely. Let $\bar{r}$ denote the component of $R$ that is largest in absolute value. Observe that there exists a constant $c$, that does not depend on $n$, $d$, or $W$ such that the set \begin{equation} A(b) = \left\{ \hat{y} = (\hat{y}^{(j)})_{j=0}^d \in \mathbb{R}^{(d+1)\times n}: \max_{i=1\ldots n} \max_{j=0,\ldots,d} \{ \vert \hat{y}^{(j)}_i\vert\} \leq c\sqrt{b \bar{r}\log(dn)} \right\} \end{equation} satisfies \begin{equation} P\left\{ \hat{Y} \in A(b)\right\} \geq 1-\eta . \end{equation} for all $b>0$. Define the functions \begin{equation} M(b) = \sup\{ \widecheck{V}^\prime(\hat{y}, W) : \hat{y} \in A(b) \} \quad\text{and}\quad m(b) = \inf\{ \widecheck{V}^\prime(\hat{y}, W) : \hat{y} \in A(b) \} , \end{equation} respectively. Both functions are finite, as $\widecheck{V}^\prime(\cdot)$ is locally bounded. Suppose that for each $b>0$, we have both \begin{equation} b (W^\top W)^{-1} - M(b) \leq K \quad\text{and}\quad m(b) - b (W^\top W)^{-1} \leq K . \end{equation} In this case, it would hold that \begin{equation} \widecheck{V}^\prime(\hat{y}, W) \in b (W^\top W)^{-1} \pm K \end{equation} for each $\hat{y} \in A(b)$ and every $b>0$, where we recall that $V(\Sigma_b) = b (W^\top W)^{-1}$ . Now, for any fixed $\hat{y}_0$, the parameter $b_0$ can be chosen such that $\hat{y}_0 \in A(b)$ for all $b$ greater than $b_0$. But then, it would be the case that \begin{equation} \widecheck{V}^\prime(\hat{y}_0, W) \in b (W^\top W)^{-1} \pm K \end{equation} for all $b$ greater than $b_0$. However, this must fail for some sufficiently large $b$, giving a contradiction. Consequently, there exists $b>0$ such that \begin{equation} b (W^\top W)^{-1} - M(b) \leq K \quad\text{or}\quad m(b) - b (W^\top W)^{-1} \leq K . \end{equation} In the first case, for all $\hat{y} \in A(n)$, we have that $\widecheck{V}^\prime(\hat{y}_0, W) < V(\Sigma_b) - K $, and so \begin{equation} P\left\{ \widecheck{V}(Y, W) - V(\Sigma_b) < -K \right\} \geq P\left\{ A(b) \right\} \geq 1-\eta . \end{equation} In the second case, for all $\hat{y} \in A(n)$, we have that $\widecheck{V}^\prime(\hat{y}_0, W) > V(\Sigma_b) + K $, and so \begin{equation} P\left\{ \widecheck{V}(Y, W) - V(\Sigma_b) > K \right\} \geq P\left\{ A(b) \right\} \geq 1-\eta , \end{equation} which completes the proof.\qed \subsection{Proof of (ref)} Consider the loss \begin{equation} L(\delta) = (1-p_0) \int_0^\delta t\text{ d}f_1(t) + p_0 \int_\delta^\infty t\text{ d}f_0(t) . \end{equation} As we have assumed that the functions $f_0(t)$ and $f_1(t)$ are continuously differentiable, we can evaluate \begin{align} \frac{\text{d}}{\text{d}\delta} \left((1-p_0) \int_0^\delta t\text{ d}f_1(t) \right) & = (1-p_0) \delta f_1(\delta)\quad\text{and}\quad \frac{\text{d}}{\text{d}\delta} \left(p_0 \int_\delta^\infty t\text{ d}f_0(t) \right) = -p_0 \delta f_0(\delta) , \end{align} by Leibniz's rule. Thus, we find that \begin{equation} \frac{\text{d}L(\delta)}{\text{d}\delta} = \delta( (1-p_0)f_1(\delta) - p_0 f_0(\delta)) . \end{equation} For $\delta>0$, equalizing the derivative (ref) with zero gives \begin{equation} (1-p_0)f_1(\delta) = p_0 f_0(\delta) , \end{equation} or equivalently \begin{equation} \mathsf{fdr}(\delta) = \frac{p_0 f_0(\delta) }{p_0 f_0(\delta) + (1-p_0)f_1(\delta)} = \frac{1}{2} . \end{equation} Hence, it will suffice to show that $\mathsf{fdr}(\delta)$ is strictly decreasing and continuous, as, in this case, there is a unique $\delta^\star$ that solves the first-order condition (ref). Observe that the continuity of $\mathsf{fdr}(\delta)$ follows immediately from the continuity of $f_0(t)$ and $f_1(t)$. To verify that it is strictly decreasing, observe that \begin{equation} \frac{\text{d}}{\text{d}\delta} \mathsf{fdr}(\delta) = \frac{p_0(1-p_0)(f^\prime_0(\delta)f_1(\delta) - f_0(\delta)f^\prime_1(\delta))}{(p_0f_0(\delta) - (1-p_0)f_1(\delta))^2} . \end{equation} The assumption that the distribution of the non-nulls is stochastically larger than the distribution of the nulls, implies that the likelihood ratio $f_0(\delta)/f_1(\delta)$ is decreasing. But notice that \begin{equation} \frac{\text{d}}{\text{d}\delta} \log\left( \frac{f_0(\delta)}{f_1(\delta)}\right) = \frac{f^\prime_0(\delta)}{f_0(\delta)} - \frac{f^\prime_1(\delta)}{f_1(\delta)} . \end{equation} and that the sign of (ref) is determined by the term \begin{equation} f^\prime_0(\delta)f_1(\delta) - f_0(\delta)f^\prime_1(\delta)) = f_0(\delta)f_1(\delta) \left( \frac{f^\prime_0(\delta)}{f_0(\delta)} - \frac{f^\prime_1(\delta)}{f_1(\delta)} \right) . \end{equation} Hence, the function $\mathsf{fdr}(\delta)$ is decreasing, which completes the proof.\qed \subsection{Proof of (ref)} Let $Y$ and denote the $n\times d$ matrix whose $(i,j)$th element is given by $Y_i^{(j)}$. Throughout, we let $\|\cdot\|_{\text{F}}$ and $\|\cdot\|_{\text{op}}$ denote the matrix Frobenius and $\ell_2$-operator norms, respectively. (ref) requires a slightly strengthened version of (ref), stated as follows. \begin{assumption}[Sub-Gaussian Proportionality] Let $Z$ denote an $n\times d$ matrix whose components are independent, mean-zero, and sub-Gaussian with parameter $M \geq 1$. There exist positive semi-definite matrices $\Lambda_n = (\lambda_{i,i^\prime})_{i,i^\prime = 1}^n$ and $\Gamma_d = (\gamma_{j,j^\prime})_{j,j^\prime=1}^d$ such that \begin{equation} Y = \bm{1}_n \tau^\top + \Lambda_n^{1/2} Z \Gamma_n^{1/2} , \end{equation} where $A^{1/2}$ denotes the square-root of the matrix $A$. \end{assumption} In addition to the restriction \begin{equation} \Cov(Y_i^{(j)}, Y_{i^\prime}^{(j^\prime)}) = \lambda_{i,i^\prime} \gamma_{j,j^\prime} \end{equation} specified by (ref), (ref) implies that (i) the elements of $Y$ are sub-Gaussian and that (ii) $Y$ can be normalized to have fully independent, rather than just uncorrelated, entries by taking $\Lambda_n^{-1/2}(Y-\bm{1}_n \tau^\top)\Gamma_n^{-1/2}$. Analogous versions of this assumption appear elsewhere in the literature hoff2016limitations. Define the quantities \begin{equation} \bar{\gamma} = \max_{j\in[d]} \gamma_{j,j},\quad \underline{\gamma} = \min_{j\in[d]} \gamma_{j,j},\quad \bar{\lambda} = \max_{i\in[n]} \lambda_{i,i},\quad\text{and}\quad \underline{\lambda} = \max_{i\in[n]} \lambda_{i,i} . \end{equation} To ease exposition, in the statement of (ref), the constants $c$ and $C$ are allowed to depend on the constants (ref) as well as the sub-Gaussian constant $M$ defined in (ref). For the sake of completeness, we track the dependence on these terms throughout the proof. The result follows from the following two large deviation bounds. Throughout, we write $\xi_j = 1/\gamma_{j,j}$ and $\hat{\xi}_j = \hat{\gamma}_0/\hat{\gamma}_j$. Recall that we have normalized $\gamma_{0,0} = 1$. Let $\Xi_d$ and $\widehat{\Xi}_d$ denote the diagonal matrices with entries $(\xi_j)_{j=1}^d$ and $(\hat{\xi}_j)_{j=1}^d$, respectively. Let $\widehat{\Lambda}_{n,d}$ denote the matrix whose $(i,i^\prime)$th component is $\hat{\gamma}_0 \hat{\lambda}_{i,i^\prime}$. Define the analogous infeasible estimator \begin{equation} \tilde{\Lambda}_{n,d} = \frac{1}{d} Y \Xi_d Y^\top \end{equation} and let $\tilde{\lambda}_{i,i^\prime}$ denote its $(i,i^\prime)$th component. \begin{lemma} Fix a constant $0<\phi<1$.\\ (i) Suppose that (ref) hold. There exist constants $c$ and $C$ such that if \begin{equation} \frac{\bar{\lambda}^{1-q/2}}{\underline{\lambda}}\frac{\bar{\gamma}^2}{\underline{\gamma}} \sqrt{\frac{\kappa_n}{n} \log(d/\phi)} < c \end{equation} then the inequality \begin{equation} \max_{j\in [d]} \vert \hat{\xi}_j - \xi_j \vert \leq C \frac{\bar{\lambda}^{1-q/2}}{\underline{\lambda}} \sqrt{M \frac{\bar{\gamma}}{ \underline{\gamma}}\frac{\kappa_n}{n} \log(d/\phi) } \end{equation} holds with probability greater than $1- C\phi$. (ii) Suppose that (ref) holds. The inequality \begin{equation} \max_{i,i^\prime \in[n]} \vert \tilde{\lambda}_{i,i^\prime} - \lambda_{i,i^\prime} \vert \leq C \frac{1}{d} \|\Omega_d\|_\text{F} \log(n/\phi) \end{equation} holds with probability greater than $1- C\phi$. \end{lemma} We apply (ref) to bound each of the terms in the decomposition \begin{equation} \widehat{\Lambda}_n -\Lambda_n= (\widehat{\Lambda}_{n,d} - \tilde{\Lambda}_{n,d} ) - (\tilde{\Lambda}_{n,d} - \Lambda_n) . \end{equation} We begin by considering the first term. Observe that \begin{align} \max_{i,i^\prime \in [n] } \vert \hat{\gamma}_0 \hat{\lambda}_{i,i^\prime} - \tilde{\lambda}_{i,i^\prime} \vert & = \max_{i,i^\prime \in [n] } \vert \frac{1}{d} \sum_{j=1}^d (\hat{\xi}^{(j)} - \xi^{(j)} ) Y_i^{(j)} Y_{i^\prime}^{(j)} \vert \nonumber \\ & \leq \left( \max_{j\in [d]} \vert \hat{\xi}_j - \xi_j \vert \right) \left( \max_{i,i^\prime \in [n]} \vert Y_i^{(j)} Y_{i^\prime}^{(j)} \vert \right) . \end{align} The random variable $Y_i^{(j)} Y_{i^\prime}^{(j)}$ has expectation smaller than $\bar{\lambda} \bar{\gamma}$ by (ref). Moreover, it has has sub-Exponential norm smaller than $(\bar{\lambda} \bar{\gamma} M)^2$ (see e.g., Lemma 2.7.6 of vershynin2018high). Hence, it holds that \begin{equation} \max_{i,i^\prime \in [n]} \vert Y_i^{(j)} Y_{i^\prime}^{(j)} \vert \leq \bar{\lambda} \bar{\gamma} + 2( \bar{\lambda} \bar{\gamma} M)^2 \log(n/\phi) \end{equation} with probability greater than $1-\phi$. Thus, (ref), Part (i), and the condition $\varphi_{n,d}<c$ (allowing the constant $c$ to depend on the terms (ref) and $M$) imply that \begin{align} \max_{i,i^\prime \in [n]} \vert \hat{\gamma}_0 \hat{\lambda}_{i,i^\prime} - \tilde{\lambda}_{i,i^\prime} \vert & \leq C \left( \frac{\bar{\lambda}^{2-q/2}\bar{\gamma}^{3/2}}{\underline{\lambda}\underline{\gamma}^{1/2}} + \frac{\bar{\lambda}^{3-q/2}\bar{\gamma}^{5/2}}{\underline{\lambda}\underline{\gamma}^{1/2}} \right) M^{5/2} \sqrt{\frac{\kappa_n}{n}} \log^{3/2}(dn/\phi) . \end{align} with probability greater than $1-C\phi$. To control the second term in the decomposition (ref), we can apply (ref), Part (ii), directly. Putting the pieces together, we find that \begin{align} & \max_{i,i^\prime \in [n]} \vert \hat{\gamma}_0 \hat{\lambda}_{i,i^\prime} - \lambda_{i,i^\prime} \vert \nonumber\\ & \leq C \left( \frac{\bar{\lambda}^{2-q/2}\bar{\gamma}^{3/2}}{\underline{\lambda}\underline{\gamma}^{1/2}} + \frac{\bar{\lambda}^{3-q/2}\bar{\gamma}^{5/2}}{\underline{\lambda}\underline{\gamma}^{1/2}} \right) M^{5/2} \left( \frac{1}{d} \|\Omega_d\|_\text{F} + \sqrt{\frac{\kappa_n}{n}} \right)\log^{3/2}(dn/\phi) . \end{align} with probability greater than $1-C\phi$. Absorbing the terms (ref) and $M$ into $C$ gives the desired result.\qed \subsection{Proof of (ref)} Recall the definition of the matrix Frobenius norm $\|\cdot\|_{\text{F}}$. Observe that \begin{equation} v^\star = \frac{\lambda^2}{d^2} \sum_{j=1}^d \eta^2_j = \frac{\lambda^2}{d^2} \| \Omega_d\|_{\text{F}}^2. \end{equation} We are interested in establishing the bound \begin{align} &\bigg\vert \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \mathbb{I}\{\hat{\gamma}_0\hat{\lambda}_{i,i^\prime} \leq \delta \} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right) \bigg\vert \nonumber \\ & \quad\quad \leq C\left(\frac{\sum_{j=1}^d \eta^3_j}{(\sum_{j=1}^d \eta_j^2)^{3/2}} + \sqrt{\frac{\log^{3}(dn/\phi)}{n}}\right) . \end{align} To ease exposition, we provide the details for the proof of the upper bound encoded in (ref). The analogous lower bound will follow from a similar argument. Recall the definition of the infeasible estimator $\tilde{\lambda}_{i,i^\prime}$ introduced in (ref). Observe that, on the event that \begin{equation} \max_{i,i^\prime \in [n]} \vert \hat{\gamma}_{0}\hat{\lambda}_{i,i^\prime} - \tilde{\lambda}_{i,i^\prime} \vert \leq t , \end{equation} we have that \begin{align} & \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \mathbb{I}\{\hat{\gamma}_{0} \hat{\lambda}_{i,i^\prime} \leq \delta\} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right) \nonumber \\ & \quad \quad\leq \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \mathbb{I}\{\tilde{\lambda}_{i,i^\prime} \leq \delta + t\} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta + t\right) \nonumber\\ &\quad \quad \quad + \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta + t\right) - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right) \nonumber \\ & \quad \quad\leq \left(\frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \mathbb{I}\{\tilde{\lambda}_{i,i^\prime} \leq \delta + t\} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta + t\right)\right) + t , \end{align} where the second inequality follows from the fact that the Gaussian cumulative distribution function has Lipschitz constant less than one. Thus, the result will follow by bounding the first term in (ref) and choosing a suitable value of $t$ such that the the event (ref) holds with high probability. To this end, we obtain a bound for the term (ref) from the following Lemma. \begin{lemma} Suppose that (ref) holds, that $\Lambda_n = \lambda I_n$ for some constant $\lambda$, and that the data $Y_i^{(j)}$ are Gaussian. Let $\eta_1,\ldots,\eta_d$ denote the eigenvalues of $\Omega_d$. Fix a constant $0<\phi<1$. It holds that \begin{equation} \bigg\vert \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \mathbb{I}\{\tilde{\lambda}_{i,i^\prime} \leq \delta \} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right) \bigg\vert \leq C\left(\frac{\sum_{j=1}^d \eta^3_j}{(\sum_{j=1}^d \eta_j^2)^{3/2}} + \frac{\log(1/\phi)}{\sqrt{n}}\right) \end{equation} with probability greater than $1-\phi$. \end{lemma} Moreover, observe that as $\Lambda_n = \lambda I_n$ and the sub-Gaussian norm of a standard Gaussian random variable is less than 2, the inequality (ref) implies that, by setting \begin{equation} t = C \left( \frac{\lambda^{1-q/2}\bar{\gamma}^{3/2}}{\underline{\gamma}^{1/2}} + \frac{\lambda^{2-q/2}\bar{\gamma}^{5/2}}{\underline{\gamma}^{1/2}} \right) \sqrt{\frac{\log^{3}(dn/\phi)}{n}} , \end{equation} the event (ref) holds with probability greater than $1-C\phi$. Hence, the decomposition (ref) and (ref) imply that the bound (ref) holds with probability greater than $1-C\phi$, where we have subsumed $\lambda$ and $\bar{\gamma}$ into the constant $C$.\qed \section{Proofs for Auxiliary Results} \subsection{Proof of (ref). Part (i)} The result follows from an application of the hanson1971bound inequality for sub-Gaussian quadratic forms. See e.g., rudelson2013hanson for a modern treatment. \begin{lemma}[Theorem 1.1, rudelson2013hanson] Let $X = (X_1,\ldots,X_m)$ be a random vector with independent, centered, and $M$-sub-Gaussian components. Let $A$ be an $m\times m$ deterministic matrix. There exists a constant $C>0$ such that \begin{equation} P\left\{ \vert X^\top A X - \mathbb{E}\left[ X^\top A X \right] \vert \geq t \right\} \leq 2 \exp\left( -C \min \left\{ \frac{t^2}{M^4 \|A\|_{\text{F}}^2} , \frac{t}{M^2 \| A \|_{\text{op}}} \right\} \right) \end{equation} for every $t\geq 0$. \end{lemma} In particular, we apply (ref) to the statistic $\hat{\gamma}_j$. To this end, collect the observations of the $j$th outcome into the vector $Y^{(j)} = (Y^{(j)}_i)_{i=1}^n$. Recall the matrix $Z$ defined in (ref). Observe that the matrix \begin{equation} \tilde{Z} = Z \Gamma_d^{1/2} \end{equation} has independent rows and has components whose sub-Gaussian norm is at most $M \bar{\gamma}$. Let $\tilde{Z}^{(j)}$ denote its $j$th column. Observe that we can write \begin{equation} \hat{\gamma}_j = \frac{1}{n} (Y^{(j)})^\top Y^{(j)} = \frac{1}{n} (\tilde{Z}^{(j)} )^\top \Lambda_n \tilde{Z}^{(j)} , \end{equation} by (ref). Moreover, we can evaluate \begin{equation} \mathbb{E}[\hat{\gamma}_j] = \gamma_{j,j} L_n , \quad\text{where}\quad L_n = \frac{1}{n}\sum_{i=1}^n \lambda_{i,i} . \end{equation} Thus, (ref) implies that \begin{equation} P\left\{ \vert \hat{\gamma}_j - \gamma_{j,j} L_n \vert \geq t \right\} \leq 2 \exp\left( -C \min \left\{ \frac{t^2}{\bar{\gamma}^4 M^4 \|\frac{1}{n}\Lambda_n\|_{\text{F}}^2} , \frac{t}{\bar{\gamma}^2 M^2 \|\frac{1}{n}\Lambda_n \|_{\text{op}}} \right\} \right) \end{equation} for all $t \geq 0$. Now, observe that (ref) implies that \begin{align} \| \Lambda_n \|_{\text{op}} &\leq \max_{i\in[n]} \frac{1}{n} \sum_{i^\prime=1}^n \vert \lambda_{i,i^\prime}\vert \leq \frac{\bar{\lambda}^{1-q}}{n} \kappa_n \quad\text{and} \\ \| \Lambda_n \|_{\text{F}}^2 & \leq \frac{1}{n^2} \sum_{i=1}^n \sum_{i^\prime =1}^n \lambda_{i,i^\prime}^2 \leq \bar\lambda^{2-q} \frac{1}{n^2} \sum_{i=1}^n \sum_{i^\prime =1}^n \vert \lambda_{i,i^\prime}\vert^q \leq \frac{\bar{\lambda}^{2-q}}{n} \kappa_n, \end{align} respectively. Consequently, the inequalities (ref), (ref), and (ref) imply that \begin{equation} P\left\{ \vert \hat{\gamma}_j - \gamma_{j,j} L_n \vert \geq t \right\} \leq 2 \exp\left( -\frac{C}{\bar{\gamma}^4 M^4 \bar{\lambda}^{2-q}} \frac{t^2 n}{\kappa_n} \right) \end{equation} for all $t\leq \bar{\gamma}^2 M^2 \bar{\lambda}$. Thus, so long as $n^{-1}\kappa_n \log(d/\phi) <1$, by choosing \begin{equation} t = c \bar{\lambda}^{1-q/2} \bar{\gamma}^2 M^2\sqrt{\frac{\kappa_n}{n} \log(d/\phi) } \end{equation} for a sufficiently small constant $c$, the inequality (ref) implies that \begin{equation} \max_{j \in [d]} \Big\vert \hat{\gamma}_j - \gamma_{j,j} L_n \Big\vert \leq c \bar{\lambda}^{1-q/2} \bar{\gamma}^2 M^2 \sqrt{\frac{\kappa_n}{n} \log(d/\phi) } \end{equation} with probability $1-\phi$, by a union bound. To conclude the proof, we translate our large deviation bound for $\hat{\gamma}_j$ into a large deviation bound for $\hat{\gamma}_j$. To this end, observe that, on the event \begin{equation} \max_{j\in[d]} \left\{ \bigg\vert \frac{\hat{\gamma}_j}{\gamma_{j,j}L_n}- 1\bigg\vert \right\}\leq \frac{1}{2} \end{equation} it holds that \begin{align} \bigg\vert \frac{\hat{\gamma}_1}{\hat{\gamma}_j} - \frac{1}{\gamma_{j,j}}\bigg\vert & \leq \frac{L_n}{\hat{\gamma}_j} \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert + \frac{1}{\gamma_{j,j}} \bigg\vert \frac{\gamma_{j,j} L_n}{\hat{\gamma}_j} - 1\bigg\vert \nonumber \\ & \leq 2\left( \frac{L_n}{\hat{\gamma}_j} \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert + \frac{1}{\gamma_{j,j}} \bigg\vert \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} - 1\bigg\vert\right) \nonumber \\ & \leq 2\left( \left( \frac{1}{\gamma_{j,j}} + \bigg\vert \frac{1}{\gamma_{j,j}} - \frac{L_n}{\hat{\gamma}_j} \bigg\vert \right) \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert + \frac{1}{\gamma_{j,j}} \bigg\vert \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} - 1\bigg\vert\right) \nonumber \\ & \leq 2\left( \left( \frac{1}{\gamma_{j,j}} + \frac{1}{\gamma_{j,j}} \bigg\vert 1 - \frac{\gamma_{j,j} L_n}{\hat{\gamma}_j} \bigg\vert \right) \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert + \frac{1}{\gamma_{j,j}} \bigg\vert \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} - 1\bigg\vert\right) \nonumber \\ & \leq 4\left( \left( \frac{1}{\gamma_{j,j}} + \frac{1}{\gamma_{j,j}} \bigg\vert 1 - \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} \bigg\vert \right) \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert + \frac{1}{\gamma_{j,j}} \bigg\vert \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} - 1\bigg\vert\right) \nonumber \\ & \leq \frac{4}{\gamma_{j,j}}\left( \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert + \bigg\vert 1 - \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} \bigg\vert + \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert \bigg\vert 1 - \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} \bigg\vert \right) , \end{align} where we have applied the elementary inequality \begin{equation} \vert z^{-1} - 1 \vert \leq 2\vert z-1\vert ,\quad\text{for}\quad\vert z-1\vert\leq{1/2} , \end{equation} to obtain the second and fifth inequalities. Now, observe that \begin{align} \bigg\vert \frac{\hat{\gamma}_j}{\gamma_{j,j}L_n}- 1\bigg\vert & \leq c \frac{ \bar{\lambda}^{1-q/2} \bar{\gamma}^2}{\gamma_{j,j} L_n} \sqrt{M \frac{\kappa_n}{n} \log(d\phi) } \nonumber \\ & \leq c \frac{\bar{\lambda}^{1-q/2}}{\underline{\lambda}} \frac{\bar{\gamma}^2}{\underline{\gamma}} \sqrt{M \frac{\kappa_n}{n} \log(d/\phi) } \end{align} with probability greater than $1-\phi/d$, by the inequality (ref). The condition (ref) implies that \begin{equation} c \frac{\bar{\lambda}^{1-q/2}}{\underline{\lambda}} \frac{\bar{\gamma}^2}{\underline{\gamma}} \sqrt{M \frac{\kappa_n}{n} \log(d/\phi) } \leq \frac{1}{2} \end{equation} with probability greater than $1-\phi/d$. Consequently, the union bound implies that the event (ref) holds for every $j$ in $[d]$ with probability greater than $1-\phi$. Thus, the inequality (ref) implies that \begin{equation} \max_{j\in [d]} \vert \hat{\xi}_j - \xi_j \vert \leq \max_{j\in [d]} \left\{C{\gamma_{j,j}}\left( \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert + \bigg\vert 1 - \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} \bigg\vert + \bigg\vert \frac{\hat{\gamma}_1}{L_n} - 1\bigg\vert \bigg\vert 1 - \frac{\hat{\gamma}_j}{\gamma_{j,j} L_n} \bigg\vert \right)\right\} , \end{equation} with probability greater than $1-\phi$. Appealing, again, to the fact that the event (ref) holds for every $j$ in $[d]$ with probability greater than $1-\phi$, and applying the bound (ref), we can conclude that the inequality \begin{equation} \max_{j\in [d]} \vert \hat{\xi}_j - \xi_j \vert \leq C \frac{\bar{\lambda}^{1-q/2}}{\underline{\lambda}} \sqrt{M \frac{\bar{\gamma}}{ \underline{\gamma}}\frac{\kappa_n}{n} \log(d/\phi) } \end{equation} holds with probability greater than $1-C\phi$, as desired.\qed \subsection{Proof of (ref). Part (ii)} The result again follows from an application of (ref). Recall the matrix $Z$ defined in (ref). Observe that the matrix \begin{equation} \widecheck{Z} = \Lambda_n^{1/2} Z \end{equation} has independent columns and has components whose sub-Gaussian norm is at most $M \bar{\lambda}$. Let $\widecheck{Z}_i$ denote its $i$th column. By (ref), we have that \begin{equation} \tilde{\lambda}_{i,i^\prime} = \frac{1}{d} Y_i \Xi_d Y_i^\top = \frac{1}{d} \widecheck{Z}_i \Lambda_d^{1/2} \Xi_d \Lambda_d^{1/2} \widecheck{Z}_i^\top . \end{equation} Thus, (ref) and a union bound imply that \begin{equation} P\left\{ \max_{i,i^\prime \in [n]} \vert \tilde{\lambda}_{i,i^\prime} - \lambda_{i,i^\prime} \vert \geq t \right\} \leq 2 n^2 \exp\left( -C t d \min\left\{\frac{ td}{ \|\Omega_d\|^2_{\text{F}}}, \frac{1}{\|\Omega_d\|_{\text{op}}}\right\} \right) . \end{equation} Re-expressing this bound, we find that the inequality \begin{align} \max_{i,i^\prime \in [n]} \vert \tilde{\lambda}_{i,i^\prime} - \lambda_{i,i^\prime} \vert & \leq C \frac{1}{d} \left(\|\Omega_d\|_\text{F} \sqrt{\log(n/\phi)} + \|\Omega_d\|_\text{op} \log(n/\phi) \right) \nonumber \\ & \leq C \frac{1}{d} \|\Omega_d\|_\text{F} \log(n/\phi) \end{align} holds with probability greater than $1-C\phi$, as desired.\qed \subsection{Proof of (ref)} Throughout, we let the operator $\sim$ denote equality in distribution. The result relies on the following distributional decomposition of Gaussian quadratic forms. Let $\mathcal{W}_q(\Sigma,p)$ denote the $q$-variate Wishart distribution, with parameter $\Sigma$ and degrees of freedom $p$. \begin{lemma}[Corollary 2.3, singull2012distribution] Fix an $n\times n$ positive semi-definite matrix $\Lambda_n = (\lambda_{i,i^\prime})_{i,i^\prime=1}^n$ and a $d \times d$ full-rank matrix $\Gamma_d = (\gamma_{j,j^\prime})_{j,j^\prime=1}^d$. Let $X$ denote a mean-zero, Gaussian random matrix, whose $(i,j)$th component is given by $X_i^{(j)}$ and that satisfies \begin{equation} \Cov(X_i^{(j)}, X_{i^\prime}^{(j^\prime)}) = \lambda_{i,i^\prime} \gamma_{j,j^\prime} . \end{equation} Let $A$ denote a deterministic, symmetric, full-rank $d \times d$ matrix. The statistic $Q = X A X^\top$ has the distribution \begin{equation} Q \sim \sum_{j=1}^d \eta_j S_j , \end{equation} where $\eta_1,\ldots,\eta_d$ are the eigenvalues of $\Gamma^{1/2}_d A \Gamma^{1/2}_d$ and $S_1,\ldots,S_d$ are independent random matrices with distribution $\mathcal{W}_d(\Lambda_n,1)$. \end{lemma} As the data are Gaussian, and satisfy (ref), (ref) implies that \begin{equation} \tilde{\Lambda}_{n,d} = \frac{1}{d} Y \Xi_d Y^\top \sim \frac{1}{d} \sum_{j=1}^d \eta_j S_j , \end{equation} where $\eta_1,\ldots,\eta_d$ are the eigenvalues of $\Omega_d = \Gamma^{1/2}_d \Xi_d \Gamma^{1/2}_d$ and $S_1,\ldots,S_d$ are independent random matrices with distribution $\mathcal{W}_d(\Lambda_n,1)$. Let $Z$ denote a mean-zero, Gaussian random matrix whose components are mutually independent. Let $Z^{(j)}$ denote the $j$th column of $Z$ and set $\Var(Z^{(j)}) = \lambda I_n$. Likewise, let $Z_i^{(j)}$ and $Z_i$ denote the $(i,j)$th component and $i$th row of $Z$, respectively. By the definition of the Wishart distribution, it holds that \begin{equation} S_j \sim Z^{(j)} (Z^{(j)})^\top \end{equation} for each $j$ in $[d]$. Thus, the characterization (ref) implies that \begin{equation} \tilde{\lambda}_{i,i^\prime} \sim \frac{1}{d} \sum_{j=1}^d \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} \end{equation} for each $i,i^\prime$ in $[n]$. Therefore, the statistic \begin{equation} \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \left(\mathbb{I}\{\tilde{\lambda}_{i,i^\prime} \leq \delta \} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right)\right) \end{equation} is equi-distributed with \begin{align} \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \psi_{\delta}(Z_i, Z_{i^\prime}) \end{align} where \begin{equation} \psi_{\delta}(Z_i, Z_{i^\prime}) = \mathbb{I}\left\{\frac{1}{d} \sum_{j=1}^d \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} \leq \delta \right\} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right) . \end{equation} The random variable (ref) can be recognized as a $U$-statistic of order 2. Thus, we apply the following large-deviation bound, due to hoeffding1963probability. \begin{lemma}[Lemma A.5, song2019approximating] Let $X_1,\ldots,X_n$ denote a collection of independent and identically distributed random variables. Suppose that $\psi(\cdot,\cdot)$ is a real-valued function that is symmetric in its two components and satisfies \begin{equation} \mathbb{E}\left[\psi(X_1,X_2)\right] = 0 , \quad \Var(\psi(X_1,X_2)) = \nu , \quad\text{and}\quad \vert \psi(X_1,X_2) \vert \leq \theta \end{equation} almost surely. Fix a constant $0<\phi<1$. The inequality \begin{align} \bigg\vert \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \psi(X_i, X_{i^\prime}) \bigg\vert \leq C \left(\sqrt{\frac{\nu \log(1/\phi)}{n}} + \frac{\theta \log(1/\phi)}{n}\right) \end{align} holds with probability greater than $1-\phi$. \end{lemma} Observe that, as $\Lambda_n = \lambda I_n$ and the components of $Z$ are independent, the random variables $Z_1,\ldots,Z_n$ are independent and identically distributed. Moreover, we have that \begin{equation} \Var(\psi_{\delta}(Z_i, Z_{i^\prime})) \leq 1 \quad\text{and}\quad \vert\psi_{\delta}(Z_i, Z_{i^\prime}) - \mathbb{E}\left[\psi_{\delta}(Z_i, Z_{i^\prime})\right]\vert \leq 1 \end{equation} almost surely. Thus, (ref) implies that the inequality \begin{equation} \bigg\vert \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} (\psi_{\delta}(Z_i, Z_{i^\prime}) - \mathbb{E}\left[\psi_{\delta}(Z_i, Z_{i^\prime}) \right]) \bigg\vert \leq C \frac{\log(1/\phi)}{\sqrt{n}} \end{equation} holds with probability greater than $1-\phi$. Hence, it suffices to bound the expectation \begin{equation} \mathbb{E}\left[\psi_{\delta}(Z_i, Z_{i^\prime}) \right] = P\left\{\frac{1}{d} \sum_{j=1}^d \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} \leq \delta \right\} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right) . \end{equation} To do this, we apply the following Berry-Esseen inequality for non-identically distributed sums \begin{lemma}[Theorem XVI.5.2, feller1991introduction] Let $X_1,\ldots,X_n$ be independent variables such that \begin{equation} \mathbb{E}\left[X_k\right] = 0,\quad \mathbb{E}\left[\vert X_k\vert^2 \right] = \sigma_k^2,\quad\text{and}\quad \mathbb{E}\left[\vert X_k\vert^3 \right] = \theta^k . \end{equation} Define the sequences \begin{equation} s_n^2 = \sum_{i=1}^n \sigma_i^2 \quad\text{and}\quad r_n = \sum_{i=1}^n \theta_i \end{equation} For all $\delta$ and $n$, the inequality \begin{equation} \bigg\vert P\left\{ \frac{1}{s_n} \sum_{i=1}^n X_k \leq \delta \right\} - \Phi\left(\delta\right)\bigg\vert \leq C\frac{r_n}{s_n^3} \end{equation} holds. \end{lemma} Consider the sum \begin{equation} \sum_{j=1}^d \frac{1}{d} \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} . \end{equation} Observe that \begin{align} \mathbb{E}\left[ \frac{1}{d} \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)}\right] &= 0 ,\nonumber \\ \mathbb{E}\left[ \vert \frac{1}{d} \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} \vert ^2\right] &= \frac{\eta^2_j}{d^2} \lambda^2 ,\quad\text{and}\nonumber \\ \mathbb{E}\left[ \vert \frac{1}{d} \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} \vert ^3\right] & \geq C \frac{\eta^3_j}{d^3} \lambda^3 , \end{align} where we have used facts about the second and third moments of Gaussian distributions for the second and third relations. Consequently, (ref) implies that \begin{align} & \bigg\vert P\left\{\frac{1}{d} \sum_{j=1}^d \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} \leq \delta \right\} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right) \bigg\vert \nonumber \\ & \leq \bigg\vert P\left\{ \frac{d}{\lambda \sqrt{\sum_{j=1}^d \eta_j^2}} \frac{1}{d}\sum_{j=1}^d \eta_j Z_i^{(j)} Z_{i^\prime}^{(j)} \leq \frac{d}{\lambda \sqrt{\sum_{j=1}^d \eta_j^2}} \delta \right\} - \Phi\left(\frac{d}{\lambda \sqrt{\sum_{j=1}^d \eta_j^2}} \delta\right) \bigg\vert \nonumber \\ & \leq \frac{\sum_{j=1}^d \eta^3_j}{(\sum_{j=1}^d \eta_j^2)^{3/2}} . \end{align} Hence, by combining the inequalities (ref) and (ref), we can conclude that \begin{equation} \bigg\vert \frac{2}{n(n-1)} \sum_{i=1}^n \sum_{i^\prime < i} \left(\mathbb{I}\{\tilde{\lambda}_{i,i^\prime} \leq \delta \} - \Phi\left(\frac{d}{\lambda}\frac{1}{\|\Omega_d\|_{\text{F}}} \delta\right)\right) \bigg\vert \leq C\left(\frac{\sum_{j=1}^d \eta^3_j}{(\sum_{j=1}^d \eta_j^2)^{3/2}} + \frac{\log(1/\phi)}{\sqrt{n}}\right) \end{equation} holds with probability greater than $1-\phi$, as required.\qed \section{Further Details for Empirical Applications} \subsection{bazzi2023other} bazzi2023other examine the impact of the migration of millions of Southern whites in the twentieth century. They use cross-sectional data concerning the demographic and political characteristics of 1,888 U.S.\ counties. The main finding is that the Southern white migration in the early twentieth century is associated with significantly higher Republican vote shares in the twenty-first century. bazzi2023other consider the model \begin{equation} \text{Vote}_c =\beta \cdot \text{Southern Whites}_{c,1940} + X_{c} + \alpha_s + \varepsilon_{c} , \end{equation} where $c$ indexes counties and $s$ indexes states. The main outcome, $\text{Vote}_c$, is the vote share received by Donald Trump in 2016. The regressor of interest, $\text{Southern Whites}_{c,1940}$, is the Southern-born white population share in 1940. The regression is instrumented by a shift-share variable constructed with predetermined Southern white migration networks as of 1900 and predicted aggregate migration flows out of the South for each decade from 1900 to 1940. The control variables $X_{c}$ are historical economic factors (population density, manufacturing employment, average farm values), ideological factors (Union Army enlistment, mortality rates from the U.S. Civil War), and the vote share for Woodrow Wilson in 1912. The variable $\alpha_s$ denotes a state fixed effect. Standard errors are clustered across counties in $60 \times 60$ mile grid cells. The main result are displayed in their Table 2, Panel A, Column (4).\footnote{The authors do not explicitly indicate the preferred specification. However, on page 1591, they write “Our main estimating equation takes the following form." Thus, we infer that Table 2 displays their “main” results. On Page 1595, they write "Secondary specifications control for additional potential sorting correlates," and so columns 5 and 7 are not the main specifications. Therefore, we use Panel A column 4 as the main specification.} The authors implement several adjustments for spatial correlation. They report alternative standard errors, constructed with the methods proposed by conley1999gmm, colella2023acreg, and adao2019shift, as well as through a wild bootstrap with clustering by state. The corrections of the main specification are shown in their Column (4) of Table A.3. The standard error based on the conley1999gmm, with a cutoff at 500 km, is the largest and is 68% larger than the baseline result. For the set of auxiliary outcomes in the TMO procedure, we collect all feasible variables available in the replication package. We exclude outcomes that are highly correlated with the main outcome, the regressor of interest, or the primary controls variables. We also exclude outcomes that are missing for over 50% of observations or are collinear with right hand side of the model. The auxiliary outcomes are all from the replication package and include variables such as the Republican vote share in earlier years and the percent of Southern blacks. \subsection{bernini2023race} bernini2023race study the impact of the Voting Rights Act (VRA) on the racial makeup of local governments in the U.S.\ South. They consider cross-sectional data consisting of 971 counties in 11 states. They find that federal scrutiny over Southern states post VRA had a sizable impact on the extent to which enfranchisement led to Black office holding. bernini2023race consider the long differences regression model \begin{align*} \Delta \text{Share Black Elected}_{cs} &=\gamma \text{Percent Black}_{1960} + \theta \text{Percent Black}_{1960}*\text{Covered}_{cs} \\ &+ X_{cs}+X_{cs}*\text{Covered}_{cs} + I_s + \epsilon_{cs} , \end{align*} where $c$ denote county and $s$ denotes state. The main outcome, $\Delta \text{Share Black Elected}_{cs}$, is the change in the share of Black elected officials between 1964 and 1980. The regressor of interest, $\text{Percent Black}_{1960}*\text{Covered}_{cs}$, is the interaction of the Black share in 1960 and a dummy for federal scrutiny. Control variables $X_{cs}$ are pre-VRA county characteristics (unemployment rate, the percent of families below the poverty line, the percent of unskilled workers, agricultural productivity, population (ln), percent urban population, cotton share, pro- and anti-Black activism, Republican share), $I_s$ is state dummies. Standard errors are clustered by judicial division. Results are shown in Table 2 Column 4, which is stated as the paper's preferred specification\footnote{On page 1019-1020, bernini2023race write “Our preferred specification in column 4 of table 2"}. For the set of auxiliary outcomes in the TMO procedure, we collect all feasible variables available in the replication package. We exclude outcomes that are highly correlated with the main outcome, the regressor of interest, or the primary controls variables. We also exclude outcomes that are missing for over 50% of observations or are collinear with right hand side of the model. The auxiliary outcomes are all from the replication package and include variables such as the change in the log turnout for governor elections from 1940 to 1960 and the change in NAACP branches from 1942 to 1964. \subsection{caprettini2023new} caprettini2023new study the complementarity between patriotism in World War II and public-good provision in the New Deal. They consider cross-sectional data considering 2329 counties. They find that higher government spending at the county level in the 1930s is positively correlated with more patriotic actions during World War II. caprettini2023new consider the OLS regression model \begin{equation*} \text{WWII Patriotism}_i = \beta \text{New Deal Grants}_i + \gamma \text{WWI Patriotism}_i+\delta X_i +\xi_s +u_i , \end{equation*} where $i$ denotes county. The main outcome, $\text{WWII Patriotism}_i$ denotes the average purchases of war bonds per capita in 1944. The regressor of interest, $\text{New Deal Grants}_i$ denotes the New Deal grants per 1930 population. Control variables $\text{WWI Patriotism}_i$ and $X_i$ are WWI volunteering rate, WWI medals per 1000 inhabitants, socioeconomic characteristics in 1930 (log of population, unemployment rate, share of veterans, share of African Americans, and share of people born in major Axis countries), share of people living on a farm in 1930, 1929 farm income, an urban dummy, 1898-1928 average Democratic vote share, value of WWII war contracts per capita, average 1939 wage of employees, and 1930 share of men. The variables $\xi_s$ denote state fixed effects. Heteroskedasticity-robust standard errors are used. Results are shown in Table 2 Panel A Column 1.\footnote{On Page 481, caprettini2023new write “Figure I summarizes our main result". Similarly, on Page 483 they write “to go beyond the graphical evidence", suggest results in this section (Table 2) are the main results. The first three columns report results for three equally important measures of patriotism, so we randomly choose the first as the main specification} caprettini2023new discuss spatial error correlation. They report conley1999gmm standard errors in Table 2 Panel B and Table 3. The conley1999gmm standard errors are between 0% and 171% larger than the heteroskedasticity-robust standard errors. For the set of auxiliary outcomes in the TMO procedure, we collect all feasible variables available in the replication package. We exclude outcomes that are highly correlated with the main outcome, the regressor of interest, and the primary controls variables. We also exclude outcomes that are missing for over 50% of observations or are collinear with right hand side of the model. The auxiliary outcomes are all from the replication package and include variables such as the number of WW2 volunteers per 100 people and quantiles of 1940 FDR vote share. \subsection{moscona2023does} moscona2023does study how innovation shapes the economic impact of climate change. They consider panel data consisting of 3000 counties between 1950 and 2010. They find that that counties' exposure to climate-change-induced innovation significantly decreases the local economic damage from extreme temperatures. moscona2023does consider the regression model \begin{align*} \text{log Agr Land Price}_{i,t} &= \delta_i + \alpha_{s(i),t}+\beta \text{Extreme Exposure}_{i,t}+\gamma \text{Innovation Exposure}_{i,t} \\ &+\phi (\text{Extreme Exposure}_{i,t}* \text{Innovation Exposure}_{i,t}) + \epsilon_{i,t} \end{align*} where $i$ denotes county, $s$ denotes state, and $t$ denotes time (1959 or 2017). The main outcome, $\text{log Agr}$ $\text{Land Price}_{i,t}$, is the price per acre of cultivated land. The regressor of interest, $\text{Extreme Exposure}_{i,t}*\text{Innovation Exposure}_{i,t}$, is the interaction of measures of a county's exposures to extreme temperature and innovation, respectively. The variable $\delta_i$ denote county fixed effects. The variables $\alpha_{s(i),t}$ denote state by year fixed effects. Standard errors are double clustered at the county and state-by-year levels. Results are shown in Table 3 Column 1 in the paper.\footnote{The authors do not explicitly indicate the preferred specification. However, on page 645, they write “Sections IV and V present our main results." Thus, we conclude that Table 3 displays the paper's main county level results. Within Table 3, we pick Column 1 as the main specification, because the authors choose this specification to perform down-stream analysis. For example, on Page 680, they write “Figure VI reports the marginal effect of exposure to extreme heat (y-axis) for quantiles of the innovation exposure distribution (x-axis), using the specification from column (1)".} moscona2023does report conley1999gmm standard errors in Table A21. The conley1999gmm standard errors are smaller than the clustered standard errors reported in the main text. As there are a limited number of variables in the replication package, we construct the set of auxiliary outcomes using external sources. Given the agricultural context of this study, we take county-level outcomes from the United States Censuses of Agriculture compiled in hainesetal2018icpsr, such as the number of farms. We also collect outcomes that are in both the historical County Data Book dataset\footnote{\url{https://doi.org/10.3886/ICPSR07736.v2}} and the contemporary USA Counties dataset\footnote{\url{https://www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}} such as employment in the agricultural sector. We use the outcomes for which there are data in both 1950 and 2010, which are the years of interest in this study. \subsection{esposito2023reconciliation} esposito2023reconciliation study the effect of the Lost Cause narrative---a revisionist and racist retelling of the U.S.\ Civil War---on national reunification and racial discrimination. They consider a panel dataset consisting of 786 counties between 1910 to 1920. They find that screenings of \textit{The Birth of a Nation} shifted public discourse toward more patriotic and less divisive language, increased military enlistment, and fostered cultural convergence. This paper considers an instrumental variables regression model \begin{equation*} \text{Reconciliation}_{ct}=\beta \text{BON}_{ct} + X_{ct}+ \alpha_c +\alpha_t+\epsilon_{ct} \end{equation*} where $c$ is county and $t$ is year-month. The main outcome, $\text{Reconciliation}_{ct}$, represents the relative log frequencies of patriotic keywords to divisive ones in local newspapers. The regressor of interest, $\text{BON}_{ct}$, is equal to 1 after the screening of \textit{The Birth of a Nation}. The regression is instrumented with the screening of an another movie, \textit{The Million Dollar Mystery}. The control, $X_{ct}$, is total number of newspaper pages available in data. Standard errors are clustered at the county level. Results are shown in Table 3 Panel A Column 3.\footnote{On page 1479, the authors write “In our baseline analysis...." Likewise, on page 1479, they write “we compute our main dependent variable, $\text{Reconciliation}_{ct}$."} We choose the specification reported in Panel A Column 3 because the authors specify that this specification is their preferred model.\footnote{On page 1472, esposito2023reconciliation write “our preferred solution is to instrument the treatment in a two-stage least squares (2SLS) version of equation (1)."} The paper also clusters the standard errors at the state level and applies the estimator proposed in colella2023acreg (see Table B34 Panel C Column 1). Clustering at the state level increases the standard error by 34%, while the correction based on colella2023acreg decreases the standard error. For the set of auxiliary outcomes in the TMO procedure, we begin with all feasible variables available in the replication package, such as whether "White-only" is observed in a newspaper job advertisement in a given month. With only the variables in the replication package, however, the set of auxiliary outcomes does not have enough power for the TMO procedure. We therefore supplement the set of auxiliary outcomes using general economic and demographic variables aggregated to the county level from the U.S. Cenus Bureau via IPUMS\footnote{\url{https://usa.ipums.org/usa/}} such as the proportion of foreign-born population in the county. For years between 1910 and 1920, we interpolate the Census outcomes. We exclude outcomes that are highly correlated with the main outcome, the regressor of interest, and the primary controls variables. We also exclude outcomes that are missing for over 50% of observations or are collinear with right hand side of the model. \subsection{calderon2023racial} calderon2023racial study the political effects of the migration of four million African Americans from the South to the North. They consider with a panel data consisting of 1263 non-southern counties from 1940 to 1960. calderon2023racial find that the great migration increased support for the Democratic Party, increased Congress members' propensity to promote civil rights legislation, and encouraged pro-civil rights activism outside the US South. This paper considers the instrumental variables regression model \begin{equation*} \Delta \text{demsh}_{c \tau}=\delta_{s \tau} + \beta \Delta \text{Bl}_{c \tau} + \gamma X_{c \tau} +u_{c \tau} , \end{equation*} where $c$ is county and $\tau$ is decade. The main outcome, $\Delta \text{demsh}_{c \tau}$, is the change in the Democratic vote share during decade $\tau$. The regressor of interest, $\Delta \text{Bl}_{c \tau}$, is the change in the Black population share. The regression is instrumented with a shift-share predictor of Black inflow. Control variables $X_{c \tau}$ are interactions between decade dummies and 1940 county characteristics (Black population share and Democratic incumbency in Congressional elections). The variable $\delta_{s \tau}$ collects interactions between decade and state dummies. The regression is weighed by 1940 county population. Standard errors are clustered at the county level. Results are displayed in Table 2 Column 6. The authors note that this is the paper's preferred specification \footnote{On page 179, calderon2023racial write "especially for our preferred specification (Column 6)"}. calderon2023racial additionally clusters their standard errors at the commuting zone level and reports the confidence interval constructed with the method of adao2019shift. These results are displayed in Table D.21 and D.22. Clustering at the commuting zone level increases the standard errors by 1%. For the set of auxiliary outcomes in the TMO procedure, we begin with all feasible variables available in the replication package, such as the share of manufacturing. With only the variables in the replication package, however, the set of auxiliary outcomes does not have enough power for the TMO procedure. We therefore supplement the set of auxiliary outcomes using outcomes from the other historical county-level papers in (ref), such as total manufacturing wages. We also exclude outcomes that are missing for over 50% of observations, are collinear with right hand side of the model, or are outside the time period in the paper. \subsection{cook2023evolution} cook2023evolution study factors correlated with the provision of nondiscriminatory services. They consider a panel dataset consisting of 3092 counties between 1939 to 1955. They find that declines white population led to increases in the number of nondiscriminatory businesses. This paper considers the linear regression model \begin{equation*} \text{Asinh}(N_{0ct})=\beta_0 +\beta_1 \text{Asinh}(\text{casualties}_c) * \text{postWWII}_t +\phi_c +\xi_t +\epsilon_{ct} \end{equation*} where $c$ denotes county and $t$ denotes year. The main outcome, $N_{0ct}$, is the number of Green Book establishments in county $c$ at time $t$. The regressor of interest, $\text{Asinh}(\text{casualties}_c) * \text{postWWII}_t$, is the interaction of the number of white casualties in World War II and a post-war indicator. The variables $\phi_c$ and $\xi_t$ denote county and year fixed effects. Standard errors are clustered at the county level. Results are shown in Table 3 Panel A Column 4 in the paper, which is stated as the paper's preferred specification.\footnote{On page 77, the authors write “our preferred specification uses county and year fixed effects".} For the set of auxiliary outcomes in the TMO procedure, we begin with all feasible variables available in the replication package, such as the number of other establishments including barber shops, restaurants, and hotels. With only the variables in the replication package, however, the set of auxiliary outcomes does not have enough power for the TMO procedure. We therefore supplement the set of auxiliary outcomes using outcomes from the historical County Data Book dataset\footnote{\url{https://doi.org/10.3886/ICPSR07736.v2}} that are available from 1939 to 1955, such as employment in services and retail sales. We interpolate the County Data Book outcomes between the survey years. \subsection{chetty2014land} chetty2014land study features of inter-generational mobility in the United States. One of the findings of this paper is that the absolute upward mobility (the expected income rank of children from families at the 25 percentile of the national parent income distribution) is negatively correlated with the proportion in a commuting zone that is African American. This paper considers the linear regression model \begin{equation*} \text{absmob}_c=\beta_0 +\beta_1 \text{frac black}_c +\delta_s + \epsilon_{c} \end{equation*} where $s$ denotes state and $c$ denotes commuting zone. The main outcome, $\text{absmob}_c$, is the measure of absolute upward mobility in commuting zone $c$. The regressor of interest, $\text{frac black}_c$ is the share of African American people in commuting zone $c$. Standard errors are clustered at the state level. Results are shown in the first row of Figure VIII in the paper\footnote{We pick this specification because the authors write “Perhaps the most obvious pattern from the maps in Figure VI is that inter-generational mobility is lower in areas with larger African American populations, such as the Southeast" on Page 1605.}. For the set of auxiliary outcomes in the TMO procedure, we begin with all feasible variables available in the replication package, such as the CZ-level Gini coefficient and labor force participation rate. With only the variables in the replication package, however, the set of auxiliary outcomes does not have enough power for the TMO procedure. We therefore supplement the set of auxiliary outcomes using contemporary county-level outcomes from the USA Counties dataset\footnote{\url{https://www.census.gov/library/publications/2011/compendia/usa-counties-2011.html}} aggregated to the CZ level, such as the enrollment in public schools and the divorce rate. \subsection{acemoglu2019democracy} acemoglu2019democracy estimates the impact of democracy on economic growth. They consider panel data of 175 countries between 1960 to 2010. They find that democratization increases GDP per capita in the long run. This paper considers the linear regression model \begin{equation*} \text{lgdppc}_{c,t} = \beta \text{demo}_{c,t}+\sum_{j=1}^4 \gamma_j \text{lgdppc}_{c,t-j}+\alpha_c+\delta_t+\epsilon_{c,t} \end{equation*} where $c$ denotes country and $t$ denotes year. Results are shown in Table 2 Column 3, which is stated as the paper's preferred specification\footnote{On page 59, acemoglu2019democracy write “Column 3, which is our preferred specification, includes four lags of GDP per capita".}. The main outcome, $\text{lgdppc}_{c,t}$, is the log of GDP per capita. The regressor of interest, $\text{demo}_{c,t}$, is a dichotomous measure of democracy. Control variables are the lags of log GDP per capita. The quantities $\delta_t$ and $\alpha_c$ are year and country fixed effects. Standard errors are clustered at the country level. For the set of auxiliary outcomes in the TMO procedure, we begin with all feasible variables available in the replication package, such as the percentage of population with at most primary education and an index measure of market reforms. With only the variables in the replication package, however, the set of auxiliary outcomes does not have enough power for the TMO procedure. We therefore supplement the set of auxiliary outcomes using publicly available sources such as the UN and the World Bank, including variables such as the Gini index and central government debt as a percentage of GDP. \end{spacing} \newgeometry{top=0.75in, bottom=0.75in} \begin{figure}[p] \caption{Distribution of Correlations Between Locations in Applications } \begin{tabular}{c} \textit{cook2023evolution} \\ \\ \textit{caprettini2023new} \\ \\ \textit{esposito2023reconciliation} \\ \end{tabular} \justifying \end{figure} \begin{figure}[p]\ContinuedFloat \caption{Distribution of Correlations Between Locations in Applications} \begin{tabular}{c} \textit{bernini2023race} \\ \\ \textit{bazzi2023other} \\ \\ \textit{calderon2023racial} \\ \\ \end{tabular} \justifying \end{figure} \begin{figure}[p]\ContinuedFloat \caption{Distribution of Correlations Between Locations in Applications} \begin{tabular}{c} \textit{moscona2023does} \\ \\ \textit{chetty2014land} \\ \\ \textit{acemoglu2019democracy} \\ \end{tabular} \justifying { Notes: (ref) displays a histogram of the Fisher-transformed correlations $\tilde{\rho}_{i,i^\prime}$ between each pair of locations in each of the empirical applications in (ref). The density curve denotes the null density estimate, obtained with the procedure outlined in (ref). The vertical line marks the optimal threshold $\hat{\delta}^{*}$. } \end{figure} \begin{figure}[p] \caption{TMO Relative to Original Standard Error across Thresholds in Applications } \begin{tabular}{c} \textit{cook2023evolution} \\ \\ \textit{caprettini2023new} \\ \\ \textit{esposito2023reconciliation} \\ \end{tabular} \justifying \end{figure} \begin{figure}[p]\ContinuedFloat \caption{TMO Relative to Original Standard Error across Thresholds in Applications} \begin{tabular}{c} \textit{bernini2023race} \\ \\ \textit{bazzi2023other} \\ \\ \textit{calderon2023racial} \\ \\ \end{tabular} \justifying \end{figure} \begin{figure}[p]\ContinuedFloat \caption{TMO Relative to Original Standard Error across Thresholds in Applications} \begin{tabular}{c} \textit{moscona2023does} \\ \\ \textit{chetty2014land} \\ \\ \textit{acemoglu2019democracy} \\ \end{tabular} \justifying { Notes: (ref) plots the ratio of the TMO standard error relative to the original standard error in the papers as the threshold $\delta$ varies. The vertical line marks the estimate of the optimal threshold $\hat{\delta}^{*}$. } \end{figure} \restoregeometry