{'fraction_non_alphanumeric': 0.04268015692752426, 'fraction_numerical': 0.01614701631220318, 'mean_word_length': 4.570623418449506, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
DeepKriging: Spatially Dependent Deep Neural Networks for Spatial Prediction

In spatial statistics, a common objective is to predict values of a spatial process at unobserved locations by exploiting spatial dependence. Kriging provides the best linear unbiased predictor using covariance functions and is often associated with Gaussian processes. However, when considering non-linear prediction for non-Gaussian and categorical data, the Kriging prediction is no longer optimal, and the associated variance is often overly optimistic. Although deep neural networks (DNNs) are widely used for general classification and prediction, they have not been studied thoroughly for data with spatial dependence. In this work, we propose a novel DNN structure for spatial prediction, where the spatial dependence is captured by adding an embedding layer of spatial coordinates with basis functions. We show in theory and simulation studies that the proposed DeepKriging method has a direct link to Kriging in the Gaussian case, and it has multiple advantages over Kriging for non-Gaussian and non-stationary data, i.e., it provides non-linear predictions and thus has smaller approximation errors, it does not require operations on covariance matrices and thus is scalable for large datasets, and with sufficiently many hidden neurons, it provides the optimal prediction in terms of model capacity. We further explore the possibility of quantifying prediction uncertainties based on density prediction without assuming any data distribution. Finally, we apply the method to predicting PM2.5 concentrations across the continental United States.


Introduction
Spatial prediction is at the heart of spatial and spatio-temporal statistics. It is aimed at predicting values of a spatial process at unobserved locations by accounting for the spatial dependence in the region of interest. Traditional applications of spatial prediction are in the fields of geological and environmental science (Cressie, 2015), and they have been extended to other fields, such as biological sciences, computer vision, economics and public health (Anselin, 2001;Austin, 2002;Waller and Gotway, 2004;Franchi et al., 2018).
The primary collection of spatial prediction methods are based on the best linear unbiased prediction (BLUP), also referred to as Kriging (Matheron, 1963). Kriging prediction is a weighted average of observed data, where the weights are determined by the spatial covariance function or var-iogram of the random process. Under the Gaussian assumption, Kriging also provides the full predictive distribution. Applying Kriging requires estimating the spatial covariance function, which is commonly assumed to be stationary. However, physical processes tend to be non-Gaussian and non-stationary. For instance, the data on wind speed and fine particles (PM 2.5 ) exposures are positive, right-skewed, and sometimes heavy-tailed (Hennessey Jr, 1977;Adgate et al., 2002), and the spatial covariance typically varies across space, e.g., in urban versus rural areas (Sampson et al., 2013). It is possible to derive the best linear prediction for certain parametric non-Gaussian processes (Xu and Genton, 2017;Rimstad and Omre, 2014) and certain non-stationary covariance structures (Fuentes, 2002;Paciorek and Schervish, 2004;Li and Sun, 2019), but Kriging for more general spatial processes remains an open problem. Another drawback of Kriging is that it is computationally prohibitive for large spatial datasets, since it involves computing the inversion of an N × N covariance matrix, where N is the number of observed locations (Heaton et al., 2019), and the compu- Recently, deep learning or deep neural networks (DNNs) have become the most powerful prediction tools for a wide range of applications, espe-cially in computer vision and natural language processing (LeCun et al., 2015). DNNs are effective for predictions with complex features such as non-linearity and non-stationarity, and they are computationally efficient in analyzing massive datasets using GPUs (Najafabadi et al., 2015). Therefore, it would be promising to apply DNNs to spatial predictions. However, classical DNNs cannot incorporate the spatial dependence appropriately.
Applications in spatial prediction with neural networks usually simply include spatial coordinates as features (Cracknell and Reading, 2014), which may not be sufficient. Recently, convolutional neural networks (CNNs, Krizhevsky et al. 2012) have been claimed to successfully capture the spatial and temporal dependencies in image processing through the relevant filters. However, the framework is designed for applications with a large feature space, and often requires large training labels as the ground truth, which does not fit for many spatial prediction problems, where only in-situ and sparse observations are available.
To tackle the above-mentioned problems, we develop an effective deep neural network for spatial prediction that 1) builds a direct link between DNNs and Kriging in spatial prediction; 2) models the spatial dependence through a set of basis functions; 3) does not require matrix operations and is scalable for large datasets; 4) provides a non-linear predictor in covariates or generally in observations; 5) has a Gaussian process representation, and brings more flexible spatial covariance structures than simply using the coordinates as features; 6) suits for different data types, e.g., non-Gaussian or non-stationary data; 7) potentially measures the uncertainty through predictive density functions without assuming any data distribution.
We call our method "DeepKriging" with the aim of achieving the optimal spatial prediction, similar to the original use of Kriging (Cressie, 1990), but by using deep neural networks. We also conduct simulation studies and apply our approach to the PM 2.5 concentration data across the continental United States to show the performance of DeepKriging compared to Kriging and other naive DNN methods. The rest of our paper is organized as follows. Section 2 introduces the construction of our DeepKriging method. Section 3 provides its theoretical properties. Section 4 presents some simulation studies to show the performance of DeepKriging. Section 5 applies DeepKriging to predict PM 2.5 concentration in the U.S. Section 6 summarizes the main results and suggests directions for future work.

Deep learning in spatial prediction
Suppose z = {z(s 1 ), . . . , z(s N )} T are measurements observed at N spatial locations from a real-valued spatial process {Y (s) : s ∈ D}, D ⊆ R d . The goal of spatial prediction is to find the optimal predictor Y opt (s 0 ) of the true process at an unobserved location s 0 , as a function of z. In decision theory, Y opt (s 0 ) is the minimizer of an expected loss function or risk function (DeGroot, 2005). That is, where L(·, ·) is a loss function and R(·, ·) is a risk function. Under the mean squared error (MSE) loss, the optimal predictor is Y opt (s 0 ) = E{Y (s 0 )|z} if it is finite. This predictor has multiple good properties such as unbiasedness and asymptotic normality under regularity assumptions (Lehmann and Casella, 2006). In particular, if Y (s 0 ) and z are jointly Gaussian, the conditional mean is a linear combination of z; if Y (s 0 ) and z are not jointly Gaussian, the conditional mean obtained with Gaussian assumption remains the best linear unbiased prediction (BLUP), which is called Kriging. However, as mentioned before, the Kriging predictor is sub-optimal for non-Gaussian data, and it is not scalable for large data size.

Deep learning in spatial prediction
In this work, we use deep learning to approximate the optimal predictor Y opt (s 0 ) in (2.1) by the output of the neural network. The optimal neural network predictor is given by f opt where f NN (·) ∈ F can be any function in the function space F expressible by a family of neural networks, and f opt NN (·) is the best function in F in terms of minimizing a certain risk R(·, ·). The inputs of the neural network can be relevant covariates x(s 0 ) and other features at s 0 . Typically, we write f NN (s; θ) as a parametric model with unknown parameters θ, which include the weights and biases in the neural network. Note that the opti- is unknown. In practice, we approximate the predictor by minimizing the empirical loss function over the training set z (Goodfellow et al., 2016); i.e., L{f NN (s n ; θ), z(s n )}. (2.2) Applying this framework of classical neural network directly to spatial prediction is problematic in at least two aspects: classical DNNs does not account for the spatial dependence, and spatial prediction typically has limited observed features rather than excessive features in common applications of neural networks. In particular, assume that the spatial process , where x(s) ∈ R P is a vec-2.1 Deep learning in spatial prediction tor process of P known covariates, β is a vector of coefficients, and ν(s) is a spatially dependent and zero-mean random process with a generally non-stationary covariance function: Cov(ν(s), ν(s )) = C(s, s ). In neural networks, we usually assume that Y (s) are mutually independent conditional on the features x(s). However, this assumption is not reasonable in spatial prediction because the covariates x(s) only contribute to the mean structure of Y (s) and ν(s) remains a spatially correlated process. Hence, more features apart from x(s) are needed to model the spatial dependence in applying the neural networks.
To account for the spatial information, the most natural way is to add d coordinates (e.g., longitude and latitude) to the features, in the hope that the neural networks can learn the dependent term ν(s) as a function of s (Cracknell and Reading, 2014). By doing that, the adjusted features become x adj (s) = (x(s) T , s) T . However, this does not help much in enlarging the feature space since usually the dimension of coordinates has d ≤ 3.
Moreover, the associated neural network may not be efficient, since if the true function is far from linear, it may take huge effort for the neural network to achieve a good approximation. For instance, the optimal predictor under the Gaussian assumption and MSE loss is the Kriging predictor, which is linear in x(s) but obviously non-linear in coordinates s; this is a special 2.1 Deep learning in spatial prediction case where the natural structure of neural networks may not work.
Going deeper into the form of Kriging prediction may give us a hint about the appropriate way to incorporate the spatial dependence in the DNN. Suppose z is observed from a generalized additive model: where X = (x(s 1 ), . . . , x(s N )) T is an N × P matrix, c(s 0 ) = Cov(Z, Z(s 0 )), Σ = Cov(Z, Z T ), and β = (X T Σ −1 X) −1 X T Σ −1 z. The spatial dependence is incorporated in Y UK (s 0 ) via a linear function of the covariance vector c(s 0 ), but its coefficient Σ −1 (z − X β) is unknown. This motivates us to use a set of known nonlinear functions as the embedding of s in the features to characterize the spatial process ν(s) in the neural network. This can be done by the lights of the Karhunen-Loève (KL) theorem (Adler, 2010), which establishes that ν(s) admits a decomposition ν(s) = ∞ k=1 w k φ k (s), where w k 's are pairwise uncorrelated random variables and φ k (s)'s are pairwise orthogonal basis functions in the domain of ν(s). Hence, ν(s) can be linearly quantified by nonlinear basis functions of s.

DeepKriging: a spatially dependent neural network
In practice, the prediction of ν(s) is typically the truncated KL expansion based on the property that given any orthonormal basis functions φ k (s), we can find some large integer K, so that ν(s) can be approximated by the finite weighted sum of basis functions, i.e., ν(s) = K k=1 w k φ k (s). Based on the KL theorem, the form of basis functions is not as important as the number of basis functions to approximate the spatial random effect ν(s).
This can also be supported by the additional simulations we conduct in Section S4.1 of the Supplementary Material. Multiple types of basis functions can be used, such as the smoothing spline basis functions (Wahba, 1990), the wavelet basis functions (Vidakovic, 2009), and the radial basis functions (Friedman et al., 2001). By adding an embedding layer with sufficiently large K, the width of the neural network is greatly increased so that the network incorporates more spatial information than using the coordinates alone. A similar idea has been used in the recommendation systems First, we need to choose the value for K and basis functions to approximate the spatial process ν(s). We adopt the idea in Nychka et al. (2015), who developed a multi-resolution model for spatial prediction for large datasets. The radial basis functions at each level of resolution are constructed using a Wendland compactly supported correlation function with the nodes arranged on a rectangular grid. In particular, at a certain level of resolution, let {u j }, j = 1, . . . , m, be a rectangular grid of points (or node points in the radial basis function terminology) and let θ be a scale parameter. The basis functions are given by φ * Hence, the embedding layer uses mutual distance locally to each knot location, implying that the spatial patterns are location invariant locally.
As a result, the proposed DeepKriging is able to model the spatial nonstationarity; as we will show in Section 3.3, the induced covariance functions of an infinitely wide DeepKriging network are in general non-stationary.
The scale parameter θ is set to be 2.5 times the associated knots spacing according to Nychka et al. (2015). The grid at each finer level increases by a factor of two and the basis functions are scaled to have a constant 2.2 DeepKriging: a spatially dependent neural network overlap. In particular, in the h-th level, the number of knots is chosen to be where d is the spatial dimension. For a massive dataset and to obtain K ≥ N , we need H = 1 + log 2 ( d √ N /10) levels. Therefore, for a four-level model for instance, we need K = 10 + 19 + 37 + 73 = 139 basis functions in one dimensional space and K = 10 2 +19 2 +37 2 +73 2 = 7159 basis functions in two dimensional space. This scheme gives a good approximation to standard covariance functions and also has the flexibility to fit more complicated shapes. The approach of multi-resolution approximation for massive spatial datasets has also been adopted in other research works; see Katzfuss (2017) and the references therein.
Then, for any coordinate s, we compute the K basis functions to get the embedded vectors φ(s) = (φ 1 (s), . . . , φ K (s)) T . The basis functions are recommended to be orthogonal based on the KL expansion. Then, let x φ (s) = (x(s) T , φ(s) T ) T be the embedded input of length P + K, and specify an L-layer DNN as u 2 (s) = W 2 a 1 (s) + b 2 , a 2 (s) = ψ 2 (u 2 (s)); . . . (2.4) For the l-th layer with N l neurons, W l is the N l × N l−1 weight matrix, b l is 2.2 DeepKriging: a spatially dependent neural network the bias vector of length N l , a l is the neuron vector of length N l , and ψ l (·) is the activation function. The output of this neural network is f DK (s), which is a function of the weights and biases. Let θ be the vector of unknown weights and biases, and θ be the estimate via Equation (2.2) based on the training sample. The final DeepKriging prediction at an unobserved location s 0 is defined as Y DK (s 0 ) = f DK (s 0 ; θ).
One major advantage of our DeepKriging method is that we can adjust the number of neurons, activation functions and loss functions to fit for different data types and model interpretations. For example, for predicting continuous variables as in a regression problem, we choose N L = 1, ψ L (·) to be an identity function, and the loss function to be the MSE. Figure 1 provides a visualization of a DeepKriging structure in two dimensional prediction for continuous data. For predicting categorical variables as in a classification problem, we choose N L to be the number of categories, ψ L (·) to be a softmax function, and the loss function to be the cross entropy loss.
For the activation functions in the hidden layers, we choose the rectified linear unit (ReLU) in default, which allows us to keep the linear relationship in the KL expansion but add some deactivated neurons to select the best number of basis functions. The DeepKriging structure also allows the covariate effects to be spatially varying.

Theoretical Properties of DeepKriging
DeepKriging provides a novel spatial prediction framework using deep learning. It differs from classical Kriging methods in several aspects. First, Kriging prediction is a linear combination of observations; in contrast, DeepKriging prediction is linked to the observations via the weights and biases through model training and is typically nonlinear in observations (see Section S3.1). Second, DeepKriging does not assume a Gaussian process with a certain covariance function but models spatial dependence by basis functions. Last, unlike Kriging which predicts the random process Y (s) at an unobserved location, DeepKriging approximates the process using a deterministic continuous function.
In this section, we provide important theoretical properties of DeepKriging including 1) the underlying relationship between DeepKriging and Kriging; 2) how accurate Deepkriging can be in terms of the prediction error compared to Kriging; and 3) how the spatial dependence is measured in the DeepKriging framework. These three aspects are critical for understanding our DeepKriging method, and will be illustrated in the following subsections, respectively.

The link between DeepKriging and Kriging-based methods
DeepKriging is closely related to Kriging and its associated variants, which can be classified as multi-resolution processes (Nychka et al., 2015;Kleiber and Nychka, 2015;Katzfuss, 2017) and Gaussian predictive processes (Banerjee et al., 2008(Banerjee et al., , 2010, all leading to spatial predictions that can be treated as linear functions of embedded features x φ (s 0 ), and thus can be potentially approximated by DeepKriging.
One example is the fixed rank Kriging (FRK) proposed by Cressie and Johannesson (2008), who used one of the low-rank approximations of the covariance matrix in order to speed up the computation of universal Kriging.
Similar to DeepKriging, they represent the spatial random effects ν(s) by K basis functions, i.e., ν(s) = φ(s) T η, where η is a K dimensional Gaussian random vector with Cov(η) = Σ K . They also assume that the model for The FRK prediction as a linear function of z is given by and Σ −1 has a computationally simple form which involves inverting the fixed rank K×K positive definite matrix Σ K and the N ×N diagonal matrix P covariates x(s 0 ) and K basis functions φ(s 0 ). This is a special case of DeepKriging when we set all of the activation functions to be linear.
FRK usually chooses K to be much smaller than N in order to speed up the computation for large datasets. Since the covariance ΦΣ K Φ T has at most rank K, such a low-rank approximation of the covariance matrix may fail to capture the high-frequency variation or small-scale spatial dependence in the spatial process (Stein, 2014). In contrast, for DeepKriging, K needs to be sufficiently large (K > N ) in order to have a good approximation of the spatial random effect ν(s), so that our method captures more spatial information in the prediction.
By setting K = N in the FRK, we can see that the (universal) Kriging A detailed proof is provided in Section S1.1 in the Supplementary Materials. This result implies that the Kriging prediction with any covariance function can be linearly expressed by the embedding features x φ (s 0 ). In this 3.2 DeepKriging in decision theory sense, DeepKriging generalizes Kriging by allowing for nonlinear functions of x φ (s 0 ) in the prediction.

DeepKriging in decision theory
Our DeepKriging prediction procedure conventionally follows an approximation- The approximation error relates to the model capacity and is defined as the risk between the true process Y (s 0 ) and the optimal predictor Y opt  (2001)) claims that every continuous function of the features x φ (s), denoted as C(x φ ), can be arbitrarily well approximated with a feed-forward neural network with a single hidden layer that contains finite number of hidden neurons and with arbitrary activation function. This indicates that the optimal DeepKriging prediction with a single hidden layer and finite loss function has the largest model capacity . This result holds for any type of data (i.e., continuous or discrete) and for any type of task (i.e., regression or classification). Therefore, the optimal DeepKriging prediction has larger capacity than the Kriging prediction in terms of minimizing the approx- The detailed proof is provided in Section S1.2 in the Supplementary Materials.
Similarly, the optimal DeepKriging prediction also has larger model capacity than the FRK prediction. FRK can be seen as DeepKriging with a single hidden layer containing finite number of neurons and a linear activation function. By allowing for a large number of basis functions, multiple layers, more flexible activation functions and a wide network, DeepKriging yields non-linear predictions that can appropriately capture the spatial 3.3 DeepKriging as a Gaussian Process dependence in the spatial process.

DeepKriging as a Gaussian Process
Neal (1994) showed that a single-layer fully-connected neural network with an i.i.d. prior over its parameters (i.e., weights and biases) is equivalent to a Gaussian process (GP), in the limit of infinite network width (i.e., infinite number of hidden neurons). Later, Lee et al. (2017) derived the exact equivalence between infinitely wide deep networks and GPs. Consequently, a similar correspondence to GPs also holds for our DeepKriging network.
We start from a regression-type DeepKriging model with a singe hidden Weights (w 1 j , w 0 ji ) and biases (b 1 , b 0 j ) are independent and randomly drawn to have zero mean and variances σ 2 w /N 1 and σ 2 b , respectively. Consequently, the post-activations a 1 j and a 1 j are independent for j = j . Moreover, sincê Y DK (s) is a sum of i.i.d terms, it follows from the Central Limit Theorem that in the limit of infinite width N 1 → ∞,Ŷ DK (s) will be Gaussian distributed. Likewise, from the multi-dimensional Central Limit Theorem, any finite collection of {Ŷ DK (s 1 ),Ŷ DK (s 2 ), · · · ,Ŷ DK (s n )} will have a joint 3.3 DeepKriging as a Gaussian Process multivariate Gaussian distribution, which is exactly the definition of a GP.
Therefore, we conclude that with sufficiently large N 1 ,Ŷ DK is a GP with zero mean and covariance function where C(s, s ) is obtained by integrating against the distribution of w 0 , b 0 as in .
For DeepKriging with deeper layers, the induced covariance function can be obtained in a recursive way according to Lee et al. (2017): where F ψ (·) is a deterministic function that depends only on the activation function ψ. An iterative series of computations can be performed to obtain the covariance C L for the GP describing the network's final output,Ŷ DK (s).
For the base case, The aforementioned results require the assumption of infinitely many hidden neurons in each layer. However, when the prior distribution of weights and biases is Gaussian, this condition is not needed. b and w are combined biases and weights, respectively. In this case, the induced covariance function ofŶ DK is given by which is the basis approximation of a spatial covariance function.
In the case of ReLU non-linearity, Equation (3.6) has a closed form of the well-known arc-cosine kernel (Cho and Saul, 2009): When no analytic form of the resulted covariance function exists, it can be computed numerically, as described in Lee et al. (2017).
Consider a regression-type DeepKriging model with a single hidden layer and no covariates effects. It can be shown that with infinitely many hidden neurons, the covariance function of the outputŶ DK (s) for any two nearby locations has the form where φ(s) is the basis vector at location s, v(s) > 0 is related to the variance when s = s , and c is the scaling parameter. The proof is provided in Section S1.3 in the Supplementary Materials. As a special case, if only the coordinates are used in the features, then φ(s) − φ(s ) 2 = s − s 2 ,

DeepKriging as a Gaussian Process
v(s) = v(s ) = v and thus C(s, s ) = v − c s − s 2 , which contains less information than in Equation (3.7). Therefore, the embedding layer in DeepKriging brings more flexible spatial covariance structures than simply using the coordinates.
Further, we can show how the DeepKriging induced covariance function can approximate the common stationary covariance functions in spatial statistics. Let the basis functions be φ l (s) = k(s, u l ) based on a certain kernel function k(·, ·) and knot u l , l = 1, · · · , K. If the u l 's form a fine grid of knots covering the spatial domain, then Note that the last term is the kernel convolution approximation to a covariance function. Higdon (2002) shows that by selecting an appropriate kernel function, we can approximate any stationary covariance function based on the kernel convolution. Further, the induced covariance function of DeepKriging also possesses favorably physical interpretations. For example, DeepKriging can yield the Matérn covariance function, also commonly used in Kriging since it is related to a stochastic partial differential equation (SPDE) of Laplace type (Whittle, 1954). In addition, DeepKriging can induce a GP that approximates a fractional Brownian motion based on the example of DNN provided in .

DeepKriging on a 1-D Gaussian process
We first consider the performance of DeepKriging when data are simulated from a 1-D stationary GP, where the Kriging prediction is optimal. We also compare DeepKriging to two naive DNNs: a DNN with only the intercept

DeepKriging on 2-D non-stationary data
In this section, we evaluate the performance of DeepKriging on 2-D nonstationary data so that the procedure is designed to resemble the real data application in Section 5. The simulation details are included in Section S3.2 of the Supplementary Material.
We use the 10-fold cross-validation method to show the performance of 4.2 DeepKriging on 2-D non-stationary data DeepKriging, Kriging with an estimated stationary covariance function and the baseline DNN with only coordinates s in the features. We calculate the RMSEs and MAPEs on the testing dataset, and show the results in Figure   S2(b) and is essential to obtain a high-resolution map of PM 2.5 exposure in order to assess its impact. The measurements from monitoring networks are the best characterization of PM 2.5 concentration at a given time and location.
However, data from monitoring locations are often sparsely distributed so that they are out of spatial and temporal alignment with health outcomes.
Meanwhile, it is known that PM 2.5 concentration is associated with meteorological conditions such as temperature and relative humidity (Jacob and Winner, 2009), where the meteorological data or data products are often easy to access with good spatial coverage and resolutions. Hence, the interpolation of PM 2.5 concentration by making use of data from monitoring networks and other meteorological data has been a promising field of research (Di et al., 2016), where spatial prediction plays a central role.
The modeling and prediction of PM 2.5 concentration are challenging.
First, PM 2.5 concentration data are obviously non-Gaussian, and thus classical Kriging methods are inappropriate here. Second, PM 2.5 data from monitoring stations are irregular and sparse, but many interpolation methods require lattice data. Third, it is more important but challenging to understand the risk of high pollution and predict pollution levels, such as 5.2 Data and preprocessing being low, medium and high; statistically, these two questions are related to estimating the probability over a threshold and a classification problem, respectively. Quantile regression and convolutional neural networks have been employed to overcome some of the above issues (Reich et al., 2011;Porter et al., 2015;Di et al., 2016). However, a unified method to handle all of the aforementioned tasks has not yet been sufficiently developed.

Data and preprocessing
To tackle the above-mentioned problems, we apply the proposed DeepKriging method to the spatial prediction of PM 2.5 concentrations based on meteorological variables. Meteorological data are obtained from the NCEP North American Regional Reanalysis (NARR) product. Reanalysis is a gridded dataset that represents the state of the atmosphere, incorporating observations and outputs of numerical weather prediction models from past to present-day. Reanalysis data are often used to represent the "true state" of the atmosphere according to observations, and thus we use the reanalysis data as the "observed data" for the covariates.

Model fitting and results
Our aim is to predict the PM 2.5 concentration values at unobserved grid cells where the six meteorological variables are provided. We use the 10-fold concentration, the level of pollution and the risk of high pollution level over the threshold 12 µg/m 3 at unobserved locations based on the NARR data.  ure 2(b) shows a smooth map of the predicted PM 2.5 concentration from DeepKriging. We also provide the distribution prediction (details and algorithms are included in Section S5 in the Supplementary Material) in order to obtain the predicted risk defined as P{PM 2.5 > 12 µg/m 3 }, shown in Figure 2(c). This map implies that high PM 2.5 pollution risks exist over a vast area of Eastern US. We further compare the results to the Kriging prediction in Figure 2(d), which implies that DeepKriging provides more local features/patterns than Kriging.

Discussion
In this work, we have proposed a new spatial prediction model using deep neural networks which incorporates the spatial dependence by a set of ba-  1983,1984). This work has provided a new perspective on deep learning in spatial prediction with a large number of basis functions. We have shown that the proposed method is superior to Kriging in many aspects both theoretically and numerically in our simulation and real application. For instance, DeepKriging is more scalable for large datasets and suits for more data types than Kriging. DeepKriging also has a GP representation with flexible spatial covariance structures, which enables Bayesian inference on regression tasks by evaluating the corresponding GP. More importantly, the proposed DeepKriging framework connects the regression-based prediction and spatial prediction so that many other machine learning algorithms can be applied.
In general applications, it is possible that the covariates at the new location s 0 are not observed. One promising approach for coping with this problem is to find the true values of the missing covariates for a subset of the observations and then train a machine learning algorithm to predict the values of those covariates for the rest (see, e.g., Imai and Khanna (2016)).
However, Fong and Tyler (2021) showed that plugging in these predictions without regard for prediction error renders regression analyses biased, inconsistent, and overconfident. They described a procedure to avoid these inconsistencies. This approach combines a new sample splitting scheme and a general method of moments (GMM) estimator to make an efficient and consistent estimator. Overall, it is non-trivial to address the problem of missing covariates: intuitive strategies such as plugging in machine learning predictions lead to bias and inconsistency, while the implementation of a more complicated method such as that in Fong and Tyler (2021) requires extra assumptions (e.g., the exclusion restriction condition) and increases computational burden. If the goal is to predict both the response and the covariates instead, a multivariate version of DeepKriging could be developed. These are left as our future work.

Supplementary Materials
The Supplementary Material contains details referenced in the main manuscript, including the proofs of the lemmas and theorems (Section S1), the settings for the DeepKriging network structure (Section S2), details of the simulation studies (Section S3), additional simulation studies (Section S4), dis-

Supplementary Material
This supplementary material contains the proofs (Section S1), the settings for the DeepKriging network structure (Section S2), details of the simulation studies (Section S3), additional simulation studies (Section S4), distribution prediction and uncertainty quantification (Section S5), and the source codes and data for reproducible research (Section S6).

S1.1 Kriging prediction is linear in x φ (s) = (x(s), φ(s)) T
Proof : Let the covariance matrix associated with the random process ν(s) be Σ ν , where Σ ν i,j = C(s i , s j ). We can build a spatial random ef- {φ(s 1 ), . . . , φ(s N )} is an N × K basis matrix with rank N , and Φ −1 R is the right inverse of Φ. From the result of Banerjee (1973), the right inverse exists. Hence, Cov(η) is also a valid covariance matrix and Cov(µ) = Therefore, the Kriging prediction of µ(s) is the same as that of ν(s), i.e., Y UK µ (s 0 ) = Y UK ν (s 0 ). On the other hand, based on the spatial random effect representation of µ(s), we can get the equivalent fixed rank Kriging prediction shown in the Equation (5) in the manuscript, such that Csáji (2001)) shows that ∀f ∈ C(x φ ), ∀ε > 0: ∃n ∈ N and certain choice of weights and biases, such that |f (s 0 ) − A n f (s 0 )| < ε, where A n f is the mapping from the standard multi-layer feed-forward networks with a single hidden layer that contains n hidden neurons. Then we have where A opt n f (s 0 ) is the optimal neural network that approximating f (s 0 ).
In the context of spatial prediction, suppose the optimal prediction in and E{| Y opt F DK (s 0 )−Y (s 0 )|} 2 < ∞ by assumptions. Using the dominated convergence theorem, we have lim

S1.3 Proof of Equation (3.7)
S1.3 Proof of Equation (3.7) Proof : Recall that the induced covariance function ofŶ DK (s) is According to , When s is close to s , we have ψ 1 (x) − ψ 1 (y) = α(x − y) for a smooth activation function, where α is a scaling coefficient. When the covariate vector x(s) is not available, we have a 1 Finally, the covariance function for nearby locations is

S2. Settings for the DeepKriging network structure
The model settings for the DeepKriging network structure include the following considerations. First, when the sample size is small, we add a dropout layer to avoid overfitting. For a large dataset with Gaussian priors, dropout and other types of regularization are not very helpful. Second, batch-normalization is another useful strategy since the covariates typically have different units, and the scales of basis function may vary too much for irregularly spaced spatial data. Third, we normalize the covariates before we run the automated batch-normalization since the covariates and basis functions required different types of normalization. Last, when the data are irregularly spaced, the value of the compactly supported basis functions for some knots that are far apart from any other locations will be zero; in this case, we remove these knots and the relevant columns in the basis matrix.
To summarize, the default setting of DeepKriging network is as follows: (1) Normalize (min-max normalization) the observed covariates x(s); (2) Build the embedding layer using a 3 to 5 level multi-resolution radial basis functions φ(s) with the corresponding basis matrix Φ; the form of the basis function and the choice of the number of basis functions K are illustrated at the beginning of Section 2.2; (3) Remove the all-zero columns of the basis matrix Φ; (4) Add the 1st dense layer with 100 hidden neurons and ReLu activation; (5) Add the 1st dropout layer with 0.5 dropout rate; (6) Add the 1st batch-normalization layer; (7) Add the 2nd dense layer with 100 hidden neurons and ReLu activation; (8) Add the 2nd dropout layer with 0.5 dropout rate; (9) Add the 3rd dense layer with 100 hidden neurons and ReLu activation; (10) Add the 2nd batch-normalization layer; (11) Add the output layer.
The default number of epochs for DeepKriging is 200. For regression tasks, the loss function is set to be the mean squared error (MSE), while for classification tasks, the loss function is set to be the cross-entropy loss. For optimization, we use the Adam optimizer (Kingma and Ba, 2014). One can tune the hyper-parameters in the network (e.g., the value of K, the number of layers, the width of each layer, the dropout rate, the batch size and the number of epochs) in order for the optimization to converge and to achieve a better performance compared to other basic methods such as Kriging.

S3.1 DeepKriging on a 1-D Gaussian process
We first consider the performance of DeepKriging when data are simulated from a 1-D stationary GP, where the Kriging prediction is optimal. We also compare DeepKriging (with x(s) and an embedding layer with basis functions of coordinate s as the input) to a DNN that does not account for the spatial dependence (with only x(s) as the input), and a DNN that incorporates the spatial dependence by including directly the coordinates in the features (with x(s) and coordinate s as the input). Specifically, the data are generated from a GP with a constant mean: z(s) = µ+ν(s)+ε(s), where µ = 1, ν(s) is zero mean GP with an exponential covariance function C(s, s ) = σ 2 exp{−|s − s |/ρ}, where the variance σ 2 = 1 and the range parameter ρ = 0.1, and ε(s) is a Gaussian white noise with the nugget vari-S3.1 DeepKriging on a 1-D Gaussian process ance τ 2 = 0.01. We generate 100 replicates for {z(s 1 ), · · · , z(s N )} from the GP with N = 1, 000 equally spaced locations over [0, 1], with 800 locations randomly selected as training data and the remaining treated as testing data. In this example, there are no observed covariates except for the intercept, i.e., x(s) = 1 for any s ∈ R. We also compare these models to Kriging prediction with the true exponential covariance function and that with an DeepKriging, thus K = 10 + 19 + 37 + 73 = 139. Figure S1 shows the prediction for one of the sample datasets using each of the five prediction methods. Table S1 shows the root mean squared error (RMSE) and mean absolute percentage error (MAPE) on the both the training and testing sets over the 100 replicated samples.

S3.2 DeepKriging on 2-D non-stationary data
In this section, we evaluate the performance of DeepKriging on 2-D nonstationary data so that the procedure is designed to resemble the real data application in Section 5 of the main manuscript. The goal is to predict values from the true process: Y (s) = sin{30(s − 0.9) 4 } cos{2(s − 0.9)} + (s − 0.9)/2, where s = (s x , s y ) T ∈ R 2 ands = (s x + s y )/2. A similar example is evaluated in many computer experiments (Ba et al., 2012), where both Kriging and neural networks are popularly applied. In our simulation, S3.2 DeepKriging on 2-D non-stationary data Table S1: Mean and standard deviations (SD) of the root mean squared errors (RMSEs) and mean absolute percentage errors (MAPEs) on both the training and testing sets across the 100 datasets from the five predictions of a Gaussian process. Kriging I and II are the Kriging prediction with the true exponential covariance function and that with an estimated Matérn covariance function, respectively. DeepKriging prediction is based on the MSE loss. DNN I and DNN II are the DNN prediction with intercept only and that with intercept and coordinates, respectively.  Figure S2(a). This process presents obvious spatial non-stationarity; for example, the smoothness over the region [0, 0.4] 2 is significantly smaller than that over [0.4, 1] 2 .
The network structures of DeepKriging as well as the baseline DNN are the default setting as stated in Section S2, except that the batch size is set to be 64, and the number of hidden layers is set to be 4 in DeepKriging, which yields the best performance (in terms of RMSE) for both the training and testing datasets among networks with 1 to 6 hidden layers (results not shown). A three-level multi-resolution model is used to generate the basis functions for DeepKriging, thus K = 10 2 + 19 2 + 37 2 = 1, 830. We use the 10-fold cross-validation method to show the performance of DeepKriging, Kriging with an estimated stationary covariance function and the baseline DNN with only coordinates s in the features. The boxplot of RMSEs are shown in Figure S2(b). The RMSEs and MAPEs on both the training and testing sets are shown in Table S2.

S4.1 Choice of basis functions
As we have stated in the main manuscript, based on the KL theorem, the form of basis functions is not as important as the number of basis functions to approximate the spatial random effect ν(s). This can be supported by the additional simulations we conduct below.   In Nychka et al. (2015), they stated that "the choice of the Wendland family of RBFs is not crucial and other compactly supported, positive definite functions will work." They used compactly supported kernels so that the key matrices are sparse to improve the computational efficiency. Our perimental designs for the 1-D simulation as before. Figure S3 shows the prediction for one of the sample datasets using each of the six prediction methods. DeepKriging prediction with both basis function choices almost overlap with the optimal Kriging prediction. We then compare the RMSEs and MAPEs on both the training and testing datasets, as shown in Ta Similarly for the 2-D simulation study in Section 4.2, replacing the radial basis functions using Wendland kernel with those using Gaussian kernel does not change the conclusion, if we compare Figure S4(b) with Figure S2(b), and compare Table S4 with Table S2 Table S4: Model performance for both the training and testing sets based on 10-fold cross-validation. Data are generated from Y (s) = sin{30(s − 0.9) 4 } cos{2(s − 0.9)} + (s − 0.9)/2, where s = (s x , s y ) T ∈ [0, 1] 2 ands = (s x +s y )/2. RMSEs and MAPEs are computed, and the mean and standard deviations (SD) for the 10 sets of validation errors are given. DeepKriging I and II are the DeepKriging predictions with Wenland kernel and Gaussian kernel, respectively, based on the MSE loss. The baseline DNN includes only coordinates s in the features. The Kriging prediction uses an estimated exponential covariance function.

S4.2 DeepKriging is non-linear
It is important to investigate how the DeepKriging prediction is linked with the observations. We generate 100 observations generated by z(s) =  Figure S6 shows the results. Although Kriging is faster for small sample sizes (N < 1, 500), DeepKriging is much more scalable when the sample size increases. For example, when N = 12, 800, which is the largest sample size we have considered, it takes more than 1.5 hours (5, 663 seconds) to train a Kriging model, but DeepKriging only costs 3.5 minutes (214 seconds) without GPU acceleration and 1.5 (94 seconds) minutes with a Tesla P100 GPU. Note that the GPU we use is freely accessed in Kaggle.
For a more powerful GPU, the computational cost will be further reduced. Bayesian learning . Based on the correspondence between an infinitely wide DeepKriging and GPs as illustrated in Section 3.3 in the main manuscript, we can potentially compute the covariance function for the induced GPs and then use the resulting GPs to perform Bayesian inference for wide DeepKriging networks. However, it is infeasible to train an infinitely wide DNN in practice, and Lee et al. (2017) found that the test performance increases as finite-width trained networks are made wider and S5.1 Deep Distribution Spatial Regression more similar to a GP, and thus that GP predictions typically outperform those of finite-width networks. The related computation burden is another issue for applications in large datasets.
The Bayesian methods mentioned above require configurations on the prior distributions for hyper-parameters.  proposed a completely distribution-free method for density estimation called deep distribution regression (DDR), where the density was approximated using histograms with discrete bins, and the discretized density was estimated by multi-class classification using neural networks. Although the method was designed for a regression problem, we can extend the idea to our DeepKriging framework for predictive distribution estimation. The idea is straightforward and the method is easy to implement. We call the method deep distribution spatial prediction (DDSP), as illustrated below.

S5.1 Deep Distribution Spatial Regression
We quantify the uncertainty of DeepKriging prediction by the conditional distribution of Y (s 0 )|z, x φ (s 0 ) at an unobserved location s 0 . We denote the probability density function (PDF) as f UQ (y|x), whose support is [l, u], |T im | 1{y ∈ T im } Note that although the density regression provides an uncertainty quantification without any distribution assumption, it may bring new uncertainties from the neural networks. Therefore, the GP process representation may be a better way for uncertainties if it is quantified by the standard deviation, even though the computational burden will increase.
In addition, the density prediction requires further computations if we apply ensembles, where each step in the ensemble will implement classification separately with a new choice of random cut-points. However, the main objective of the ensemble in the density prediction is to obtain a smooth density. If the research goal is to get uncertainties using quantiles without the density curve, then the ensemble is not always needed.
Another hyperparameter we need to specify in the density prediction is the bin size M + 1. Typically, large M will provide a more complex pat-S5.2 Uncertainty Quantification: A Simulation Study tern of the density, and small M will give a more stable estimation. In practice, we can follow the Freedman-Diaconis rule (Freedman and Diaconis, 1981) as a robust estimation in the histogram approximation such that M = {max(z) − min(z)} 3 √ N 2×IQR(z) , where IQR(z) is the interquartile range of the observations.

S5.2 Uncertainty Quantification: A Simulation Study
The simulated data is generated by a 1-D Gaussian mixture model: z(s) = {sin(5s) + 0.7 + τ 1 (s)}π(s) + {2 sin(8s) + τ 2 (s))}(1 − π(s)), where τ 1 (s) iid ∼ N (0, 0.2 2 ), τ 2 (s) iid ∼ N (0, 0.3 2 ), π(s) iid ∼ Bernoulli(0.5), and τ 2 (s), τ 2 (s), and π(s) are mutually independent. Observations are drawn from 2, 500 regularly spaced locations in [0, 1], out of which 2, 000 are training samples and 500 are testing samples. Figure S7 shows the prediction results on both training and testing samples. The advantage of DeepKriging is obvious for the uncertainty quantification. The process is heteroscedastic with an obviously large variance when s is larger than 0.8. The heteroscedasticity is well captured by the DeepKriging but missed by the Kriging. Moreover, the prediction band of the Kriging is too narrow for large s and two large for small s. The results show that the RMSE of the Kriging (0.490) is comparable to DeepKriging (0.485). But, in terms of the quantile loss, the AQTL of the Kriging (60.47) is much larger than that of DeepKriging (0.12), which further implies that Kriging is only a good estimation in terms of the mean squared error which reflects the mean and variance, but cannot reflect other important properties such as quantiles.

S6. Source codes and data
• The source codes can be found in the following Github repository: https://github.com/aleksada/DeepKriging.
• The methods are implemented in Keras, a library of both Python and R. Most of the codes for this work are written in Python. We also provide some simple examples of DeepKriging using Keras in R.
• The PM2.5 station data from AQS used in the application can be found in https://aqs.epa.gov/aqsweb/airdata/download files.html.
• The meteorological data from NARR used in the application can be found in https://psl.noaa.gov/data/gridded/data.narr.html.