{'fraction_non_alphanumeric': 0.04624967725277562, 'fraction_numerical': 0.014297702039762458, 'mean_word_length': 4.327458734525447, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Valid model-free spatial prediction

Predicting the response at an unobserved location is a fundamental problem in spatial statistics. Given the difficulty in modeling spatial dependence, especially in non-stationary cases, model-based prediction intervals are at risk of misspecification bias that can negatively affect their validity. Here we present a new approach for model-free nonparametric spatial prediction based on the conformal prediction machinery. Our key observation is that spatial data can be treated as exactly or approximately exchangeable in a wide range of settings. In particular, under an infill asymptotic regime, we prove that the response values are, in a certain sense, locally approximately exchangeable for a broad class of spatial processes, and we develop a local spatial conformal prediction algorithm that yields valid prediction intervals without strong model assumptions like stationarity. Numerical examples with both real and simulated data confirm that the proposed conformal prediction intervals are valid and generally more efficient than existing model-based procedures for large datasets across a range of non-stationary and non-Gaussian settings.


Introduction
Providing valid predictions of the response at an unobserved location is a fundamental problem in spatial statistics. For example, epidemiologists may wish to extrapolate air pollution concentrations from a network of stationary monitors to the residential locations of the study participants. There are a number of challenges one faces in carrying out valid prediction at a new spatial location, but one of the most pressing is that existing methods are model-based, so the reliability of the predictions depends crucially on the soundness of the posited model. For example, prediction intervals based on Kriging-see Cressie (1992) and Section 2.1-often rely on normality, and stationarity is often assumed to facilitate estimating the spatial covariance function required for Kriging. It is now common to perform geostatistical analysis for massive datasets collected over a vast and diverse spatial domain (e.g., Heaton et al., 2019). For complex processes observed over a large domain, the normality and stationarity assumptions can be questionable. Failing to account for nonstationarity can affect prediction accuracy, but typically has a larger effect on uncertainty quantification such as prediction intervals (Fuglstad et al., 2015). While there are now many methods available for dealing with nonstationary (see Risser, 2016, for a recent review) and non-Gaussianity (e.g., Gelfand et al., 2005;Duan et al., 2007;Reich and Fuentes, 2007;Rodriguez and Dunson, 2011), these typically involve heavy computations. This exacerbates the already imposing computational challenges posed by massive datasets. Further, fitting the entire stochastic process may be unnecessary if only prediction intervals are desired. Nonparametric machinelearning methods can be used for prediction (Kim et al., 2016a,b;Lim et al., 2017;Tai et al., 2017;Hengl et al., 2018;Franchi et al., 2018;Wang et al., 2019;Li et al., 2020), but these methods typically focus on uncertainty estimation. In this paper we propose a method with provably valid prediction intervals-exact in some cases, asymptotically approximate in others-for the response at a single location without requiring specification of a statistical model, and hence not inheriting the risk of model misspecification bias.
In recent years, the use of machine learning techniques in statistics has become increasingly more common. While there are numerous examples of this phenomenon, the one most relevant here is conformal prediction. This method originated in Vovk et al. (2005) and the references therein (see, also, Shafer and Vovk, 2008), but has appeared frequently in the recent statistics literature (e.g., Lei and Wasserman, 2014;Lei et al., 2018;Guan, 2019;Romano et al., 2019;Tibshirani et al., 2019). What makes this method especially attractive is that it provides provably valid prediction intervals without specification of a statistical model. More precisely, the conformal prediction intervals achieve the nominal frequentist prediction coverage probability, uniformly over all data distributions; see Section 2.2. The crucial assumption behind the validity of conformal prediction is that the data are exchangeable.
Whether it is reasonable to assume exchangeability in a spatial application depends on how the data are sampled. On the one hand, if the locations are randomly sampled in the spatial domain, then exchangeability holds automatically; see Lemma 1. In such cases, standard conformal prediction can be used basically off the shelf. On the other hand, if the locations are fixed in the spatial domain, then exchangeability does not hold in general. We show, however, that for a wide range of spatial processes, the response variables at tightly concentrated locations are approximately exchangeable; see Theorem 1. Therefore, a version of the basic conformal prediction method applied to these tightly concentrated observations ought to be approximately valid.
Using this insight about the connection between exchangeability and the sampling design, we propose two related spatial conformal prediction methods. The first, a so-called global spatial conformal prediction (GSCP) method, described in Section 3, is designed specifically for cases where the spatial locations are sampled at random. In particular, this global method produces a prediction interval which is marginally valid, i.e., valid on average with respect to the distribution of the target location at which prediction is desired; asymptotic efficiency of this global method is also investigated. The second, a local spatial conformal prediction (LSCP) method, described in Section 4, is designed specifically for the case when the spatial locations are fixed. Since our goal is to proceed without strong assumptions about the spatial dependence structure, it is only possible to establish approximate or local exchangeability. Therefore, the proposed local spatial conformal prediction method can only provide approximately valid predictions; see Theorem 2. But our goal in the fixed-location case resembles the "conditional validity" target in the conformal prediction literature (e.g., Barber et al., 2019;Chernozukov et al., 2021) so, given the impossibility theorems in the latter context, approximate validity is all that can be expected.
For both the global and local formulations, our proposed method is computationally feasible for large datasets and model-free in the sense that its validity does not depend on a correctly-specified model. In Sections 5 and 6, we show using real and simulated data that the proposed methods outperform both standard global Kriging and local approximate Gaussian process regression (laGP; Gramacy and Apley, 2015) for non-stationary and non-Gaussian data. In addition to be useful for spatial applications, it is also an advancement in conformal prediction to the case of dependent data, and establishes the conditions on the spatial sampling design and data-generating mechanism that ensure (approximate) validity of the conformal prediction intervals.
The remainder of this paper is organized as follows. Section 2 reviews spatial and conformal prediction. Sections 3 and 4 introduce the proposed methods, which are examined using simulations in Section 5 and a real data analysis in Section 6. Additional numerical and theoretical results, along with all proofs, are given in the Supplemental Materials.

Spatial prediction methods
Let Y i ∈ R and X i ∈ R d , with d ≥ 1, be the observable pairs at spatial location s i ∈ D ⊆ R 2 . Note that X i can include covariates that are deterministic functions of the spatial location, such as elevation, genuinely stochastic covariates like wind speed, or even non-spatial covariates such as the smoking status of the resident at location. Write the data points as triples Z i = (s i , X i , Y i ), for i = 1, 2, . . .. We assume only a single observation is made at each location and thus often adopt the notation Y i = Y (s i ) and X i = X(s i ).
Geostatistical analysis often assumes that the data follow a Gaussian process model, Y i = X i β + θ i + ε i , for i = 1, . . . , n, where β is the vector of regression coefficients, ε 1 , . . . , ε n are independent Normal(0, τ 2 ) errors, θ i = θ(s i ), and θ is a mean-zero Gaussian process with isotropic covariance function C(θ i , θ j ) = σ 2 ρ(d ij ), a function of the distance d ij between locations s i and s j . A common example is the Matérn correlation function ρ(d; φ, κ) parameterized by correlation range φ and smoothness parameter κ. Denote the spatial covariance parameters as Θ = {σ 2 , τ 2 , φ, κ}. The main assumptions of this model are that the data are Gaussian and the covariance function is stationary and isotropic, i.e., it is a function only of the distance between spatial locations and is thus the same across the spatial domain.
Consider data Z n+1 = (Z 1 , . . . , Z n , Z n+1 ). In our applications, Z n will be observed and (s n+1 , X n+1 ) will be given, and the goal is to predict the corresponding Y n+1 . However, the ordering of the data is irrelevant so one can imagine different orderings that correspond to data point i in the last position, where i = 1, . . . , n, n + 1. That is, imagine we have the observed data z n+1 (i) = {z 1 , . . . , z n+1 } \ {z i }, along with (s i , x i ) and parameter estimates β and Θ; then the predictive distribution of Y i is normal with meanμ n+1,i (s i , x i ) and variance σ 2 n+1,i (s i , x i ), where both the mean and variance depend on Θ and the configuration of the spatial locations; see the Supplementary Material for the specific expressions. The standardized residuals are and the corresponding 100 where q α is the upper α/2 quantile of a standard normal distribution.

Conformal prediction
Here we take a step back and review conformal prediction for non-spatial problems; for a detailed treatment, see Vovk et al. (2005) and Shafer and Vovk (2008). Suppose we have a data sequence Z 1 , . . . , Z n , Z n+1 , . . ., assumed to be exchangeable with joint distribution P, that is, Z 1 , Z 2 , Z 3 , . . . and Z ξ(1) , Z ξ(2) , Z ξ(3) , . . . have the same joint distribution for any permutation ξ defined on the positive integers. This data may be response-only, i.e., Z i = Y i , or may be response-covariate pairs, i.e., Z i = (X i , Y i ); we will focus on the latter more general case. No assumptions about P are made here, beyond that it is exchangeable. We observe Z n = z n , and the goal is to predict Y n+1 at a new value X n+1 of the covariate. More specifically, we seek a procedure that returns, for any α ∈ (0, 1), a prediction interval Γ α (Z n ; X n+1 ) that is valid in the sense that where P n+1 is the distribution of (Z 1 , . . . , Z n , Z n+1 ) under P. That we require the inequality (2) to hold for all exchangeable distributions P rules out the use of model-based procedures, such as likelihood or Bayesian methods. The original conformal prediction method proceeds as follows. Define a non-conformity measure ∆(B, z), a function that takes two arguments: the first is a "bag" B that consists of a finite collection of data points; the second is a single data point z. Then ∆(B, z) measures how closely z represents the data points in bag B. For example, if π is a prediction rule and d is some measure of distance, then we might take ∆(B, z) = d(π B , z), the distance between z and the value π B returned by the prediction rule π applied to B. The choice of ∆ depends on the context, though often there is a natural choice. Throughout this paper we will assume ∆ is symmetric in its first argument, so that shuffling the data in bag B does not change the value of ∆(B, z).
Given the non-conformity measure ∆, the next step is to appropriately transform the data via ∆. Specifically, augment the observed data z n = (z 1 , . . . , z n ) with a provisional value z n+1 of Z n+1 = (X n+1 , Y n+1 ); this z n+1 value is generic and free to vary. Define . . , n, n + 1. Note that δ n+1 is special because it compares the actual observed data with this provisional value of the unobserved future observation. Next, compute the plausibility (e.g., Cella and Martin, 2022) of z n+1 as a value for Z n+1 according to the formula where 1{A} denotes the indicator of event A. Note that this process can be carried out for any provisional value z n+1 , so the result is actually a mappingỹ → p(ỹ | z n ,x), for a givenx, which we will refer to as the plausibility contour returned by the conformal algorithm. This function can be plotted to visualize the uncertainty about Y n+1 based on the givenx, data z n , the choice of non-conformity measure, etc. Moreover, a prediction set Γ α (z n ;x) can be obtained as 3 Global spatial conformal prediction 3.1 GSCP algorithm The first approach we consider is a direct application of the original conformal algorithm to spatial prediction but with spatial dependence encoded in the non-conformity measure. In contrast to the local algorithm presented in Section 4, this method equally weights the non-conformity across all spatial locations in the plausibility contour evaluation. Therefore, we refer to this as global spatial conformal prediction, or GSCP for short. From Section 2.1, recall that z i = (s i , x i , y i ) and z n+1 = {z 1 , . . . , z n+1 }; also, z n+1 (i) denotes z n+1 \ {z i }, the full data set with z i excluded. Now define the non-conformity measure for the GSCP algorithm as whereμ n+1,i (s i , x i ) andσ n+1,i (s i , x i ) are, respectively, the mean response and standard error estimates at spatial location s i with covariate x i based on data in z n+1 (i) . Then we define a plausibility contour exactly like in Section 2.2, with obvious notational changes. That is, for provisional values (s n+1 , x n+1 , y n+1 ) of (S n+1 , X n+1 , Y n+1 ), we have The corresponding 100(1 − α)% prediction interval for Y n+1 , denoted by Γ α (z n ; s n+1 , x n+1 ), is just an upper level set of the plausibility contour, consisting of all those provisional y n+1 values with plausibility exceeding α, analogous to (4). Any reasonable choice of (μ,σ) estimates can serve the purpose here, including inverse distance weighting predictions (Henley, 2012), deep learning predictions (Franchi et al., 2018), Kriging predictions, etc. In our numerical results presented below, we use the Kriging estimates as defined in Section 2.1, so that δ i is a standardized Kriging residual, |e i |, from (1). Conformal prediction is invariant to monotone transformations of its δ i 's, and we found that similar results are obtained with other related measures, such as unstandardized Kriging residuals. A particular advantage of our recommended choice of δ i 's is that we can quickly compute the plausibility contour and prediction interval by exploiting the inherent quadratic structure of the Kriging-based non-conformity measure; see the Supplementary Materials. Moreover, note that validity of the GSCP-based prediction intervals does not require the Gaussian model associated with the Kriging method be correctly specified, nor does it depend on our choice of the δ i 's.

Theoretical validity of GSCP
Given the importance of exchangeability to the validity of conformal prediction and the fact that the spatial dependence generally is incompatible with exchangeability, we might have some concerns about the validity of GSCP. However, there are practically relevant cases in which exchangeability does hold, in particular, when the spatial locations are sampled independently and identically distributed (iid). The following elementary lemma explains this.
Since randomly sampled spatial locations makes the data exchangeable, a validity property for GSCP follows immediately from the general theory in, e.g., Shafer and Vovk (2008).
The upper bound in (8), which follows from the same arguments as in Lei et al. (2018), implies that the GSCP method is not only valid but also efficient in the sense that the coverage probability is not too much larger than the nominal level. That is, the coverage condition is not being achieved simply giving excessively wide intervals. Some further details on the efficiency of the global spatial conformal prediction procedure are investigated in the Supplementary Materials. Note, also, that Theorem 1 makes no assumptions about the distribution of (X, Y ), so it surely covers non-Gaussian and non-stationary processes.
Theorem 1 gives a marginal validity result in the sense that it accurately predicts the response Y (S n+1 ) at X(S n+1 ), for a randomly sampled spatial location S n+1 . However, it does not ensure conditional validity, i.e., the case where S n+1 = s with s being a fixed spatial location. There are negative results in the literature (e.g., Lei and Wasserman, 2014) which state that strong conditional validity-for all P and almost all targets s -is impossible with conformal prediction. Considerable effort has been expended recently trying to achieve "approximate" conditional validity in some sense; see, e.g., Lei and Wasserman (2014), Tibshirani et al. (2019), Barber et al. 6 (2019), and Chernozukov et al. (2021). Remarkably, there is at least one scenario in which a strong conditional validity result can be achieved in our context. In particular, as we show in the Appendix, the GSCP-based intervals are both marginally and conditionally valid for the special case of an isotropic process sampled uniformly on a sphere. Admittedly, these are rather strong conditions, so one would hope for (approximately) valid prediction under much less. In Section 4, we show that asymptotically valid prediction intervals at a fixed location can be obtained under only mild conditions on the sampling scheme and the unknown process.
4 Local spatial conformal prediction

LSCP algorithm
For valid prediction at a fixed location s , we propose a local spatial conformal prediction (LSCP) approach that is based on only those data points in the neighborhood of s . Fix an integer m > 0 and select a neighborhood around s that contains m many locations, s i j , for j = 1, . . . , m. Note that {s i j : j = 1, . . . , m} is a subset of the full set of spatial locations s 1 , . . . , s n . Without structural assumptions about the response process, such as stationarity, the data at locations far from s are not obviously relevant to prediction at s , so removing-or down-weighting (Section 4.3)-them from the local analysis is reasonable. Plus, in applications where the infill asymptotic regime is appropriate, there are many observations nearby s , so m could be taken to be large.
From here, we can proceed very much like in Section 3. For notational simplicity, assume that indices i = 1, . . . , m correspond to those m spatial locations closest to s . Now let Z i = (s i , X i , Y i ), for i = 1, . . . , m, denote the observations at these m closest locations to s . With a slight abuse of that notation, set s m+1 = s and (x m+1 , y m+1 ) as the provisional values of X and Y at s . Then define the non-conformity scores exactly as before: (i) , z i ), i = 1, . . . , m + 1. With this, we can readily obtain the plausibility contour function: Specific details are presented in Algorithm 2 in the Supplementary Materials. The output of this algorithm is a 100(1 − α)% prediction interval for Y m+1 = Y (s ), depending on Z m and the observed X m+1 = X(s ), which we denote by Γ α s (Z m ; X m+1 ).

Theoretical validity of LSCP
Our theoretical results hinge on a definition of local exchangeability. Let D ⊂ R 2 be a compact spatial domain, e.g., [0, 1] 2 . For a generic R d -valued stochastic process T defined on D, with d ≥ 1, define the localized version of T , relative to a location s ∈ D, as 7 indexed by the unit disk U and the radius r > 0. Now suppose that T can be decomposed as where L and E are independent R d -valued stochastic process, L is a continuous spatial process, E is a non-spatial process, and ψ is a deterministic, continuous, R d -valued function. More specifically, suppose that L and E, respectively, satisfy the following conditions: • L is L 2 -continuous at s in the sense that its localized version L r satisfies E L r (u) − L 0 (u) 2 → 0 as r → 0 for any u ∈ U; • E is locally iid at s , that is, its localized version E r converges in distribution to an iid process as r → 0.
This formulation is too abstract to be useful, but formulating a general result here is appropriate. Appendix A.2 describes several common spatial models that satisfy (11), and generalizes the above formulation to the case where the covariates are also considered stochastic processes. These assumptions yield a certain kind of local exchangeability which will be used below to show the LSCP algorithm achieves a desired validity property. We first establish this local exchangeability result, which may be of independent interest.
Proposition 1. Suppose that T can be decomposed as in (11), where L is L 2 -continuous at s , E is locally iid at s , and L and E are independent. Then the localized process T r in (10) converges in distribution as r → 0, and the limit is an exchangeable process in the sense that its finitedimensional distributions are exchangeable.
Using Proposition 1, we can establish the (asymptotically approximate) theoretical validity of the LSCP method. To set the scene, those m spatial locations closest to s fall in a neighborhood of some radius r. As is common in the spatial statistics literature (e.g., Stein, 1990;Cressie, 1992), we adopt an infill asymptotic regime in which the region D remains fixed while the number of observations n goes to infinity, hence filling the space. The relevant point for our analysis is that under this regime the number of observations made in any neighborhood of s ∈ D will go to infinity. Such a regime is natural-and necessary-in cases without structural assumptions about Y , e.g., stationarity, where it is simply not possible to learn the local features of a process at s if data are not concentrated in a neighborhood around s . Under the infill asymptotic framework , if m is fixed and the number of locations n is increasing to fill the bounded space D, then the radius of the neighborhood in which those m points fall is vanishing. For example, if the spatial locations are (roughly) uniformly distributed in D, then the number of points in a neighborhood of radius r would be proportional to nr 2 ; setting this equal to m gives r = r n = (m/n) 1/2 → 0 as n → ∞.
It follows from Proposition 1 that the joint distribution of the response Y at these m-many spatial locations around s (corresponding to m-many vectors in U) would be approximately exchangeable and, consequently, a conformal prediction algorithm that creates non-conformity scores using only these m observations would be valid for predicting Y (s ).
Theorem 2. Consider an infill asymptotic regime with n spatial locations in the bounded domain D, with n → ∞. Fix an integer m > 0 and let r = r n → 0 be such that the m closest locations to s fall in a neighborhood of radius r. Under the assumptions of Proposition 1, the non-conformity measure ∆ is a continuous function of its inputs, and if the limiting distribution in Proposition 1 is continuous, then the LSCP prediction intervals are asymptotically valid at s in the sense that where P m+1 is the joint distribution of Z 1 , . . . , Z m , Z m+1 at the m spatial locations and at s .

The smoothed LSCP algorithm
Theorem 2 implies that the local spatial conformal prediction with m nearest neighbors is approximately valid under the infill asymptotic regime. However, in practice completely disregarding the contribution of the observations outside the m nearest neighbors may be unsatisfactory, so we propose a smoothed version of the LSCP algorithm (sLSCP).
The GSCP algorithm weights all n + 1 non-conformity measures δ i equally in the plausibility contour computation in (3), but this is questionable for non-stationary processes with stochastic properties that vary throughout the spatial domain. To allow for non-stationarity, the sLSCP algorithm weights the non-conformity measures δ i by how far the corresponding s i is from the prediction location. Let f be a non-increasing function, and define weights , and the proportionality constant ensures that n+1 i=1 w i = 1. Different f functions can be applied, but we recommend the normalized Gaussian kernel function with bandwidth η, , i = 1, . . . , n + 1.
Note that, if η → ∞, then w i → (n + 1) −1 for each i, which corresponds to the GSCP algorithm. Finally, with these new weights, the plausibility contour at a provisional value (s n+1 , x n+1 , y n+1 ) of Y (s ) is given by As before, we recommend the Kriging-based strategy with δ i the standardized residual in (1).
Since we are interested only in the local structure of Y , it is natural that locations far from s have negligible weight, as in (12). But including all n observations requires some non-trivial calculations, e.g., inverting a large n × n covariance matrix. Therefore, to avoid cumbersome and ultimately irrelevant computation, we recommend using only the M n observations closest to Algorithm 1: Smoothed local spatial conformal prediction (sLSCP). Input: observations z i = (s i , x i , y i ), i = 1, . . . , n; predict location s ; non-conformity measure ∆; significance level α; and a fine grid of candidate response values Tuning parameter: weight parameter η ∈ (0, ∞); number of neighbors to consider M ≤ n Output: (12); 5 for y M +1 in the specified grid do s for both the Kriging predictions that determine δ i and in the plausibility scores in (13). The resulting method is both locally adaptive and computationally efficient even for large data sets.
The tuning parameter η can be selected using cross validation, as illustrated in Sections 5 and 6. The value of M is determined by the bandwidth η so that all observations with substantial w i are included, as are observations that are required for the Kriging prediction of these observations. Typically the number of nearby observations required to approximate the Kriging prediction is a small subset of the total number of observations (Stein, 2002). As a rule of thumb, M could be selected to roughly include all observations within 2η + r radius of s , where 2η captures observations with substantial weights, and r is selected so that all the M observations include the nearest 15 neighbors of the observation within 2η of s . We summarize the details of Algorithm sLSCP in Algorithm 1. For simplicity, we use sLSCP and LSCP indistinguishably.
It is important for our proposed methods to be computationally feasible for large datasets. Conformal prediction itself is relatively expensive since it requires fitting the underlying model once for each held-out data point being predicted. In particular, the Kriging residuals and the associated non-conformity score computations require us to computeμ n+1,i (s i , x i ) andσ n+1,i (s i , x i ) for each i = 1, . . . , n, which involves n many evaluations of ( β, Θ). To overcome this computational bottleneck, various adjustments have been considered in the literature. One is the split conformal prediction strategy-also called inductive conformal prediction in Vovk et al. (2005)-which is common; see, e.g., Lei et al. (2018, Sec. 2.2). The idea is to split the data into two parts: one for fitting the underlying model and the other for running conformal prediction with the fitted model from the first part fixed. The theoretical validity of split conformal prediction is now well-known, e.g., Section 3 of the Supplementary Materials. Alternatively, as is common in parametric Kriging, one could use the entire dataset to estimate ( β, Θ) and then use the entire dataset again for prediction with the parameter estimates plugged in as if they were the "true values." Given that the number of parameters in the working spatial model is relatively small, both approaches should perform well for moderate to large n. The simulation results presented in the Supplemental Materials suggest that this plug-in conformal is more efficient than split conformal in terms of width of the corresponding prediction intervals (or, more precisely, in terms of the interval score as defined in Section 5.2), so the numerical results in Sections 5 and 6 below are based on plug-in versions of the proposed GSCP and LSCP algorithms.

Data generation
We consider one mean-zero Gaussian stationary process (Scenario 1) and seven non-Gaussian and/or nonstationary data-generating scenarios (Scenarios 2-8). Data are generated based on transformations of a latent Gaussian process Z(s) and a white noise process E(s) with standard normal distribution, where s = (s x , s y ) ∈ [0, 1] 2 . The mean-zero stationary Gaussian process process Z(s) has a Matérn covariance function with variance σ 2 = 3, range φ = 0.1, and smoothness κ = 0.7. Data are sampled on the N ×N grid of n = N 2 points in the unit square, s ∈ {N −1 , 2N −1 , . . . , 1} 2 , with N = 20 or N = 40. The scenarios are: Scenario 1 is Gaussian and stationary, Scenarios 2-4 are stationary but non-Gaussian, and Scenarios 5-8 are nonstationary either in the spatial variance (Scenarios 5 and 6), error variance (Scenario 7), or mean (Scenario 8). Scenario 3 generates skewed data to assess the method's performance when the symmetry of the base Kriging model is violated.

Prediction methods and metrics
For each dataset we apply the global and local (with η = 0.1) conformal spatial prediction algorithms. For the parametric Kriging method and the initial Kriging predictions of our proposed conformal prediction, we estimate the spatial covariance parameters using empirical variogram methods (Cressie, 1992). The empirical variograms are calculated using the variog function in the R package geoR, and the covariance parameters are chosen to minimize the weighted (by number of observations) squared error between the empirical and model-based variograms. We compare the proposed conformal prediction methods with standard global Kriging prediction and the local Kriging (laGP) method of Gramacy and Apley (2015) that dynamically defines the support of a Gaussian process predictor based on a local subset of the data. For laGP, we use the function provided by the laGP package in R the local sequential design scheme starting from 6 points to 50 points through an empirical Bayes mean-square prediction error criterion.
Methods are trained using a completely random set of 90% of the observations and tested on the remaining 10%. Each scenario is repeated 100 times, and performance is evaluated using average coverage of (1 − α)100% prediction intervals, average interval width, and average interval score (Gneiting and Raftery, 2007), defined as is the 100(1 − α)% prediction interval, y n contains the observations y 1 , . . . , y n , and z + = z ∨ 0 denotes the "positive part." A smaller interval score is desirable as this rewards both high coverage and narrow intervals. We use α = 0.1 in this simulation study.

Results
We present results averaged over data sets and all spatial locations in Table 1. For the nonstationary scenarios varying across s x (Scenarios 5-7), we present the results by the first spatial coordinate (s x ) averaged over the data sets and the second coordinate (s y ) in Figure 1, e.g., the value of coverage plotted at s x = N −1 is the average of the coverage over the N points of the form In Scenario 1, the Gaussian and stationary process, the performance of GSCP, LSCP, and Kriging are comparable (Table 1). Kriging performs well in this case since the data generating mechanism aligns with its underlying assumption, but the conformal methods are competitive with the parametric model in terms of both coverage and interval width. In Scenarios 2, 3, and 4, the non-Gaussian but stationary processes, GSCP, LSCP, and Kriging perform more or less the same in terms of interval score and outperform laGP. However, the coverage of the conformal methods, especially the GSCP algorithm, is closer to the nominal level than Kriging and the Kriging intervals are generally wider than the conformal intervals. Figure 1 shows the results for nonstationary Scenarios 5 and 6 when N = 40. LSCP performs the best among the four methods for these nonstationary scenarios. For Scenario 5, the process is Gaussian to the west and non-Gaussian to the east. The global prediction methods GSCP and Table 1: Performance comparison for simulation scenarios ("Scen") without a covariate. The metrics are the empirical coverage of 90% prediction intervals ("Cov90"), the width of prediction intervals ("Width") and the interval score ("IntScore"), each averaged over location and dataset. The methods are global (GSCP) and local (LSCP) conformal prediction, stationary and Gaussian Kriging ("Kriging") and local approximate Gaussian process ("laGP") regression.  Table 2: Performance comparison for the canopy height data. The metrics are the width, coverage ("Cov90") and interval score ("IntScore") of 90% prediction intervals, each averaged over 10,000 randomly chosen test locations. The methods are local conformal prediction (LSCP), stationary and Gaussian Kriging ("Kriging") and local approximate Gaussian process ("laGP") regression. Kriging generate prediction intervals with similar width for all s x (ignoring edge effects), while LSCP and laGP provide wider intervals on the east (s x near 1) than the west (s x near 0). LSCP has coverage around 90% for all s x , and the lowest interval scores, especially in the east where the process is more non-Gaussian. Similarly, in Scenario 6, the correlation is stronger in the east than the west, and the LSCP performs the best by providing adaptive prediction interval width and valid coverage across space.
We also conducted a simulation study when spatial locations are sampled uniformly on [0, 1] 2 . The performance is very similar to that when locations are fixed at equally-spaced grid points, so we only show the latter in the paper. Additional results for the scenarios with covariates (thus a comparison with universal Kriging) and a sensitivity analysis confirming our method's robustness to the estimates of the spatial covariance parameters are included in the Supplemental Materials.

Real data analysis
This section demonstrates the performance of conformal prediction method using the canopy height data in Figure 2a. The data were originally presented in Cook et al. (2013) and were analyzed using a nearest-neighbor Gaussian process model in Datta et al. (2016). The data are available in the R package spNNGP (Finley et al., 2017). There are n = 1, 723, 137 observations and clear nonstationarity and non-normality. For example, there are several heterogeneous areas with small height canopies around the location with longitude and latitude being 729,000 and 470,000, respectively.
We compare methods using 90% prediction intervals for 10,000 test locations chosen randomly from the full dataset. Since the data clearly exhibit nonstationarity we do not apply GSCP. We select the kernel function and bandwidth parameter using cross-validation over the validation locations. The average interval score is consistently smaller for the Gaussian kernel (ranging between 4.3 and 4.4 by η) than the uniform kernel (ranging between 4.9 and 5.2 by η) and minimized by the Gaussian kernel with η = 6 × 10 −4 . Table 2 compares the performance of LSCP on the 10,000 test locations with Kriging and laGP. LSCP outperforms the other methods as the empirical coverage of LSCP is the closet to the desired 90% and the LSCP minimizes the interval score. Figure 2 plots the interval widths for each method. Unlike Kriging, the LSCP and laGP interval widths are locally adaptive with wider intervals in heterogeneous areas and more narrow intervals  in homogeneous areas. Comparing LSCP and laGP, LSCP generally provides narrower intervals than laGP, which means the proposed method is more efficient than laGP. In addition, the locations of the observations that fall outside the prediction intervals are uniformly distributed for LSCP and laGP, but clustered in the heterogeneous areas for Kriging. In short, the proposed spatial conformal prediction algorithm shows its superiority in this real data analysis.

Discussion
In this paper, we proposed a spatial conformal prediction algorithm to provide valid, robust, and model-free prediction intervals by combining spatial methods and the classical conformal prediction. We provided both global and local versions to accommodate different stationarity cases and sampling designs. We proved their validity under various sampling designs and data-generating mechanisms. To the authors' knowledge, this work is among the first in making the classical conformal algorithm work for non-exchangeable data. Our simulation studies and real data analyses demonstrate the advantage of the proposed spatial conformal prediction algorithms. We also developed an R package entitled scp (https://github.com/mhuiying/scp) to compute the plausibility contours and generate spatial prediction intervals using either Kriging residual or any other user-defined non-conformity measure.
An attractive feature of the proposed algorithms is that they are model-free in the sense that their theoretical validity does not depend on correct specification of a model. In our implementations we use the squared residuals from a simple parametric model to define the non-conformity terms. However, as anticipated by our theoretical results, our simulation study shows that the methods work well even if the parametric family or the mean and covariance functions are misspecified, and that the results are insensitive to inaccurate estimation of the parameters in the parametric model. This robustness allows the methods to be applied broadly and with confidence.
The local conformal prediction method relies on a dense grid of points around each prediction location and thus a large dataset. Computation time often prohibits application of spatial methods to large datasets. Fortunately, we are able to apply our method to large datasets by exploiting local algorithms and an explicit formula for the plausibility contour using the Kriging-based nonconformity score. Similar derivations are needed for prediction procedures other than Kriging in order to maintain computational efficiency.
Future research directions include extending the work to spatial processes with discrete observations, e.g., when Y (s) is binary or a count. Generalized spatial linear models are compatible with our current framework (see the Appendix), but continuity in the distribution function is required. Therefore, further studies would be required to establish the validity of conformal prediction for discrete data. Another limitation of the proposed algorithms is that they only produce intervals for a single location. Generalizing the algorithms to produce joint intervals for multiple locations would be useful in some applications. One option is to use a Bonferroni correction; of course, this may be inefficient for many simultaneous predictions, but greater efficiency would require model assumptions to link the multiple locations and facilitate information sharing. It would also be of interest to extend the proposed spatial conformal prediction methods to spatiotemporal data, perhaps building on recent work for time series data (Xu and Xie, 2020;Zaffran et al., 2022).

A.1: Conditional validity on a sphere
An obstacle that prevents a conditional validity result in the existing literature is an "edge effect." That is, conditional validity is typically achieved at targets in the middle of the domain, but fails at targets in the extremes; see Figure 1(b) in Lei and Wasserman (2014). So if it were possible to eliminate the edge effect-even if in a trivial way, by eliminating the edge itself-then there is hope for establishing a conditional validity result. In our spatial context, but perhaps not in other cases, it may not be unreasonable to assume that the spatial locations are sampled iid from a uniform distribution on a sphere. Since the sphere has no edges and a uniform distribution has no extremes, there is no "edge effect" preventing conditional validity. Some additional structure in the (X, Y ) process is also needed here, in particular, it should be isotropic in the sense that the correlation structure only depends on the distance between spatial locations. Note that, if the mean of the parametric base model is correctly specified, so that the conditional distribution of Y − Xβ, given X, is free of X, then the stationarity assumption about X can be removed.
To our knowledge, Proposition 2 below gives the first finite-sample conditional validity result for conformal prediction in the literature, albeit under rather strong conditions. Proposition 2. Let (X, Y ) be an isotropic stationary process over the sphere D = {s ∈ R 3 : s = 1}, and suppose that the locations S 1 , . . . , S n , S n+1 are independent and uniformly distributed on D. For Γ α as described above, define the conditional coverage probability function Then the GSCP-based predictions are conditionally valid, i.e., c(s | α, n, P) ≥ 1 − α for all (α, n, s ) and all P under which (X, Y ) is stationary and isotropic and S are iid uniform on D.

A.2: Locally-exchangeable processes
To better understand how the locally-exchangeable processes in Section 4.2 relate to our prediction problem, consider a simple case with no covariates, where {Y (s) : s ∈ D} is the only random process under consideration. In that case, we want to show that T (s) = Y (s) has this local exchangeability property. Then the sufficient condition (11) above amounts to assuming there exists a suitable real-valued function ψ Y , along with appropriate processes L Y and E Y , such that There are a number of common models for continuous responses that meet this condition, including the additive model in Section 2.1, certain generalized spatial linear models Diggle et al. (1998), spatial copula models (Krupskii and Genton, 2019), and max-stable processes (Reich and Shaby, 2012). For example, in a generalized spatial linear model, with suitable spatial process L Y and Gaussian white noise E Y , take where H ξ is the distribution function for an exponential family with natural parameter ξ, g is the link function, and Φ is the standard normal distribution function. Of course, to meet the continuity requirement, the exponential family must have a density with respect to Lebesgue measure. For the practically relevant case with both a response Y and covariate X process, the idea is similar but the notation is more complicated. The goal is to find conditions under which the joint process T (s) = (X(s), Y (s)) has a representation as in (11). Admittedly, it is challenging to consider the joint process directly, which is why we aim to give a simpler sufficient condition based on the marginal distribution of X and the conditional distribution of Y , given X. Consider the following decomposition: Roughly, this amounts to assuming that each of X and Y has a decomposition like that described above for Y alone. Individual assessments of the distributional properties of X and Y are more manageable than directly considering their joint distribution. And as the following simple lemma states, separate considerations of its marginal and conditional structure suffice to establish a decomposition of the joint structure.
Lemma 2. If (X, Y ) can be decomposed as in (15), if the pairs of processes (L X , L Y ) and (E X , E Y ) in are individually L 2 -continuous and locally iid, respectively, and if ( , e) → ψ X ( , e) and ( , e, x) → ψ Y ( , e | x) are both continuous, then T (s) = (X(s), Y (s)) satisfies the conditions of Proposition 1.
As we discussed above, the decomposition of the X marginal as in (15) is quite flexible. For example, it is quite common that X could be expressed as X = L X + E X for a spatial and nonspatial components, L X and E X , respectively, but the additive form is not necessary. And just like in our discussion above Proposition 2 above, if the mean of the parametric base model is correctly specified, then these assumptions about X here can be dropped. Similarly, if the decomposition of Y in the response-only model was flexible, then the corresponding conditional decomposition in (15) must be equally flexible. For example, the same generalized spatial linear model can be considered, but now it is allowed to depend smoothly on the covariate.
Applying the predictive distributions in (16) and provisionally setting Y n+1 = y, for i = n + 1 we can write Therefore, the y satisfying δ i −δ n+1 ≥ 0 are within two roots, denoted as a i ≤ b i , of the quadratic equation U i +V i y +W i y 2 = 0. Then the plausibility calculation (6) becomes which is a step function with jumping points being a i 's and b i 's, and the plausibility function for sLSCP in (13) simply becomes n i=1 w i 1{a i ≤ y ≤ b i }. With the help of the step function, we can solve for the prediction interval directly with no need to enumerate for possible y satisfying p(y|Z n+1 (n+1) ) ≥ t n (α). The GSCP and LSCP steps are summarized in Algorithms 2 and 3, and smoothed LSPC algorithm is given in the main text.

Proofs from Sections 3-4 Proof of Proposition 1
Write the localized version of the T process as T r (u) = ψ L r (u), E r (u) , u ∈ U, Algorithm 2: Global spatial conformal prediction (GSCP).
Input: observations z i = (s i , x i , y i ), i = 1, . . . , n; prediction location and covariate (s n+1 , x n+1 ); non-conformity measure ∆; significance level α; and a fine grid of candidate y n+1 values Output: 100(1 − α)% prediction interval, Γ α = Γ α (z n ; s n+1 , x n+1 ), for the response Y n+1 = Y (s n+1 ) at x n+1 = X(s n+1 ). 1 for provisional values y n+1 in the specified grid do 2 for i = 1, . . . , n + 1 do compute plausibility for y n+1 as p(y n+1 | z n , s n+1 , x n+1 ) in (6); where ψ is a continuous function of two arguments and the L and E components are L 2 -continuous and locally iid, respectively. We are interested in the finite-dimensional distribution of the localized process, so take a fixed set of m vectors u 1 , . . . , u m in U. Define the vectors First, a bit of notation. If w is a generic object, let w ⊗m denote a copy of m versions of w in in rows. For example, if w is a scalar, then w ⊗m = w1 m , where 1 m is an m-vector of unity. Alternatively, if w is a (column) vector, then w ⊗m is a matrix with m identical rows, each containing w . Also, let · denote either the usual Euclidean norm on R m or the Frobenius norm on matrices with m rows, depending on the dimension of its argument. For the L part in the above representation, Markov's inequality implies The expectation in the upper bound can be rewritten as and, since m is fixed, the assumed L 2 -continuity of the L process implies that the right-hand side vanishes as r → 0. Therefore, we have that L m r → L(s ) ⊗m in probability and, hence, in distribution.

Algorithm 3: Local conformal spatial prediction (LSCP).
Input: observations z i = (s i , x i , y i ), i = 1, . . . , n; prediction location and covariate (s , x m+1 ); non-conformity measure ∆; significance level α; and a fine grid of candidate y values Tuning parameter: number of neighbors to consider m ≤ n Output: 100(1 − α)% prediction interval, Γ α s = Γ α s (z m ; x m+1 ), for Y (s ) with X(s ) = x m+1 : 1 form z i , i = 1, . . . , m, based on m locations closest to s ; 2 s m+1 ← s ; 3 for provisional values y m+1 in the specified grid do 4 for i = 1 to m + 1 do 5 z m+1 For the E part in the above representation, locally iid implies where E m 0 is an iid vector. This and the continuous mapping theorem gives where ψ is being applied component-wise. The right-hand side is clearly conditionally iid, given L(s ), which implies exchangeability. Finally, since this holds for all m and for all (u 1 , . . . , u m ), we know that the process T r has a distributional limit, and that limit is an exchangeable process.

Proof of Proposition 2
Let ρ be a generic 3 × 3 rotation matrix that takes a location s on the sphere D to a new location ρs, also on the sphere. Write ρP for the distribution of the spatial process post-rotation by ρ, (ρS, X(ρS), Y (ρS)). From our stated distributional assumptions, namely, that S is uniform on D and (X, Y ) is isotropic and stationary, it follows that the joint distribution is invariant to rotations of the sphere, i.e., ρP = P. Consequently, the conditional coverage probability function satisfies c(ρs | α, n, ρP) = c(ρs | α, n, P).
Next, by the Kriging-based construction of the conformal prediction interval, it is similarly easy to see that the coverage event is invariant to rotations too. That is, if ρZ n is the data after rotating the spatial locations by ρ, then Γ α (ρZ n ; ρS n+1 , X(ρS n+1 )) Y (ρS n+1 ) ⇐⇒ Γ α (Z n ; S n+1 , X(S n+1 )) Y (S n+1 ).
Putting together the equalities in (18) and (19), we conclude that s → c(s | α, n, P) is constant. But the result in Theorem 1 applies in present case; therefore, if the marginal coverage probability exceeds 1 − α, then the constant conditional coverage probability must exceed 1 − α too, which completes the proof.

Proof of Theorem 2
Let T (s) = (X(s), Y (s)) be the joint response-covariate process. We have m spatial locations in a neighborhood of s , which can be expressed as s i = s + ru i , i = 1, . . . , m, where u i ∈ U and r is the neighborhood's radius. Set s m+1 = s and define X i = X(s i ) and Y i = Y (s i ), for i = 1, . . . , m, m + 1. Let T m+1 r denote the collection of all m + 1 triples (s i , X i , Y i ), including the (m + 1) st entry that corresponds to the values at s . The non-conformity scores, δ 1 , . . . , δ m+1 , are functions of T m+1 r , and we will collect these into an (m + 1)-vector δ m+1 r , whose dependence on the radius r is now being made explicit. Then where rank(m + 1; δ m+1 r ) denotes the rank of the (m + 1) st entry within the collection δ m+1 r . By Proposition 1, we have that T m+1 r converges in distribution, as n → ∞ or, equivalently, as r → 0, to an exchangeable T m+1 0 . Since the non-conformity score vector, δ m+1 r , is a continuous function of T m+1 r , it follows from the continuous mapping theorem that, as r → 0, δ m+1 r converges in distribution to δ m+1 0 , say, which is just the non-conformity measure applied to the limit T m+1 0 . Now, turning to the ranks, each permutation of 1, . . . , m + 1 is a possible value of the ranks, and it corresponds to an event A in the space of δ m+1 0 . The boundary of that event consists of cases in which there are ties in δ m+1 0 . By our continuity assumptions, this boundary has probability 0, which makes A a continuity set. Then it follows from the Portmanteau lemma (e.g., van der Vaart, 1998, Lemma 2.2) that P(δ m+1 This holds for the A corresponding to each configurations of the ranks which, in turn, implies rank(m + 1; δ m+1 r ) → rank(m + 1; δ m+1 0 ) in distribution, as r → 0.

Theoretical efficiency of GSCP Setup and statement of the result
Here we establish an asymptotic, theoretical efficiency result for (a version of) our global spatial conformal prediction framework. Since the focus is on asymptotics, we can consider a simplified version, a variation on the so-called split conformal prediction strategy (e.g., Lei et al., 2018;Lei and Wasserman, 2014;Vovk et al., 2005). Without any real loss of generality, assume we have a sample of 2n location-covariate-response triples, i.e., . . , n} and i ∈ I n 2 = {n + 1, . . . , 2n}.
We will use the I n 2 data to estimate the prediction mean and standard deviation functions, and the I n 1 data to do conformal prediction. This simplifies the analysis considerably. To see this, let (μ n ,σ n ) be the prediction mean and standard deviation functions estimated based on data of size n in I n 2 . If we are going to apply the non-conformity measure only to data in I n 1 , then it is just applying a fixed function to each Z i , and we do not need a subscript "(i)" to indicate that observation i was excluded in training the prediction rule. That is, in split conformal, the data (in I n 2 ) for training the prediction rule does not mix with the data (in I n 1 ) for constructing the prediction interval as it does the in full conformal prediction algorithm.
Under this setup, the (split) conformal prediction interval can take a relatively simple form. Letq n,α denote the upper α/2 quantile of {∆ n (Z i ) : i ∈ I n 1 }. If we assume symmetry in the distribution of the signed ∆ n (Z i )'s-see below-then a 100(1 − α)% (split) conformal prediction interval for a new response valueỸ associated with a new location-covariate pair (S,X), iŝ µ n (S,X) ±σ n (S, X)q n,α .
Clearly, the width of this interval is 2σ n (S,X)q n,α . Our goal is to show that this is nearly the width of the "optimal" prediction interval, as n → ∞, and hence that (split) conformal is not only valid but asymptotically efficient.
For the analysis that follows, we assume that the working model presented in Section 2.1 is correct, i.e., that the structural mean is linear in covariates X, that there is an underlying Gaussian process θ with isotropic Matérn covariance function depending on parameters (σ 2 , φ, κ), and a Gaussian nugget with variance τ 2 . This assumption is necessary because very little is known about the convergence properties of the Kriging estimates under misspecification; in other words, we would not be able to characterize an "optimal" prediction interval under weaker assumptions.
When we consider "n → ∞," we are assuming an infill asymptotics regime wherein the range of spatial locations D is a fixed, compact set but the number n (or 2n) of observed triples (S, X, Y ) is going to infinity. Also assume that the range of X is compact. Let µ and σ be the true prediction mean and standard deviation functions under the working model. More specifically, µ = µ θ depends on the Gaussian process θ as is given by and σ (s, x) ≡ τ is a constant, the nugget standard deviation. If all these starred quantities were known, then one could construct a prediction interval for a new responseỸ associated with a new location-covariate pair (S,X) as follows. Define which is just the non-conformity measure applied toZ using the true prediction mean and standard deviation functions. Under the true working model, the standardized residual in the above display is standard normal. Let q α denote the upper α/2 quantile of Normal(0, 1). Then the "optimal" prediction interval forỸ at (S,X) would be µ (S,X) ± q α σ (S,X).
Clearly, the length of this optimal prediction interval is 2q α σ (S,X). Define the absolute difference (ignoring the factor 2) between the widths of the split conformal and optimal prediction intervals as Ξ n (S,X) = |q n,ασn (S,X) − q α σ (S,X)|, which is a function of the data {Z i : i ∈ I n 2 } and the new location-covariate pair (S,X). Our claim is that Ξ n (S,X) vanishes in probability at a certain rate as n → ∞. This is similar to what Lei et al. (2018) prove in their Theorem 6.
Finally, given a vanishing sequence η n → 0, define the event A n as where the mean and standard deviations are functions of location-covariate pairs in a compact domain, and · ∞ is the supremum norm. Note that the best possible rate is η n ∼ c n n −1/2 , where c n → ∞ arbitrarily slowly.
Proof. See below.
Here we make two remarks about the theorem-one about the condition and the other about the conclusion. First, Lei et al. (2018) consider an event similar to our A n . Ours appears more complicated because we work with standardized residuals and, therefore, rely on both prediction mean and standard deviation functions, whereas they use unstandardized residuals and only analyze a prediction mean function. It is well known that under infill asymptotics, the estimators of the spatial correlation parameters generally are not consistent (e.g., Zhang, 2004), but the prediction mean and standard deviation functions can be estimated consistently (e.g., Tang et al., 2019, and references therein). Moreover, under certain regularity conditions as described in Tang et al. (2019), we expect that the rate η n in the above theorem would be n −r , where r = κ 2(κ+1) < 1/2 and κ > 0 is the Matérn covariance function's smoothness parameter. Whether the assumption of Theorem 3 holds for this η n has yet to be confirmed.
Second, the logarithmic term in the rate is a result of the proof technique and so likely is not needed. Therefore, we conjecture that the result could be improved to Ξ n (S,X) = O P (η n ), but this slight improvement in the rate would require a different or substantially more involved proof. So we leave establishing the improved convergence rate result as a topic of future research.
Then for a generic sequence ζ n ,

Split versus plug-in conformal prediction results
In split conformal prediction (e.g., Lei et al., 2018;Lei and Wasserman, 2014;Vovk et al., 2005), the training data are split (equally and randomly) into two groups. The first group is used to estimate the model parameters Θ and is then discarded. The algorithm then proceeds with only the remaining training observations for conformal prediction, including evaluating the base model (given parameters estimates) and plausibilities. Table 3 compares the results using the settings defined in the main document with and without split conformal. Although split conformal maintains good coverage, the intervals are a little wider on average and thus the interval scores are higher.
where the covariate process X(s) is a mean-zero stationary Gaussian process processes with Matérn covariance function with variance σ 2 = 3, range φ = 0.1 and smoothness κ = 0.7. Table 4 summarizes the performance of the methods with and without using the covariate in the analysis. In Scenario 9, the covariate enters the model additively as assumed by all methods and thus including the covariate reduces the interval score. In Scenario 10, the covariate is not additive and the Kriging models are misspecified. However, all conformal methods maintain coverage near the nominal level when N = 40.

Sensitivity analysis
Here we study the proposed method's sensitivity to tuning parameter choice and to the initial Matérn covariance parameter estimates. Because these results are not addressing parameter estimation, in this subsection we estimate the model parameters using the entire dataset (or treat them as fixed, as described in the next paragraph) and then perform leave-one-out cross validation. Figure 4 gives the average interval score of LSCP with different tuning parameter η for Scenarios 5 and 6, along with the average interval score of GSCP, Kriging, and laGP. We find that there is a V-shaped relationship between resulting interval score and tuning parameter η. The average performance is the best when η = 0.2, and LSCP outperforms the other methods regardless the selection of η. To achieve the best performance, in practice, we would select η using cross validation. The proposed spatial conformal prediction methods rely on estimation of the Matérn covariance parameter, Θ. To test for sensitivity to this step, we repeat simulation Scenario 1 since in this case the optimal Θ is known to be the true Matérn covariance parameters used to generate the data. When generating prediction intervals, rather than estimating the Matérn covariance parameterΘ from the simulated data, we use the values ofΘ listed in the first four columns of Table 5 for conformal prediction and Kriging, where each of the four parameters is modified by plus or minus 50% of the true value. We use the GSCP method, α = 0.1 and N = 20 for this sensitivity analysis. Table 3: Performance comparison for simulation scenarios ("Scen") with ("-split") and without split conformal. The metrics are the empirical coverage of 90% prediction intervals ("Cov90"), the width of prediction intervals ("Width") and the interval score ("IntScore"), each averaged over location and dataset. The methods are global (GSCP) and local (LSCP) conformal prediction.  : Sensitivity analysis for simulation Scenarios 5 and 6. The interval score for the LSCP method is computed for several kernel bandwidth parameters, η, and compared to other methods. Table 4: Performance comparison for simulation scenarios with a covariate. The metrics are the empirical coverage of 90% prediction intervals ("Cov90"), the width of prediction intervals ("Width") and the interval score ("IntScore"), each averaged over location and dataset. The methods are global (GSCP) and local (LSCP) conformal prediction, stationary and Gaussian Kriging ("Kriging"). Methods that use the covariate are denoted "(X)". The corresponding average coverage and average width of the conformal prediction intervals and the Kriging intervals are listed in the last four columns. While using the true Θ leads to slightly more narrow intervals, the average coverage of the conformal prediction intervals are between 89.67% and 89.70% across Θ. In contrast, the coverage of the Kriging intervals range from 82.76% to 95.13% across Θ. Therefore, the conformal methods are more robust to the estimation of Matérn covariance parameters than the parametric Kriging model. Table 5: Sensitivity analysis for the simulation study. The first four columns represent the Matérn covariance parameter (Θ) used to generate the conformal and Kriging prediction intervals; the first row gives the true value. The last four columns give the average coverage and interval width of the 90% intervals produced by the conformal and Kriging methods when using the corresponding Θ.

Nugget
Partial