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.
42,467 characters · 9 sections · 37 citation commands
Spatially Robust Inference with Predicted and Missing at Random Labels
\onehalfspacing
{\it Inference with predicted data} arises when data are sparsely labeled but predicted labels are available and used for statistical inference. This is increasingly common in global health surveillance fan2024narratives, where limited audits are combined with full-coverage prediction maps bhatt2015malaria,golding2017child,di2017air, land-use monitoring, which relies on dense remote-sensing predictions with sparse field verification bullock2020amazon,sexton2013treecover, and citizen-science projects that pair scalable automated classifiers with limited expert labels willett2013galaxyzoo2. Here, the scientific goal is typically not individual-level predictions but valid uncertainty quantification for population-level statistics salerno2025we.
Recent works develop methods for drawing inference with predicted labels wang2020inference,angelopoulos2023ppi, gronsbell2024another,zrnic2024crossppi, but these approaches are almost uniformly derived under independent and identically distributed (iid) labeling. Real-world settings routinely violate that assumption in two important ways. First, label acquisition is rarely done at random (i.e., outcome labels are rarely missing completely at random, or MCAR). Labeling commonly depends on observed features and geography, inducing a missing at random (MAR) mechanism. Second, observations are often spatially dependent. Both issues matter for inference. Under MAR, naive estimates can be biased when labeling probabilities are correlated with certain predictors. Further, under spatial dependence, classical variance estimators break down even when point estimators remain approximately unbiased. Addressing both sources of error, selection bias from covariate-dependent labeling and dependence in the data-generating process, is essential for valid inference.
Our framing separates the problem into two components: a {\it base prediction} and a subsequent stage of {\it local nuisance correction}. The base prediction model, $f: \mathcal{X}\to\mathcal{Y}$, produces predicted labels $\widehat{Y}_i = f(X_i)$ for every unit based on predictors, $X_i$. In many applications, this model is trained externally and enters the analysis as a deterministic function. Further, $f$ may be arbitrarily complex, and we do not require any knowledge of its training procedure. For inference, the key point is that $\widehat{Y}$ is generally not a selection-robust substitute for $Y$ under MAR. Correcting both prediction bias and locally varying missingness therefore requires estimating two nuisances using our analytic sample: an outcome model, $\hat m(W_i,s_i)$, and a label propensity model, $\hat\pi(W_i,s_i)$, where $W_i=(X_i,\widehat{Y}_i)$ and $s_i$ are spatial coordinates.
Estimating these nuisance functions in practice is challenging because only a small fraction of units are labeled, making cross-fitting necessary to prevent overfitting chernozhukov2018dml. In spatially dependent settings, however, cross-fitting introduces an unintended consequence: units evaluated using the same nuisance estimates inherit shared estimation noise, creating artificial correlations in the resulting estimator. When standard spatial variance estimators are applied directly, this induced correlation can be mistaken for genuine spatial dependence, leading to unstable or miscalibrated confidence intervals.
We propose a jackknife-based variance correction that separates genuine spatial dependence from artificial correlation introduced by cross-fitting. The method removes fold-induced components from our doubly robust (DR) estimator before applying a Conley-style spatial heteroscedasticity and autocorrelation consistent (HAC) covariance estimator and then restores between-fold variability through an ANOVA-style adjustment conley1999gmm. This jackknife-HAC construction leaves the DR point estimator unchanged, modifying only the variance calculation, making the approach modular and easy to adapt to existing workflows.
Existing methods address individual aspects of this problem. Many approaches for inference with predicted data assume independent sampling salerno2025moment,angelopoulos2023ppi,zrnic2024crossppi,wang2020inference. Doubly robust estimators accommodate covariate-dependent labeling in classical missing data and semiparametric settings bang2005doubly,zeng2010adjustment,butera2022doubly,little2019missing. Spatial econometrics provides methods for inference under dependence, including Conley-type covariance estimators and random-field central limit theorems conley1999gmm,bester2011cluster,jenish2009clt,chandrasekhar2023affinity. Our contribution is to integrate these ideas into a valid approach and to address a specific issue that arises when nuisance functions are cross-fitted on spatially dependent data. We focus on mean estimation for clarity, the same principles extend to other smooth functionals via influence functions. While our approach addresses spatial dependence with missing outcomes, chen2025unified provide a complementary framework for independent data, utilizing pattern-stratified Z-estimation to guarantee valid inference with machine learning imputations under general missing-at-random (MAR) mechanisms across both outcomes and covariates.
Note that {\it externally supplied predictions} correspond to settings where the base prediction model is trained on auxiliary data and then treated as fixed in the analytic sample. When predictions are instead {\it learned from the current dataset}, they are generated out-of-fold to avoid reusing labeled outcomes in both prediction and subsequent correction (see Assumption (ref)). The proposed jackknife-HAC correction targets a different phenomenon: fold-shared training noise arising from cross-fitted local nuisance estimation (Section (ref)), rather than dependence induced by fitting the base prediction model itself. Empirical illustrations of this distinction appear in our simulation mechanism (Section (ref)) and our mechanism ablation study (Appendix (ref)).
In the following, Section (ref) formalizes the estimand and identification under MAR labeling and overlap. Section (ref) introduces our estimator, describes cross-fitting the nuisances, and presents the jackknife-HAC variance formula. Section (ref) states high-level asymptotic normality results under dependent-array CLT conditions. Section (ref) reports simulated and real data results, along with diagnostic analyses.
For units $i=1,\ldots,n$, let $Y_i \in \mathbb{R}$ denote our outcome of interest and $R_i\in\{0,1\}$ indicate whether $Y_i$ is observed. Let $X_i$ collect observed predictors and $s_i\in\mathbb{R}^2$ be a location or spatial dependence coordinate. We assume a predicted outcome, $\widehat{Y}_i$, for every unit is either by {\it externally supplied} or {\it learned from the current dataset} (see Section (ref)). For brevity, let $W_i=(X_i,\widehat{Y}_i)$. The observed record is \[ O_i = (W_i,s_i,R_i,Y_i^{\rm obs}), \] where $Y_i^{\rm obs} = Y_i$ if $R_i = 1$ and is missing if $R_i = 0$.
For simplicity, consider the target as the population mean, \[ \theta_0 = \mathbb{E}[Y_i], \] though this extends to other smooth functionals. We assume $\mathbb{E}[Y_i]=\theta_0\ \forall i$ and $\mathbb{E}[Y_i^2]<\infty$. Define nuisance functions \[
\]
Assumption (ref) is standard in the missing data literature and allows label availability to depend on observed covariates and spatial location. The overlap condition ensures that every covariate-location stratum contributing to the target mean has a nonzero probability of being labeled. Without overlap, valid inference for the global mean requires either changing the estimand (e.g., restricting to regions with sufficient labeling) or imposing additional assumptions.
Based on scharfstein1999adjusting, for any measurable pair $(m,\pi)$, we define the doubly robust estimating function \[ \psi_i(\theta; m,\pi) = m(W_i,s_i)+\frac{R_i}{\pi(W_i,s_i)}\{Y_i-m(W_i,s_i)\} - \theta. \]
Proof of Proposition (ref) can be found in Appendix (ref).
The intuition for $\psi_i$ is that the outcome model term $m_i$ is a calibrated function of the predictions for each unit. The correction term $R_i(Y_i- m_i)/\pi_i$ adds back observed residual information with an inverse-probability of labeling weight so that, under MAR, the labeling bias is removed when either the outcome model is correct or the propensity is correct. Dependence correction enters only through the variance, as point identification is a MAR problem, while calibrated uncertainty is a dependence problem.
The bias of iid methods can be written directly. For exposition only, assume a common marginal label rate, $\bar\pi=\mathbb{E}[R_i]$, and consider one such method zrnic2024crossppi \[
\] under MAR labeling. The second term captures selection bias induced by nonrandom labeling: unless the prediction residual, $Y_i-m(W_i,s_i)$, is uncorrelated with the labeling probability, $\pi_0(W_i,s_i)$, there is bias. This representation clarifies why such methods perform well in MCAR settings but fail under MAR, particularly when labeling probabilities vary with features or locations associated with prediction error. Our MAR experiments explicitly construct such dependence by letting label propensities depend on both coordinates and observed covariates.
Our doubly robust estimator solves the sample analog of Proposition (ref), but the variance is modified to respect the fold-level noise induced by the cross-fit nuisance functions.
\paragraph{Base Prediction Model and Nuisance Functions.} We assume a predicted outcome, $\widehat{Y}_i$, for every unit is produced by a base prediction model that may be trained externally and treated as fixed in the analytic sample. We absorb $\widehat{Y}_i$ into the covariate vector, $W_i=(X_i,\widehat{Y}_i)$. The nuisance functions \[
\] must be estimated within the analytic sample to correct local miscalibration and local missingness patterns.
\paragraph{Cross-Fit Local Nuisance Estimation.} We partition the data into $K$ folds, $\mathcal{I}_1,\ldots,\mathcal{I}_K$. For fold $k$, we estimate $\hat m_{-k}$ and $\hat\pi_{-k}$ on the other $K-1$ folds and evaluate them on $k$: \[ \hat m_i=\hat m_{-k}(W_i,s_i),\quad \hat\pi_i=\hat\pi_{-k}(W_i,s_i), \quad i\in \mathcal{I}_k. \] To reduce train-test leakage under spatial dependence, we use a buffered set based on distances $d_{ij}=\lVert s_i-s_j\rVert$. Let \[
\] If $|\mathcal{T}_k|$ is too small, we use standard jackknife folds. All propensity predictions are clipped to $[\pi_{\min},1-\pi_{\min}]$ for a fixed $\pi_{\min}\in(0,1/2)$ (in experiments, $\pi_{\min}=0.10$).
\paragraph{DR Point Estimator.} Define the uncentered DR estimate \[ \hat\psi_i^\circ = \hat m_i+\frac{R_i}{\hat\pi_i}(Y_i-\hat m_i), \quad \hat\theta=\frac{1}{n}\sum_{i=1}^n \hat\psi_i^\circ. \]
\paragraph{Naive Spatial HAC Instability.} Cross-fitting prevents overfitting the labeled subset, but it has a side effect under spatial dependence: all units in a fold, $k$, share the same fitted nuisance functions, $(\hat m_{-k},\hat\pi_{-k})$, so their estimates $\hat\psi_i^\circ$ share {\it fold-level noise}. A stylized decomposition is \[ \hat\psi_i^\circ=\theta_0+\phi_{0i}+a_{k(i)}+\varepsilon_{n,i}, \] where $\phi_{0i}=\psi_i(\theta_0;m_0,\pi_0)$ is the centered score function, $a_{k(i)}$ is a fold-level random effect induced by nuisance estimation in settings of spatial dependence, and $\varepsilon_{n,i}$ is a remainder. For $i$ and $j$ in the same fold, $\mathrm{Cov}(\hat\psi_i^\circ,\hat\psi_j^\circ)$ includes $\mathrm{Var}(a_{k(i)})$ regardless of distance. Since a nontrivial fraction of near-neighbor pairs lie in the same fold, a distance-weighted Conley HAC estimator can misattribute this shared nuisance noise to genuine short-range spatial dependence, inflating standard errors in finite samples.
\paragraph{Jackknife-HAC Variance.} Let $\bar\psi_k=n_k^{-1}\sum_{i\in\mathcal{I}_k}\hat\psi_i^\circ$, $n_k=|\mathcal{I}_k|$, and define fold-centered scores $\tilde\psi_i=\hat\psi_i^\circ-\bar\psi_{k(i)}$. Let $\kappa(\cdot)$ be a triangular kernel and set the HAC bandwidth $h_n$ by a distance quantile: \[
\] In experiments we choose $h_n$ (and the buffer radius $r_n$) by distance quantiles, $h_q$, for simplicity and reproducibility. In the asymptotic statement, we treat these as deterministic sequences that satisfy standard spatial HAC and buffered-splitting regularity conditions. Define weights $w_{ij}=\kappa(d_{ij}/h_n)$. The within-fold covariance is a Conley HAC quadratic form applied to $\tilde\psi$: \[ \hat V_{\mathrm{within}} =\frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^n w_{ij}\tilde\psi_i\tilde\psi_j. \] We keep only the off-diagonal contribution and add the between-fold term computed from the fold means: \[ \hat V_{\mathrm{off}} =\hat V_{\mathrm{within}}-\hat V_{\mathrm{diag}}, \quad \hat V_{\mathrm{diag}}=\frac{1}{n^2}\sum_{i=1}^n \tilde\psi_i^2, \] \[ \hat V_{\mathrm{between}} =\frac{K}{K-1}\sum_{k=1}^K\left(\frac{n_k}{n}\right)^2(\bar\psi_k-\hat\theta)^2. \] The jackknife-HAC variance estimator is \[ \hat V_{\mathrm{JK}} =\hat V_{\mathrm{off}}+\hat V_{\mathrm{between}}. \] Intuitively, fold centering removes fold-shared components such as $a_{k(i)}$ from the spatial covariance calculation, while the between-fold term reintroduces fold-mean variation only through the fold averages. The diagonal subtraction prevents counting idiosyncratic (diagonal) variation twice. In our implementation, we further let $\hat V_{\mathrm{JK}}^{+}=\max(\hat V_{\mathrm{JK}},\epsilon)$ with small $\epsilon>0$ for numerical stability. A two-sided $1-\alpha$ confidence interval is then given by \[ \hat\theta\pm c_{1-\alpha/2}\sqrt{\hat V_{\mathrm{JK}}^{+}}, \] where $c_{1-\alpha/2}=z_{1-\alpha/2}$ in the spatial branch.
\paragraph{Optional Moran Gate.} In implementation, we optionally gate the within-fold covariate by a labeled residual Moran test: compute Moran's $I$ moran1950 on $\{Y_i-\hat m_i:R_i=1\}$. If the test does not reject a spatial signal, we replace the within-fold HAC covariance by an iid covariance matrix and use a Student's $t_{K-1}$ critical value for $c_{1-\alpha/2}$. This gate is a pragmatic stability heuristic. Our analytic results treat $\hat V_{\mathrm{JK}}^{+}$ abstractly through a high-level consistency assumption.
\paragraph{Algorithm Summary (One Replicate).} Our procedure is as follows: (i) obtain the base predicted outcomes, $\widehat{Y}_i$; (ii) choose folds and buffer, $q_b$, then cross-fit local nuisance models, $\hat m(W,s)$ and $\hat\pi(W,s)$, on the analytic sample; (iii) compute $\hat\theta=n^{-1}\sum_i \hat\psi_i^\circ$; (iv) compute $\hat V_{\mathrm{JK}}^{+}$; and (v) report \[ \hat\theta\pm c\sqrt{\hat V_{\mathrm{JK}}^{+}}. \]
\paragraph{Asymptotic Fold Scaling.} For the high-level asymptotic statement in Section (ref) we allow the number of folds $K=K_n$ to increase with $n$ (Assumption (ref)), which makes the ANOVA term stable in large samples. In finite samples we use small $K$ for feasibility with limited labels.
Two estimand conventions should be distinguished in finite-sample reporting. Our analytic results target the superpopulation mean, $\theta_0=\mathbb{E}[Y_i]$, which treats $(W_i,s_i,Y_i)$ as random draws from a dependent array. Here, our main results use superpopulation-style coverage. In synthetic designs, this is for a known data generation process mean, $\mu$, but in observed benchmark resampling, this is a fixed full-pool mean used across repeated draws. Finite-target coverage against each realized sample mean is reported only as a comparison in the Appendix, where it is typically more attenuated. This convention follows the finite-versus-superpopulation distinction in modern causal inference jin2024tailored and is conceptually parallel to the fixed-versus-random spatial structure distinction emphasized in spatial confounding gilbert2025consistency.
Let $\phi_{0i}=\psi_i(\theta_0;m_0,\pi_0)$ and $Z_{n,i}=\phi_{0i}-\mathbb{E}[\phi_{0i}]$. Our asymptotic interpretation follows increasing-domain spatial asymptotics: as $n\to\infty$, the spatial domain expands, local dependence neighborhoods remain bounded in physical scale while long-range covariance accumulation remains weak enough for HAC stabilization.
Let $\mathcal{A}_n(i)\subset\{1,\ldots,n\}$ be affinity sets with $i\in\mathcal{A}_n(i)$, instantiated as distance neighborhoods $\mathcal{A}_n(i)=\{j:d_{ij}\le h_n\}$, and define \[ \Omega_n=\sum_{i=1}^n\sum_{j\in\mathcal{A}_n(i)}\mathrm{Cov}(Z_{n,i},Z_{n,j}). \] We take $h_n$ to be the same bandwidth sequence used by the spatial covariance in Section (ref).
Assumption (ref) allows the number of folds $K_n$ to grow while each fold remains large, which stabilizes the between-fold term in $\hat V_{\mathrm{JK}}$. Assumption (ref) separates external- versus within-sample prediction regimes; when predicted outcomes are generated from current labels, cross-fitting is needed to avoid own-observation reuse. Assumption (ref) is a high-level orthogonality remainder condition standard in debiased machine learning with cross-fitting chernozhukov2018dml. In iid settings, cross-fitting plus bounded moments and the DML product-rate condition \[ \Big(\frac{1}{n}\sum_{i=1}^n(\hat m_i-m_{0i})^2\Big)^\frac{1}{2}\Big(\frac{1}{n}\sum_{i=1}^n(\hat\pi_i-\pi_{0i})^2\Big)^\frac{1}{2}=o_p(n^{-\frac{1}{2}}) \] is sufficient for Assumption (ref) chernozhukov2018dml. Under dependence, we impose Assumption (ref) directly rather than deriving it from primitive mixing conditions. Importantly, the identification-level double robustness in Proposition (ref) does not by itself imply Assumption (ref): asymptotic normality requires nuisance errors to be sufficiently small in the sense of the remainder bound above. Assumption (ref) concerns the jackknife-HAC variance estimator; sufficient conditions combine standard kernel, bandwidth, and covariance-decay conditions for Conley HAC conley1999gmm,bester2011cluster with fold regularity and the requirement that centering the folds removes only asymptotically negligible fold-shared nuisance noise.
For the complete proof, see Appendix (ref).
This section addresses one question: under MAR or MCAR designs with sufficient overlap, does dependence-aware inference with predicted labels improve coverage relative to iid baselines, and how does that pattern change as dependence and sampling concentration increase?
We compared our approach (referred to as Spatial DR-JK-HAC) to the following methods: Cross-PPI zrnic2024crossppi, PPI++ angelopoulos2023ppi++, and Bootstrap-PPI efron2026boostrap. These serve as example baselines, as they assume independent labeling and iid sampling. Spatial DR-JK-HAC is our doubly robust point estimator paired with a jackknife-HAC variance designed for MAR labeling and spatial dependence. Additional DR baselines (DR-iid and spatial DR-HAC) are reported in the Appendix as robustness comparisons; they follow the same DR point estimator but pair it with an iid variance or a direct Conley-style HAC variance, respectively.
We generate a synthetic spatial population on an $n_{\mathrm{side}}\times n_{\mathrm{side}}$ grid (here $n_{\mathrm{side}}=250$, so $N=62{,}500$ population units). Each population replicate draws smoothed Gaussian random fields on this grid. The dependence level $\sigma$ is the smoothing radius (in grid cells); larger $\sigma$ induces longer-range spatial correlation. To keep the marginal scale comparable across $\sigma$, we center and rescale each $\sigma$-smoothed field to have approximately zero mean and unit standard deviation.
Each unit has observed features $W=(X,U_{\mathrm{obs}},\mathrm{coords})$, where $X$ is a weakly dependent field (fixed smoothing radius $2$) and $\mathrm{coords}$ are normalized coordinates. Outcomes follow \[ Y = \mu + \beta X + \lambda U_{\mathrm{obs}} + \lambda U_{\mathrm{unobs}} + \varepsilon, \] with $\mu=0$, $\beta=0.8$, $\lambda=1.0$, and $\varepsilon\sim\mathcal{N}(0,0.6^2)$. $U_{\mathrm{unobs}}$ is a second $\sigma$-smoothed field not included in $W$ to create spatial dependence that cannot be removed by conditioning on $W$, so predictions trained on $W$ are intentionally imperfect at high dependence.
Predictions, $\hat Y$, are based on an auxiliary training set xu2025unifiedframeworksemiparametricallyefficient. Throughout, we use an {\it auxiliary-split} construction: we split the population into an auxiliary training pool (35%) and an analysis pool (65%), fit a gradient-boosted tree predictor of $Y$ from standardized $W$ on the auxiliary pool, and evaluate $\hat Y$ only on the disjoint analysis pool. We then draw repeated analysis samples of size $n=600$ from the analysis pool under either iid sampling or soft-block sampling. In soft-block sampling, a small “core” of size $0.05n$ is taken as nearest neighbors of a random anchor point, and the remainder is filled iid from the remaining pool, making dependence more visible to the sampling distribution, but without making global superpopulation inference impossible.
Given a drawn sample, we take a fixed label budget of 20% either MCAR or via an adversarial MAR mechanism that concentrates labels away from high-uncertainty regions while respecting $\pi\in[0.10,0.90]$ (see Appendix (ref)). To keep the simulation focused on dependence-driven calibration, we use the true $\pi$ when forming the DR estimator.
Figure (ref) reports the 90% coverage probabilities for $\sigma\in\{0,40,80,120\}$ (100 population replicates per $\sigma$, 200 sample draws per population). Under iid sampling and MCAR labeling, the baseline methods have near nominal coverage, even as $\sigma$ grows (e.g., at $\sigma=120$, Cross-PPI and PPI++ cover at $0.90$ and $0.90$). Under iid sampling and MAR, these methods under-cover (at $\sigma=120$, Cross-PPI: $0.46$, PPI++: $0.57$, Bootstrap-PPI: $0.33$), reflecting MAR-induced bias. Under soft-block sampling, these methods under-cover as dependence grows, even in MCAR (at $\sigma=120$, Cross-PPI: $0.81$, PPI++: $0.80$, Bootstrap-PPI: $0.65$), because spatial dependence is visible in the sampled data. Spatial DR-JK-HAC preserves near-nominal coverage across all four cells (at $\sigma=120$, block-MCAR: $0.91$, block-MAR: $0.90$).
Appendix (ref) summarizes tradeoffs in CI widths. At $\sigma=120$ under soft-block sampling, Spatial DR-JK-HAC intervals are $1.36\times$ wider than Cross-PPI in MCAR ($0.39$ vs $0.28$ mean length) and $1.57\times$ wider in MAR ($0.43$ vs $0.27$).
We evaluated five observed benchmark datasets, three from zrnic2024crossppi: (i) {\it forest disturbance}, a geolocated Amazon deforestation dataset with remote-sensing covariates bullock2020amazon,sexton2013treecover; (ii) {\it galaxy zoo morphology}, based on Galaxy Zoo 2 labels with associated SDSS imaging metadata willett2013galaxyzoo2,york2000sdss; and (iii) {\it census income}, constructed from the American Community Survey via Folktables ding2021folktables, and two harder spatial benchmarks: (iv) {\it malaria} (Colombian municipalities; outcome is a transformed malaria burden variable), and (v) {\it Health+} (a tract-level life-expectancy outcome with socioeconomic covariates). In each dataset, we repeatedly draw $n=320$ observations under iid and soft-block sampling and mask labels under MCAR or MAR at $20\%$ budget. Coverage uses the fixed full-pool benchmark mean as the repeated-sampling target (the observed-data analog of a superpopulation-style target in this design). All observed-data results use the same auxiliary-split prediction construction described above, with 40 replicates per dataset-arm-sampling cell for $K=5$ methods.
Table (ref) summarizes coverage for the observed data under this auxiliary-split prediction design, averaged over iid and soft-block sampling. Averaged across datasets in MAR, coverage is $0.665$ (Cross-PPI), $0.726$ (PPI++), $0.571$ (Bootstrap-PPI), $0.810$ (DR-iid), $0.816$ (Spatial DR-HAC), and $0.874$ (Spatial DR-JK-HAC, $K=5$). Forest and Malaria show the largest MAR gains for Spatial DR-JK-HAC ($K=5$, $t$) relative to Spatial DR-HAC. In MCAR, DR-based methods are generally near or above nominal, while Bootstrap-PPI remains low in several datasets.
CI-width differences are directionally consistent with the synthetic design. In MAR settings, Spatial DR-JK-HAC (K=5, $t$) has about $1.38\times$ the mean width of Spatial DR-HAC when averaging widths across datasets. Appendix (ref) reports K=10 fold-sensitivity and width tables: K=10 ($t$) reduces width inflation relative to K=5 while retaining higher coverage than Spatial DR-HAC in most datasets, and K=10 ($z$) is close to Spatial DR-HAC in both coverage and width.
We also ran targeted diagnostic checks to isolate mechanisms and rule out obvious analysis artifacts. First, we validate coordinate proxies by residual Moran and nearest-neighbor diagnostics rather than narrative caveats, and we report that proxy quality is strongest for forest and weaker for galaxies and census (Appendix (ref)). Second, we report overlap diagnostics directly, including design and estimated propensity quantiles, clipping rates, and inverse-propensity effective sample size (Appendix (ref)). Third, we provide a mechanism attribution check that isolates the roles of DR correction, dependence-aware variance estimation, and the jackknife-HAC correction in the hardest synthetic cell (Appendix (ref)). Replicate-level interval visualizations that make MAR centering failures visible are shown in Appendix (ref).
Under Assumption (ref), the point estimator in this paper is deliberately classical: it is the standard doubly robust (DR) estimator for missing outcomes bang2005doubly,little2019missing. The novelty is in the {\it covariance correction} required once the DR nuisance functions are estimated locally in the analytic sample with cross-fitting. Existing methods are built for independent designs and typically do not address MAR labeling plus spatial dependence in one end-to-end pipeline angelopoulos2023ppi++,zrnic2024crossppi. Spatial HAC estimators are well developed conley1999gmm,bester2011cluster, but applying them directly to cross-fitted DR scores can be unstable because fold members share the same fitted nuisance functions.
This paper isolates and addresses these issues. Cross-fitting is essential with small label budgets, but it induces fold-shared training noise in $\hat m$ and $\hat\pi$, hence, in the estimated score functions, $\hat\psi_i^\circ$. When nearby units share a fold, this noise can masquerade as short-range spatial correlation and inflate a Conley HAC variance estimate. The jackknife-HAC variance in Section (ref) attenuates this fold-shared component by centering within folds and recombining within- and between-fold variation. This correction leaves the identification score unchanged and targets the variance estimation.
Our assumptions are intended to match the logic used in Conley-style spatial econometrics. Researchers typically do not require global independence; they require that dependence be predominantly local and that long-range covariance accumulation be weak enough for CLTs and HAC estimators to stabilize conley1999gmm,bester2011cluster,jenish2009clt. The overlap requirement is the standard positivity condition from missing data analysis, interpreted here as a design requirement for audit allocation: every covariate-location stratum that contributes to the target mean must have nonzero probability of manual labeling.
The limitations are also clear. First, if overlap for the global mean is badly violated, valid confidence intervals are generally impossible without changing the estimand or adding structure; overlap-restricted estimands, adaptive trimming, or partial identification become appropriate alternatives crump2009overlap. Second, when dependence coordinates are proxy embeddings rather than true spatial coordinates, both HAC weighting and Moran-style diagnostics can be misaligned with the actual dependence structure. Third, the jackknife-HAC correction is designed for the practical regime of small labeled fractions and cross-fitted nuisance estimation; a complete theory for its finite-sample behavior under spatial dependence is beyond the current scope and is a natural target for future work.
Finally, this framework is modular: the same DR score can be paired with other dependence-appropriate covariance structures beyond spatial HAC. Appendix (ref) illustrates this with a jackknife-corrected two-way clustered variance estimator for overlapping dependence structures cameron2011robust,chandrasekhar2023affinity.
We propose a method for spatially robust inference with predicted and missing at random labels. With cross-fitted nuisance estimation in the analytic sample, the standard DR estimator inherits fold-shared training noise, and naive Conley HAC variance estimation can misinterpret this noise as spatial dependence. We proposed a jackknife-HAC variance estimator with standard DR estimator. Under standard identification, dependent-array CLT conditions, and high-level nuisance and variance consistency assumptions, the resulting intervals are asymptotically valid.
Empirically, the jackknife-HAC correction can be evaluated across label budgets, fold choices, and coordinate quality. Theoretically, a natural direction is deriving the conditions under which this jackknife-HAC approach is asymptotically equivalent to HAC, while improving finite-sample robustness when labels are limited. More broadly, this approach is modular and can be paired with other dependence-matched covariance structures (e.g., Newey-West for time series newey1987simple or multiway clustering for overlapping dependence cameron2011robust,chandrasekhar2023affinity). These results suggest a broader framework for dependence-robust inference when machine learning predictions supplement sparse labels.