Smooth Sailing: Lipschitz-Driven Uncertainty Quantification for Spatial Associations
Authors: David R. Burt
Renato Berlinghieri
Stephen Bates
Tamara Broderick
Venue: NeurIPS 2025 poster ()
Keywords: Spatial
Confidence Intervals
Linear Regression
Lipschitz
Uncertainty Quantification
Primary area: probabilistic_methods
TL;DR: We provide a method for reliable uncertainty quantification for spatial associations in the face of model misspecification and nonrandom spatial locations.

Abstract: Estimating associations between spatial covariates and responses — rather than merely predicting responses — is central to environmental science, epidemiology, and economics. For instance, public health officials might be interested in whether air pollution has a strictly positive association with a health outcome, and the magnitude of any effect. Standard machine learning methods often provide accurate predictions but offer limited insight into covariate-response relationships. And we show that existing methods for constructing confidence (or credible) intervals for associations can fail to provide nominal coverage in the face of model misspecification and nonrandom locations — despite both being essentially always present in spatial problems. We introduce a method that constructs valid frequentist confidence intervals for associations in spatial settings. Our method requires minimal assumptions beyond a form of spatial smoothness and a homoskedastic Gaussian error assumption. In particular, we do not require model correctness or covariate overlap between training and target locations. Our approach is the first to guarantee nominal coverage in this setting and outperforms existing techniques in both real and simulated experiments. Our confidence intervals are valid in finite samples when the noise of the Gaussian error is known, and we provide an asymptotically consistent estimation procedure for this noise variance when it is unknown.

Smooth Sailing: Lipschitz-Driven Uncertainty Quantification for Spatial Association
David R. Burt"
Renato Berlinghieri'
Stephen Bates
Tamara Broderick
MIT LIDS
(dburt,reni) s_bates, tamarab) @mit. edu
Abstract
Estimating associations between spatial covariates and responses - rather than merely predicting responses - is central to environmental science, epidemiology, and economics. For instance, public health officials might be interested in whether air pollution has a strictly positive association with a health outcome, and the mag- nitude of any effect. Standard machine learning methods often provide accurate predictions but offer limited insight into covariate-response relationships. And we show that existing methods for constructing confidence (or credible) intervals for associations can fail to provide nominal coverage in the face of model misspecifi- cation and nonrandom locations despite both being essentially always present in spatial problems. We introduce a method that constructs valid frequentist confi- dence intervals for associations in spatial settings. Our method requires minimal assumptions beyond a form of spatial smoothness and homoskedastic Gaussian error assumption. In particular, we do not require model correctness or covariate overlap between training and target locations. Our approach is the first to guarantee nominal coverage in this setting and outperforms existing techniques in both real and simulated experiments. Our confidence intervals are valid in finite samples when the noise of the Gaussian error is known, and we provide an asymptotically consistent estimation procedure for this noise variance when it is unknown.
Introduction
Scientists and social scientists often seek to understand direction and magnitude of associations in settings where variables vary spatially. And since these associations are often used to inform policy, communicating uncertainties crucial. Example associations interest include the relationship between acrosol concentrations and regional precipitation changes [76], the link between proximity to major highways and the prevalence of dementia [35], and the link between air pollution exposure and birth weight [34]. In each case, the covariates (e.g., aerosol concentrations, proximity highways, and air pollution) may be viewed as functions of spatial location. Moreover, scientists often have data at some spatial locations but want infer associations at others. For instance, country might measure birth weight and air pollution at the municipal level for some municipalities and wish to understand their relationship in municipalities without data.
Our goal in the present work is to provide valid and useful confidence intervals for an estimator of these associations when (i) the underlying model may be misspecified and (ii) inference is needed spatial locations that may differ from those in the observed data. First we argue that existing methods do not already solve this problem. In particular, we make this argument in turn for modern flexible machine learning methods, linear methods, spatial regression methods, and debiasing approaches.
* 'Equal Contribution
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Flexible Machine Learning Methods. A priori, we might expect a nonlinear relationship between air pollution (as a covariate) and birth weight (as the response). A natural idea is to fit a flexible model such as (deep) Gaussian process [52, 15], transformer [70, 33], or XGBoost [10]. These methods can achieve high predictive accuracy, but the methods on their own often lack interpretability and do not immediately yield conclusions about covariate-respons. associations [57, 17]. Researchers have proposed a number of post hoc interpretability methods, such as Shapley values [61, 38). LIME [53], partial dependence plots [20] and Accumulated Local Effects plots [3]. These methods work directly with the fitted model and can be used to describe associations between individual covariates and model predictions. And because they are interpreting the fitted model (rather than the response itself), there is no need (or mechanism) to quantify uncertainties due to insufficient data. When the fitted model closely approximates the true response, one can view the output of these interpretability methods as describing the relationship between the covariates and response.
In many applications (e.g., in engineering, advertising, and marketing), the available data often contains enough information to be confident that the highly flexible machine learning model ap- proximates the response everywhere of interest. Conversely, in many applications in the sciences and social sciences (such as those cited above), there is often insufficient data to confidently fit a nonparametric or high-dimensional model closely to the true (latent) response everywhere relevant. In particular, in many spatial problems, there is often not enough information for a flexible method to reconstruct the response well in many spatial locations of interest. Nonetheless, might think there could still be sufficient data to capture associations - and, by quantifying uncertainty, we can check our confidence in any conclusions. The discrepancy between applications that are data-rich and data-poor in this sense can help explain the phenomenon observed by Rudin [57]: that post hoc approximations of black-box models may be worse than fitting interpretable models directly. For example, as illustrated in a real-life recidivism case, a linear model used to approximate a black-box model can seriously misrepresent key relationships including that between race and recidivism.
Ordinary Least Squares. So it seems natural to choose an appropriate interpretable model (rather than post hoc method) and quantify its uncertainty. Buja et al. [7] demonstrated that (directly-fit not post hoc) linear models can be used to interpretably summarize associations even in the face of potentially nonlinear relationships (i.e., under misspecification). This observation helps explain why all of the applied studies cited above use linear regression. ? Moreover, ordinary least squares (OLS) comes equipped with classical confidence intervals. However, when the linear model is misspecified, Buja et al. [7] notes that the OLS estimator depends on the covariate values, and it follows that classical confidence intervals valid only at the observed (training) values. Since we're we're interested in valid intervals at other locations, which will have different covariate values, we need another approach. The sandwich estimator [30, 77, 78. 7] offers valid intervals under misspecification, but only if all covariates are drawn from single distribution, which isn't It the case in our setting where we want to draw inferences at different locations from where we observe data. Thus, even though OLS provides interpretable estimates, it does not solve our key challenge: constructing valid confidence intervals in the common scenario where both (i) the linear model is misspecified and (ii) inference is needed at spatial locations that may differ from those in the observed data.
Spatial Regression Methods. One might hope that method specifically tailored for spatial regression could address these issues. Generalized least squares (GLS) regression [2] is designed to handle spatial correlation in residuals [13, pp. 22-24], but does not address the bias introduced by misspecification and nonrandom locations. Bayesian spatial models, such as those based on Gaussian processes, are also common, but their credible intervals tend underestimate uncertainty in the presence of both misspecification and nonrandom locations [74, 42).
Debiasing Approaches. Another natural idea is to construct debiased estimator of the association, and then account for the variance of this debiasing procedure. Importance weighting methods from the covariate shift literature [62] pursue this goal by reweighting each observation according to an estimated density ratio between target and source covariate distributions, aiming to remove bias from distribution shift. Semiparametrie inference in partially linear models (e.g., [55, 54, 11]) takes related approach: it first fits a flexible model to capture variation explained by a "nuisance" parameter (in our case, the spatial location) and then constructs debiased estimator of the target association by regressing out the variation explained by this nuisance, so that the remaining signal reflects the
In fact, Castro Torres and Akbaritabar [9, Figure lJ surveyed quantitative methods in papers up to 2022 and found that over half of those in fields such as agricultural science, social sciences, and health sciences reported results from linear model.

---

direct relationship between the covariate and the response. However, these methods face two key limitations. First, they assume that the debiasing step can perfectly remove bias, which need not be true (especially in finite samples). Second, their validity relies on the assumption that observations are drawn independently and identically. And this assumption rarely holds spatial applications, where locations are fixed and often spread unevenly due to physical, logistical, or policy constraints. We discuss these and related approaches in more detail in Appendix A.
Our Contribution. In what follows, we show through real and simulated experiments (Section 4) that existing approaches can yield confidence intervals with coverage far below the nominal level. Our principal contribution is to introduce the first method for constructing confidence intervals in spatial associations that guarantees frequentist coverage the nominal level even when the underlying model is misspecified and inference is required at fixed, nonrandom locations that differ from those observed. We are able to account for misspecification and nonrandom locations simultaneously by making more spatially appropriate assumptions than prior work. In particular, prior work relies on a (spatially inappropriate) assumption of independent and identically distribution data; we instead assume the response is smooth function of space observed with homoskedastic Gaussian noise. When the variance of this noise is known, our confidence intervals are valid in finite samples. To address the common case where the variance is unknown, we provide an asymptotically consistent estimator for it. In our experiments, our method is the only one that consistently attains nominal coverage (or even comes close).
2 Problem Setup
We start by describing the available data, After reviewing well-specified linear regression, we set up the misspecified case with different target and source data, and we establish our estimand in this case.
Data. Following the covariate-shift literature [4, 47, 14], we refer to our fully observed (training) data as the source data; likewise, we let targer data denote the (test) locations and covariates where we do not observe the response but would like to understand the association between covariates and response. In partieular, the source data consists of N triplets (S,,, Xn, Y, Y,) with spatial location S, E S, covariate X,, E IR and response" Y, E R. Here S represents geographic space; we assume S a metric space with metric dg. We collect the source covariates in the matrix X E R XxP and the source responses in the N-long column vector Y. The target data consists of M pairs (S%, X )MM with S S. X E R.. The corresponding responses Y R. 1 S < M are unobserved. We collect the target covariates in X E MP and responses in column vector y E T...
Review: Well-specified I Linear Model. Though will focus the misspecified case, start by reviewing the classic well-specified case for comparison purposes, In the classic well-specified setup, we have Y, Xn + En and Y = 0 + en with column-vector parameter dolS E R² and 6n; Gm A/(0, o²) for some (unknown) 0² > 0. For any fixed set of source data points (and assuming invertibility holds as needed), we can recover the parameter exactly as
N Σ (Y, -Vx-- 1=
folS arg min
E PERP
(X"X) (1)
N
An analogous formula holds for target points; fOlS is constant across covariate values in any case. Since the population expectation is unknown, analysts typically estimate dOlS via maximum likeli- hood. The standard estimator (fOLs,p) and confidence interval at level a (foLs,p) p) for the pth coefficient (p) are
e (x -TTY, Tors.p = OoLS.p + zaP,
(2)
where o2=
Ze. the a-quantile of a standard normal distribution, and Ep denotes the P-dimensional vector with a in entry and Os elsewhere. Under correct specification, the confidence interval is valid: that is, it provides nominal coverage. E.g., a 95% confidence interval contains the true parameter at least 95% of the time under resampling. The noise variance o² in used in Eq. (2) is typically replaced with an estimate à².
'We propose possible extensions to multivariate responses in Appendix C. But we leave the multivariate- response case largely as an area for future work.
3

---

Misspecified Spatial Setup: Data-Gener Process. In what follows, we assume the data is generated as Y, g(Xn, S,) + cn and Y =((" S)) + e for a function g that need not have a parametric form and with en fm iid _N(00, o²) for some (unknown)) o² >
We assume that source and target covariates are fixed functions of spatial location. Recall that all of the covariates in our examples (acrosol concentrations, proximity to highways, and air pollution) can be expected to vary spatially and be measured with minimal error. We similarly expect meteoro- logical variables such precipitation, humidity, and temperature to be reasonably captured by this assumption.
Assumption 1. There exists a fanction x : S - R.P. such that X ((5)))orr S S M and X, = x(S., for S n S N.
Under Assumption 1, our data- aata-generating process simplifies. Assumption 2. There exists a function f : S + R such that Vm E (1,...,,M N.,,,,,,,,, Ym = f(S))) + e. and Vn E (I....,N): Y, - S(S.) + En. where En en ~ NN(0 o²).
sid
We emphasize that we are still allowing the response to be a function of both the covariates and the spatial location. But we need not state the dependence on the covariates explicitly in Assumption 2 due to Assumption 1. Formally, let g be the function satisfying E|y X. S = [XX S). By Assumption 1, X XS).. and so E S] = (x((S). S). Define f - R by f(s) = ( ((S)).)) S). In other words, the response can be function of the covariates, but because the covariates are fixed function of space, we can also write the response as a function of space alone.
Our Estimand. At a high level, our goal is to capture the relationship between a covariate and the response variable target locations using data from source locations, while taking into account that these two sets of locations may differ. From this perspective, we can define our estimand as the parameter of the best linear approximation to the response, where "best" is defined by minimizing squared error over the target locations.
M
Ors = arg min BCRP
E
Σ (Y
As in [78, 7], since the data-generating process may be non-linear. OolS is no longer constant like the well-specified case; instead it is function of the target locations.
Before solving the minimization in Eq. (3), we argue that covariate shift, as commonly defined, does not solve the problem of interest in this paper. Instead of using our estimand above, one might instead treat the source and target locations as random draws from separate distributions, drop the conditioning S (in Eq. (3)), and rely on covariate-shift methods. Yet to the best of our knowledge these methods currently offer no finite-sample confidence intervals, leaving the question of valid confidence intervals unresolved. Moreover, in the spatial applications we study, locations are rarely ii.d.: for example, the United States Environmental Protection Agency (EPA) places monitoring stations strategically, and the targets of interest may be only few municipalities one - with missing data. In such settings, it is unclear that a meaningful population-level distribution of locations exists. We therefore condition the training and target locations in Eq. (3), i.e. we treat the locations as deterministic. Alternatively, we can view our setting as a special case of covariate shift in which the training and target distributions are the respective fixed sets of locations, and therefore have disjoint support. However, under this view, standard covariate-shift estimators are inapplicable.
or even single
To solve the minimization in Eq. (3), will be convenient to assume invertibility, as for OLS. Assumption 3. X X* is invertible
With Assumptions 1 and 3, we can solve the minimization in Eq. (3) to find
(4)
See Appendix F.1 for derivation. To the best of our knowledge, the target-conditional estimand in Eq. (4) has not been previously proposed or analyzed in spatial linear regression.
'Conversely demographic covariates may more reasonably be thought of as noisy functions of space, and further work is needed to handle the noisy case.
4

---

In order to estimate E|y* [*] sufficiently well to in turn estimate OLS well, we need to make regularity assumptions. In classical and covariate-shift settings, these take the form of i.id. and hounded-density
assumptions. Since we' 'v seen that the classic assumptions are inappropriate in this spatial setting, we instead assume f is shared across source and target (Assumption 2) and not varying so quickly in space as to be difficult to learn from limited data.
Assumption 4. The conditional expectation, f. is L-Lipschitz as a function from (S, (S,ds) (R. I.). That is, for any 8, E S, [[(() - [[+')] & (s,))),')..
As an illustration of how this assumption might be satisfied, consider X an Li-Lipschitz function of the spatial domain, and f(s) 3 x(s) + h(s) with h fixed, Lz-Lipschitz function of S. Then is a (113)|22) + L2 L2>-Liischitz function of the spatial domain. A similar assumption was recently considered in Burt et al. [8] to derive consistent estimator for prediction error in spatial settings.
Recall (from the discussion of y vs. I after Assumption 2) that we ultimately allow our response to be a function of both the covariates and the spatial location. But, using our assumption that covers both dependencies (Assumption 2), we impose our smoothness restrictions on the average response as function of space alone. We think it is often easier to reason about smoothness in physical space, rather than in a (potentially high-dimensional) space of covariates.
3 Lipschitz-Driven Inference
We next provide a confidence interval for OLLS, the pth coefficient of the target-conditional least squares estimand (Eq. (4)). We support its validity with theory (in the present section) and experiments (in Section 4). To that end, we start by providing an fffipuatiy-coomputaalle point estimate. We end by discussing the role and choice of the Lipschitz constant L.
Lipschitz-Dree Driven Point Estimation. Since the target covariates are known, the key challenge in estimating Eq. (4) is estimating the unknown quantity [SY*
For our first approximation, recall that, by Assumption 4, f varies smoothly in space. Since the conditional distribution of the responses given the spatial locations is the same function for both target and source data (Assumption 2), we can approximate R [S*] by a weighted average of TE[Y [|] values for locations 5 near S*. Concretely, let v E RMXN be (non-negative) matrix of weights. If V assigns weight mostly to source locations near each corresponding target location, then by the Lipschitz assumption (Assumption 4), [S** s SE E[ [S]] is also unobserved, so next approximate it by observed values of Y (at each source location in S): WE(Y [S] x V.Y Together, these two approximations yield the estimator: 6% = =T (x* x VY. In our experiments, we construct v as follows.
Definition 5 (Nearest-Neighhor Weight Matrix). Define the 1-nearest neighbor weight matrix by
V mn {
Sn = closest source location S We break ties uniformly at random.
(5)
otherwise
While this simple construction works well in our present experiments, we discuss the potential benefits of other constructions in Appendix B. This point estimation approach closely related to KNN imputation [66], used for missing data. But KNN imputation does not account for repeated use of training responses or bias in estimation due to imputation, problems we address in the next section. We provide more complete discussion of KNN imputation in Appendix A.
Lipschitz-Driven Confidence Intervals. We detail how to efficiently compute our proposed con- fidence interval for Sts., in Algorithm 1 (Appendix B). We prove its validity in Theorem 7 be- low. Before stating our result, we establish relevant notation and intuition for how our method works. First, we show the difference between our estimand and the estimator is normally dis- tributed. Toward that goal, we start by writing 0% POLS,P for ep (X*Tx*) -1 XT E M and v w RN By Assumption 2 and the previous equation,
6% M
Wmf En= Y
OLLS,r
0 Σ Wm (S
Σ
M t ((S.)
n cn
V
(6)
bias
randommess
5

---

That is, Eq. (6) expresses OLS 6% as the sum of (i) a bias term due to differing locations between the source and target data and (ii) a meean-zero Gaussian randonness term due to observation noise.
Since the spatial locations are fixed, the bias term is not random and can be written as E IR. follows that 00 0% ~ (bb, 1211"111) since the variance of the randomness term is the sum of the variances of its (independent) summands.
Our strategy from here will be to (1) bound b, (2) establish a valid confidence interval using our bound on b while assuming fixed o², and (3) estimate a² consistently (as N oc).
To bound the bias b, we use Assumption 4 to write
|
sup gEFf
Σ
Σ Un (Sm)
where FL is the space of L--Lipshitz functions from S + R. In Appendix F.3, we show that is possible to use -antorovich-Ruuinsteinn duality to restate the right side of Eq. (7) Wasserstein- distance between discrete measures. This alternative formulation is useful since it can be cast as a linear program [49, Chapter 3]; see Appendix F.3. Let B denote the output of this linear program.
Given B |6]. the following lemma (with proof Appendix F.4) allows us to construct a confidence interval for LLLS, centered on 6% We discuss the benefits of this construction over alternative approaches in Appendix B.5.
Lemma 6. Let b E (-B B), ē > O, and a E (0.1). Then the narrowest 1 confidence interval that is symmetric and valid for all (b, 2²) of th (A) +b (-2B /2 A) = 11 with & the cumulative density function of a standard normal distribution. Also, the satisfying this inequality is E [ [ - (- -]}.
form --B- B + زدة where A is the solution
The resulting confidence interval appears in Algorithm 1. We next establish its validity. So far, we have covered only the known 2² case. We handle the unknown r² case after the following theorem.
Theorem 7. Suppose (Sm X mm)m- and (Sn, Xns Yn) satisfy Assumptions I to 4 with known o². Define the (random) interval I as in Algorithm I using the known value of a². Then with probability at least a, fors.p 1 That is, I has coverage (conditional on the test locations) at least -aa a.
In Appendix F.2, we prove validity of our confidence interval for generic choice of weight matrix V. Theorem an immediate corollary of that result.
Consistent Estimation of the Noise Variance o². Generally, the noise variance a² is unknown, so we will substitute an estimate for o² in the calculation of the confidence interval I (Step 5 in Algorithm 1). In Corollary 9 below, we show that the resulting confidence interval has asymptotically valid coverage. To that end, we first show that the estimator in Eq. (10) consistent for o².
Proposition 8. Suppose the spatial domain S A- A] for some A > O, E N. Take Assumptions 2 and 4. For any sequence of source spatial locations (S,) take N'y as in Eq. (10).
Then x + o² in probability N 0O.
See Appendix F.5 for proof. For intuition, recall that the conditional expectation minimizes expected squared error over all functions. Since the conditional expectation is L-Lipschitz (Assumption 4),
(Y - [(S)))
inf gETL
(8)
R=1
The empirical version of Eq. (8) is Eq. (10), which we show a good estimate for x² for large N.
When S is subset of Euclidean space or a subset of the sphere, the minimization in Eq. (10) yields quadratic program for computing o-.. We provide implementation details in Appendix B.7.
Given Proposition 8, it follows from Slutsky's Lemma that the resulting confidence interval is asymptotically valid. We provide a formal statement below, and a proof in Appendix F.5
Corollary 9. For S as in Proposition 8, with the assumptions and notation of Theorem 7, but with o2 unknown, the confidence interval T from Algorithm has asymptotic (with fixed M and as N - o0) coverage at level 1 a.

---

Coverage
Widsss
1.0
Coverage
Methods
0.5
Ouzs
GP BCIs
0.0
Nominal
KDEIW
1.0
Sandwich
Covarage
CI Width
O1.5
0.5
GLS
0.0
0.0
Shilt
Shift
Figure 1: Coverages (left) and confidence interval widths (right) for our method as well as 5 other methods (3 methods in the lower experiment). In the upper experiment, method and GP BCIs consistently achieve the nominal coverage (95%); the GP BCIs line (dashed blue) overlaps with ours (solid black) for most shifts. Of the two methods with correct coverage, method yields much narrower intervals. the lower experiment, only our method achieves the nominal coverage. The shaded region for coverage is a (conservative) 95% confidence interval while the shaded region for CI width is +2 standard deviations; for more detail, see Appendix D.1.
Choice of the Lipschitz Constant L. The Lipschitz assumption allows us make inferences about the target data from the source data. The Lipschitz constant subsequently enters intervals in two principal ways: via B and via o- Intuitively, larger values (allowing for less smooth responses) lead to our algorithm constructing confidence intervals with larger bounds (B) the bias but smaller estimated residual variance (o*)) ). We depict this trade-off concrete example in Section 4 (Fig. 2).
Ultimately the choice of Lipschitz constant must be guided by domain knowledge. We give one conerete example describing our choice of the Lipschitz constant our real-data experiment on tree cover (Section 4). We give a second concrete example of how to select the Lipschitz constant in Appendix B.3; this case, the response annual average PM2.5 over California. experiments, we know the minimum value for which Assumption 4 holds; call it Lg. So we first choose Lo- Then we perform ablation studies in both simulated and real data showing that we essentially maintain coverage while varying L over an order of magnitude around our initial choices. show that further decreasing can decrease coverage and discuss why it is useful to err on side of conservatism (i.e., a larger L). However, we expect even small values of L to improve upon classical confidence intervals in terms of coverage, since classical confidence intervals do not account for bias at all: our method similarly ignores bias when L 0.
our simulated
Experiments
In simulated and real data experiments, find that method consistently achieves nominal coverage. whereas all the alternatives dramatically fail to do so. We also provide ablation studies to evaluate the effect of varying the Lipschitz constant in both simulated and real settings.
Baselines. We compare to five alternative constructions.
Ordinary Least Squares (OLS). We treat the noise variance as unknown and estimate it as the average squared residual, with a correction for the number of degrees of freedom and t-statistic instead of z-statistic [24, pp. 50-52].
Sandwich Estimator The sandwich estimator [30, 77, 78] uses the same point estimate OLS but a different variance estimate. We take the variance estimate from MacKinnon and White [39, Equation 6): - +_ (Xx (X ep. where R is diagonal matrix with the squared residuals as entries.

---

CI Width Shift
Shif=000
Shift=0.8
Lipechitz Constants
10
Width
Widtr
7.5
Wiith
Raudommess
Raudommess
5
Bias
Blas
0.1
3.5
0.5
1
0.55
0.0
0.5
0
10
0
2
Shift
Tipschitz Coustant
Tápschitz Coustant
Figure 2: Left: the confidence interval width of our method as a function of shift for each Lipschitz constant L. All L yield coverage of 1.0. Middle and right: the confidence interval width (solid line with dot marker) as function of the Lipschitz constant shift (middle) and shift 0.8 (right). The vertical axis is shared across all three plots. The bias contribution to the width (dashed line, marker) is monotonically increasing in L. The randomness contribution (dashed line, square marker) is monotonically decreasing.
Importance Weighted Least Squares (KDE/W). As suggested by Shimodaira [62, Section 9], calculate importance weights via kernel density estimation (KDE). We select the bandwidth parameter with 5-fold cross validation; see Appendix B.9. Given the KDE weights, use the point estimate and confidence interval from weighted least squares.
Generalized Least Squares (GLS). We maximize the likelihood of Y ~ ~N(0T x, E), with E specified by an isotropic Matérn 3/2 covariance function and nugget, to select the parameters of the covariance function and nugget variance. Then we use the restricted spatial-regresson framework [29]; since project the spatially-correlace error term onto the orthogonal complement of the covariates, the point estimate coincides with OLS.
Gaussian Process Bayesian Credible Intervals (GP BCIs). We use the model (SS) 8 XS) h(S) cc, with T ~ (O, ,, h ~ GP(D, by an isotropic Matérn 3/2 kernel function with hyperparameters ", and 6 ~ (0, )2). We select (A,7,5) 5) by maximum likelihood. We report posterior credible intervals for 0p-
Single Covariate Simulation. In our first simulation, the source locations are uniform on [-1, 1J2 (blue points in Fig. 4 left plots). The target locations are uniform on | -1+shift
1+shift
T+|shift] [sifti
where shift controls the degree of distribution shift between source and target (orange points in Fig. left plots). In this experiment, the single covariate X x(S) = 5(1) S(2) (Fig. 4, third plot). And the response is Y X + 1/2 (() + (S()))3) + c, with (V.0,..13). 0.13). Fig. 4, fourth plot, shows the conditional expectation of the response given location. We can compute the ground truth parameter in closed form because have access to the conditional expectation of the response (Eq. (4)). We vary shift. E [0, +0022 +0.44 +-066, +0.88 and run 250 seeds for each shift.
Figure I, top left, shows that only our method and the GP consistently achieve nominal coverage. Given correct coverage, narrower (i.c., more precise) confidence intervals are desirable; Fig. 1, top right, shows that our method yields narrower intervals than the GP. KDEIW comes close to achieving nominal coverage when there is no shift. But under shift with limited data, it is not able to fully debias the estimate, and coverage drops. For large M (here M 100), expect the sandwich estimator to achieve nominal coverage at shift = since it is guaranteed to cover the population under first-order misspecification without distribution shift [30, 78). But under any of the depicted non-zero shifts, the sandwich, OL.,, and GL.S achieve zero coverage. The extreme narrowness of the OLS, sandwich, KDEIW, and GLS intervals (Fig. 1, top right) suggests that the problem with these methods is exactly their overconfidence. Because these approaches assume that the estimator unbiased (or can be debiased), and that errors are independent and Gaussian, their intervals contract far too quickly, even with small amounts of data (here, N 300). Essentially, these methods' reliance on strong modeling assumptions leads to non-robust coverage.
We note that the confidence intervals and coverages for each method have (approximately) the same values for either =shift in this experiment (Fig. 1, top right) because the covariate is symmetric around the line ((1) s(2); see the middle right plot of Fig. 4 in Appendix D. So positive and
8

---

Aridity Index
Elevation
Slope
Averoge Whith Conerago
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
0g
GLSS
GP BCIs
KDEI
OLS
Sandwich
Oure
GLS
GP BCls
KDEIW
OLS
Saudekhh
Oure
GLS
GP BCIs
KDEW
OLS
Sandwish
Ours
Figure 3: Coverages (upper) and confidence interval widths (lower) for our method as well other methods. Each column represents a parameter in the tree cover experiment. Only our method consistently achieves the nominal coverage.
negative shift values move the source and target locations symmetrically across the 5(1) 5 line. We do not expect such symmetry for general covariates and will not see it in our next simulation (Fig. 1, bottom right).
Simulation with Several Covariates. Our second simulation generates locations as in the previous experiment. Now we use N = 10,000 and = 100. And we generate 3 covariates, X
(Sn() + (os())))) ((2 cos(5))) - (in() and x S(1) + S2. The response is + e with ~ (O, 0.1²).
Y= x + coefficient. Calculation 0% scales poorly with N due to needing to solve a quadratic program. So here instead estimate o² using the squared eиrоr of leave-one-out 1-nearest neighbor regression fit on the source data. See Appendix B.8 for details. We compare against the same set of methods except do not include GLS GP BCIs since these require further approximations to scale for this N.
focus on inference for the first
We again find that our method achieves coverage while the other methods do not (Fig. 1, bottom left). In this experiment, no other method achieves coverage over 30% across any shift value (even 0). As before, competing methods are overconfident, with very small CI widths (Fig. 1, bottom right).
For methods besides our own, coverage levels at O are generally lower in this experiment at shift = than in the previous experiment. The difference is that N much larger here (making confidence intervals narrower and exacerbating overconfidence) while M is the same. The sandwich estimator covers the analogue of (ors.p where spatial locations are treated random. So it has good coverage shift 0 when M >> N, but when M << N.
Effect Lipschitz Constant on Confidence Intervals in Simulation Experiment. Above, we know and use the minimum value (T.) for which Assumption holds; that is, I, To (L.p 2 and 3 3 respectively). Now we repeat the first simulation but vary L E {0.1, 0.5, 1.0, 2.0, 3.5, 7.5, 10). All of these L values yield coverage of 1.0 for our method, above the nominal coverage of 0.95, even though coverage
guaranteed by our theory for L S 2 2.8.
We next show that the confidence interval width reflects bias-randomness trade-off as varies. If the noise were known, the confidence interval would be monotonically increasing in L. Since the noise is unknown, only the bias contribution to the interval width (2B, Step 4 in Algorithm 1) increases (Fig. 2, middle and left, x). Conversely, smaller values of L. yield larger values for d3y- so the randomness contribution to the interval width (2c, Step 5 in Algorithm I) increases (Fig. 2, middle and left, D). The full confidence interval width, 2B 2cA (ar),
monotonic in L.
=p 9p.Xp to quantify how tree cover percentage the contiguous United States (CONUS) in the year 2021 relates to three variables, P (Aridity Index. Elevation, Slope}. We use the 983 data points from Lu et al.
Tree Cover Linear Regression. We use a linear regression model Yiee Cover =

---

[37], who in turn draw on [68, 65, 43]. We define our target region in the West portion of CONUS as locations with latitude in the range (25, 50) and longitude in the range (-125, -10). Out of all points in this region, we designate 50% totaling 133 sites data by taking a uniform random sample of 20% of the remaining spatial locations, repeated over 250 random seeds to assess coverage performance. Each seed yields 170 source locations. Fig. 7 illustrates the spatial split between source and target data for a representative seed. We discuss the data and our pre-processing in detail in Appendix E.
target data. Next, we select the source
We compare coverage for confidence intervals of the three parameters, 9 Ariiiity_iid Indexx Eleeaaiion and Slope We discuss how we evaluate coverage in Appendix E.3. In the top row Fig. 3, we see that our method is the only one to achieve the 95% nominal coverage for all three parameters. Conversely, for the Slope parameter, every other method achieves coverage at most 54%. In the bottom row of Fig. 10, we again see that alternative methods fail to provide coverage due to their overconfidence (small widths); see also Fig. 10, which shows all methods' constructed confidence intervals across all three parameters for of the 250 seeds. In Appendix B.4, we further discuss how our intervals are the narrowest intervals among those that maintain validityy. In Appendix E.5, we conduct a similar analysis but with target locations in the Southeast, rather than West, of CONUS. The results align with our discussion here.
Choice of Lipschitz Constant in the Tree Cover Experiment. For the tree cover experiment, we leverage domain knowledge to set the Lipschitz constant to L = 0.2, in units of percent tree cover per kilometer (km). This choice implies that 1% change in tree cover corresponds to moving 1/0.2 5 km. To arrive at this choice, we observe that in certain regions of the U.S., such as the Midwest, tree coverage remains relatively uniform over several kilometers, so smaller Lipschitz constants (e.g., I. = 0.02, corresponding to 1% change over 50 km) would be appropriate. However, in other regions- - such as the western U.S., where elevation changes are more pronounced (e.g.. the Rockies, California, and the Pacific Northwest) distances. To account for these variations conservatively, we choose L = 0.2. More generally, for real-world applications, we recommend the following strategy: (i) use domain knowledge to select a reasonable Lipschitz constant for the response variable, and (ii) inflate the Lipschitz constant to ensure conservative estimate (which is more likely to satisfy Assumption 4).
tree cover can change sharply over short
Effect of Lipschitz Constant on Confidence Intervals in Tree Cover Experiment. In Fig. 9, we show the coverage and width of our confidence intervals across three orders of magnitude of Lipschitz constants (L from 0.001 to 1). For varying between 0.1 and (a single order of magnitude variation around our chosen value of 0.2), we find that coverage is always met except in one case (Aridity Index with L = 0.1), where it is very close to nominal (89% instead of 95%). For Slope and Elevation, coverage is met or very nearly met for all L values. For Aridity Index, coverage is low for L < 0.1. Meanwhile, confidence intervals become noticeably wider for L > 0.5 while remaining relatively stable for smaller values. These results support intuition to err on the side of larger L values to be conservative and maintain coverage.
5 Discussion
We propose a new method for constructing confidence intervals for spatial linear models. We show via theory and experiments that our intervals accurately quantify uncertainty under misspecification and covariate shift. In experiments, our method is the only method that consistently achieves (or even comes close to) nominal coverage. We observe that, very commonly in spatial data analyses, covariates and responses may be observed at different, nonrandom locations in space. Since method does not actually use the source covariate values in inference for G it can be applied in this common missing-dal scenario. Though it requires additional work, we believe the ideas here will extend naturally to the widely used class of generalized linear models.
Acknowledgments
This work was supported in part by Social and Ethical Responsibilities of Computing (SERC) seed grant, an Office of Naval Research Early Career Grant, Generali, a Microsoft Trustworthy Al Grant, and NSF grant 22 2417177.
10

---

References
[1] Agrawal, A., Verschueren, R., Diamond, S., and Boyd, S. (2018). A rewriting system for convex optimization problems. Journal of Control and Decision, 5(1):42-60.
[2] Aitken, A. C. (1936). IV least squares and linear combination of observations. Proceedings of the Royal Society of Edinburgh, 55:42-48.
[3] Apley, D. W. and Zhu, J. (2020). Visualizing the effects of predictor variables in black box super- vised learning models. Journal of the Royal Statisfical Society Series B: Statistical Methodology, 82(4))1059-100)
[4] Ben-David, S., Blitzer, J., Crammer, K., and Pereira, F. (2006). Analysis of representations for domain adaptation. In Advances in Neural Information Processing Systems (NeurIPS).
[5] Brent, R. (1971). An algorithm with guaranteed convergence for finding a zero of a function. The Computer Journal, 14(4):422-45 425.
[6] Brunsdon, C., Fotheringham, A. S., and Charlton, M. E. (1996). Geographically weighted regression: a method for exploring spatial nonstationarity. Geographical analysis, 28(4):281-29)
[7] Buja, A., Brown, L.,, Berk, R., George, E., Pitkin, E., Traskin, M., Zhang, K., and Zhao, L. (2019). Models as approximations I: Consequences illustrated with linear regression. Statistical Science, 34(4):523-544
[8] Burt, D. R., Shen, Y., and Broderick, T. (2025). Consistent validation for predictive methods in spatial settings. In Artificial Intelligence and Statistics (AISTATS).
[9] Castro Torres, A. F. and Akbaritabar, A. (2024). The use of linear models in quantitative research. Quantitative Science Studies, 5(2):426 446
[10] Chen, T. and Guestrin, C. (2016). XGBoost: scalable tree boosting system. In Proceedings of the 22nd ACM SIGKDD international conference knowledge discovery and data mining, pages 785-794.
[11] Cheroozuukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1):C1-C68.
[12] Chow, J. C., Chen, L.-W. A., Watson, J. G., Lowenthal, D. H., Magliano, K. A., Turkiewicz, K., and Lehrman, D. E. (2006). PM2.5 chemical composition and spatiotemporal variability during the California regional PM 1 air quality study (CRPAQS). Journal of Geophysical Reseurch: Atmospheres, 11((1))).
[13] Cressie, N. (2015). Statistics for spatial data. John Wiley & Sons.
[14] Csurka, G. (2017). A Comprehensive Survey on Domain Adaptation for Visual Applications, pages 1-35. Springer International Publishing.
[15] Damianou, A. and Lawrence, N. D. (2013). Deep Gaussian processes. In Artificial Intelligence and Statistics (AISTATS).
[16] Diamond, S. and Boyd, S. (2016). CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1-5.
[17] Velei-Velez, F. and Kim, B. (2017). Towards rigorous science of interpretable machine learning, arXiv preprint
[18] Flamary, R., Courty, N., Gramfort, A., Alaya, M. Z.. Boisbunon, A., Chambon, S., Chapel, L., Corenflos, A., Fatras, K., Fournier, N., Gautheron, L.., Gayraud, N. T., Janati, H., Rakotomamonjy. A., Redko, I., Rolet, A., Schutz, A., Seguy, V., Sutherland, D. J., Tavenard, R., Tong, A., and Vayer, T. (2021). POT: Python Optimal Transport. Journal of Machine Learning Research, 22(78):1
[19] Fotheringham, S., Brunsdon, C., and Charlton, M. (2002). Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Wiley.
11

---

[20] Friedman, J. H. (2001). Greedy function approximation: A gradient boosting machine. The Annals of Statisties,, 2965)1119--122.
[21] Gilbert, B., Ogburn, E. L., and Datta, A. (2024). Consistency of common spatial estimators under spatial confounding. Biometrika.
[22] Goulart, P. J. and Chen, Y. (2024). Clarabel: An interior-point solver for conic programs with quadratic objectives. arXiv preprini arXiv:2405.12762.
[23] Gramacy, R. B. (2020). Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences. Chapman and Hall/CRC.
[24] Greene, W. H. (2012). Econometric Analysis: 7th Edition. USA: Pearson Prentice Hall.
[25] Gretton, A., Smola, A., Huang. J., Schmittfull, M., Borgwardt, K., and Schölkopf. B. (2008). Covariate shift by kernel mean matching. In Dataset Shift in Machine Learning. MIT Press.
[26] Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi. R., Gerber, F., Gramacy, B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F., and Zammit- Mangion, A. (2019). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological and Environmental Statistics, 24:398-425.
[27] Heckman, N. E. (1986). Spline smoothing in partly linear model. Journal of the Royal Statistical Society Series B: Statistical Methodology, 48(2):244 248.
[28] Hinkley, D. V. (1977). Jackknifing in unbalanced situations. Technometrics, 1933):285-292.
[29] Hodges, J. S. and Reich, B. J. (2010). Adding spatially-correlated errors can mess up the fixed effect you love. The American Statistician, 64(4):325-33 3344.
[30] Huber, P. J. (1967). The behavior of maximum likelihood estimates under nonstandard condi- tions. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics, pages 221-234 University of California Press.
[31] Kanamori, T., Suzuki, T., and Sugiyama, M. (2012). Statistical analysis of kernel-based least-squares density-ratio estimation. Machine Learning, 86:335-367.
[32] Kirszbraun, D. (1934). Über die rusammenziehende und lipschitzsche transformationen. Fundamenta Mathematicae, 22:77-108.
[33] Klemmer, K., Rolf, E.,, Robinson, C., Mackey, L., and Rußwurr, M. (2025). Satclip: Global, general-purpose location embeddings with satellite imagery. Proceedings of the AAAI Conference on Artificial Intelligence, 39(4))4347-435.
[34] Lee, J., Costello, S., Balmes, J. R.. and Holm, S. M. (2022). The Association between Ambient PM2.5 and Low Birth Weight in California. International Journal of Environmental Research and Public Health, 192204133554
[35] Li, C., Gao, D., Cai. Y. S., Liang, J., Wang, Y., Pan, Y., Zhang, W., Zheng, F., and Xie, W. (2023). Relationships of Residential Distance to Major Traffic Roads with Dementia Incidence and Brain Structure Measures: Mediation Role of Air Pollution. Health Data Science, 3:0091.
[36] Liu, D. C. and Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical Programming, 45(1-3):503-528.
[37] Lu, K., Kluger, D. M., Bates, S., and Wang. S. (2025). Regression coefficient estimation from remote sensing maps. Remote Sensing of Environment, 330:1 14949.
[38] Lundberg, S. M. and Lee, S.-L (2017). A unified approach to interpreting model predictions. Advances in Neural Information Processing Systems (NeurIPS).
[39] MacKinnon, J. and White, H. (1985). Some heteroskedastiiitycconiccc covariance matrix estimators with improved finite sample properties. Journal of Econometrics, 29(3):305-325.
12

---

[40] Matthews, A. G. d. G., van der Wilk, M., Nickson, T., Fujüi, K., Boukouvalas, A., León- Villagrá, P., Ghahramani, Z., and Hensman, J. (2017). GPflow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40):1-6.
[41] Mayor, J. R., Sanders, N. J., Classen, A. T., Bardgett, R. D., Clément, Fajardo, A., Lavorel, S., Sundqvist, M. K., Bahn, M., Chisholm, C., Cieraad, E., Gedalof, Z., Grigulis, K., Kudo, G., Oberski, D. L., and Wardle, D. A. (2017). Elevation alters ecosystem properties across temperate treelines globally. Nature, 542(7639):91-95.
[42] Müller, U. K. (2013). Risk of Bayesian inference in misspecified models, and the sandwich covariance matrix. Econometrica, (1)5): ((5):1805-184)
[43] NASA Jet Propulsion Lab (2020). NASADEM merged DEM global arc second v001 [data sct]. NASA EOSDIS Land Processes DAAC.
[44] Nobre, W. S., Schmidt, A. M., and Pereira, J. B. M. (2021). On the effects of spatial confounding hierarchical models. International Statistical Review, (9)2)3002-322.
[45] Paciorek, C. J. (2010). The importance of scale for spatial-confounding bias and precision of spatial regression estimators. Statistical Science, 25(1):107-125.
[46] Page, G. L., Liu, Y., He, Z., and Sun, D. (2017). Estimation and prediction in the presence of spatial confounding for spatial linear models. Scandinavian Journal of Statistics, 44(3):780-797.
[47] Pan, S. J. and Yang, Q. (2010). A survey on transfer learning. IEEE Transactions on Knowledge and Data Engineering, 22(10):1345-135)
[48] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. (2011). scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 2:2825-2883..
[49] Peyré, G. and Cuturi, M. (2019). Computational Optimal Transport: With Applications to Data Science. Foundations and Trends@ in Machine Learning, 11(5-6)335-6077
[50] Portier, F., Truquet, L., and Yamane, I. (2024). Nearest neighbor sampling for covariate shift adaptation. Journal of Machine Learning Research, 25(410):
[51] Que, Q and Belkin, M. (2013). Inverse density an inverse problem: the Fredholm equation approach. In Advances in Neural Information Processing Systems (NeurIPS).
[52] Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. The MIT Press.
[53] Ribciro, M. T., Singh, S., and Guestrin, C. (2016). "Why should I trust you?" Explaining the predictions of any classifier. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1135-1144.
[54] Robins, J. M., Mark, D., and Newey, W. (1992). Estimating exposure effects by modelling the expectation of exposure conditional on confounders. Biometrics, pages 479-495,
[55] Robinson, P. M. (1988), Root-N-consistent semiparametric regression. Econometrica: journal of the Econometric Society, pages 931-954.
[56] Rolf, E.,, Klemmer, K., Robinson, C., and Kerner, H. (2024). Position: Mission critical satellite data is a distinct modality in machine learning. In Proceedings of the 41st International Conference on Machine Learning.
[57] Rudin, C. (2019). Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. Nature machine intelligence, 1(5):206-215.
[58] Rudin, (1976). Principles of Mathematical Analysis. International series in pure and applied mathematics. McGraw-Hill.
13

---

[59] Sandel, B. and Svenning, J.-C. (2013). Human impacts drive a global topographic signature in tree cover, Nature Communications, 4(1):2474.
[60] Seabold, S. and Perktold, (2010), statsmodels: Econometric and statistical modeling with python. In 9th Python in Science Conference.
[61] Shapley, L. S. (1953). A value for n-person games. In Contributions the Theory of Games II. Princeton University Press Princeton.
[62] Shimodaira, H. (2000). Improving predictive inference under covariate shift by weighting the log- log-likelihood function. Journal of Statistical Planning and Inference, 90(2):27-2244.
[63] Sugiyama, M., Suzuki, T., and Kanamori, (2012). Density Ratio Estimation Machine Learning. Cambridge University Press.
[64] Sugiyama, M., Suzuki, T., Nakajima, S., Kashima, H., von Bünau, P., and Kawanabe, M. (2008). Direct importance estimation for covariate shift adaptation. Annals of the Institute of Statistical Mathematics, 60(4):699-74)
[65] Trabucco, A. (2019). Global aridity index and potential evapotranspiration (ETO) climate database v2. CGIAR Consortium for Spatial Information.
[66] Troyanskaya, O., Cantor, M., Sherlock, G., Brown, P., Hastie, T., Tibshirani, R., Botstein, D., and Altman, B. (2001). Missing value estimation methods for DNA microarrays. Bioinformatics. 17(6):520-525.
[67] Tuia, D., Persello, C., and Bruzzone, L. (2016). Domain adaptation for the classification of remote sensing data: An overview of recent advances. IEEE Geoscience and Remote Sensing Magazine, 4(2):41-57.
[68] USFS (2023). USFS tree canopy cover v2021.4 (Conterminous United States and Southeastern Alaska).
[69] Valentine, F. A. (1945). A Lipschitz condition preserving extension for a vector function. American Journal of Mathematics, 67(1):83--93
[70] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, u., and Polosukhin, I. (2017). Attention is all you need. In Advances in Neural Information Processing Systems (NeurIPS).
[71] Villani, C. (2009). Optimal Transport: Old and New. Springer.
[72] Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, van der Walt, J., Brett, M., Wilson, J., Millman, K. J.. Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, I., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H.. Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272.
[73] Wainwright, M. J. (2019). High-Dimensional Statistics: A Non-Asymptotic Viewpoint. Cam- bridge Series Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge.
[74] Walker, S. G. (2013). Bayesian inference with misspecified models. Journal of Statistical Planning and Inference, 143(10):1211-1633
[75] Weber, L. M., Saha, A., Datta, A., Hansen, K. D., and Hicks, S. C. (2023). nnSVG for the scalable identification of spatially variable genes using neirest-n Gaussian processes. Nature Communications, 14(1):1-12.
[76] Westervelt, D. M., Conley, A. J., Fiore, A. M., Lamarque, J.-F,, Shindell, D. T., Previdi, M., Mascioli, N. R., Faluvegi, G., Correa, G., and Horowitz, L. W. (2018), Connecting regional aerosol emissions reductions to local and remote precipitation responses. Atmospheric Chemistry and Physics, 12616):124....2455
14

---

[77] White, H. (1980a). A heteenskedasttitty-ossssee covariance matrix estimator and direct test for heteroskedasticity. Econometrica, 48(4):817-838.
[78] White, H. (1980b). Using least squares to approximate unknown regression functions. Interna- tional Economic Review.
[79] Yu, H., Fotheringham, A. S., Li, Z., Oshan, T., and Wolf, L. J. (2020). On the measurement of bias in geographically weighted regression models. Spatial Statistics, 38: 100453.
15

---

A Extended Related Work
In this section, we discuss related work on bias in linear regression in spatial settings, local approaches to regression that often rely on data- cata-borro strategies similar to our nearest neighbor approach. and covariate shift.
Bias in Spatial Regression. Linear models with Gaussian process random effect that is, models of the form
+ ((S)= +
with a Gaussian process and independent and identically distributed noise are classic tool in spatial regression and remain widely used in applications [75. 23, 26]. Whether to treat the covariates as fixed functions of spatial location or random is a topic of significant debate. In order for the model to be identifiable, it generally necessary for the covariates to be thought of as random [21, Proposition 1]. However, Paciorek [45] observed it may be more reasonable to think of covariates as fixed (the perspective we also take). If X and g are not independent (in the random case) or close to orthogonal in the fixed case, Paciorek [45] illustrates that bias is introduced into the estimation of A. And that the degree of this bias depends on the relative scales over which X and g vary. He gives a closed form for this bias under strong linear- Gaussian assumptions that depend on parameters that would generally be unknown. Page al. [46] builds on the result of Paciorek [45], and shows that bias in predictions made by the spatial model due to confounding may be small, despite biases in parameter estimation. Hodges and Reich [29] proposed using restricted spatial regression, essentially explaining as much of the variation in the response as possible by the covariates to reduce bias in estimation due to spatially-correlatte our experiments. Nobre et al. [44] extended earlier work about bias due to spatial confounding by considering the case when multiple observations are available at the same spatial location with independent noise. They showed that spatial regression models can still be biased with repeated observations due to confounding.
error. We take this approach in the GLS baseline used in
In our work, we focus on bias in estimation of the target conditional OLS. Because we do not assume the model is well-specified, there is no 'true' value of the parameter to estimate. Instead, our goal is to estimate the 'best' (in least squares sense) lincar approximation to the response. We focus on the case when X is a fixed function of location (Assumption 1) and make weak non-parametric, smoothness assumptions (Assumption 4). This contrasts with prior works that generally assume the covariates are random and make Anear-Gaussian assumptions to calculate the bias [45, 46). Our work shows how to directly incorporate bias into confidence intervals for & under weak assumptions, but does not address issues of inon-Jidentifaabiiity in the estimation problem.
Gilbert et al. [21] showed that despite finite-sample bias in estimation of 0 under spatial confounding. consistent estimation of 0 is possible in the identifiable case when X random under infill asymp- toties. This is essentially orthogonal to our work, we focus on stmite valid) confidence intervals in the case when X is fixed.
Local Versus Global Regression Estimates. An alternative perspective on regression in spatial set- tings, referred to as geagraphically weighted regression (GWR), focuses on estimating the association between the response and covariates at each location in space [6, 19]. For this regression problem to be well-defined, the covariates at each spatial location must be treated as random variables. Like our method, this approach uses nearest neighbors or local smoothing approaches to estimate the local coefficients at each location in space. And - as in our work - because of this borrowbo from nearby locations, bias is introduced into the estimation of these local regression coefficients. The confidence intervals reported for GWR typically do not account for this bias [79]. While Yu et al. [79] provide a formula for the bias introduced by "data- bittaaborrowing,,. their formula relies on the assumption that the model is well-specified and depends on unobservable parameters.
We focus on the estimation bias in the estimation of global parameters. Although local parameter estimates are often of interest in spatial analyses [19, Chapter 1], global estimates are also useful as summary statistics for making decisions that impact regions or multiple localities. And we characterize the bias in these estimates under weak, non-parametric assumptions and incorporate it into our confidence intervals. Extensions of our approach to GWR coefficients are promising direction for future research. However, this extension not simple, as it involves accounting for randomness in the covariates.
16

---

Covariate Shift and Importance Weighting. Differences between the source locations (with observed responses) and the target locations (at which we want to study the relationship between variables) can be seen as an instance of covariate shift [56, 67, 75J. Covariate shift is often dealt with by importance weighting, reweighting each training example to account for differences in the density of the source and target distributions. Prior work in spatial machine learning has generally focused on addressing covariate shift in the context of estimatinn a method's prediction error. However, these approaches can, in principle, be extended to parameter estimation in OLS, mirroring the broader covariate shift literature [62].
While in the main text we focused on density ratio estimation using kernel density estimation [62], many other approaches for density ratio estimation can be used [63. 31, 25, 51, 64]. Importance weighting de- biiaesses the estimator, but can only be used if the distribution of the source data is supported on a region containing the test data, And the confidence intervals obtained by importance weighting with estimated weights do not account for errors in estimation of the density ratio (if this exists). These confidence intervals are therefore not guaranteed to achieve the nominal coverage rate. And they are not applicable in cases where extrapolation is required, the estimator cannot necessarily be de-biased in these cases.
In contrast, our approach ensures nominal coverage even when de-biasing the estimate via importance weighting is not possible. This advantage comes at the cost of an additional regularity assumption (Assumption 4), as well as often empirically wider confidence intervals.
Comparison to Semiparametrie Inference for Partially Linear Model. Partially linear models take the form
E Y X, S = 8 + ((S) [S] = X(S)
where 3 is the parameter of interest, S a nuisance (or control) variable, X is the covariates, and y is an unknown and possibly complicated function. These models are widely studied in the semiparametric literature. Among many others, Robins et al. [54], Robinson [55]. Chernozhukov et al. [11] focus estimation of 3, under the assumption that the triples (Sn, Xn, Yv,) are independent and identically distributed across data indices n.
In many spatial applications, it isn't 't reasonable to think of the nuisance variable (geographic space) as sampled independently and identically, or even at regularly spaced locations. Observational data are often collected in highly non-uniform way - densely in some regions, sparsely or not at all in others - due to physical constraints, accessibility, or policy decisions. This non-uniformity introduces distribution shifts when attempting to generalize inferences from one region to another. In the present work, we focus on inference with fixed spatial locations and do not impose regularity conditions on the sampling design. This setup allows to quantify uncertainty in associations in cases where extrapolation to poorly-sampled geographic areas is required, or in cases with heavily clustered training locations and more uniform target locations.
notable exception to the assumption of fully i.i.d. data in the semiparametric literature is Heckman [27]. Heckman [27] considered time as nuisance variable, and assumed this nuisance variable was one-dimensional and sampled densely and in sufficiently regular way. In contrast, we allow for multiple spatial dimensions and do not require regularity assumptions about the sampling design.
Comparison to KNN Imputation. Our estimator is closely related to KNN imputation [66]. The 1-Nearest Neighbor point estimate we consider is equivalent to INN imputation for the response, but using only the spatial coordinates - not the covariate values- as features when performing imputation. Directly applying KNN imputation to fill in the target response and then calculating standard confidence intervals would not lead correct estimates of the variance because of repeated use of training data when imputing missing values- we calculate the variances accounting for potential repeated use of these responses. Additionally, using our smoothness assumptions, we are able to quantify additional uncertainty due to potential bias introduced when imputing the missing response values at the target locations. Because our confidence intervals account for potential bias and calculate the variance accounting for the weight assigned to each training example, our confidence intervals are guaranteed to be conservative, whereas confidence intervals calculated after data imputation need not be.
17

---

Algorithm Lipschitz-Driven CI for OL.S.
input ((S, X,, Y,)] ((5m X x Lipschitz constant I,, confidence level a, o² (op-
tional)
output A (1 -cconideeee interval I for OL...
1: v + I-NN mn
L.W
[5m 1m/m=1) (Definition 5) X
2: 6% ofp [ X
Y
3: t op (x+1 x*)
X
w
4: t supgEFL ΣM um g(5m)
En of (S,)) Compute with linear program (Ap-
pendix B.2).
5: of unknown, solve for the following estimator vig, quadratae program (Appendix B.7);
a² =
inf gEFL
E 1Vy=-(S5)) -(,)))
(10)
N
6: c c a 1|
7: Find A(a) satisfying $ (A(())
(--2B/c- A(a)) a, by root finding algorithm
(Appendix B.6).
8: g
0 B - c()) 0 B++
B Implementation Details
In this section, we describe the implementation details of our method. Particularly, this involves computing the upper bound on the bias, computation of A Lemma 6 and computation of 0% from Eq. (10). To summarize our method, we state it an algorithm.
B.1 Selecting the Matrix V.
In the main selected v by using the nearest source location to each target location. We now provide additional discussion about alternative choices of y/ and possible trade-offs.
A generalization of the 1-nearest neighbor approach is to consider V matrix determined by K- nearest neighbors.
Definition 10 (Nearest-Neighbor Weight Matrix), Define the K-nearest neighbor weight matrix by
K
Sn E {K closest locations S..
V mn
(11)
otherwise
For definiteness, assume that, if multiple sources are the same distance from target, ties are broken uniformly at random.
When using a K-nearest neighbor matrix, there is an inherent birs-variance trade-off in choosing K. Increasing broadens the geographic range of source observations used, hence the bias will increase; however, it also lowers the variance by averaging over more responses. The degree to which increasing K introduces additional bias depends on the density data relative to the Lipschitz constant the benefit of smoother weighting schemes (e.g. K > 1) generally becomes more pronounced as the density of spatial sampling increases since smoother weighting schemes naturally leverage multiple nearby points, reducing variance while controlling bias effectively. We do not investigate this trade-off in detail. the experiments we ran, we found 1-nearest neighbor to work well. This is consistent with results in the mean estimation literature suggesting that if source and target distributions overlap substantially, the variance term remains manageable, even when using 1-nearest neighbor approaches [50].
B.2 Implementation of the Wasserstein Bound Calculation
We show in Appendix F.3 how Eq. (7) reduces to computing a Wasserstein-1 distance between empirical measures. To implement the Wasserstein distance calculation, rely on the Python
18

---

Optimal Transport library [18]. For the simulated experiments, we compute the cost matrix using Euclidean distances between spatial locations. For the real-world experiment, we use the Haversine distance to account for Earth's 's curvature.
B.3 Choice of Lipschitz Constant Based on Prior Research: Example for Air Pollution
As concrete example of selecting the Lipschitz constant based on domain-expertise we consider the problem of selecting a Lipschitz constant in an analysis where the response is annual average PM2.5 over California. Chow et al. [12] claim that "Zones of representation for PM2.5 varied from 5-10km for the urban Fresno and Bakersfield sites, and increased 15-20km for the boundary and rural sites" where "[t]he zone of representation is defined as the radius of circular area in which species concentration varies by no more than +00 as extends outward from the center monitoring site." The annual PM2.s concentrations in the study area do not exceed 30pg/m² The combination of zone representation between 5-20km, and variation of not more than 30 ug/mm within this zone of representation suggests a range of possible Lipschitz constants: 0.25-1.2(pg/m) h2./gg/gm////// //km The authors also point out that topographical and meteorological phenomena contribute to this scale of variation. So we would not expect this proposed constant to be "universal" for problems related to PM2.5. but we might expect that this range of Lipschitz constants a reasonable starting point for other studies involving annual average PM2.5 with similar weather and topography California. We showed our real-data analysis that a range of Lipschitz constants can still produce qualitatively similar results (Figure 9) and correct coverage. To err towards the side of conservative analysis, /km)
would recommend a user selects the largest Lipschitz constant in this range (i.e. I. 1.2(pg/m²) /km)
B.4 Confidence-Intervs Widths
In some experiments, most notably the tree-cover analysis in Fig. 3, our confidence intervals can ap- pear wide. This increased width arises naturally from our explicit control of bias under extrapolation: when test locations differ substantially from training locations, the Wasserstein- enlarges, yielding intervals that faithfully represent genuine uncertainty rather than methodological conservatism. By contrast, cxisting methods (c.g., OLS, sandwich, and importance-weighter ap- proaches) produce much narrower intervals but fail to achieve nominal coverage, often approaching zero coverage in Fig. I and Fig. 2. Our method therefore yields the narrowest intervals among those that do maintain validity. Producing sill-narrower intervals would be trivial if one were willing to sacrifice coverage- after all: a zero-width interval achieves minimal size but no inferential mean- ing. Finally. note that intervals typically informative in practice: in many settings they remain sufficiently tight to exclude zero, providing clear conclusions about the sign and magnitude of associations.
bias bound
In practice, interval width also reflects the spatial density of available training data. When observations are closely spaced relative to the smoothness scale implied by the Lipschitz constant, the bias bound remains small and the resulting intervals are narrow. As the spacing grows, the bound naturally increases, producing wider intervals that correctly reflect the added uncertainty from extrapolation. Thus, even when data are sparse, the width of our intervals provides a meaningful indication of the reliability of inferences under the assumed smoothness level.
B.5
Use of Confidence Interval Lemma 6
In all experiments, we construct confidence intervals in Lemma 6. A simpler alternative, which guarantees coverage for all Gaussian distributions (bb e²) with b E [-B, two-sided confidence intervals for each b E [-B, B) and then take their union. This produces an interval of the form
/BBB). B). is to form
[ B - c$ 10-1(1-1))
2) 2) B + cΦ c [1- :]
However, the confidence intervals from Lemma 6 are never longer than this union-based approach, and the root-finding step required to compute them adds negligible overhead. Consequently, we opt for the intervals in Lemma 6 in all our experiments.
19

---

B.6 Computation of A for Lemma 6
Define ( (b)) ( 2B A) + a. Our goal is find a root of g in the interval [+-1(1 (1)a)), a), 2 (1 a/2)]. We first show that g is monotonic, that a root exists in this interval, and that the root is unique.
e
2B A), where & denotes the Gaussian probability density function. This is strictly positive, and so is strictly monotone increasing.
Differentiating,
(4) ((
Also,
a)) 1 (-e
2B
--1(1-
-tta<l-a--t-a- (12)
by non-negativity of the CDF. And,
2B Φ Φ- - aa/2)
g( - -1/22)) xx/2 +b
-14 +a
(13)
= aa2 (1 -1- +
2B
+ (1 a/2
(14)
c
We used symmetry of the Gaussian in the second equality. Then,
2B (?
1- ------11----
21 a/2) (1-2))
a/20/2
(15)
and so 1 (1- a/2)) < O. We conclude that g has a root in the interval. By strict monotonicity
of g, this root is unique.
We Brent's method [5] as implemented in Scipy [72] compute this root numerically.
B.7 Computation of 8 via quadratic programming
To estimate the noise parameter, we need solve the minimization problem
o inf SEJL
Σ OYn -(S,))
(16)
The first obstacle is that the infimum is taken over an infinite-dimensionn space. However, by the Kirszbraun theorem [32, 69], every L-Lipschitz function defined on subset of RD or S can be extended to an L-Lipschitz function on the whole domain. Since our objective function depends only on the values of g at the source spatial locations, only need to enforce the Lipschitz condition between all pairs of source spatial locations, not over the entire domain. Enforcing the Lipschitz condition these N source locations amounts to N(N 1)/2 linear inequality constraints.
Particularly, if define G (g(S1), (S2), ... ... gS.))) then we have the constraints
AG < Lvee(T)
(17)
where L is the Lipschitz constant, I' E N is the matrix of pairwise distances between distinct points in S, and A E NNN a sparse matrix with exactly one and one -1 in distinct rows of each column, representing all such pairs in its rows.
The objective a symmetric, positive-definite quadratic form since it a sum of squares. Therefore, the optimization problem is quadratic program.
In practice, the Scipy sparse matrix algebra [72] and the CLARABEL solver [22] through the CVXPY optimization interface [16, 1] to solve this quadratic program.
B.8 Scalable Estimation of o
The quadratic programming approach for estimating o² outlined in the main text (Eq, (10) and described in detail the previous section does not scale well to large numbers of source locations. Therefore, for synthetic experiment with N 10,000, we take different approach.
20

---

For In N, let 7(n)= = arg min #n ds(Sn, (S(Sn,,)) Sa). That is s7(n) is the index of the nearest point to Sn.. Define the estimator,
N
(Xn E
xn(n))
(18)
2N
n=1
Then as long as (S Si,) 0, by the Lipschitz assumption, Y, Yy(n)) = Cn (n))) and so
N
N
N
NN
~ 2N
Σ
En - En())
(n + Σ
(v(1)
(19)
2N
==1
W=1
n=l
n=1
Then,
E E N
0².
(20)
In general, we expect the estimate concentrate around its expectation, provided that no single point in the source data is the nearest neighbor of too many other points in the source data.
B.9 Implementation of Baseline Methods
We now describe the details of the implementation of the baseline methods.
Ordinary Least Squares. We use the ordinary least squares implementation from statsmodel [60]. We use the default implementation, which calculates the variance as the average sum of squared residuals with offfreddm correction, as in Eq. (2). We use a I-statistic to compute the corresponding confidence interval, which is again the default in statsmodel [60].
Sandwich Estimator. We use the sandwich estimation procedure included in ordinary least squares in statsmodels [60]. We use the HC1 function, which implements the sandwich estimator with the degrees-of- of-fredoom correction from Hinkley [28], MacKinnon and White [39]. That is, the (XX)(x"x)-+ ep. We use the default settings in
variance is estimated as N- Xx_x) (X"XX"X statsmodels, which use a z-statistic with the sandwich estimator to compute the corresponding
confidence interval.
Importance Weighted Least Squares. We use the scikit-learn [48] implementation of kernel density estimation estimate the density of test and train point separately. We use Gaussian kernll (the default) and perform 5-fold cross validation to select the bandwidth parameter, maximizing the average log likelihood of held-out points in each fold. For the simulation experiments, cross-validation to select the bandwidth parameters over the set {0.01, 0.025, 0.05, 0.1, 0.25, 0.5). We selected this set of possible bandwidths span several orders of magnitude from very short bandwidths, to bandwidths on the same order as the entire domain. We select the bandwidths for the source and target density estimation problems separately. For the real data experiments, the bandwidths was selected from the set {0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0,5, 13. This range was selected since the maximum Haversine distance between points in the spatial domain of interest is approximately 1. Once density estimates are obtained, we evaluate the ratio of the density function on the training locations, and use these as weights to perform weighted least squares. Weighted least squares is performed using the default settings in statsmodels [60].
performed
Generalized Least Squares. We perform generalized least squares in two-stage manner. We first approximate the covariance strueture with a Gaussian process regression model. Then we use this approximation fit generalized least squares model with restricted spatial regression [29].
More precisely, first we optimize the maximum likelihood of a Gaussian process regression model with linear mean function depending the covariates including intercept and Matérn 3/2 kernel defined on the spatial source locations. We use the GPFlow [40] implementation of the likelihood, and L-BFGS [36] for numerical optimization of the likelihood. The optimization is initialized using the GPFlow default parameters for the mean and covariance functions.
Once the maximum likelihood parameters have been found, we use the found prior covariance function for defining the covariance between datapoints to be used in the generalized least squares routine. We use restricted spatial regression [29], and so the covariance matrix is defined as PK P + > where 2 the noise variance selected by maximum likelihood in the GP model, K the N N matrix formed by evaluating the Matémn 3/2 kernel with the maximum likelihood kernel parameters on the
21

---

source locations, and P is the orthogonal projection onto the orthogonal complement of the covariates (including intercept). i.c. IN
This covariance matrix PKP + I is passed to the GLS method in statsmodel, and confidence intervals as well as point estimates are computed using the default settings.
Gaussian Process Bayesian Credible Intervals.
We first optimize the maximum likelihood of a Gaussian process regression model with Matérn 3/2 kernel defined on the spatial source locations summed with a linear kernel defined on the covariates. This has the same likelihood as having linear mean function in the covariates with Gaussian prior over the weights, see Rasmussen and Williams [52. Page 28]. We use the GPFlow [40] implementation of the likelihood, and L-BFGS [36] for numerical optimization of the likelihood. The optimization initialized using the GPFlow default parameters for the mean and covariance functions. Once we have caleulated the maximum likelihood parameters, The posterior over is Gaussian and has the closed form,
compute the posterior credible interval for 0.
past (xppp X E-E- Spoot))
(21)
where
(X 11)
Σ Epost
(22)
where 22 is the prior variance of 0, = K 82IN. where 52 is the variance of the noise, and K is the N x kernel matrix formed by evaluating the Matérn 3/2 kernel on the source spatial locations. The posterior for post is therefore also Gaussian,
T post
(e 2 Ppost X Ep).
(23)
Credible intervals are then computed using z-scores together with this mean and variance.
C Extension To Multivariate Responses
We focus on a scalar response Yn E TR. We expect similar machinery can be adapted to confidence intervals for each coordinate of Y, E RR. Namely, separate univariate problem, and apply our Lipschitz-baas based bias bound and variance calibration to each coordinate. If one desires simultaneous coverage over all D outputs, straightforward Bonferroni correction (i.e. defining a' = a and applying our method to construct for each coordinate) another family-wise error control can be used. An exploration of improvements of this Bonferroni correction approach for multivariate responses would be an interesting direction for future work.
treat each component Y:(4) R
a' confidence intervals
D Simulations Additional Details
In this section, we present figures to visualize the data generating process for simulation experiments. In all experiments, an intercept is included in the regression well the covariates described. All simulation experiments were run on Intel(R) Xeon(R) W-2295 CPU @ 3.00GHI using 36 threads. The total time run all simulation experiments was under two hours. The single covariate experiment took 9-10 minutes to run; the three covariate experiment took around 80 minutes to run, and the Lipschitz ablation study took around 70 minutes run.
D.1 Reported Uncertainty in Simulation Experiments
In Fig. I, we provide error bars for the empirical coverage as well as point estimate. The upper side of the confidence interval indicates the largest value of the (true) coverage such that with probability 97.5%, the empirical coverage would be less than or equal to the observed value. Conversely, the lower edge of the confidence interval indicates the smallest value of the (true) coverage such that with probability 97.5%, the empirical coverage would be greater than or equal to the observed value. This is therefore (conservative) 95% confidence interval for the true coverage.
To calculate the upper and lower bounds, we observe that the empirical coverage is binomial random variable, with parameters equal to the number of seeds and the true coverage. We then numerically
22

---

Shift 0.2
Shift 0.6
Covariate
Expected Response
1.0
0.5
0.0
-1.0
-2
-1 0 -1
Figure 4: Spatial sites for the source (blue) and target (orange) data are shown in the left most plots for different values shift used in generating the data. More extreme values of the shift parameter lead to larger biases in parameter estimation from the training data without adjustment. The third plot from the left shows the covariate surface, while the fourth shows the expected response each spatial location.
invert the binomial cumulative mass function to calculate the upper and lower bounds by performing bisection search the probability parameter.
We also provide =2 standard deviation error bars on the confidence interval width. As the confidence interval width is not necessarily normal this may not be 95% confidence interval for the confidence interval width. But is meant as an indicator of spread of confidence interval widths each method obtains.
D.2 Data Generation for the Single Covariate Experiment
We show example datasets used in the simulation experiment with a single covariate and the simulation experiment in which investigated the impact of Lipschitz constant confidence interval width and coverage in Fig. 4. The left most two panels show the distribution of source (blue) and target (orange) locations for values of the shift parameter. Large positive values of the shift parameter lead to target distributions clustered to the top right, values close to zero lead to the target locations being approximately uniformly distributed and large negative values lead to target locations clustered to the bottom left. The third panel from the left in Fig. shows the covariate plotted as function of spatial location.
((S0),s()) = 5(1) + S(2)
(24)
The right most panel shows the conditional expectation of the response plotted as a function of spatial location,
E E
5() + 5(2) + ((( + (5(2))3).
(25)
The gradient of this conditional expectation is,
1 + 5(1) 1 + S(2)
(26)
which has norm
(1 + 8(1)2 (1 + 5(2))2.
(27)
This obtains maximum of 2 on [-1, 1J2. and so this is the Lipschitz constant of f(S) |
D.3 Data Generation for the Three Covariate Experiment
In Fig. we show data generated for the three covariate shift experiment. Source and target locations are generated as the one covariate simulation, but with N 10,000 (M 100 is still used). These not shown in Fig. 5.
23

---

Covariate
Covvriate
Expected Response
10
-05
-10
Figure 5: The first 3 plots from the left show the covariate surfaces, while the fourth shows the expected response at each spatial location for the second simulated experiment. The source and target locations (not shown) are the same as Fig. 4, though with N 10,000.
The covariates are, from left to right in Fig.
(1) sin(s(1)) cos( (5(2))
(28)
(2) cos(5(1)) - sin sin (5(2))
(29)
(3) s(1) + 5(2)
(30)
The conditional expectation the response (right most panel in Fig. 5)
E| = 5(22) = (1) (2) + マ (() + (5(2))3)
(31)
(sin(5(1)) + cos(( (cosss))))) sin(s(2)) + (15 (8(2))
(32)
The gradient of this conditional expectation is
cos(25(1)) - sin(((()) + 5(2) 5(1) - cos(2 + sin( 5 5(2) 5(2)
(33)
We see that both arguments of this gradient are less than 3 in absolute value, and therefore the norm of this Lipschitz constant is less than or equal to v32 + 3 3/2.
E
Tree Cover Experiment Additional Details
this section we provide additional details the real data tree cover experiment, as well as figures to visualize the data and additional experimental results. The tree cover experiments was run Intel(R) Xeon(R) (R) W-2295 CPU @ 3.00GHz using threads. The total time to run the experiment was under 5 minutes. The Lipschitz ablation study also took less than minutes to run.
E.1 Tree Cover Data
Our analysis draws on data from Lu al. [37], who manually labeled 983 high-resolution Google Maps satellite images for tree cover percentage. While no license is specified, we used this data with permission the authors. These images were selected from random locations within the 2021 USFS Tree Canopy Cover (TCC) product [68]. This dataset is public domain. We follow Lu et al. [37] and use three covariates:
1. Global Aridity Index (1970-2000)) Averaged at 30 are- arc-seconds resolution [65]. This index is calculated as the ratio of precipitation to evapotranspiration, with lower values indicating more arid conditions.
2. Elevation: Provided by NASA's 30-meter resolution dataset [43].
Slope: NASA's 30 m Digital Elevation Model. Also provided by NASA's 30-meter resolution dataset [43]. While we could not find specific license information, as the Slope and Elevation datasets are produced by US government agency (NASA), we understand this data be public domain following section 105 of the Copyright Act of 1976.

---

Figure 6 provides a visual overview of both the tree cover and covariates, As preprocessing step, we convert the (latitude, longitude) coordinates of each data point into radians. This allows to use the Haversine formula to compute distances in kilometers for the Wasserstein-1 cost and the nearest. neighbor weighting procedure.
Elevation and slope are important factors influencing tree cover worldwide. Generally, areas at lower elevations tend have more tree cover, often due higher temperatures [41], and sloped terrains also support greater tree coverage [59]. As done in Lu et al. [37], we focus on these three covariates and did include additional factors that might affect tree cover. This decision was made because our primary objective is to demonstrate uncertainty quantification rather than to provide a comprehensive explanation of tree cover dynamics.
Figure provides a visual overview of the distribution for both the tree cover and covariates.
Tarpe Veriablee (Thon Dower X0
Artiiiy inces
50
Lartntn
300
30
Larrice
2000
10
125
100
231
Loogituse
Lougtoue
Figure 6: Tree cover response and covariates. The dots represent the 983 locations considered. Top left: distribution of tree cover percentage. Top right: Average Aridity Index, measured as the ratio of precipitation to evapotranspiration. Bottom left: Elevation, measured in meters. Bottom right: Slope, measured in degrees.
E.2 Source and Target Data Split and Spatial Preprocessing
We define target region be the Western portion of the Continental United States. In particular, we consider locations that have latitudes between 25t and 50% and longitudes between - -125 and -110", Within spatial locations in this defined region, we pick 50% of all spatial locations totaling 133 sites - as larget data.
To select the source data, perform the following steps:
Consider all the remaining spatial locations, i.e. exclude the 133 target points from the pool of 983 spatial points. This leaves us with 850 points.
2. Since we are interested in evaluating whether our method and the baselines achieve nominal coverage, uniformly randomly sample 20% of the remaining locations across 250 different random seeds. By doing this, for each random seed we have 170 source locations.
Fig. 7 visually represents the spatial distribution of the source and target locations for representa- tive random seed.
As preprocessing step, we also convert the geographical coordinates (latitude and longitude) of each data point from degrees radians. This conversion is essential because it allows to apply the Haversine formula, which calculates the great-cirele distance between two points on the Earth's
25

---

50
45-
Lattude
40-
35*
30-
Source
Target
120
-70
Longitude
Figure 7: Split of the tree cover dataset in target distribution in the West of the United States. Target locations are shown in orange, while source locations are shown in blue.
surface in kilometers, to compute distances in kilometers for the Wasserstein-1 cost and the nearest neighbor weighting procedure.
E.3 Estimating the Ground Truth Parameters to Evaluate Coverage
In order to evaluate coverage, we repeat the data subsampling process described above 250 times. Ideally, would estimate the coverage as the proportion of these seeds in which the estimand OLS. falls inside the confidence interval we construct for each method. However, we cannot evaluate OLS.p directly, even though we have access the target responses. To account for sampling variability, we use our method and each baseline to construct confidence interval for the difference.
0z - 6p
(34)
where = -xxxxxxxxx--xxxx---
the estimated parameter using the target data (which our method and baselines don't have access to) and 0g is the estimated parameter we compute with each method using the source responses. If this confidence interval contains 0, we count the method as having covered the true parameter, while if it doesn't 't count the method not having covered the parameter. We estimate the variance of 6; using the model-trusting standard errors 8² N E r where are the residuals of the model fit on the training data. We expect this to inflate our estimate of the standard variance of the target OLS estimate if the response surface is nonlinear, the residuals will be larger due to bias. By possibly overestimating and incorporating this sampling variability into confidence intervals, we expect the calculated coverages to overestimate the true coverages. The resulting coverages are shown in the top row of Fig. 3.
In Fig, 8 we also report confidence intervals by calculating the proportion of times that 07 is contained in each confidence interval. 0p is an unbiased estimate for OLss,p whether or not the model is well- specified. However, we might expect that the coverages reported with this approach underestimate the actual coverage of each method's confidence intervals due to not accounting for sampling variability in 07-
E.4 Experiment with Varying Lipschitz Constant
In this section, we provide an additional experiment with these real data where the focus to assess how varying the underlying Lipschitz constant in Assumption 4 changes coverage and interval width. For this experiment we consider 9 different values for the Lipschitz constant, L E {0.001, 0.005, 0.01, 0.02, 0.05, 0.2, 0.5, 1). This range of values is a reasonable range of values that a practitioner with domain knowledge in the tree cover field might be willing to specify. Indeed, the smallest I. Lipschitz constant considered, 0.001, corresponds to assuming that to have 1 1% increase in tree cover we need move 1/0.001 1000 kilometers. This is quite an extreme value, but in some parts of arid regions such New Mexico and Arizona it can be true. On the other hand,
26

---

Aridity Index
Elevation
Slope
1.0
Coneaae
0.5
0.5
0.5
0.0
Average CJ Witth Coverage tof point ostimate)
Aridity Index
Elevation
Slope
1.0
1.0
1.0
0.5-
0.5
0.5-
0.0
0.0
2.0
0.05
0.0
0.0
0,0
GLS
GP BCls
KDETW
0LS
Sardeidi
Oue
GLS
GP BC
KDEIW
OLS
Saudledch
Ouri
GLS
CP nCc
KDERW
OLS
Saidwich
Oure
Figure Coverages for the difference (upper), coverages for the point estimate (middle), and confidence interval widths (lower) for method well as 5 other methods for the West US data. Each column represents a paramcter in the tree cover experiment. Only our method consistently achieves the nominal coverage.
the largest Lipschitz constant considered, 1, corresponds to assuming that to have 1% increase in tree cover we need move kilometer. This also extreme, but in regions in the US where elevation changes very pronounced (c.g. in the Colorado Rockies), tree cover vary sharply over very short distances.
In Fig. 9, we present the results of this experiment. We find that varying the Lipschitz constant within this range does not significantly impact coverage. Specifically, for slope and elevation, coverage remains consistent across all constants except L = 0.001. For aridity index, coverage is low for 0.1 but exceeds 95% for 2 0.2. Meanwhile, confidence intervals become noticeably wider for L > 0.5 while remaining relatively stable for smaller values.
E.5 Additional Experiment: Target South-East US
In this experiment, define target region in the Southeastern portion of CONUS locations with latitude in the range (25, 38) and longitude in the range (-100, -75). Out of all spatial points in target data. Next, select the source
this region, 50% totaling 118 sites data by taking a uniform random sample of 20% the remaining spatial locations, repeated over 250 random seeds to assess coverage performance. Each seed yields 173 source locations. Fig. 11 illustrates the spatial split between source and target data for a representative seed. As before, as preprocessing step we convert the (latitude, longitude) coordinates of each data point into radians.
designated
Results. We report the results for the confidence interval coverage and width Fig. 12. As before, our method consistently achieves or exceeds 95% nominal coverage for all parameters. All competing methods but KDEIW here achieves 95% nominal coverage (or close 95% nominal coverage) for the aridity index and slope parameter when considering the coverage for the 07. O (top row). All competing methods fall short of the nominal threshold for the clevation parameter. The wider
27

---

Aridity Index
Elevation
Slope
1.00
1.00
1.00
0.75
0.75
0.75
Covergee
0.50
0.50
0.50
0.25
0.25
0.25
0.00
0.00
0.00
Aridity Index
Elevation
Slope
1.00
1.00
Coverage (somerrsaiion
0.75
0.75
0.75
0.50
0.50
0.50
0.25
0.25
0.25
0.00
0.00
0.00
20
Avenge CI Waith
0.2
200
10-
100
0.1
0.001
0.000
0.01
0.02
0.05
0.1
0.2
0.5
0.001
0.005
0.1
0.02
0.05
0.1
0.2
0.0
0.001
0.01
0.02
0.05
0.1
0.2
0.0
Figure 9: Coverage and average confidence interval widths over 250 seeds for different values for the Lipschitz constant L E {0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1}. The horizontal axis is shared across all plots.
intervals produced by our method (last row) reflect the trade-off between achieving reliable coverage and maintaining narrower intervals.
In the middle row we show the coverage for the point estimate 4z. Here we see how our method is the only one that achieves nominal coverage for all the parameters. In particular, all the other methods do not achieve nominal coverage for any of the parameters.
Varying Lipschitz constant. Finally, we report also for this experiment results when varying the assumed Lipschitz constant. As explained Appendix F,4, we consider 9 different values for the Lipschitz constant, I, {0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1}. We report the results in Fig. 14. As before, we find that varying the Lipschitz constant within this range does not significantly impact coverage. Specifically, here see that for aridity index and elevation, coverage remains consistent across all constants above = 0.01. For slope, 95% nominal coverage is achieved for
L > 0.05. And as before confidence intervals become noticeably wider for L > 0.5 while remaining relatively stable for smaller values.
28

---

Atidity Todlet
Elombtoo
Sloge
GP RCI
KDEIW
Sandwich
OLS
Ours
GP RCI
KDRIW
Sandwich
OLS
Ours
GP DCLs
G1.8
KDEIW
Sandwich
Ours
BCI
GLS
KDEIW
Sundwich ich
Ouz
BCI
GLS
Sundwich
0LS
Ous
GLS
KDEIW
Sundwich
0LS
Ours
GL.S
KDEIW
Sandwich
GLS
100
0.05
0.00
0.05
2.5
2.5
Figure 10: Confidence intervals for different seeds for the West US. Each row shows confidence intervals for the various methods over the three parameters for given seed. The dashed vertical lines represent the true parameters (black is the point estimate, orange is 95% confidence interval). The blue dots the point estimates for the different methods, and the blue lines are the confidence intervals.
29

---

Latinde
Source Turgel
25
- 120
-110
-100
90
80
70
Longitude
Figure 11: Spatial sites for the source (blue) and target (orange) data. The target data are chosen from the south-eastern part of the CONUS, whereas the source data cover the whole region.
Aridity Index
Elevation
Slope
Cuverage
0.5
0.5
0.5
0.0
0.0
0.0
Coverage (of point estiude)
Aridity Index
Elevation
Slope
1.0
0.5
0.5
0.5
0.0
0.0
Average CI Width
5.0
200.0
0.2
0.0
GLS
GP BCIs
kDEnw
OLS
Sandwih
Our
GLS
GP BCIa
KDEIW
OLS
Sundeich
Oum
GLS
GP BCIs
KDETW
OLS
Sundwich
Oum
Figure 12: Coverages for the difference (upper), coverages for the point estimate (middle), and confidence interval widths (lower) for method well other methods for the Southeast US data. Each column represents a parameter in the tree cover experiment. Only our method consistently achieves the nominal coverage.
30

---

Arility Index
Elontion
Stope
GP RCI
KDEIW
Sandwich
OLS
Ours
GP RCI
KDRIW
Sandwich
OLS
Ours
GP DC1s
G1.8
KDEIW
Sandwich
Ours
BCI
GLS
KDEIW
Sundwich ich
Ouz
BCI
GL.
Sundwich
0LS
Ouz
GLS
KDEIW
Sundwich
0LS
Ours
GL.S
KDEIW
Sandwich
GLS
200
0
200
0.2
0.0
0.2
10
Figure 13: Confidence intervals for different seeds for the South East US. Fach row shows confidence intervals for the various methods over the three parameters for given seed. The dashed vertical lines represent the true parameters (black is the point estimate, orange is 95% confidence interval). The blue dots the point estimates for the different methods, and the blue lines are the confidence intervals.
31

---

Aridity Index
Elevation
Slope
1.00
1.00
Coverage
0.75
0.75
0.75
0.50
0.25
0.25
0.25
0.00
0.00
0.00
Aridity Index
Elevation
Slope
Coverage: (consevative)
1.00
1.00
1.00
0.75
0.75
0.50
0.50
0.50
0.25
0.25
0.25
0.00
Average CI Width
1000
40
0.5
0.0
0
0.001
0.055
0.01
0.02
0
0.1
0.2
0.5
0.001
0.005
0.01
11:12
0.05
0.1
0.2
0.5
0.001
0.005
0.01
0.12
11.05
0.1
0.2
0..5
Figure 14: Coverage and average confidence interval widths 250 seeds for different values for the Lipschitz constant I. E (0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1).
32

---

Proofs
F.1
Derivation of Target-Condittonal Ordinary Least Squares
We derive the OLS (Ordinary Least Squares) estimand for the given optimization problem:
M
fos arg min E DERP
Σ
0
x
Sm
(35)
where Y* is unobserved, and X is observed and fixed function of S.
First, expand the squared term inside the expectation:
(Ym
x
(2m))²
X 0
+ (0T X x²)
(36)
Substitute this back into the expectation:
YY 0 xx
S
E
Σ (Y
2Y 0 X (0
(37)
Since the expectation is linear, separate the terms:
M
E
(Ym)
S
2E
S
Σ
(0 x
(38)
Now we can simplify each of the terms as follows:
The first term, Σm E ((Ym) S does not depend on 0, so it can be treated as constant with respect the optimization problem
The second term, 2 ΣM expectation, the fact that is not random, and Assumption 1:
O X Sin can be rewritten using the linearity of
M
20
X
S
(39)
m=1
The third term, E EM (0T rewritten as:
X
S
is non-random by Assumption and can be
M
00 x
(40)
m=1
The optimization problem then becomes
M
Go = arg min PERP
constant
TT
x
y
Σ
(
x2
(41)
m7
m=
m=1
And since the constant term does not affect the optimization, we can drop it to get
M
M
OLLS = arg min BE3" -200
X
Y
5
+
Σ
(42)
m=1
m=l
To find the minimizer, we take the derivative of the objective function with respect to and set it to zero:
0
M
M + Σ
20
E
y+ m
X+ )2 0 m
(43)
00
m
m
m=1
12
M5
-
Σ
X E
y+
S
+ 2 Σ
X+ (x* m
T
00
(44)
33

---

And by inverting this and using the fact that we assumed that X X is invertible
M
M
OIs
(x
X
Y
S
(45)
m=1
which in matrix form can be written as
OLs
(46)
as in Eq. (4).
F.2 Proof of Theorem 7
In the main text, we stated our confidence interval for the INN choice of V. We now state version for general non- non-negative matrices with columns summing to 1.
) and (Sn, Xn, Yn Tn)n satisfy Assumptions I to 4. Define matrix with non-negative entries with columns summing to 1, define " = wV. As in the main text, take 6% = (x+ X+TWY. Define c = alli"1 Il with o² the variance of the additive noise from Assumptiin 2. Define the (random) interval I as Algorithm 1 with known o². Then with probability at least a, OoLS,p I That is, 1 has coverage (conditional on the test locations) at least 1 aa
Theorem 11. Suppose (S) X
M IR EN
as in Algorithm 1. For any V m
Proof. We being in the main text. We show the difference between estimand and estimator is normally distributed. To that end, we decompose the difference between our estimand and estimator into bias term and mean-zero noise term.
N
OLSS, - 0
Σ m=1
Wm f(S
Σ ==1
L
Y..
(47)
for w := of (x* X+T of Eq. (47) can be written
wi E IRN, By Assumption 2, the righthand side RM and V
M
M Σ m=1 ((5
Σ n=1
Y,
Wm ((S
Σ Un ((S.)
(48)
hias
randommess
Since the spatial locations are fixed, the bias term is not random and can be written as E We can calculate the variance directly,
VOOLs
0.1
Up 6n
2 [(n] 021|0""1
(49)
We used the en are independent and identically distributed with variance σ² (Assumption 2). It follows that
ols., 8 ~ o
(50)
To bound the bias b, we use Assumption 4 to write
|b|
< sup g
Σ
(n.((S)))) Σ EoN9(Sm)
(51)
We therefore apply Lemma 6 to complete the proof.
34

---

F33 Computing an Upper Bound on the Bias of Our Estimand with Wasserstein-1 Distance
First of all observe that if Em= Wm = T.En=) "n 1, rM 2 0 for 1 < S N, and Wm 2 0 for S m I M then the supremum in Eq. (7) would be equal to a Wasserstein-1 distance by Suatorovvic-Ruubistte
duality [71, Theorem 5.10, Case 5.16].
Next, consider what happens if Em Wm - vx + O. We show that the right-hand side of Eq. (7) infinite. Assume by contradiction that there exists > 0 that upper bounds this supremum. Because Em=1 Win > U/ * 0 and all constant functions are L-Lipschitz, for any
y 0, taking, for all (S) G _m=1 Wmm -n
C+7
M
N Σ (Sw))
M Σ
N Σ
sup OEFL
Σ
g
Un
Wm
G n2
m=11
m=
n=1
+
(52)
This contradicts the assumption that C' upper bound on the supremum. Because C was arbitrary, the right hand side of Eq. (7) is infinite, as desired.
Our assumption that V is non-negative matrix whose columns sum to avoids this situation. We formalize our upper bound on the bias in the following proposition.
Proposition 12. Suppose that Em=1
L mmn
for all Let RM and v (WW) E B
M
N
sup gEFL
(m)()mm)
Σ
2
(Sn))
m=1
n=1
= ALW1
Wm
8g:
+
8s.: Σ Wmm
85s + Σ
U ds
(53)
A
mE!
NEI
mEJ
nEJ
where 1 = (1 M : w. 2 0}. 1' -lsis
N 10 < <0},
[1 j S M w/ 0} and
J' {1= v 0} and A =
|mmml + Σ l0
Proof. First, observe that Σ Σ n
En
M8 Σ m=
Wm V m.n
Em= Um Ex=1 Vm,n = Σm=1 Wm-
Next, we normalize the weights to sum absolute value to 2, and rescale the function class to consist of 1-Lipschitz function. Define A =
M
|wm| + Σm 1 [xx
l).
Then
M
N Σ Lg(Ss)
sup Σ gEFL
Aing(S) (St)
Σ
v un
sup (Sa)
Σ m=1
Wm Lg(
(54)
g
m=1
n=1
n=1
M
N
w
Wm
AL sup
g
(S
g(nn)
(55)
9EFi
where have used that function is I.-Lipschitz if and only if it can be written by scaling 1-Lipschitz function by L.
Define 1 (1 < i < M w, 0}. 1' (1sisNN: x < = (1 <isM my 0} and
J' = (1 وj < N of 2 0}. Then
M Σ
N
AL sup gEFi m=1
Wm
((S))
Σ
n
gSS))
(56)
A
n=1 A
AL sup gEFi
Σ
Wm
g(S
+ Σ
2
(SS)) -
Σ
W
((S )m + Σ b
g(Sn)
A
A
A
A
mE!
n
mEJ
neJ
(57)
35

---

Because E 00 = =E m= tmm and and J/ partition the index sets,
Σ
Wm +
Wm +
Σ
(58)
mEI
REN
mEJ
MEJ
And because the set I,1', ,, J' sort the indices into positive and negative parts
J,,
E Wm
Σ
Wm
n
+
U
+
1
24.
(59)
mel
MEJ
REI
REJ
Therefore,
Wim
Wm
V
8s.
and
(60)
A
A
mE!
nEI
MEJ
nEJ!
are probability measures. We can apply cantorovich-Rubinsteer duality to write,
B ALW1
Σ
Wm
85
ō5.:
Wyp 8s
+ Σ
(61)
TEE
WEI
MEJ
WEJ
where W1 denotes the I-Wasserstein distance.
We compute the Wasserstein-1 distance by linear programming, see discussion in Appendix B.2. Scalable upper bounds could also be computed by exhibiting coupling between the measures (for example by solving an entropy regularized optimal transport problem). See [49, Chapters and 4] for details on computation of Wasserstein distances.
F.4 Proof of Lemma 6
Proof of Lemma 6. We aim prove that the interval [-B A, B+ is the narrowest 1 confi- dence interval that is valid for all b E [- B, B) where A is the solution of T (A) - + (-2B/ë 18/ē-1 A)
1 &.
Ensuring Coverage Probability. Suppose that X6 ~ N(b, (22) for E [-B - A, B + AJ. Then,
Pr(X, A,A))- = +b (B+A-D +b
B -B- (BB-A-b) Č
To construct valid confidence interval for any E [- B), we require that
Pr(X, E [-B -, B A 2 a, vb E I- B). This ensures 1 - a coverage over all possible values of b in [-B, B, B}.
Reduce the problem to Worst-Case Coverage. To find the narrowest interval, we identify the worst-case value of that minimizes the coverage probability. Let
B +
A
(-B-A-D)
-
Cbb; A) dp
denote the coverage probability of the interval [-B B, B + A for X6 ~ (bb,²²). In order to ensure the interval is valid for all coverage, want to bound below.
inf C(b:A)
bE[-B.B]
The interval [-B A, B A is symmetric about 0, and the Probability Density Function for Gaussian of mean is symmetric about b. Thus, the coverage probabilities b -B and b = are equal. Consequently, it suffices to consider b E [0, 13).
Moreover, observe that C(b; A) is a strictly decreasing function of on [0, B] since (i) Φ (B+A-b) decreases b increases (because B + A b decreases and ((2) is monotonic) and (ii) & -B-
36

---

also decreases as increases (because -B A b becomes more negative). Thus, C(b; A) is strictly decreasing function of [0, B]. The minimum value C(;; A) occurs b B.
Ensuring coverage in the worst case. At the worst-case value h = B, the coverage probability is:
2B A
C(B:A) =
To ensure that the interval [-B - -A, B + A] achieves at least 1 coverage for all b E [-- [-B, B), we solve:
A + Φ
2B - (-26--
Č -l-a.
(62)
This guarantees that the interval is valid for all b and achieves the desired coverage level.
Narrowest interval. The narrowest interval corresponds to the smallest A that satisfies the Eq. (62).
By construction, any smaller A would fail to achieve the required coverage for b = =, violating the validity condition.
F5 Proof of Proposition 8
In this section, we prove Proposition 8. For simplicity of exposition, we prove the result for = [0, 1) The result generalizes to any spatial domain which a compact metric space.
Proof of Proposition 8. Using Assumption and expanding the quadratic form (en + (f(S,) g(S,)), we have
o - o== -N + CN.
(63)
where CN inf gefz N 1
- g(S,))² +NN En g(S)) and ZN
0². Since NN is an average of independent and identically distributed variable, and since R[[x] = O, the law of large numbers (LLNN implies Z - in probability. Because ZN 0, if ÇN - in probability
N 2 c2
conclude by Slutsky's Lemma that 0 - o² in probability. Therefore, the remainder the proof involves showing ÇN + 0 in probability.
Define SN to be the empirically centered version of S, that JN = \ ΣM f(S.). Then since the space of Lipschitz functions is invariant to shifts by constant functions
CN inf gEft N
(SN o(s)))
-
(64)
Define the process.
N Σ (fN(S))) - -(())))
N Σ n=1
TN (g)
- (f.
(65)
N
N
n=1
so that CN is the infimum of TN. (N(fN) = =0. Therefore, for ÇN < 0 almost surely. It remains to show that for any O, lim Pr((N < <) 0.
The essential challenge to showing that for any 5 < O, limy- Pr(CN < 5) + is the infimum the space of Lipschitz functions. Our proof has three steps. First, we show that it suffices to consider a subset of the space of Lipschitz functions with bounded infinity norm. Second, we show that this space is compact as subset of L Third, we show that because the infimum then compact set, it can be well-approximated by a minimum over a finite set of functions. And then a union bound suffices to prove the claim.
Step 1: It's Enough to Consider Bounded Subset of Lipschitz Functions.
Because fN has empirical mean 0 and is continuous, by the intermediate value theorem it takes on the value 0 somewhere on [o. 1 Because fN is L-Lipschitz and defined on a set of diameter vD, and 0 somewhere inside this set, || < LVD.
37

---

Define the set xz (2LVDD222 where (x) denotes the space of functions uniformly bounded by constant on [0, 1| By subadditivity of measure, for any
Pr(CN < <) < Pr
inf TN(g) <6) + Pr
inf
TN(g) <8
(66)
gE
gEFLNFE
$ Pr
inf TN(g) <
+ Pr
inf
[g(g)
(67)
DEFL
gEFL
We first consider the second term in this sum and show it tends to O. We apply a crude Cauchy-Schwarz bound the second term in Eq. (65) so that for any g.
N Σ n=1
N
(g) N
Σ (f(S.) -( n=1
1
Σ
2
(68)
N
N
N
1=1
Therefore, (gg) < 0 implies
E (
(69)
For any E Fr. Fr. because y takes on a value at least 2Lv 2o and is L-Lipschitz, is larger
than LVD 2o² over the entire unit cube. And because [ S LvD, (S,) - g(S.) 2 20² for all n. Therefore for all FL FL
N Σ
- 20
(70)
We conclude that
lim Pr inf N-0
Tx(2) <0
S lim Pr
Σ n=l
2
20 =0.
(71)
DEFLVFF
N
where the final inequality is the law of large numbers.
Fr. is Compact Subset of the Space of Bounded Functions with Sup Norm. All that remains to show is that for any
Pr inf TN (g) <8) + O. DEFL
(72)
The idea (following standard arguments made with empirical processes) that we can take a cover of Fr. of finite size, such that each element of (g(gg) almost constant over elements of this cover. This essentially lets us approximate the infimum with a minimum over a finite set, up to small error. And once the problem reduced to a minimum we can apply union bound and the law of large number. We will formalize this in the next paragraph, but we first show that Fz compact.
Because every Lipschitz continuous function is equicontinuous, functions in Fz. are pointwise bounded by 2LVD + 20², and [0, 1] is compact, we may apply Arzela-Ascoli [58, Theorem 7.25] to conclude that Fr, with the sup norm is sequentially compact. It is therefore compact, as sequentially compact metric space compact.
Step 3: Reduction to a Minimum over a Finite Set and Union Bound. For any 8' < and any g Ft.
1
lim PI(TN( (g) <0) <8) N
lim Pr
Mec.oh C..(SS)) < 8') =
<8¹
(73)
N
where the final equality is the law of large numbers.
38

---

Therefore, for any finite set C C FL.
N
lim Pr(min (g( <8) -
lim
Σ DEC
Pr
(ng(
<8)
(74)
N-00
SEC
N 40000
n=1
lim Pr
(ng < 0.
(75)
N
N
SEC
We used countable subaddirivity in the inequality.
Now for any ? > O, there exists a finite set of functions C2 C FL such that for any 9 E FL, there exists g' E C, with lig o' xc y. And since (9) is pathwise uniformly (in N) continuous Fr. equipped with sup norm,
inf TN(9)= min TN p(()
(76)
gEFL
SEC,
where is nonnegative function such that lim p(7) 0. Therefore, for any 8 < O, can find y such that p(7) - For this y, applying Eq. (74) with 8' 1 allows conclude that lim
Pr( inf (g( - 8) = 0.
(77)
N gEFL
This is uniform law of large number for the class of Lipschitz and bounded functions. More quantitative results are likely possible using empirical process theory, see Wainwright [73, Chapter 5].
Proof of Asymptotic Coverage (Corollary 9). We now prove Corollary 9. We begin by recalling the definition of an asymptotically valid confidence interval.
Definition 13. We say sequence of (rundom) intervals (IN) has asymptotically valid coverage of at level (1 - a. if
lim Pr(0 E IN) 1-
(78)
N
Corollary follows from Theorem 7 and Proposition 8 by the following lemma, which is a special case of Slutsky's ' lemma.
Lemma 14. Let o² > O. Suppose that t.s.,p 8% - bn ~ NN, where (bx (bw) is fixed
sequence. Suppose 0x converges in probability to 0². Then,
1
(FOLS. 0g bN) +0000
(79)
0%
where convergence is in probability.
Proof. The result is special case of Slutsky's 's lemma, using 03x + 0²200 0.
Proof of Corollary 9. Define by = =M= wmf(S() Em U f(St). Because OL. 6
(bb o² 2 IB))
(BOLs 6% -b) ~ N'
(80)
If 2 0, then 0% 0 for all N, because the conditional expectation is L-Lipschitz function leading to squared error the minimization algorithm used to calculate Ø2-- Therefore, the resulting confidence interval has coverage
-a) for all N by Theorem 7.
For 22 > 0 we apply Lemma to conclude that
1
(0rss 6/ bN) + (0,
(81)
where convergence is in probability.
Convergence in probability implies convergence in distribution, and therefore convergence of quantiles at continuity points. The Gaussian CDF is continuous. Therefore, the quantile computation in Lemma 6 using y in place of x² produces asymptotically valid confidence interval in the limit as N + co.
39

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claims of the paper are that 1.) Existing methods for constructing
uncertainty estimates for associations do not provide correct coverage under misspecification and covariate shift. 2.) We provide method that does. The first claim is supported in two parts. The first part is our experiments, in which we show that existing methods for confidence intervals for linear models fail in this setting. The second part is that, argue in the introduction, post hoc interpretability methods for black-box machine learning methods don't account for uncertainty in the response and essentially don't have clear mechanism for doing so. The second claim also supported two parts. First, it is supported by theoretical results (in which prove that under smoothness assumptions and Gaussian noise, our confidence intervals achieve nominal coverage). Second, the second claim is also supported by our experiments, in which we verify both simulated and real data our method produces valid (conservative) confidence intervals.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Two limitations of our method discussed are: 1.) The need to choose Lipschitz constant. We discuss this in detail in a simulation experiment and in the real-world experiment, looking the effect of this choice the constructed confidence intervals. The choice L = O leads to an assumption that there is no bias, and so leads to classical confidence intervals (for the weighted estimator consider). 2.) Our method, described, only works for linear models. This is clearly outlined in Section 2. As detailed Section 1, linear models a useful tool for describing associations between variables, even in instances where non-linear relationships exist. This is especially the case when there is limited information in the data, and a.) quantifying uncertainty crucial to the validity of scientific inference and b.) complicated machine learning model might not be able to precisely and fully capture the response. 3.) Computationally, our estimate of the noise variance is slow for large datasets. We instead derived and used a different estimator for 10000, but did not prove that it consistent. We discuss this point in Appendices B.7 and B.8.
Justification:
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well vell-specification asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
40

---

The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: The complete set of assumptions for our method are stated in Section 2. Complete and detailed proofs of theoretical results are stated in Appendix F.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: Details of experiments are described in Appendix Anonymized code is available at https: http:////nonoouuu/// 4open. wwweensscennee/eee Meteieee////DDI-eeIrSSSSS
Guidelines:
The answer NA means that the paper does not include experiments.
the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same
41

---

dataset, provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c)
If the contribution is new model (e.g.. a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g.. with an open-source dataset instructions for how construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We do not generate new datasets. Citations are provided for the data used, and code to download the data and run experiments is provided at https: / /anonymous. 4open sopen.scienn www.sseeeeeeeeieeeeeeeeeeeeeeeee /tteeeer//DIIII-e Neur
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/ Mncides/CodeSubbnn
connolicy) for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / //nips. wctppclccpccc/ccccccccccccccccce hptppccgguudds/////ppppppo gwwggg/iee///// /ewggdd///edoeuenniiiiieeeeeeeee (oodeSumiisssiiicccccooco.oo...e asssonoliccy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question:
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
42

---

Justification:
All experimental details are described in Section and Appendix B.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
For the simulation experiments, we provide confidence intervals for the coverage and confidence interval width. Interpretation of these error bars is in Appendix D. In the real- wolld experiment, we took into account sampling variability in estimation of the ground truth parameter, but did not report the confidence intervals for the empirical coverage confidence interval widths.
Justification:
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
We describe computation times for the simulation experiments and the com- puting environment used in Appendix D. For the real-world experiment, computational resources are deseribed in Appendix E.
Justification:
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
43

---

The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics attps://neurre
Answer: [Yes]
Justification: We do use or introduce datasets involving human subjects. Datasets used are publicly available. The method developed does not present significant and specific social concerns.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e-g.. there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: The introduction motivates the importance of (reliable) uncertainty quantifica- tion in making scientific inference trustworthy. We see this as the primary societal impact of our work.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: The paper focuses on quantification associations in spatial data through linear modeling framework. This does not pose high risk for misuse.
44

---

Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: A description of the dataset used, including citations to the relevant data sources and license information, are provided in Appendix E.I. Relevant code packages cited within the implementation details section (Appendix B).
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The code released contains a readme file with instructions on how to recreate experiments, well setup file indicating dependencies.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details the dataset/code/moodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
45

---

Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA].
Justification:
No human subjects were used.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Question:
Answer: [NA].
Justification: No human subjects were used.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA].
Justification:
LLMs were not used in the core method development for this paper.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our L.MM policy /https://neut //tu:r////eusss/e ce/ for what should or should be described.
46


Reviews:
[{"review_id": "FCfwSwCNEE", "reviewer": "Reviewer_W9eu", "summary": "The paper considers constructing confidence intervals for a linear regression problem under covariate shift, where both training and test data is associated with a \"spatial information\". \n\nConcretely, with $\\mathcal{S}$ being some metric space (which should be thought of as low-dimensional - e.g. earth surface), the authors consider a model where there is an (unkown) Lipschitz function $\\chi : \\mathcal{S} \\to \\mathbb{R}^D$, mapping the spatial location into covariates, and the (noisy) observations are given as $Y_i = f(S_i) + \\varepsilon$ for some Lipschitz function $f: \\mathcal{S} \\to \\mathbb{R}$. The algorithm observes $N$ tuples $(S_i, X_i, Y_i)$ for $S_i \\sim \\mathcal{D}_1$, such that $X_i = \\chi(S_i)$, followed by $M$ tuples $(S_i\\^*, X\\_i\\^\\*)$ where $S_i$ is drawn from a potentially different distribution $\\mathcal{D}_2$. \n\nThe relation between $Y_i$ and $X_i$ is not assumed to be linear -- the authors nevertheless suggest fitting a linear model in this misspecified case, attempting to approximate $Y_i$ as $\\langle \\theta, X_i\\rangle$ as well as possible - as they claim even if it does not necessairly yield the best possible predictions for $Y^*$, it is more interpretable than alternatives, and often might provide qualitative understanding of which covariates $X_i$ have impact on $Y$.\n\nAs such, they attempt to provide a confidence intervals for coordinates of the least square linear regression parameter at the \"target locations\" $S_i^\\*$\n$\\theta_{OLS} := \\mathrm{arg min}\\_{\\theta} \\sum\\_{i \\leq M} \\mathbb{E} \\left[(Y\\_i^\\* - \\theta^T X\\_i^\\*)^2 | S\\_i^\\*\\right],$ \nwhere, again $Y_i^\\*$ are not observed by the algorithm, (and averaged over in the formula above)\n\nTo this end, they approximate $f(S\\_i^\\*) = \\mathbb{E}[Y\\_i^\\* | S\\_i^\\*]$ as a weighted average $Y\\_i$ for nearby points $S\\_i$ in the training set (for experiments, they used simply putting weight $1$ on the nearest neighbor and $0$ on all remaining entries), producing estimated vector of responses $\\hat{Y}\\_i^\\*$. They use this to compute $\\hat{\\theta}\\_{OLS} = (X^{\\*T} X^\\*)^{-1} X^{\\*T} \\hat{Y}^\\*$, and calculate the confidence interval (valid under the relatively weak assumption that the response $f(S)$ is Lipschitz in $S$) for coordinates $\\hat{\\theta}\\_{OLS, p}$, by bounding bias introduced by approximating $Y\\_i^\\*$ with $\\hat{Y}\\_i^\\*$. Thanks to the fact that the $\\hat{Y}\\_i^\\*$ was taken as weighted averages of nearby points, obtaining bound on a bias reduces to a specific Wasserstein distance calculation.", "questions": "See \"Strengths and weaknesses\"; Are there some guidelines for choosing the Lipschitz constant for $f$? How likely is it in practice that the space is covered densely enough by the training data to provide useful confidence intervals? What is the effect on the mechanism for the choice of weights for the approximation of $\\mathbb{E}[Y^* | S^*]$, on the results (and especially, how it depends on the density of the sampled points $S_i$ in the space $\\mathcal{S}$?)", "limitations": "Some limitations have been addressed.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "The paper is extremely well motivated, with both natural and relevant setting - authors explicitly provide a number of applications across various sciences for which the data is associated with spatial information, and the distributions of the spatial variable differs on the training and test data. The linear regression is already ubiquitous in sciences, in particular due to being more interpretable than more complex models. As such, having reliable confidence interval for regression coefficients in the misspecified case is particularly helpful.\n\nAssumption that the covariates are just a (deterministic) Lipschitz function of the spatial variable is a relatively strong limitation, but it has been explicitly acknowledged by authors - this work is a step in the right direction, and as they mention, there is a room for further research. It seems reasonable to assume, on the other hand, that in many scenarios the conditional expectation of the response variable on the spatial variable is Lipschitz in the spatial variable.\n\nIt is unclear to me, how to pick a specific Lipschitz constant for the dependence of the expected response on the spatial variable, in a given practical situation, and moreover - how often one should expect that the space $\\mathcal{S}$ is covered by a training data within the radius much smaller than the chosen Lipschitz constant - this already shows up in Figure 3 (the comparison of this method with other methods for the tree cover experiment): while the proposed method, is the only one of the compared ones that achieves a nominal coverage (i.e. confidence intervals are reliable), the width of the confidence intervals is an order of magnitude larger than the one produced by, for example, GP-BCI, and the coverage is in fact way larger than one would intuitively deem necessary - which might lead to unnecesairly rejecting many potential scientific results as statistically insignificant - arguably, the problem less prevalent than the opposite problem of false-positives in science.\n\nThe proposed framework allows for approximating $\\mathbb{E}[Y^\\* | S^\\*]$ as a weighted linear combination of nearby training points, but in most calculations a very simple weight distribution is used, specifically the \"nearest neighbor\" weight $1$ put on the nearest neighbor of $S^\\*$ in the dataset. Slightly more discussion on why this choice was made, and how it affects the results, particularly depending on how densely space $\\mathcal{S}$ is covered by the training data. This seems like a somewhat unnatural choice, as it approximates a Lipschitz function $f$ by essentially piecewise constant function on the Voroni diagram given by the data points $S_i$ -- which is far from Lipschitz. The authors briefly mention $K$-nearest neighbors as a different choice; another natural choice would be to use weights proportional to, say, gaussian kernel i. e. $w\\_i \\propto \\exp(-d(S^\\*, S\\_i)^2 / L^2)$; it is slightly more discussion (and potentially a follow-up work) on how this choice affects the results would be great.", "was_revised": true, "final_justification": "After clarification, I continue to recommend acceptance.", "rebuttal": "We are grateful to the reviewer for their careful reading, insightful comments, and detailed, thoughtful feedback. We are excited to hear that the reviewer found our paper well-motivated and our setting natural and relevant.\n\n**Assumptions on the covariates:** (*“Assumption that the covariates are just a (deterministic) Lipschitz function of the spatial variable is a relatively strong limitation”*)\n\nWe see two potential assumptions that may be of concern here: (A) assuming the covariates are Lipschitz in the spatial variables, and (B) assuming that the covariates are a deterministic function of the spatial variables. To address (A), we want to start by clarifying that we do not require a Lipschitz assumption on the covariates as a function of the spatial variables. The (only) smoothness assumption we make (namely, Assumption 4) is on the conditional expectation of the response given the spatial variables. For point (B), we agree with the reviewer that requiring the covariates to be a deterministic function of space (namely, Assumption 1) means our method may not apply straightforwardly to certain types of data (e.g., demographic data with multiple individuals in a single household). However, we think there remain many practical settings of interest where Assumption 1 is not problematic; for instance, consider covariates such as precipitation, humidity, temperature, or particulate matter concentration — and either treat time as a spatial variable, average across a particular time range of interest, or consider a fixed time point. \n\n**Choice of Lipschitz constant:** Choosing the Lipschitz constant requires user knowledge about the specific domain of the problem. Generally, especially in scenarios requiring extrapolation, the Lipschitz constant cannot be reliably estimated solely from the observed data without additional assumptions — since function values are not available far from observed locations. We emphasize that we selected the Lipschitz assumption specifically because it is relatively interpretable in physical and spatial contexts compared to common statistical assumptions (e.g., that the function lies within the unit ball of a reproducing kernel Hilbert space). We believe that some form of smoothness assumption is often made, implicitly or explicitly, when working with spatial data. And making this assumption explicit is preferable to leaving it implicit. \n\nFor a concrete example, consider the problem of selecting a Lipschitz constant in an analysis where the response is annual average PM2.5 over California. Chow et al. 2006 claim that “Zones of representation for PM2.5 varied from 5 to 10 km for the urban Fresno and Bakersfield sites, and increased to 15–20 km for the boundary and rural sites” where “[t]he zone of representation is defined as the radius of a circular area in which a species concentration varies by no more than ±20% as it extends outward from the center monitoring site”. The annual PM2.5 concentrations in the study area do not exceed 30 ug/m^3. The combination of a zone of representation between 5 and 20 km, and a variation of not more than 30 ug/m^3 within this zone of representation suggests a range of possible Lipschitz constants: 0.25–1.2 (ug/m^3) / km. The authors also point out that topographical and meteorological phenomena contribute to this scale of variation. So we would not expect this proposed constant to be “universal” for problems related to PM2.5, but we might expect that this range of Lipschitz constants is a reasonable starting point for other studies involving annual average PM2.5 with similar weather and topography to California. We showed in our real-data analysis that a range of Lipschitz constants can still produce qualitatively similar results (Figure 9) and correct coverage. To err towards the side of a conservative analysis, we would recommend a user selects the largest Lipschitz constant in this range (i.e. 1.2 (ug/m^3)/km).\n\n**Width of confidence intervals:** We appreciate the reviewer’s point regarding the substantial width of confidence intervals produced by our method in certain experiments, such as the tree cover experiment (Figure 3). We agree that excessively wide intervals could limit practical utility and potentially lead to overly conservative conclusions, rejecting otherwise meaningful scientific findings as insignificant. We highlight that the increased width arises naturally from our explicit control of bias under extrapolation conditions, especially when the test locations significantly differ from training locations. The large interval widths observed in Figure 3 result directly from this extrapolation scenario, reflecting genuine uncertainty rather than methodological conservatism. We also want to point out that our experiments clearly illustrate that existing confidence intervals typically fail to achieve the nominal coverage — and often have near-zero coverage (as in Figs. 1 & 2). Thus, our method yields the narrowest intervals among all evaluated methods that consistently achieve nominal coverage. We emphasize that producing narrower intervals by sacrificing validity is trivial: a zero-width confidence interval achieves minimal width but does not yield meaningful coverage guarantees. Hence, the widths of invalid intervals are not directly comparable to those of valid intervals. We also note explicitly that our intervals are practically meaningful: in many cases shown in our experiments, intervals are narrow enough to exclude zero, providing clear scientific conclusions about the sign and magnitude of effects. We will explicitly discuss and illustrate these points in the revised manuscript.\n\n**Practical coverage and training data density:** The reviewer rightly raises a practical question regarding how densely the spatial domain must be covered by training data to yield useful confidence intervals. In practice, we expect that our intervals will be relatively narrow whenever the data spacing is substantially smaller than the scale implied by the Lipschitz constant (the distance scale over which meaningful response variations occur). Conversely, if sampling is sparse relative to the Lipschitz scale, intervals will naturally widen. Since this increasing width represents legitimate uncertainty due to extrapolation, we would argue that our intervals are still useful in this case, just (appropriately) wider than the dense-data case. We will explicitly discuss this point in the revised manuscript, clarifying when we can expect our method to provide narrow intervals and how this relates to spatial data density and the Lipschitz constant. \n\n**Effect of weighting mechanism choice:** The reviewer correctly points out that our method allows general weighting schemes for approximating the conditional expectation. Our current choice — nearest-neighbor weighting (weight = 1 for the closest point) — is intentionally simple, transparent, and yields the tightest possible Lipschitz bound on bias. It simplifies derivations and clearly demonstrates our theoretical guarantees. We agree that alternative weighting schemes (such as k-nearest neighbors or the Gaussian kernel suggested by the reviewer) would also be natural choices. These methods produce smoother approximations at the expense of slightly looser bias bounds, potentially reducing variance and practical interval widths. The benefit of smoother weighting schemes generally becomes more pronounced as the density of spatial sampling increases since smoother weighting schemes naturally leverage multiple nearby points, reducing variance while controlling bias effectively. We believe exploring this tradeoff is an interesting and practically relevant direction. We will include additional discussion and comparisons between these weighting schemes in the revised manuscript, explicitly highlighting how data density impacts these trade-offs. We also agree with the reviewer that follow-up work exploring these weighting schemes would be very valuable.\n\nWe hope these clarifications address your questions.  Thank you again for your constructive feedback!\n\n**References:**\n\nChow, Chen, Watson, Lowenthal, Magliano, Turkiewicz, Lehrman. PM2.5 Chemical Composition and Spatiotemporal Variation During the California Regional PM10/PM2.5 Air Quality Study. Journal of Geophysical Research. 2006."}, {"review_id": "sE7TfMlHEX", "reviewer": "Reviewer_pd1Y", "summary": "This paper proposes confidence intervals for association in spatial setting where the linear model is misspecified. Additive Gaussian noise assumption and covariates are fixed function of spatial location assumption are imposed on the data-generating process. Further, under the assumption that the response is a Lipschitz function of the spatial location is assumed, and the authors provide confidence intervals with valid coverage based on these conditions.", "questions": "In addition to the questions raised in above, I found the Assumption 1 very restrictive. It effectively washes away the need for covariates in this setting. Assumption 1 states that for fixed spatial location parameter in the model, the corresponding covariate value is always the same. I found this assumption unreasonable across many applications that the authors mentioned (e.g., air pollution). The explanation that the authors provided in the paper does not directly justify why this assumption holds in applications.", "limitations": "Yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is well-organized, with clear logic between different sections. I didn’t check the mathematical details of the paper, but it seems solid at first glance. I do not work in the spatial association area, so it is hard for me to judge the contribution of the paper given the assumptions made in the paper. \n\nI found the introduction poorly written. The second to the fourth paragraph talks about possible relationships between source and spatial location, and lacks focus. It is unclear what the authors are trying to convey. I would recommend the authors to focus directly on the linear model and discuss other possible models either in related work or discussion section. \n\nI found the claim on the contribution of the paper in accurate: the authors claim that their result hold under non parametric assumptions on the data-generating distributions. However, as stated in their assumptions 1 and 2, additive Gaussian noise assumption is clearly imposed on their data-generating process. Additionally, they require smoothness assumptions on the DGP. The statement is very misleading.", "was_revised": true, "final_justification": "The provided justification by the authors further clarified the contributions of the paper. I was not able to check the mathematical details of the paper, so i will keep my score at 4.", "rebuttal": "We thank the reviewer for their thoughtful feedback. We are glad to hear they found our paper well-organized and mathematically rigorous. Before responding to individual points, we recall that our goal in the present paper is to provide valid confidence intervals for the signed strength of the association between covariates and a response in regression.\n\n**Purpose of paragraphs 2–4 in the introduction.** The reviewer felt that paragraphs 2–4 of the introduction lacked focus. An important premise for our work is that linear regression is still a useful tool, even when the linear model is misspecified (which is always the case). We think making an explicit argument for the relevance of linear regression is particularly important given that the NeurIPS community often emphasizes algorithmic and computational approaches to prediction problems, over model-based approaches. The goal of these paragraphs, and the introduction, is to describe why the reader should care about the problem we study, before getting into the assumptions we make and our proposed solution.\n\n**Use of the term nonparametric.** The reviewer suggests that, because we assume additive Gaussian noise and make a smoothness assumption on the latent function, we should not call our method nonparametric. We believe our use of the term *nonparametric* in the present work follows standard usage in the machine learning literature: namely, we use it to describe a model that is not completely specified in terms of a finite-dimensional parameter vector. For example, Gaussian process regression is generally referred to as nonparametric, even when the latent function is smooth and the observation noise is i.i.d. Gaussian (Giordano et al, 2022, Zhang et al, 2022; El Harzli et al 2024). And papers analyzing kernel ridge regression refer to the method as nonparametric, even if the analysis is done assuming i.i.d. Gaussian noise and a smooth latent function (Cui et al, 2021). As Michael Jordan (2010) wrote, “The word ‘nonparametrics’ needs a bit of explanation. The word does not mean ‘no parameters’; indeed, many stochastic processes can be usefully viewed in terms of parameters (often, infinite collections of parameters). Rather, it means ‘not parametric,’ in the sense that [...] nonparametric inference is not restricted to objects whose dimensionality stays fixed as more data is observed.” We contrast this usage with the term *distribution-free inference* (Cherrian, 2024; Correia 2024), which is frequently used for inference that requires exchangeability assumptions, but no other distribution assumptions. Like classical confidence intervals, our confidence intervals are not valid in a distribution-free setting. Our assumptions are neither strictly stronger nor weaker than distribution-free approaches: we do not require the exchangeability assumption (because the spatial locations are fixed for us), but we do assume a parametric (Gaussian) model for the noise. We see in our experiments that approaches that are known to be valid in the distribution-free setting, like the sandwich estimator, do not provide correct coverage in the setting we study.\n\n**Extensions beyond Gaussian noise.** An asymptotic analysis of our method under weaker noise settings (e.g. bounded variance) would be an interesting area for further work. If one allows the number of neighbors to tend to infinity, then asymptotic normality of our estimator follows by the classical central limit theorem together with the continuous mapping theorem. However, a detailed and complete asymptotic study of our method would require significant additional mathematical development.\n\n**Restrictiveness of Assumption 1.** The reviewer suggests that Assumption 1 requires that, for a fixed spatial location, the covariate will be always the same. In fact, we do not need to assume that covariates are fixed across time in our framework. In particular, as we detail next, there are at least three ways that time can be accounted for within our framework, and which way is most appropriate depends on the scientific question of interest:\n1. Time can be treated as another “spatial” dimension. The test points and train points both then have a temporal location as well as a spatial location. We might take this approach if we were interested in the association between hourly air temperature and hourly air pollution, including information for all hours from 2010–2020. In this case, we would need to assume that hourly average air pollution satisfies a Lipschitz assumption in both space and time.\n2. For each spatial location, each covariate and response could be aggregated over time (for example, via taking a mean or median at each spatial location, particularly if at many points in time the same spatial location is observed). The regression analysis can then be run on this temporally-aggregated data. We might take this approach if we were interested in the association between decadal average temperature and decadal average air pollution in major US cities, averaged over 2010–2020.\n3. Data from a single point in time can be analyzed. In this case, many covariates (e.g. temperature) are essentially fixed functions of space, since time has been fixed. We might take this approach if we were interested in the association between current temperature and current air pollution in major US cities at a particular fixed time\nWe think that each type of question above is of practical interest. We also agree with the reviewer that there exist questions of practical interest for which our method is not appropriate (e.g., analyzing demographic data with multiple individuals in a single household). We will add more discussion of the various ways to handle time within our framework, as well as limitations, in an updated manuscript.\n\nWe hope these clarifications address your questions.  Thank you again for your constructive feedback!\n\n**References:**\n\nGiordano, Ray, Schmidt-Hieber. On the inability of Gaussian process regression to optimally learn compositional functions. NeurIPS 2022.\n\nZhang, Yuan, Zhu. Byzantine-tolerant federated Gaussian process regression for streaming data. NeurIPS 2022.\n\nEl Harzli, Grau, Valle-Perez, Louis. Double-Descent Curves in Neural Networks: A New Perspective Using Gaussian Processes. AAAI 2024.\n\nJordan. \"Bayesian nonparametric learning: Expressive priors for intelligent systems.\" Heuristics, probability and causality: A tribute to Judea Pearl, 2010.\n\nCui, Loureiro, Krzakala, Zdeborova. Generalization Error Rates in Kernel Regression: The Crossover from the Noiseless to Noisy Regime. NeurIPS 2021.\n\nCherian, Gibbs, Candes. Large Language Model Validity via Enhanced Conformal Prediction methods. NeurIPS 2024.\n\nCorreia, Massoli, Louizos,  Behboodi. An Information Theoretic Perspective on Conformal Prediction. NeurIPS 2024."}, {"review_id": "pHvhRUOEJQ", "reviewer": "Reviewer_jtMd", "summary": "The paper concerns estimating associations between covariates and responses with the correct frequentist coverage of the confidence interval. With this aim, the authors proposed a new method for associations in spatial settings.", "questions": "please see the strength and weakness section.", "limitations": "yes", "rating": 4, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The authors described the ideas and demonstrated them using simulation studies. I have some concerns about the novelty, although I might overlook it.\nI struggled to understand the logic behind Assumptions 1 and 2, the existence of independent maps from spatial locations to the response and covariates, respectively. i.e., Two data-generating processes. It would be more natural to consider a parametric model with covariates which are mapped from locations to account for the association between response and spatial location. In the simulation studies, the authors generated data in this way. I think it is good to clarify this point.\n\nGood estimation of a map between location and covariate influences the parameter estimation; the authors assume that they are smooth, indicated by the Lipschitz constant. To me, the smoothness assumption of this map sounds overly simplified, and there is a bias-variance trade-off and more thorough discussion and studies about its impact on parameter estimation.\n\nAlthough the authors discussed variance parameter estimation in the paper, a confidence interval with the nominal coverage challenges, as the authors noted in the simulation study. Simple linear models are considered in the paper. Although the authors plan to extend this to GLM, the current version seems to present limited studies.", "was_revised": true, "final_justification": "Thanks to the authors for clarifying the aims of the paper and responding to my questions.\nNow I have a better understanding of the paper, and I've changed my score.", "rebuttal": "We thank the reviewer for the feedback and comments. Before responding to individual points, we recall that our goal in the present paper is to provide valid confidence intervals for the signed strength of the association between covariates and a response in regression.\n\n**Lack of novelty and focus on linear regression.** Given later comments in the review, we interpret the reviewer’s concern about lack of novelty to stem from our focus on the analysis of linear regression. If the reviewer has other novelty concerns, we ask them to please clarify those. \n\nThe novelty of our paper comes from recognizing and addressing biases introduced into estimation of associations in spatial settings; in particular, the biases due to model misspecification and distribution shift are ubiquitous in real-life problems. We provide the first framework for correcting these biases when reporting associations between covariates and a response.\n\nWe would like to push back on the idea that a paper focused on linear regression cannot be novel or relevant to the ML community precisely because it is on the topic of linear regression. We first note that, as we argue in our introduction, linear regression just happens to be the current best tool (even considering more complex, recently-developed methods) for the particular data analysis task of quantifying associations (a task that is extremely widespread across the sciences and social science). So our focus is not on linear regression for its own sake but instead because of its appropriateness for the task. Indeed, there remain many tasks of interest to ML for which linear regression is an appropriate tool. For this reason, numerous papers were published even just at last year’s NeurIPS on topics related to linear regression (e.g. Xie and Huo 2024, Lin et al, 2024, Jain et al. 2024, Zhu et al. 2024, Liu and Novikov 2024).\n\nThat all being said, we agree with the reviewer that extensions of our approach beyond linear regression are interesting. We hope that our present work will lay the groundwork for that future work, which we expect to require substantial additional effort. \n\n**Motivation for Assumptions 1 and 2.** (*“I struggled to understand the logic behind Assumptions 1 and 2, the existence of independent maps from spatial locations to the response and covariates, respectively. i.e., Two data-generating processes. It would be more natural to consider a parametric model with covariates which are mapped from locations to account for the association between response and spatial location. In the simulation studies, the authors generated data in this way. I think it is good to clarify this point.”*)\n\nWe understand the reviewer to be suggesting that our assumptions (particularly Assumption 2) exclude the possibility that the response is a function of the covariate. We emphasize that, while in Assumption 2 we formally state that the response is a function of the spatial location, the response can also be a function of the covariates. We do not state this potential dependence explicitly in Assumption 2 because of Assumption 1. Formally, suppose that $E[Y|X, S] = g(X, S)$. By Assumption 1, $X = \\chi(S)$, and so $E[Y|X,S] = g(\\chi(S), S)$. Define $f: \\mathcal{S} \\to \\mathbb{R}$ by $f(s) =  g(\\chi(S), S)$. Then $E[Y|X,S] = E[Y|X] = g(\\chi(S), S)$. In particular, the response can be a function of the covariates, but because the covariates are a fixed function of space, we can also write the response as a function of space alone. We think it is often easier to reason about smoothness in physical space, than in a (potentially high-dimensional) space of covariates. And so we impose our smoothness restrictions on the average response as a function of space alone, but this formulation does not mean the response cannot depend on the covariates as outlined above. We will clarify this point in our revision.\n\nWe believe the reviewer is suggesting to instead assume the data arises according to $Y = f(X) +$ noise, with a parametric choice of f. Choosing a parametric model for f would require agreeing on an appropriate model. The linear model is a natural choice but recovers the well-specified case, which is both well-studied and unrealistic. Other parametric choices might be considered arbitrary and limiting. As we describe in the paragraph above, we are essentially already making a similar assumption, but with a nonparametric choice of f, and we allow the latent function to depend on the spatial location in additional ways beyond just via the covariates.\n\n**Smoothness of maps (functions) 1.** (*“Good estimation of a map between location and covariate influences the parameter estimation; the authors assume that they are smooth, indicated by the Lipschitz constant.”*)\n\nWe understand the reviewer to be concerned with the problem of estimating a mapping between spatial locations and covariates. We would like to emphasize that, while we assume a map between spatial locations and covariates exists, we do not need to estimate it. Nor do we make any smoothness assumptions about it. We **do** make a smoothness assumption on the expected response.\n\n**Smoothness of maps (functions) 2.** (*“To me, the smoothness assumption of this map sounds overly simplified”*)\nWhile we reiterate that we do **not** require a Lipschitz assumption on the covariates as a function of the spatial locations, the reviewer might be concerned about the smoothness assumption we do make on the expected response. We chose to make this smoothness assumption because (1) it was simple enough to be physically interpretable and (2) there are real-data analyses where it is a reasonable assumption. For a concrete example of the second point, see our response to reviewer W9eu; there we reference a study that suggests concrete scales over which the annual average PM2.5 in California varies, and these could be used to select a Lipschitz constant of, for example 1.2 $\\mu g/m^3$. We emphasize that any machine learning method must require some assumption about how much a response is allowed to vary as you move in space — so that the method can perform inference at locations where we haven’t directly observed data; if we allow arbitrarily large changes in the response within an arbitrarily small spatial distance, even seemingly simple interpolations become impossible. While one could choose a more complicated assumption, we worry that alternative assumptions could lose physical interpretability, making them difficult to use in practice.\n\n**Bias-variance trade-off.** (*“there is a bias-variance trade-off and more thorough discussion and studies about its impact on parameter estimation.”*)\nWe assume the reviewer is concerned with the bias-variance tradeoff that arises from specifying a Lipschitz constant — i.e. that with our method for estimating the Lipschitz constant, a higher variance is estimated when a small Lipschitz constant is specified, and a smaller noise variance is estimated when a larger Lipschitz constant is specified. We illustrated the impact this trade-off has on parameter estimation in Figure 2 on simulated data, as well as on real-data in Figure 9 in the appendix. On the simulated data, we saw that confidence interval width is not always monotonically increasing with the Lipschitz constant, because of smaller estimates of the noise variance (Figure 2, middle and right). But when more extrapolation is required, larger Lipschitz constants lead to wider confidence intervals (Figure 2, left). In the real-data experiment, we also saw that for the dataset considered, larger Lipschitz constants led to wider confidence intervals. We further found that our method maintained nominal coverage across a wide-range of Lipschitz constants, suggesting that choosing a small Lipschitz constant in many cases still can correct bias. These observations support that attempting to account for the bias introduced by misspecification and distribution shift is preferable to ignoring it all together.\n\n**Extensions beyond linear models and limited studies.** (*“Although the authors plan to extend this to GLM, the current version seems to present limited studies.”*)\nWhile we agree the extension of our method to GLMs is of interest (indeed we are working on this direction), it requires non-trivial work beyond what is in this paper to ensure the confidence intervals are valid. In this paper, we developed a framework for addressing bias due to distribution shift and misspecification when estimating spatial associations using linear regression; our present work already required novel assumptions and technical arguments. To address the stated concern about limited studies, we note that we included several (3) simulation studies as well as two analyses on a real data set in the current paper. We feel that these illustrate the key benefits, considerations, and limitations of our method. If the reviewer has specific other studies in mind (or substantive points that would require additional studies to address), we would be grateful for suggestions.\n\n**References:**\n\nChow, Chen, Watson, Lowenthal, Magliano, Turkiewicz, Lehrman. PM2.5 Chemical Composition and Spatiotemporal Variation During the California Regional PM10/PM2.5 Air Quality Study. Journal of Geophysical Research. 2006.\n\nXie, Huo. High-dimensional (Group) Adversarial Training in Linear Regression. NeurIPS 2024.\n\nLin, Wu, Kakade, Bartlett, Lee. Scaling Laws in Linear Regression: Compute, Parameters, and Data. NeurIPS 2024.\n\nJain, Sen, Kong, Das, Orlitsky. Linear Regression using Heterogeneous Data Batches. NeurIPS 2024.\n\nZhu, Manseur, Ding, Liu, Xu, Wang. Truthful High Dimensional Sparse Linear Regression. NeurIPS 2024.\n\nChih-Hung Liu, Gleb Novikov. Robust Sparse Regression with Non-Isotropic Designs. NeurIPS 2024"}, {"review_id": "CKMoR7Epuz", "reviewer": "Reviewer_NCSm", "summary": "This paper proposes a new method for constructing valid confidence intervals for spatial associations under model misspecification and distribution shift. By assuming only that the response is Lipschitz-continuous in space, the method avoids strong modeling assumptions and achieves nominal coverage where existing methods fail. Experiments on simulated and real data show superior coverage and efficiency compared to standard approaches.", "questions": "1. The paper considers the setting where $Y$ is an $N$-long column vector. How would the proposed method extend to cases where $Y$ is high-dimensional (i.e., multivariate responses)?\n\n2. From a mathematical perspective, what are the fundamental reasons that make the proposed method better than Gaussian Processes (GP BCI), particularly under model misspecification or distribution shift?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**:\n- **Clarity and Accessibility**: The paper is very clearly written and easy to follow. Even for readers outside the immediate subfield, the key ideas and contributions are accessible.  \n- **Strong Motivation**: The motivation is compelling and well-grounded in real-world challenges of spatial inference under model misspecification and distribution shift.  \n- **Theoretical Rigor**: The analysis is mathematically solid, with well-justified assumptions and detailed, carefully presented proofs.\n\n**Weakness**:\n- **Title Meaning**: Why is the paper titled *\"Smooth Sailing\"*? How does it relate to the method and its content?", "was_revised": true, "final_justification": "I would like to give score 4 for this paper.", "rebuttal": "We thank the reviewer for their insightful feedback and thoughtful questions. We are excited to hear that the reviewer found our paper clearly written, strongly motivated, and mathematically rigorous. Below we address each question in turn.\n\n**Why “Smooth Sailing”?** The “Smooth Sailing” part of our title is meant to be a pun, playing on two meanings of “smoothness.” (1) “Smooth sailing” is an English idiom that denotes “easy progress without impediment or difficulty” (Merriam-Webster Dictionary). (2) Under only a Lipschitz-continuity assumption (a “smoothness” assumption), our procedure guarantees valid confidence intervals — even in the presence of model misspecification and covariate shift. So our goal was to convey that, under a smoothness (Lipschitz) assumption, inference using our procedure proceeds smoothly (easily).\n\n**Extension to multivariate responses.** While our paper focuses on a scalar response $Y(s)\\in \\mathrm{R}$, very similar machinery can be adapted to confidence intervals for each coordinate of $Y(s)\\in \\mathrm{R}^D$. Namely, we can treat each component $Y^{(d)}(s)$ as a separate univariate problem, and apply our Lipschitz‐based bias bound and variance calibration to each coordinate. If one desires simultaneous coverage over all $D$ outputs, a straightforward Bonferroni correction or any family‐wise error control can be used. We haven’t tested this proposal in practice, as our focus was on scalar responses, but we think it would be an interesting future task. We will add a short paragraph in Section 3 outlining this extension.\n\n**Why not Gaussian process credible intervals?** Gaussian process–based Bayesian credible intervals are popular in spatial statistics. As is always the case with Bayesian approaches, the credible intervals represent subjective beliefs. In the presence of misspecification (which we focus on in this work) posterior credible intervals can under or overestimate uncertainties about the parameter (Müller, 2013; Walker, 2013). This is even more the case when there is also a distribution shift, as the prior is likely to have a large impact on the posterior when we are extrapolating. In contrast, our method (i) explicitly bounds worst‐case bias via a 1-Wasserstein distance and (ii) adds this bias bound to the confidence interval, yielding guaranteed frequentist coverage under arbitrary covariate shift. In contrast, our method (i) explicitly bounds worst‐case bias via a 1-Wasserstein distance and (ii) adds this bias bound to the confidence interval, yielding guaranteed frequentist coverage under arbitrary covariate shift.\n\nWe hope these clarifications address your questions.  Thank you again for your constructive feedback!\n\n**References:**\n\nMüller. “Risk of Bayesian inference in misspecified models, and the sandwich covariance matrix.” Econometrica.  2013. \n\nWalker. “Bayesian inference with misspecified models.” Journal of Statistical Planning & Inference. 2013."}]

Author rebuttal:

Meta-review: We thank the authors for their submission. 

The authors propose a confidence interval estimator in a spatial regression setting when a linear model may be misspecified.  They show that given a smoothness assumption on the response, their approach can achieve nominal coverage. The method is evaluated on synthetic and real data, favorably comparing to alternative approaches. 

Reviews were uniformly positive.  Reviewers agreed that this work is well motivated, well-organized and (mostly) clearly written, and technically solid.  The authors clearly state the limitations of their approach (i.e., the strength of the smoothness assumption).  The discussion with reviewer W9eu shed light on the practical implications of this limitation, with the remedy to incorporate domain knowledge about the response function.  The discussion with reviewer pd1Y led to multiple clarifications in the text, as well as additional discussion about the specific technical contribution of this work and the existing literature.
Decision: Accept (poster)