Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
145,424 characters · 21 sections · 46 citation commands
Assumption-Lean Shrinkage and Model Averaging for Spatial Parameters
\begingroup \onehalfspacing
\endgroup \onehalfspacing
When economic decisions rest on thousands of noisy estimates, how much should each estimate borrow from its neighbors, and which neighbors should count? The same issue arises for noisy estimates at the school kaneEstimatingTeacherImpacts2008, chettyMeasuringImpactsTeachersI2014, hospital dimickRankingHospitalsReliability2010, hullEstimatingHospitalQuality2020, firm klineSystemicDiscrimination2022, or small-area fayEstimatesIncomeSmall1979 level whenever researchers have several plausible ways to define related units. The Opportunity Atlas chettyOpportunityAtlasMapping2018 estimates of neighborhood economic mobility helped identify high-opportunity neighborhoods in the Creating Moves to Opportunity program bergmanCreatingMovesOpportunity2024, yet many of the underlying tract-level estimates are noisy due to few observations.\footnote{Title I education funding under the Elementary and Secondary Education Act is allocated using the Census Bureau's Small Area Income and Poverty Estimates for over 13,000 school districts. Medicare's Hospital Readmissions Reduction Program adjusts payments to thousands of hospitals based on noisy risk-adjusted readmission rates guptaPerformancePayHospitals2021. See waltersEmpiricalBayesMethods2024 for an overview of empirical Bayes methods in economics.} Shrinkage estimation can improve these estimates by borrowing strength from related units. But borrowing from which units? Adjacent neighborhoods tend to have similar economic mobility, while neighborhoods separated by a highway or school district boundary may differ sharply. This paper treats that choice as part of the estimation problem: we build a library of shrinkage rules that encode different notions of relatedness, then use the data to average over the resulting estimates---with weights that can, and sometimes do, concentrate on a single rule.
We observe noisy estimates $Y=(Y_1,\ldots,Y_n)^\top$ of unobserved parameters $\theta = (\theta_1,\ldots,\theta_n)^\top$, written in vector form as \[ Y = \theta + \varepsilon, \] where $\theta$ is the latent parameter vector to be estimated and $\varepsilon$ is sampling noise. The goal is to estimate $\theta$ accurately by pooling information across related units, and the question is how to pool. Empirical Bayes (EB) methods are a natural starting point. These methods estimate a prior distribution for $\theta$ and report the implied posterior means as the shrinkage estimates. Much of the nonparametric EB literature makes the prior flexible while retaining exchangeability across units:\footnote{A prior is exchangeable if it is invariant to relabeling the units: it can encode how spread out the $\theta_i$ are, but not which particular units are related. Section (ref) states this formally.} the resulting rule does not use information about which units are close in space or otherwise linked kiefer_consistency_1956, jiang_general_2009, koenkerMizeraConvexOptimization2014, soloff_multivariate_2025. Other EB approaches relax exchangeability by allowing the prior to vary with precision or covariates ignatiadisCovariatePoweredEmpiricalBayes2021, chenEmpiricalBayesWhen2024, luo_empirical_2025. We treat such prior specifications as one way to construct candidate shrinkage maps. A prior specification or covariance structure can motivate a shrinkage map $Y\mapsto f(Y)$ that returns an estimate of $\theta$, but we ask which map performs better under squared-error loss, not which specification is the right model for $\theta$. In this sense our approach is assumption-lean about the latent vector: it requires a model for the sampling noise in the estimates, but no model for the latent parameters themselves. The candidate library of shrinkage maps can therefore contain EB posterior-mean maps, spatial maps motivated by covariance models or adjacency structures, and maps whose tuning parameters are estimated from the same noisy estimates.
Candidate shrinkage maps can differ both in which other units enter each reported value and in how their tuning parameters are trained: one map may shrink through an estimated prior, another may give more weight to nearby or adjacent units, and another may estimate its smoothing weights from the observed vector $Y$. This flexibility can lower estimation error, but it creates an overfitting risk: the map that looks best for the observed vector $Y$ may be fitting the sampling noise $\varepsilon$ rather than estimating $\theta$ well. When the sampling noise is Gaussian with known covariance $\Sigma$, $\varepsilon \sim \mathcal{N}(0, \Sigma)$, Stein's Unbiased Risk Estimate (\ifmmodesure\elsesure\fi; steinEstimationMeanMultivariate1981) provides an observable, unbiased estimate of each candidate's expected squared-error loss. Selecting the candidate with the smallest \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi therefore guards against overfitting, favoring the map expected to estimate $\theta$ best rather than the one that fits the observed $Y$ most closely. Such an estimate matters because cross-validation does not measure the right quantity here: with one noisy estimate per latent parameter $\theta_i$, holding out unit $i$ does not pin down the squared error relative to $\theta_i$. \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi avoids this validation problem without placing a distributional assumption on $\theta$, and lets us compare candidate shrinkage maps on the squared-error loss scale used to evaluate estimation of $\theta$.
Our procedure has two steps. The researcher first specifies candidate classes and how each is trained. Applied to the observed data, these choices yield a finite library of trained maps $Y \mapsto f_k(Y)$, $k=1,\ldots,K$. \ifmmodesure\elsesure\fi is then evaluated for each trained map, including the correction for parameters trained on the same data, and minimizing \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi over convex weights yields the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average $Y \mapsto \sum_{k=1}^K w_k f_k(Y)$. Choosing a single map is the special case in which the weights concentrate on it, so model selection is nested within this averaging. Section (ref) formalizes this sequence; Table (ref) summarizes the workflow.
We make two methodological contributions. First, we give sufficient conditions under which \ifmmodesure\elsesure\fi minimization can be used to choose within parameterized classes of shrinkage maps $Y\mapsto f_\gamma(Y)$ that use the full vector of noisy estimates. The resulting oracle inequalities show that the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-selected map performs nearly as well as the best map in the class, with performance measured by squared-error loss against $\theta$. A close point of comparison is kwonOptimalShrinkageEstimation2025, who studies best-in-class shrinkage for panel fixed effects. That setting uses repeated observations over time and focuses on affine shrinkage rules. By contrast, the setting here is cross-sectional: one unit's estimate may enter another unit's fitted value through nonlinear rules based on geographic distance, spatial adjacency, or similarity in the observed estimates themselves. This within-class selection result also relates to work that uses estimated loss, \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi, or cross-validation to tune regularized many-parameter estimators abadieRiskMachineLearning2019, vives-i-bastidaSTRETCHINGNETMULTIDIMENSIONAL2023, adusumilliCrossValidationSURE2026.
Second, we give a \ifmmodesure\elsesure\fi-based model-averaging step for a finite library of trained shrinkage maps. This result combines the trained maps directly, taking each as a fixed building block. Each trained map must satisfy a regularity condition, checked separately for each candidate (Section (ref)). The oracle comparison then applies whether the maps come from closed-form formulas or iterative optimization. Given a finite library of estimators, the procedure chooses convex weights by minimizing \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi for the weighted average, with the weights treated as fixed inside the criterion (Section (ref)). The oracle guarantee is stated for fixed weights. When the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average is reported, its \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi value is a separate evaluation of the map $Y\mapsto f_{\hat w(Y)}(Y)$. That evaluation accounts for the data-dependence of the selected weights and trained parameters. A close model-averaging comparison is hansenLeastSquaresModel2007, who studies weights over linear least-squares fits. Here the weights range instead over trained nonlinear maps, including shrinkage rules whose tuning parameters are estimated from the data before averaging.
The empirical application uses Opportunity Atlas data to estimate tract-level economic mobility across 20 commuting zones. The tract-level estimates show strong spatial patterns: nearby tracts often have similar estimated mobility, but geography, adjacency, school-district boundaries, and historical segregation can make different forms of relatedness empirically relevant. In the main neighborhood mobility comparison, candidate rules differ by distance metric and incorporation of demographic covariates, and the empirical analysis reports the \ifmmodesure\elsesure\fi-weighted average of candidate maps as the primary estimator. A separate Cook County comparison adds the value-similarity rule of Section (ref), which lets large differences in the observed estimates reduce smoothing between nearby tracts. In this setting, the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-weighted average of candidate maps reduces \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-estimated mean squared error (MSE) by about 55% relative to the raw maximum-likelihood benchmark (\ifmmode\text{\textup{\textsc{mle}}}\else\textup{\textsc{mle}}\fi, the unshrunk tract estimates $Y$) and by about 27% relative to \ifmmode\text{\textup{\textsc{close}}}\else\textup{\textsc{close}}\fi-\ifmmode\text{\textup{\textsc{gauss}}}\else\textup{\textsc{gauss}}\fi, the closed-form non-spatial EB benchmark---the Gaussian member of the conditional location-scale (CLOSE) family of chenEmpiricalBayesWhen2024. The best individual spatial shrinkage rule varies across commuting zones. The broader lesson is that the relevant notion of relatedness is an empirical choice, and we recommend making that choice by \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted averaging rather than by committing to one form of spatial smoothing in advance.
The remainder of the paper proceeds as follows. Section (ref) sets up the Gaussian compound-decision problem for shrinkage estimation, introduces \ifmmodesure\elsesure\fi as an observable estimate of squared-error loss, gives examples of candidate shrinkage maps, and then shows how \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi can be used to choose weighted averages of trained candidates. Section (ref) states the two oracle inequalities, one for \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi selection within a class and one for \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted averaging across a library. Section (ref) applies the framework to Opportunity Atlas mobility estimates, first comparing estimators on \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-estimated MSE for the latent mobility vector and then asking how the same shrinkage estimates affect a targeting exercise that selects high-mobility tracts. Section (ref) concludes.
The researcher observes an $n$-dimensional vector $Y = \theta + \varepsilon$, where $\theta \in \mathbb{R}^n$ is a fixed but unknown parameter vector and $\varepsilon \sim \mathcal{N}(0, \Sigma)$ is Gaussian noise with known covariance matrix $\Sigma$. The Gaussian sampling model is a working approximation for microdata-derived estimates with reported precision, as in related empirical Bayes applications chenEmpiricalBayesWhen2024. In such applications, each component $Y_i$ is itself an average or regression coefficient estimated from an underlying micro-sample, and the known covariance matrix $\Sigma$ reflects the sampling precision of those estimates. In the Opportunity Atlas application, $\Sigma$ is taken to be the diagonal matrix of reported marginal variances; Appendix (ref) records what changes when \ifmmodesure\elsesure\fi is computed with an approximate covariance matrix.
This is a compound decision problem in the sense of robbins1951asymptotically: $\theta$ is fixed (no prior is placed on it), the researcher chooses a decision rule $f$ that returns a vector of actions $f(Y)$ with one action $f_i(Y)$ per unit, each allowed to depend on the entire vector $Y$, and the rule is judged by its average squared-error loss against the fixed $\theta$. Each such rule is a map $f\colon\mathbb{R}^n\to\mathbb{R}^n$, and we compare rules within a candidate class $\mathcal F$ by the realized loss \[ L_n(f)=\frac{1}{n}\|f(Y)-\theta\|_2^2. \] The infeasible oracle in this class is \[ f^* \in \operatorname*{arg\,min}_{f\in\mathcal F} L_n(f). \] This oracle uses the unknown vector $\theta$ and is therefore only a benchmark. The statistical problem is to use the observed vector $Y$ to select, from the rules $f\in\mathcal F$, an estimate whose realized loss is close to this oracle benchmark.
Throughout this paper, expectations are over the sampling noise in $Y=\theta+\varepsilon$, treating $\theta$ as fixed. We reserve the term risk for the expected realized loss, $R_n(f):=\mathbb{E}[L_n(f)]$. Because $\theta$ is unknown, neither $L_n(f)$ nor $R_n(f)$ can be evaluated directly, so effective estimation requires an observable criterion whose behavior tracks the unobserved loss. Appendix (ref) gives a simple condition under which lower squared-error estimation error also reduces errors in downstream comparisons based on the estimated vector.
As a concrete example, each $Y_i$ is a tract-level estimate of economic mobility from the Opportunity Atlas chettyOpportunityAtlasMapping2018, with known sampling variance $\Sigma_{ii} = \sigma_i^2$ reflecting the precision of the underlying microdata. The parameter vector $\theta \in \mathbb{R}^n$ represents true neighborhood-level mobility across hundreds to thousands of Census tracts in a commuting zone. Reported standard errors vary substantially across tracts because the underlying sample sizes differ. Because the units are neighborhoods, spatial structure is a natural source of pooling information. The application therefore motivates shrinkage rules that can leverage geographic relationships rather than treating all tracts as exchangeable.
For any continuously differentiable estimator $f$ such that $\mathbb{E}[\|f(Y)-Y\|_2^2]<\infty$ and $\mathbb{E}[\sum_{i,j}|\Sigma_{ij}\partial_j f_i(Y)|]<\infty$, Stein's lemma gives an observable statistic $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f)$ satisfying $\mathbb{E}[\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f)] = R_n(f)=\mathbb{E}[L_n(f)]$:\footnote{Continuous differentiability is used here as a convenient sufficient condition. Unbiasedness of \ifmmodesure\elsesure\fi also extends to weakly differentiable maps $Y\mapsto f(Y)$ satisfying the same integrability conditions. For example, $x\mapsto |x|$ has weak derivative $\operatorname{sign}(x)$. The sufficient conditions used below are stated as continuous-differentiability conditions.} \[ \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f) := \underbrace{\frac{1}{n}\|Y - f(Y)\|_2^2 - \frac{1}{n}\mathrm{tr}(\Sigma)}_{\text{noise-corrected in-sample MSE}} + \underbrace{\frac{2}{n}\mathrm{tr}\{\Sigma Df(Y)\}}_{\text{complexity correction}}. \] Here $Df(Y) = [\partial f_i(Y)/\partial Y_j]_{ij}$ is the $n \times n$ Jacobian matrix of $f$, and $\mathrm{tr}(\Sigma)=\sum_{i=1}^n \Sigma_{ii}$ is the trace of the noise covariance matrix, the total sampling variance in $Y$. The first term subtracts the irreducible noise variance $\mathrm{tr}(\Sigma)$ from the in-sample prediction error, converting it into an estimate of the estimation error $\|f(Y) - \theta\|_2^2$ rather than the prediction error $\|f(Y) - Y^{\mathrm{new}}\|_2^2$, where $Y^{\mathrm{new}}$ is an independent draw from the same model. But this estimate is generally biased downward when the rule uses the same $Y$ to form the reported values being evaluated.\footnote{Cross-fitting eliminates the downward bias issue by making the decision rule separable across sample splits ignatiadisCovariatePoweredEmpiricalBayes2021, chenCompoundSelectionDecisions2025. \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi instead accounts for the dependence directly through the complexity correction, using the full sample without splitting.} The complexity correction measures the sensitivity of $f$ to the data---a penalty that is larger for more flexible estimators and corrects for this optimism.
For a linear smoother $f_S(Y)=SY$ with a data-independent matrix $S$, the same formula reduces to \[ \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_S) = \frac{1}{n}\|(I-S)Y\|_2^2 -\frac{1}{n}\mathrm{tr}(\Sigma) +\frac{2}{n}\mathrm{tr}(\Sigma S). \] Once $S$ is fixed, every term on the right-hand side is computable from $Y$, $\Sigma$, and $S$.
The statistic $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f)$ is the observable risk criterion used below to evaluate shrinkage estimators. It can also serve as a training criterion within a candidate class. Section (ref) gives conditions under which $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f)$ tracks $L_n(f)$ closely enough to justify those choices.
Write $f_i(Y)$ for the $i$th coordinate of a shrinkage rule $f$, the value reported for unit $i$. The examples below supply the paper's candidate classes in the empirical application of Section (ref). For a parameter space $\Gamma\subset \mathbb{R}^d$, each example represents a parameterized class of maps $\{f_\gamma:\gamma\in\Gamma\}$ encoding one notion of relatedness, and they are ordered in this section by the degree to which the full vector $Y$ enters each component rule $f_i(Y)$. In the normal--normal empirical Bayes (\ifmmodenn-eb\elsenn-eb\fi) example, $f_i(Y)$ depends on $Y$ only through unit $i$'s own estimate and indirectly through trained scalars shared by all units. The Gaussian-process (\ifmmode\textup{\textsc{gp}}\else\textup{\textsc{gp}}\fi) examples let other units' estimates enter $f_i(Y)$ directly, through weighted combinations $f_i(Y)=\sum_j s_{ij}Y_j$ whose smoothing weights $s_{ij}$ are determined by spatial relationships such as geographic distance. The value-similarity example lets the smoothing weights depend on the observed estimates themselves, $s_{ij}=s_{ij}(Y)$: among nearby units, those with similar values receive more weight. Each example is motivated by a working model for the latent vector $\theta$, but the working model serves only to construct the class $\{f_\gamma: \gamma\in \Gamma\}$. The \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi comparisons in the paper evaluate the resulting maps by squared-error loss for the fixed vector $\theta$, whether or not the working model is properly specified. When the tuning parameters of a class are themselves trained on $Y$, \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi must account for that dependence; Section (ref) gives the correction.
One way to relax this working assumption keeps independence across units but drops identical distribution: conditioning on unit covariates $X_i$---which may include the standard error $\sigma_i$---replaces the two scalars with functions, so the working prior becomes $\theta_i \mid X_i \sim \mathcal{N}(\mu(X_i), \tau^2(X_i))$ and the trained rule shrinks $Y_i$ toward $\hat\mu(X_i)$ by an amount governed by $\hat\tau^2(X_i)$ and $\sigma_i^2$ ignatiadisCovariatePoweredEmpiricalBayes2021, chenEmpiricalBayesWhen2024. Units with the same covariates are still treated symmetrically. Nothing in the rule links unit $i$ to specific other units.
Figure (ref) illustrates the difference between global and spatial shrinkage on mobility estimates for Cook County tracts selected from the Chicago commuting zone. The raw estimates (panel A) are visibly noisy. \ifmmodenn-eb\elsenn-eb\fi (panel B) shrinks every tract toward the same global target $\bar Y$, the amount depending only on the tract's own precision: estimates above $\bar Y$ are pulled down and estimates below are pulled up, regardless of the values of neighboring tracts. The result both over-smooths and under-smooths: tracts in different community areas---Chicago's named groupings of Census tracts---are blurred toward one another, while a tract surrounded by similar neighbors is still dragged away from its local average toward the distant global mean $\bar Y$. The spatial \ifmmode\textup{\textsc{gp}}\else\textup{\textsc{gp}}\fi (panel C) instead shrinks each tract toward a local neighborhood average, so the smoothed map retains more of the spatial pattern in the raw estimates while reducing tract-level noise.
Figure (ref) quantifies the competing shrinkage targets. Each tract's raw estimate $Y_i$ is plotted against its leave-one-out spatial \ifmmodegp\elsegp\fi shrinkage target $\mu_i$, defined by writing the \ifmmode\textup{\textsc{gp}}\else\textup{\textsc{gp}}\fi prediction for tract $i$ as $f_i(Y)=S_{ii}Y_i+(1-S_{ii})\mu_i$. The spatial target $\mu_i$ plays the role that the global mean $\mu$ played in \ifmmode\text{\textup{\textsc{nn-eb}}}\else\textup{\textsc{nn-eb}}\fi. The figure overlays the two targets: \ifmmode\text{\textup{\textsc{nn-eb}}}\else\textup{\textsc{nn-eb}}\fi shrinks every tract toward $\bar{Y}$ (orange line), regardless of spatial context, while the spatial \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi shrinks toward $\mu_i$ (blue diagonal). Whenever a tract's raw estimate lies between the two targets---the shaded region of Figure (ref), containing 39% of the 1,304 tracts---the two procedures pull $Y_i$ in opposite directions. The tracts highlighted in Figure (ref) sit on opposite lobes of this region: the spatial \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi pulls the tract in Albany Park up toward its high-mobility neighborhood while \ifmmode\text{\textup{\textsc{nn-eb}}}\else\textup{\textsc{nn-eb}}\fi pulls it down toward $\bar Y$, and the mirror pattern holds for the tract in East Garfield Park.
The \ifmmodegp\elsegp\fi shrinkage maps above have the form $SY$ once the covariance matrix is fixed: for a given tract, the estimates that enter its reported value are determined by distance, not by the realized values. The next example keeps the spatial smoothing structure but lets differences in the observed estimates reduce the geography-based similarity between nearby tracts. The empirical motivation is that nearby tracts can be separated by highways, rivers, school-district boundaries, or boundaries associated with historical segregation. In such cases, smoothing over geographic distance may average across places whose observed mobility estimates differ sharply. The construction is the bilateral filter, originally developed as an edge-preserving smoother in image processing tomasiBilateralFilteringGray1998, adapted to the Gaussian shrinkage form above.
Figure (ref) illustrates the difference between geography-only and value-similarity shrinkage targets for a set of central-Chicago tracts in Cook County. The figure displays the leave-one-out shrinkage targets $\mu_i = (f_i(Y) - S_{ii} Y_i) / (1 - S_{ii})$. Panel (A) shows the noisy estimates, including visible local contrasts across some community-area borders. Panel (B) uses a geographic smoother, so nearby tracts enter the formula according to distance. This attenuates some of those local contrasts. Panel (C) also uses similarity in the observed estimates, so pairs of tracts with dissimilar values receive less weight and more of the visible contrast is retained. Compare panels (B) and (C) at communities \textcircled{1} and \textcircled{3}: geographic smoothing alone bleeds North Lawndale's low estimates across the boundary into the Lower West Side, while the value-similarity target preserves the contrast. Both smoothed panels use the same preliminary covariate adjustment (OLS residualization on demographic covariates), held fixed across panels, so the comparison isolates geography-only versus value-similarity smoothing.
Each way of defining which estimates are related---through geographic distance, tract adjacency, observed-value similarity, or their combinations---yields a distinct candidate class with its own tuning parameters $\gamma$. A practitioner then faces two nested choices: within each candidate class, how should $\gamma$ be trained? And when averaging across classes, what weight should each trained candidate receive? The across-candidate comparison is posed in terms of \ifmmodesure\else\textsc{sure}\fi, while within-candidate training may use \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi or another criterion chosen by the researcher. Either way, the resulting trained map is later evaluated by \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi. The value-similarity example shows one source of extra $Y$-dependence, through the matrix $S(Y)$. The next subsection turns to another: tuning parameters that are trained on $Y$.
In practice, using a candidate class $\{f_\gamma:\gamma\in\Gamma\}$ requires a training rule chosen by the researcher. This rule maps the observed vector $Y$ to parameter values $\hat\gamma(Y)$. The tuning parameter $\gamma$ may be a length scale, a variance component, a regularization strength, or a value-similarity parameter. Write \[ F(Y)=f_{\hat\gamma(Y)}(Y) \] for the trained map. For risk evaluation, \ifmmodesure\elsesure\fi is applied to the full map $F$ using the chain rule, not to $f_\gamma$ with the realized value $\hat\gamma(Y)$ plugged in and treated as fixed. With output coordinates as rows, the chain rule gives \[ DF(Y) = \underbrace{ D_y f_\gamma(Y)\big|_{\gamma=\hat\gamma(Y)} }_{\text{direct sensitivity, holding $\gamma$ fixed}} + \underbrace{ D_\gamma f_\gamma(Y)\big|_{\gamma=\hat\gamma(Y)}D_Y\hat\gamma(Y) }_{\text{sensitivity from training on $Y$}}. \] Therefore \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi for the trained map equals the fixed-parameter \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi formula evaluated at the realized $\hat\gamma(Y)$, plus an additional training correction: \[
\] We call the first brace---the \ifmmodesure\elsesure\fi formula with $\gamma$ treated as fixed---proxy \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi. The second brace is the training correction.
The training rule can be a closed-form estimator or an iterative algorithm. Method-of-moments estimates, maximum-likelihood estimates, and fixed iterative optimization routines all produce a map $Y\mapsto\hat\gamma(Y)$.\footnote{Iterative optimization includes standard stochastic-gradient methods. The Opportunity Atlas implementation uses AdamW, a decoupled-weight-decay variant of Adam, for the trainable \ifmmodegp\elsegp\fi candidates kingmaAdamMethodStochastic2017,loshchilovDecoupledWeightDecay2019.} The training choice is part of candidate construction. For \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi evaluation, the relevant object is the resulting trained map $Y\mapsto f_{\hat\gamma(Y)}(Y)$, including the sensitivity of $\hat\gamma(Y)$ to the same data.
For a \ifmmodegp\elsegp\fi candidate class indexed by covariance parameters $\gamma$, write the fixed-parameter zero-mean shrinkage map as $f_\gamma(Y)=K_\gamma(K_\gamma+\Sigma)^{-1}Y$. The conventional \ifmmode\textup{\textsc{gp}}\else\textup{\textsc{gp}}\fi training rule is maximum marginal likelihood. Given a working covariance family $\{K_\gamma:\gamma\in\Gamma\}$, this training rule chooses $\hat\gamma_{\mathrm{ML}}(Y)$ by maximizing the Gaussian marginal likelihood for $Y$ with covariance $K_\gamma+\Sigma$ rasmussenGaussianProcessesMachine2006.\footnote{Under the auxiliary \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi prior specification $\theta\sim\mathcal{N}(0,K_\gamma)$ and $Y\mid\theta\sim\mathcal{N}(\theta,\Sigma)$, integrating out $\theta$ gives the marginal likelihood $Y\sim\mathcal{N}(0,K_\gamma+\Sigma)$. If this covariance specification is correct for some $\gamma_0\in\Gamma$, the posterior mean based on $K_{\gamma_0}$ is the squared-error Bayes rule, so likelihood-based covariance training and squared-error prediction are aligned. If no such $\gamma_0$ exists, the working covariance family is misspecified, and maximum marginal likelihood targets the $\gamma$ minimizing Kullback--Leibler divergence within $\{K_\gamma:\gamma\in\Gamma\}$, which need not minimize squared-error risk of the induced shrinkage map. bachocCrossValidationMaximum2013,bachocAsymptoticAnalysisCovariance2018 show this target mismatch for \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi covariance estimation under misspecification: maximum likelihood targets Kullback--Leibler divergence, whereas cross-validation targets prediction mean squared error. \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi plays the corresponding role here by estimating the mean squared error of the shrinkage rule.} A training rule need not minimize \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi for the resulting trained map to be evaluated using \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi. A likelihood-trained \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi map could therefore be included as one trained candidate in the finite library considered for averaging in Section (ref). In the empirical application, the trainable \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi candidates in Table (ref) are trained by minimizing proxy \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi over $\gamma$.
In practice, the researcher does not need to derive $D_Y\hat\gamma(Y)$ by hand. Automatic differentiation can propagate sensitivities through an implemented training algorithm. When $\hat\gamma(Y)$ is characterized by first-order conditions, implicit-differentiation tools use those conditions to obtain the same sensitivity without deriving a new formula for each training problem blondel_efficient_2022. The trace terms in \ifmmodesure\elsesure\fi can then be computed efficiently using randomized trace estimation hutchinson_stochastic_1990,nobel_tractable_2023.\footnote{If $v$ is a random vector with $\mathbb{E}[vv^\top]=\Sigma$, then $\mathbb{E}[v^\top DF(Y)v]=\mathrm{tr}\{\Sigma DF(Y)\}$. A randomized trace estimate averages $v^\top DF(Y)v$ over several independent draws of $v$. Each draw requires the Jacobian-vector product $DF(Y)v$, which automatic differentiation computes by propagating one direction through the chain rule. This avoids constructing $DF(Y)$ explicitly and avoids the matrix-matrix products that would arise from carrying the full Jacobian through the training rule.}
The chain-rule decomposition also clarifies how the theory treats trained candidates. Once the training rule is fixed, the composite map $Y\mapsto f_{\hat\gamma(Y)}(Y)$ is the estimator evaluated by \ifmmodesure\elsesure\fi, and Appendix (ref) gives sufficient conditions under which this composite map satisfies the regularity condition used for averaging.
Suppose the preceding steps produce a finite library of $K$ candidate estimators, $f_1,\ldots,f_K$. The candidates may differ in the information they use, their preprocessing choices, or their training rules. Rather than committing in advance to one candidate class, \ifmmodesure\elsesure\fi model averaging uses \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi itself to choose convex weights across the trained maps. When candidate rules make different errors, a convex combination can beat every single candidate. In the Opportunity Atlas application, the best single candidate differs across commuting zones. Averaged across them, however, the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average has lower estimated risk than any single candidate (Section (ref)). The oracle comparison below is therefore against the best fixed convex combination---a stronger benchmark than the best single candidate, since every single candidate is a vertex of the simplex.
Specifically, for weights in the simplex we form \[ f_w(Y) = \sum_{k=1}^K w_k\, f_k(Y), \qquad w \in \Delta^{K-1} := \{w \in \mathbb{R}^K_+ : \textstyle\sum_k w_k = 1\}, \] and choose weights by minimizing the fixed-weight \ifmmodesure\elsesure\fi criterion: \[ \hat{w} \in \operatorname*{arg\,min}_{w \in \Delta^{K-1}} \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_w). \] Within this minimization problem, each proposed weight vector $w$ is treated as fixed. This fixed-weight criterion is the object used for the oracle comparison in Section (ref). After the observed data select $\hat w(Y)$, the reported map is the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-weighted average, \[ \tilde{f}(Y) := f_{\hat{w}(Y)}(Y) = \sum_{k=1}^K \hat{w}_k(Y)\, f_k(Y). \] As a function of the weights, $f_w$ is linear in $w$, and for each fixed $w$ its Jacobian satisfies $Df_w=\sum_k w_k Df_k$. Therefore the Jacobian trace $\mathrm{tr}(\Sigma\,Df_w) = \sum_k w_k \mathrm{tr}(\Sigma\,Df_k)$ is also linear in $w$, so the fixed-weight objective $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_w)$ is quadratic in $w$.\footnote{The simplex-constrained quadratic program can be solved with standard convex optimization tools. If \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi is reported for the final data-chosen average, differentiating the selected weights can be handled by differentiating the optimization conditions when the solution is locally stable.}
Before averaging over candidate estimators, each candidate's \ifmmodesure\elsesure\fi value must be computed for the trained map actually produced from $Y$. When tuning parameters are trained on $Y$, evaluating \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi as if those parameters were fixed can understate risk. This is the excess-optimism problem for \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-tuned estimators studied by tibshiraniExcessOptimismBiased2019. The same concern applies to weights chosen by minimizing \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi, so the reported \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi value for the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average evaluates the full map $\tilde f$, rather than the fixed-weight criterion used to choose $\hat w$.
\begingroup
\endgroup
Table (ref) summarizes the procedure and separates three uses of \ifmmodesure\elsesure\fi. In Step 2, \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi may be used as a training criterion within a candidate class. In Steps 3 and 4, \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi gives a common squared-error criterion for trained candidates and fixed-weight averages, so the researcher can select the smallest-\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi candidate or choose convex weights across candidates. In Step 5, \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi is used to evaluate the map actually reported, including the sensitivity induced by trained parameters and, when relevant, by selected weights.
The training and averaging steps lead to two theoretical questions. The first question concerns selection within a fixed parameterized class of shrinkage maps $\{f_\gamma:\gamma\in\Gamma\}$: if the researcher chooses $\hat\gamma$ by minimizing $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_\gamma)$ over $\Gamma$, how close is the realized loss $L_n(f_{\hat\gamma})$ to the loss $\min_{\gamma\in\Gamma}L_n(f_\gamma)$ of the best parameter choice in the class? The second concerns averaging after a finite candidate library has been assembled: if the researcher chooses convex weights by minimizing the fixed-weight \ifmmodesure\elsesure\fi criterion, how close is the resulting loss to the loss from the fixed convex combination with the smallest realized loss?
Throughout this section, the sampling experiment and the estimator classes are implicitly indexed by the dimension $n$. When the $n$-dependence needs to be explicit, we write $Y^{(n)},\theta^{(n)},\varepsilon^{(n)},\Sigma_n$ and $f_{n,\gamma}:\mathbb{R}^n\to\mathbb{R}^n$. Otherwise, we suppress the $n$-dependence and write $Y,\theta,\varepsilon,\Sigma$, $f_\gamma$, $\mathcal F=\{f_\gamma:\gamma\in\Gamma\}$, and $\Gamma$. Constants described as independent of $n$ are uniform over this sequence.
This section establishes two oracle inequalities for \ifmmodesure\elsesure\fi-based selection among shrinkage maps. The within-class guarantee, Theorem (ref), applies to a compact parameterized class $\mathcal F=\{f_\gamma:\gamma\in\Gamma\}$ of maps that report vectors of shrinkage estimates. If $\hat\gamma\in\operatorname*{arg\,min}_{\gamma\in\Gamma}\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_\gamma)$ and $\gamma^*\in\operatorname*{arg\,min}_{\gamma\in\Gamma}L_n(f_\gamma)$, then the theorem bounds the excess realized loss $L_n(f_{\hat\gamma})-L_n(f_{\gamma^*})$. The averaging guarantee, Proposition (ref), applies after a finite library of trained candidate maps $f_1,\ldots,f_K$ has been assembled. For fixed weights $w\in\Delta^{K-1}$, write $f_w(Y)=\sum_{k=1}^K w_k f_k(Y)$. If $\hat w(Y)\in\operatorname*{arg\,min}_{w\in\Delta^{K-1}}\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_w)$ and $w^*\in\operatorname*{arg\,min}_{w\in\Delta^{K-1}}L_n(f_w)$, then the proposition compares the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-weighted average $\tilde f(Y):=f_{\hat w(Y)}(Y)$ with the fixed-weight convex combination $f_{w^*}$ that has the smallest realized loss over the simplex.
The guarantees do not place a distribution on the latent vector $\theta$. For each dimension $n$, we fix $\theta\in\mathbb{R}^n$ and take expectations only over the sampling noise in $Y=\theta+\varepsilon$. Thus $\gamma^*$ and $w^*$ are realized-loss benchmarks for the candidate maps under consideration, not procedures derived from a correctly specified prior on the latent parameters. Although Bayesian or empirical Bayes specifications can motivate maps such as $f_\gamma$, the guarantees below evaluate the resulting maps directly. If every map in $\mathcal F$ or every convex average of $f_1,\ldots,f_K$ has high realized loss for the fixed vector $\theta$, the theory does not remove that approximation error; it controls the additional loss from choosing within the class or library using \ifmmodesure\elsesure\fi.
Theorem (ref) and Proposition (ref) impose different requirements because they apply at different points in the construction. Theorem (ref) studies exact minimization of \ifmmodesure\elsesure\fi over the full parameter set $\Gamma$, so its regularity condition is uniform over the parameterized class $\mathcal F$. Proposition (ref) begins after the finite candidate maps $f_1,\ldots,f_K$ have already been constructed. Those candidates may include trained parameters and may come from different training rules; Proposition (ref) does not require each $f_k$ to solve a within-class \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi minimization problem. For averaging, the optimization is over the weights $w$, and the regularity condition is imposed separately on each composite candidate map, including the dependence introduced when a candidate's tuning parameters are trained on the same data.
Theorem (ref) studies the rule selected by minimizing $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_\gamma)$ over the compact parameter set $\Gamma$. For this oracle comparison, the observable criterion $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_\gamma)$ must track the realized loss $L_n(f_\gamma)$ uniformly over the candidate class, not only at a fixed map. Existing \ifmmodesure\elsesure\fi selection guarantees cover finite collections of globally Lipschitz estimators bellecSecondOrderStein2021; the result below extends that benchmark in two ways. First, it allows $\Gamma$ to be a compact continuum of tuning parameters, such as length scales and kernel variances. Second, it replaces global Lipschitzness with a broader regularity condition: the adjustment relative to the raw estimate, its derivative, and its variation with $\gamma$ are controlled by polynomial envelopes in the input vector of raw estimates.
The value-similarity rule shows why the broader regularity condition is needed. Changing one observed estimate can change both its own reported value and the weights assigned to other estimates, so the resulting map can fail to be globally Lipschitz even in a two-dimensional fixed-parameter case.\footnote{This issue is not specific to the value-similarity example. kimLipschitzConstantSelfAttention2021 show that standard dot-product self-attention is not globally Lipschitz in its input. Self-attention is not part of the empirical library in this paper, but it is another example of a flexible data-adaptive map that would not be covered by a theory taking global Lipschitzness as a primitive condition.} Appendix (ref) gives this non-Lipschitz calculation, and Appendix (ref) verifies a fixed value-similarity building block as a single candidate for averaging. Under a bounded-row-sum condition on the fixed geographic factor of the kernel, Appendix (ref) shows that this fixed value-similarity building block is no more costly, in the regularity order used below, than fixed linear smoothers. Training of tuning parameters raises a related verification issue: regularity of the fixed maps $f_\gamma$ does not by itself establish regularity of the reported map $Y\mapsto f_{\hat\gamma(Y)}(Y)$, because the trained parameters $\hat{\gamma}$ are also functions of $Y$. Appendix (ref) gives primitive conditions under which trained maps still satisfy the per-candidate regularity condition used for averaging.
For symmetric matrices, write $A\succ0$ for positive definite and $A\succeq0$ for positive semidefinite. Let $\lambda_{\max}(A)$ denote the largest eigenvalue of a symmetric matrix. Norms $\|\cdot\|_2$, $\|\cdot\|_{\mathrm{op}}$, and $\|\cdot\|_F$ denote Euclidean, operator, and Frobenius norms, respectively. The first assumption restates the Gaussian sampling model of Section (ref) along the sequence of experiments and adds two uniform bounds: on the noise scale and on the average magnitude of the latent vector.
The bound on $\theta^{(n)}$ is an average-magnitude condition: $\|\theta^{(n)}\|_2^2/n$ remains bounded, so the fixed latent vectors do not grow in average squared size along the sequence. This condition is implied by putting each coordinate of $\theta^{(n)}$ in a fixed compact set, but it is slightly more general: individual coordinates may exceed any fixed bound as long as their squared magnitudes remain controlled on average.
For a generic input $y\in\mathbb{R}^n$, write $g_\gamma(y) := f_\gamma(y) - y$ for the adjustment relative to the raw estimate.\footnote{Under Tweedie's formula efronTweediesFormulaSelection2011, the posterior mean of the normal location model satisfies $\mathbb{E}[\theta | Y] = Y + \Sigma \nabla \log p(Y)$, where $p$ is the marginal density of $Y$ and $\nabla\log p(Y)$ is the score of that marginal density. Thus, when a candidate map is motivated by a posterior-mean formula, $g_\gamma$ can be read as an estimate of this score term. Assumption (ref) places regularity conditions on the candidate maps themselves; it does not require a correctly specified marginal density for $Y$. See ghoshSteinSUREScoreMatching2025 for a unified treatment of \ifmmodesure\elsesure\fi and Hyv\"arinen score matching.} Write $\|g(y)\|_{W} := \|g(y)\|_2 + \|Dg(y)\|_F$ for the combined function--Jacobian norm. Assumption (ref) formalizes polynomial-envelope regularity by controlling $\|g_\gamma(y)\|_W$ and the corresponding variation with $\gamma$ for every $y\in\mathbb{R}^n$, rather than only at the realized random vector $Y$.
The reference point $\gamma_0$ corresponds to the map $f_{\gamma_0}$ in the class. The first term in Assumption (ref) controls the adjustment $g_{\gamma_0}$ at this reference point, while the supremum controls how both $g_\gamma$ and $Dg_\gamma$ vary with $\gamma$. If the class contains the identity map, then $g_{\gamma_0}\equiv0$ is the natural reference choice. Otherwise, any fixed member of the class satisfying the displayed envelope can serve as the reference. The exponent $\beta$ controls how the envelope may grow with the normalized input magnitude $\|y\|_2/\sqrt n$: when $\beta=0$, the bound is uniform in $y$, while $\beta>0$ permits polynomial growth. For fixed linear smoothers $f(y)=Sy$, Lemma (ref) in Appendix (ref) shows that the pointwise envelope holds with $\beta=1/2$ and $\nu_n=O(\sqrt n)$ when the smoothing matrix $S$ has bounded operator norm and bounded maximum row Euclidean norm, $\max_{i\leq n}\|S_{i\cdot}\|_2$, where $S_{i\cdot}$ denotes the $i$th row of $S$. In spatial applications, the row bound captures a bounded-sensitivity form of local borrowing: the reported estimate for one spatial unit can average information from nearby units, but the Euclidean norm of the corresponding row of $S$ does not grow with $n$.
The continuous-differentiability requirement excludes some familiar nonsmooth shrinkage rules. For example, the one-dimensional hard-thresholding rule $f(y)=y\,1\{|y|>\tau\}$, $\tau>0$, has jumps at $y=\pm\tau$, so its derivative is not defined there and the rule is not covered by Assumption (ref). The pointwise polynomial envelope is a sufficient condition for the concentration theorem below.
With these regularity conditions in hand, we state the main result: \ifmmodesure\elsesure\fi tracks the realized loss uniformly over $\mathcal F$, and the expected excess realized loss of the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-selected rule is of order $\nu_n \max\{d_\Gamma,1\}^{4+\beta}/n$.
The two displays control different quantities: the first is an uncentered uniform approximation bound for $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f)$ over $\mathcal F$, while the second is the excess realized loss from choosing $\hat f$ by minimizing $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n$. In the oracle comparison, the component of $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f)-L_n(f)$ that does not depend on $f$ is common to $\hat f$ and $f^*$, so the $n^{-1/2}$ term from the first display does not enter. The proof is deferred to Appendix (ref), which first establishes a general Sobolev-moment concentration result and then applies Assumption (ref) to obtain the displayed bound.
The factor $\max\{d_\Gamma,1\}^{4+\beta}$ is a convention for including singleton classes in the same rate display. When $d_\Gamma\geq1$, this factor is $d_\Gamma^{4+\beta}$; when $\mathcal F$ is a singleton, with $d_\Gamma=0$, the uniform concentration bound still contains the $n^{-1/2}$ term and the $\nu_n/n$ contribution from the reference map $f_{\gamma_0}$, even though there is no variation over $\gamma$ to control. The constants hidden by $\lesssim$ depend only on the fixed bounds $C_\theta$, $\bar\sigma$, the regularity exponent $\beta$, and the diameter bound $D_\Gamma$. Apart from the displayed factors, these constants are uniform in $n$, $d_\Gamma$, and $\nu_n$.
The averaging result returns to the finite library of trained maps $f_1,\ldots,f_K$ from Section (ref) and imposes regularity on each final trained map $f_k$ separately, rather than uniformly over the full training family for each candidate. For averages of affine estimators $f_k(Y)=S_kY$ with fixed matrices $S_k$, sharp oracle inequalities based on unbiased risk estimates are available dalalyanSharpOracle2012; Proposition (ref) covers finite libraries of nonlinear trained maps whose parameters are learned from the same data.
Assumption (ref) is a per-candidate moment condition on each final trained map $f_k$, including derivative contributions from any training rule used to construct that map. Because the averaging library is finite, the assumption does not require uniform increment bounds over a parameter set. A convenient sufficient condition is a pointwise polynomial envelope: if, for all $y\in\mathbb{R}^n$ and $k=1,\ldots,K$, \[ \|g_k(y)\|_2+\|Dg_k(y)\|_F \leq \nu_n^{(k)}\left(1+\frac{\|y\|_2}{\sqrt n}\right)^{2\beta_k}, \] then, under Assumption (ref), the Gaussian moment bound in the proof of Lemma (ref) gives Assumption (ref) with the same exponent $\beta_k$ and with $\nu_n^{(k)}$ inflated by a constant depending only on $\beta_k$, $C_\theta$, and $\bar\sigma$. Appendix (ref) gives the corresponding class-level Sobolev-moment formulation, and Appendix (ref) gives sufficient conditions for the composite map $Y\mapsto f_{\hat\gamma_k(Y)}^{(k)}(Y)$ to satisfy Assumption (ref) when the parameter estimate $\hat\gamma_k(Y)$ is trained on the same data.
The proposition below applies the fixed-weight \ifmmodesure\elsesure\fi criterion to this finite library and shows that, with $\bar\beta := \max_k \beta_k$ and $\bar\nu_n := \max_k \nu_n^{(k)}$, the expected regret of the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-weighted average, relative to the best fixed convex combination, is of order $(\log(eK))^{4+\bar\beta}\,\bar\nu_n/n$. For fixed $w\in\Delta^{K-1}$, write \[ f_w(Y)=\sum_{k=1}^K w_k f_k(Y), \] where the Jacobian of $f_w$ treats the entries of $w$ as constants.
The proof, in Appendix (ref), decomposes the fixed-weight \ifmmodesure\elsesure\fi error into a noise term common to all weights plus a weighted average of candidate-specific terms. The maximum over the library yields the $(\log(eK))^{4+\bar\beta}$ factor.
The fixed-weight \ifmmodesure\elsesure\fi criterion is a quadratic function of $w$ and can be computed exactly. Since $Df_w(Y)=\sum_k w_k Df_k(Y)$ when the weights are held fixed, define $A_{k\ell}(Y):=n^{-1}f_k(Y)^\top f_\ell(Y)$, $b_k(Y):=n^{-1}Y^\top f_k(Y)$, and $c_k(Y):=n^{-1}\mathrm{tr}\{\Sigma Df_k(Y)\}$. Then the fixed-weight criterion can be written as \[ \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_w) = \underbrace{ w^\top A(Y)w }_{\text{quadratic in }w} + \underbrace{ 2\{c(Y)-b(Y)\}^\top w }_{\text{linear in }w} + \underbrace{ \frac{1}{n}\|Y\|_2^2- \frac{1}{n}\mathrm{tr}(\Sigma) }_{\text{constant in }w}. \] The matrix $A(Y)$ is positive semidefinite because $a^\top A(Y)a=n^{-1}\|\sum_k a_k f_k(Y)\|_2^2$ for any $a\in\mathbb{R}^K$. The display therefore shows that minimizing $\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f_w)$ over the simplex is a convex quadratic program (QP). The selected weights $\hat w(Y)$ are the minimizer of this QP.
Appendix (ref) verifies Assumption (ref) for the main estimator building blocks used in the Opportunity Atlas application, and gives trained-parameter and closure tools for assembling trained candidates from those pieces. These sufficient-condition checks connect Proposition (ref) to the empirical candidate library: the application uses \ifmmodesure\elsesure\fi to average over the trained library. The within-class result, Theorem (ref), is the separate guarantee for exact \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi training over parameterized shrinkage classes satisfying the stronger uniform regularity condition; Appendix (ref) records how the reported \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi values are computed. We now turn to the Opportunity Atlas application, where the candidate maps differ in how they define which tracts are neighbors and the empirical question is whether shrinkage and averaging reduce estimated squared-error loss across commuting zones.
Does spatial shrinkage improve estimates of neighborhood economic mobility, and does averaging over spatial specifications reduce sensitivity to that choice? The Opportunity Atlas chettyOpportunityAtlasMapping2018 estimates tract-level intergenerational economic mobility for over 70,000 Census tracts in the United States. The application is motivated by evidence that economic mobility varies substantially across places and that childhood exposure to neighborhoods can affect adult outcomes chettyWhereLandOpportunity2014, chettyImpactsNeighborhoodsIntergenerational2018. The target here is narrower: estimating the latent tract-level mean of the released Opportunity Atlas outcome, not re-estimating causal exposure effects. Opportunity Atlas estimates are used to rank neighborhoods in settings such as housing mobility programs bergmanCreatingMovesOpportunity2024, so reducing estimation error can change which places are identified as high-opportunity. The released tract-level estimates are noisy measurements of latent neighborhood mobility, with reported standard errors and pronounced spatial patterning across nearby tracts. Because each tract estimate is observed only once, the empirical comparison cannot be organized around holdout performance. We therefore use \ifmmodesure\elsesure\fi, with sampling variances implied by the reported standard errors, as the common risk scale for these comparisons.
The main empirical comparison yields two findings. First, spatial shrinkage substantially improves \ifmmodesure\elsesure\fi-estimated MSE relative to non-spatial empirical Bayes baselines. Second, multiple plausible spatial specifications compete across commuting zones (CZs): geographic distance has lower \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-estimated MSE in some CZs, while contiguity (tract adjacency) has lower \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-estimated MSE in others, and using OLS to residualize tract estimates on demographic covariates before spatial smoothing can change rankings within each distance family. This heterogeneity is the reason the application reports a \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average rather than choosing a single spatial specification for all CZs. Candidates are trained, evaluated, and averaged separately within each CZ. National summaries average each CZ's \emph{\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi ratio}, the CZ-level \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-estimated MSE divided by the corresponding raw-\ifmmode\text{\textup{\textsc{mle}}}\else\textup{\textsc{mle}}\fi benchmark, weighting each CZ by its tract count. Selected averaging weights are also averaged across CZs using tract-count weights. The \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average matches or improves on the best individual candidate's \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi ratio in 16 of 20 CZs.
We study 20 CZs spanning a range of sizes ($n$ from 723 to 3,859 tracts, median 1,016; 25,777 tracts in total). Together, these CZs cover more than one-third of the U.S. population. The main outcome is pooled household income rank in adulthood for children with parents at the 25th percentile of the national income distribution (kfr_pooled_pooled; we refer to this family of children's household income-rank outcomes as KFR); in the Opportunity Atlas Table 1 data, this rank is measured in 2014--2015 for the 1978--1983 birth cohorts. Each tract $i$ has an estimate $Y_i$ and reported standard error $\mathrm{se}_i$, with variance $\sigma_i^2=\mathrm{se}_i^2$. These are the raw, unshrunk tract estimates and sampling standard errors, so $f(Y)=Y$ is the maximum-likelihood benchmark.\footnote{The empirical analysis treats the reported standard errors as fixed known sampling standard errors and does not account for uncertainty in the standard-error estimates themselves.} The Gaussian location model is $Y_i = \theta_i + \varepsilon_i$, $\varepsilon_i \sim \mathcal{N}(0, \sigma_i^2)$, with variances varying by a factor of $10$--$100$ across tracts within a CZ. This heteroskedasticity reflects differences in tract-level effective sample size and in how precisely the underlying Opportunity Atlas regressions estimate outcomes at the 25th percentile of parent income.
The empirical comparison uses two metrics to encode spatial proximity. Geographic distance is the Euclidean distance between tract centroids in longitude--latitude coordinates. Contiguity distance is the shortest-path distance on the tract adjacency graph, where two tracts are adjacent if they share a boundary or vertex, so that distance 1 means direct neighbors, distance 2 means neighbors-of-neighbors, and so on. These metrics capture different notions of spatial relatedness: contiguity can reflect administrative and social boundaries that may not align with physical distance, such as tracts separated by a river or highway.
The main comparison across CZs uses a fixed library of $K = 7$ candidate maps $Y\mapsto f_k(Y)$, summarized in Table (ref); additional variants appear only in supporting analyses. The library of candidate maps is designed to vary three empirical choices: non-spatial versus spatial pooling, geographic versus contiguity distance, and spatial smoothing with versus without covariate residualization. The candidates include non-spatial baselines (\ifmmodemle\elsemle\fi, \ifmmode\textup{\textsc{nn-eb}}\else\textup{\textsc{nn-eb}}\fi, \ifmmode\text{\textup{\textsc{close}}}\else\textup{\textsc{close}}\fi-\ifmmode\text{\textup{\textsc{gauss}}}\else\textup{\textsc{gauss}}\fi) and spatial \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi candidates that vary the distance metric and preprocessing. The OLS-preprocessed spatial candidates residualize tract estimates on four tract-level demographic covariates: percent White, percent Black, percent Hispanic, and median age. This covariate-residualization step has the same motivation as the small-area-estimation use of auxiliary covariates with noisy area-level estimates fayEstimatesIncomeSmall1979. It is used only for the OLS-labeled spatial candidates. The \ifmmode\text{\textup{\textsc{close}}}\else\textup{\textsc{close}}\fi-\ifmmode\text{\textup{\textsc{gauss}}}\else\textup{\textsc{gauss}}\fi benchmark is a Gaussian, precision-dependent EB rule in the spirit of chenEmpiricalBayesWhen2024: estimates are locally centered and scaled using weights formed from log reported variance, and the standardized values are then shrunk by the same heteroskedastic normal--normal posterior-mean formula as \ifmmode\text{\textup{\textsc{nn-eb}}}\else\textup{\textsc{nn-eb}}\fi.
The preprocessing labels describe transformations that are part of the candidate map \(Y\mapsto f_k(Y)\). For rows labeled Local NW, Nadaraya--Watson weights based on log reported variance define a local mean and scale; shrinkage is applied to the standardized estimates, and predictions are then transformed back to the original rank scale. For rows also labeled OLS, the estimates are first residualized on demographic covariates; the same Nadaraya--Watson standardization and spatial \ifmmodegp\elsegp\fi shrinkage are then applied to the residualized estimates, and the fitted covariate component is added back afterward. For \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi evaluation, automatic differentiation tracks the full map \(Y\mapsto f_k(Y)\) for each candidate, including residualization, standardization, shrinkage in the transformed space, and transformations back to ranks, conditional on the fixed Nadaraya--Watson weights and the scale floors (small constants that bound the local scale estimates away from zero). The value-similarity rule in Example (ref) appears in the Cook County comparison in Section (ref). It is not part of the main seven-candidate average across CZs. Rather than select a single row of the table, we take as the primary empirical estimator the convex average whose weights are chosen by minimizing \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi over the seven candidate maps, as in Section (ref); the reported \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi for this average evaluates the full map $\tilde f$, including both trained-parameter and selected-weight dependence.
All main empirical comparisons evaluate the trained candidate maps \(Y\mapsto f_k(Y)\), including their preprocessing and training steps, and the final \ifmmodesure\elsesure\fi-weighted convex average of those maps. For a map \(f\), the loss of interest is \(L_n(f)=n^{-1}\|f(Y)-\theta\|_2^2\), and the corresponding risk is \(R_n(f)=\mathbb{E}[L_n(f)]\). Because the latent tract-level vector \(\theta\) is unobserved, the empirical tables use \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi as the observable risk estimate. Under the Gaussian location model, if the covariance matrix used in the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi formula equals the true sampling covariance of \(Y\), then \(\mathbb{E}[\ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi_n(f)]=R_n(f)\). Thus \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi puts candidate maps on a common mean-squared-error scale under the stated sampling model.
The Opportunity Atlas reports a standard error for each tract-level regression estimate. The empirical evaluation uses these standard errors to form \[ \Sigma=\operatorname{diag}(\mathrm{se}_1^2,\ldots,\mathrm{se}_n^2), \] which treats sampling errors across tract estimates as uncorrelated. This covariance assumption concerns the estimation noise in the released tract estimates, not the spatial dependence in the latent mobility vector. We interpret each \(Y_i\) as a direct estimate of tract \(i\)'s latent mobility mean \(\theta_i\); under separate tract-level estimation with disjoint underlying observations, the reported marginal standard errors are a natural working choice for the covariance input in the Gaussian location approximation. The main remaining source of off-diagonal sampling covariance would be overlap in the underlying children contributing to multiple tract estimates, for example among movers. Appendix (ref) gives the omitted-covariance bias decomposition and shows that comparisons are most affected when candidate maps differ in how they smooth across pairs of tracts with correlated estimation errors.
The reported \ifmmodesure\elsesure\fi values are also distinct from the values of proxy \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi used to train the \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi candidates. After training, \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi is recomputed for the implemented map \(Y\mapsto f_k(Y)\), accounting for trained-parameter dependence. The same issue arises for the final convex average: the fixed-weight \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi criterion selects the weights, while the reported \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi value evaluates the full map \(\tilde f\colon Y\mapsto f_{\hat w(Y)}(Y)\), including the derivative of the selected weights with respect to the data. Under the conditions in Appendix (ref), \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi remains unbiased for the risk of \(\tilde f\). Appendix (ref) compares the reported \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi values with fixed-parameter proxies that treat trained parameters as constants. For individual candidates, this gap is the training correction, and its sign shows the optimism from ignoring trained-parameter dependence. For the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average, the gap also reflects the correction for the data-selected weights. The coupled-bootstrap procedure of oliveiraUnbiasedRiskEstimation2024 is a practical derivative-free alternative: by refitting on one perturbed sample and evaluating on its coupled counterpart, it gives an unbiased estimate of the risk for the rule trained on a variance-inflated input, approaching the original-risk target as the perturbation level shrinks. Appendix (ref) reports this comparison for one CZ. The main evaluation across CZs uses \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi because the candidate maps are differentiable and the automatic differentiation runs efficiently. Appendix (ref) gives the implementation details.
The central empirical finding is that the spatial candidates have lower \ifmmodesure\elsesure\fi-estimated MSE than the non-spatial empirical Bayes baselines in every CZ, while the best spatial specification varies across CZs. Figure (ref) plots each highlighted candidate's \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi ratio in each CZ. The highlighted series are the geographic and contiguity \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi variants, with and without OLS preprocessing, together with the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average; the gray dashed line traces the best individual candidate in each CZ.
The best geographic-distance candidate wins in 8 CZs while the best contiguity-distance candidate wins in 12. The pattern does not follow CZ size or an obvious geographic rule. Neither distance metric systematically dominates. This heterogeneity extends beyond distance metrics: within a given CZ, OLS preprocessing can also change the ranking of spatial candidates. The geographic \ifmmodegp\elsegp\fi has lower \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-estimated MSE than its contiguity counterpart in 8 CZs without OLS preprocessing, while the geographic OLS-preprocessed \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi has lower \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-estimated MSE than the contiguity OLS-preprocessed version in 7 CZs.
The pattern in Figure (ref) is the empirical case for averaging. The highlighted spatial candidates do not move in parallel across CZs: a distance metric or preprocessing choice that performs well in one CZ can lie well above the lower envelope---the best individual candidate's ratio in each CZ---in another. The \ifmmodesure\elsesure\fi-weighted average uses the same \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-estimated MSE scale within each CZ to choose convex weights over the candidate maps, rather than imposing a national choice between geographic distance, contiguity distance, and OLS preprocessing. The resulting \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average closely tracks the lower envelope, matching or improving on the best individual candidate's \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi ratio in 16/20 CZs. Thus the role of \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted averaging in the application is to retain the gains from spatial shrinkage while reducing sensitivity to which spatial specification is best in a given CZ.
Table (ref) reports tract-weighted average performance across all 20 CZs. All four spatial candidates have substantially lower \ifmmodesure\elsesure\fi-estimated MSE than the non-spatial baselines: the reduction relative to raw \ifmmode\textup{\textsc{mle}}\else\textup{\textsc{mle}}\fi ranges from 48% for the geographic \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi to 53% for the OLS-preprocessed contiguity \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi, compared to 38% for \ifmmode\text{\textup{\textsc{close}}}\else\textup{\textsc{close}}\fi-\ifmmode\text{\textup{\textsc{gauss}}}\else\textup{\textsc{gauss}}\fi.
The \ifmmodesure\elsesure\fi-weighted average has a tract-weighted average \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi ratio of 0.45, lower than every individual candidate in Table (ref). The improvement illustrates the averaging logic of Section (ref): candidate rankings differ across CZs, so minimizing \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi over convex combinations can lower \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-estimated MSE without committing to a single spatial specification. The reported table value is the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi evaluation of the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average $\tilde f$. This evaluation differentiates through both the trained candidate parameters and the selected weights; Appendix (ref) gives smoothness conditions under which it is unbiased, and the combined correction relative to the fixed-weight, fixed-parameter proxy average is $+0.007{}$ on the \ifmmode\text{\textup{\textsc{mle}}}\else\textup{\textsc{mle}}\fi-normalized scale (Table (ref)). The selected averaging weights concentrate on the spatial candidates: OLS-preprocessed contiguity \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi receives 39.0%, OLS-preprocessed geographic \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi receives 27.6%, and geographic \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi receives 20.1%, with the remaining weight spread across the other candidates.
As a supporting comparison, Figure (ref) reports what happens when value-similarity smoothing is added in one setting where nearby neighborhoods differ sharply: Cook County tracts selected from the Chicago CZ. The stepwise comparison begins with non-spatial baselines, then adds the geographic \ifmmodegp\elsegp\fi, the OLS-preprocessed geographic \ifmmode\textup{\textsc{gp}}\else\textup{\textsc{gp}}\fi, and finally \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi-\ifmmode\text{\textup{\textsc{bilat}}}\else\textup{\textsc{bilat}}\fi. The \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi-\ifmmode\text{\textup{\textsc{bilat}}}\else\textup{\textsc{bilat}}\fi candidate implements the value-similarity idea from Example (ref): starting from OLS-preprocessed geographic smoothing, it gives more weight to nearby tracts whose observed mobility estimates are similar. On the same \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-ratio scale, adding \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi-\ifmmode\text{\textup{\textsc{bilat}}}\else\textup{\textsc{bilat}}\fi lowers both ratios: the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average falls from $0.519$ to $0.507$, and the best individual ratio falls from $0.544$ for the OLS-preprocessed geographic \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi to $0.522$ once \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi-\ifmmode\text{\textup{\textsc{bilat}}}\else\textup{\textsc{bilat}}\fi is added. The \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average places weight $0.655$ on \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi-\ifmmode\text{\textup{\textsc{bilat}}}\else\textup{\textsc{bilat}}\fi in the final comparison step. This focused comparison shows that, in a setting where nearby neighborhoods differ sharply, adding value similarity lowers \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-estimated MSE further.
The \ifmmodesure\elsesure\fi-estimated-MSE comparisons evaluate accuracy of the estimated mobility vector, but these estimates are often used as inputs into ranking and selection decisions, a compound-decision setting studied by guKoenkerInvidiousComparisons2023. Related work studies inference for ranks and selected high-opportunity neighborhoods from noisy Opportunity Atlas estimates mogstadInferenceRanksApplications2024, andrewsInferenceWinners2024. We use the shrinkage estimates in a targeting exercise under a top-third rule, comparing which tracts each rule selects and the selected group's average latent mobility rank. Within each CZ, each rule ranks tracts by the estimated outcome, selects the top third, and estimates that group's average latent mobility rank using an evaluation strategy motivated by the coupled-bootstrap procedure of oliveiraUnbiasedRiskEstimation2024; Appendix (ref) describes the implementation and Appendix (ref) gives the corresponding unbiasedness calculation. The targeting library contains four trained maps: \ifmmode\textup{\textsc{mle}}\else\textup{\textsc{mle}}\fi, \ifmmode\text{\textup{\textsc{nn-eb}}}\else\textup{\textsc{nn-eb}}\fi, \ifmmode\text{\textup{\textsc{close}}}\else\textup{\textsc{close}}\fi-\ifmmode\text{\textup{\textsc{gauss}}}\else\textup{\textsc{gauss}}\fi, and the geographic-distance \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi candidate from Table (ref), which uses Local NW preprocessing and no OLS residualization. Table (ref) reports the three non-\ifmmode\text{\textup{\textsc{mle}}}\else\textup{\textsc{mle}}\fi maps and the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average formed from all four, with estimated gains measured relative to the raw-\ifmmode\text{\textup{\textsc{mle}}}\else\textup{\textsc{mle}}\fi targeting rule for the pooled outcome and three subgroup outcomes. \ifmmode\text{\textup{\textsc{mle}}}\else\textup{\textsc{mle}}\fi enters as the zero benchmark rather than as a separate row. Dollar-equivalent gains use the official Opportunity Atlas 2015 percentile-dollar crosswalk, so they should be read as an interpretation of rank gains rather than as a separately estimated dollar outcome.
Appendix (ref) reports the same four-candidate comparison with the geographic-distance \ifmmodegp\elsegp\fi on the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-ratio scale for these related KFR outcomes. In Table (ref), the largest absolute targeting gain occurs for the Black-male outcome: the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average improves top-third targeting by an estimated 1.21 rank points, or about \$1,265 in dollar-equivalent terms. The table also separates overall shrinkage gains from the incremental gain of geographic smoothing over the non-spatial \ifmmode\text{\textup{\textsc{close}}}\else\textup{\textsc{close}}\fi-\ifmmode\text{\textup{\textsc{gauss}}}\else\textup{\textsc{gauss}}\fi benchmark: relative to \ifmmode\text{\textup{\textsc{close}}}\else\textup{\textsc{close}}\fi-\ifmmode\text{\textup{\textsc{gauss}}}\else\textup{\textsc{gauss}}\fi, the geographic-distance \ifmmode\text{\textup{\textsc{gp}}}\else\textup{\textsc{gp}}\fi adds more for the pooled-male and White-male outcomes than for the Black-male outcome. \FloatBarrier
We develop \ifmmodesure\elsesure\fi-based model averaging, and the selection guarantees that underpin it, for shrinkage maps that exploit spatial structure. The risk comparisons are among the resulting maps, so the prior distribution, prior covariance structure, or similarity rule used to motivate a candidate map need not be correctly specified as a model for \(\theta\). The theory gives sufficient conditions for two uses of \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi with nonlinear shrinkage maps that have cross-unit dependence. One result covers selection within a parameterized class by minimizing \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi. The model-averaging result shows that once a candidate library has been assembled, the fixed-weight \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi criterion can average its members with a finite-library oracle guarantee.
Empirically, in the main pooled-outcome comparison, every spatial candidate has a lower \ifmmodesure\elsesure\fi ratio than both non-spatial empirical Bayes baselines in every one of the 20 CZs. The best individual spatial specification varies with local geography, and the \ifmmode\textup{\textsc{sure}}\else\textup{\textsc{sure}}\fi-weighted average of candidate maps has a \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi ratio of 0.45. A supporting top-third targeting exercise finds that the \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi-weighted average selects tracts with higher estimated average income ranks than the raw-\ifmmode\text{\textup{\textsc{mle}}}\else\textup{\textsc{mle}}\fi targeting rule across the reported outcomes.
For applications with several plausible notions of similarity, the practical lesson is to build a diverse library of shrinkage estimators and let \ifmmodesure\elsesure\fi evaluate and average them, rather than committing to a single exchangeable model ex ante. The same issue arises for noisy area, school, hospital, or firm-level estimates whenever researchers have several credible ways to pool information across units. Targeting and other decision-focused analyses remain important chenCompoundSelectionDecisions2025. This paper instead applies \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi to average shrinkage maps for estimating the latent vector, with targeting treated as a downstream application of those estimates. \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi's unbiasedness leans on the Gaussian sampling model for the noise in $Y$. For tract-level regression predictions built from many observations this is an approximation we consider mild. In the application, \ifmmode\text{\textup{\textsc{sure}}}\else\textup{\textsc{sure}}\fi is computed treating the reported standard errors as known and the sampling errors across tracts as uncorrelated. The main threat to that covariance assumption is overlap in the underlying children contributing to multiple tract estimates, and Appendix (ref) characterizes what that misspecification costs.