Semiparametric Regression for Spatial Data via Deep Learning

For a spatial domain of interest $\mathcal{S}$, we consider the following semiparametric spatial regression model:

$$y(\bm{s})=f_{0}(\bm{x}(\bm{s}))+e_{1}(\bm{s})+e_{2}(\bm{s}),\,\,\,\bm{s}\in\mathcal{S}$$

(1)

where $f_{0}:[0,1]^{d}\rightarrow\mathbb{R}$ is an unknown mean function of interest, $\bm{x}(\bm{s})=(x_{1}(\bm{s}),\ldots,x_{d}(\bm{s}))^{\top}$ represents a $d$-dimensional vector of covariates at location $\bm{s}$ with $x_{i}(\bm{s})\in[0,1]$, $e_{1}(\bm{s})$ is a mean zero Gaussian random field with covariance function $\gamma({\bm{s},\bm{s}^{\prime}})$, $\bm{s},\bm{s}^{\prime}\in\mathcal{S}$, and $e_{2}(\bm{s})$ is a spatial Gaussian white noise process with mean 0 and variance $\sigma^{2}$. Furthermore, we assume that $e_{1}(\bm{s})$, $e_{2}(\bm{s})$, and $\bm{x}(\bm{s})$ are independent of each other. Thus the observation $y(\bm{s})$ comprises three components: large-scale trend $f_{0}(\bm{x}(\bm{s}))$, small-scale spatial variation $e_{1}(\bm{s})$, and measurement error $e_{2}(\bm{s})$; see, for instance, Cressie and Johannesson, (2008).

In the spatial statistics literature, it is popular to focus on predicting the hidden spatial process $y(\bm{s})^{*}=f_{0}(\bm{x}(\bm{s}))+e_{1}(\bm{s})$ using the observed information (Cressie and Johannesson,, 2008). However, the primary interest of this paper is to estimate the large-scale trend $f_{0}(\bm{x}(\bm{s}))$, where the relationship between the hidden spatial process and the covariates could be complex in nature. To capture such a complex relationship, we assume that $f_{0}$ is a composition of several functions inspired by neural networks characteristics (Schmidt-Hieber,, 2020). Hölder smoothness (see Definition 1 in Appendix) is a commonly used smoothness assumption for regression function in nonparametric and semiparametric literature. Thus it is natural to assume the true mean function $f_{0}$ is a composition of Hölder smooth functions, which is formally stated in the following assumption.

Without loss of generality, we assume $C_{i}>1$ in Assumption 1. The parameter $L_{*}$ refers to the total number of layers, i.e., the number of composite functions, $\bm{r}$ is the whole number of variables in each layer, whereas $\bm{\tilde{r}}$ is the number of “active” variables in each layer. The two parameter vectors $\bm{\beta}$ and $\bm{C}$ pertain to the Hölder smoothness in each layer, while $\bm{a}$ and $\bm{b}$ define the domain of $\bm{g}_{i}$ in the $i$th layer. For the rest of the paper, we will use $\mathcal{CS}(L_{*},\bm{r},\bm{\tilde{r}},\bm{\beta},\bm{a},\bm{b},\bm{C})$ to denote the class of functions that have a compositional structure as specified in Assumption 1 with parameters $(L_{*},\bm{r},\bm{\tilde{r}},\bm{\beta},\bm{a},\bm{b},\bm{C})$.

It is worth mentioning that Assumption 1 is commonly adopted in deep learning literature; see Bauer and Kohler, (2019), Schmidt-Hieber, (2020), Kohler and Langer, (2021), Liu et al., (2020), Li et al., (2021), among others. This compositional structure covers a wide range of function classes including the generalized additive model. In Example 1 and Figure 1, we present an illustrative example to help readers understand the concept of compositional structure within the context of the generalized additive model.

In this paper, we consider estimating the unknown function $f_{0}$ via a deep neural network owing to the complexity of $f_{0}$ and the flexibility of neural networks. So before presenting our main results, we first briefly review the neural network terminologies pertaining to this work.

An activation function is a nonlinear function used to learn the complex pattern from data. In this paper, we focus on the Rectified Linear Unit (ReLU) shifted activation function which is defined as $\sigma_{\bm{v}}(\bm{z})=(\sigma(z_{1}-v_{1}),\ldots,\sigma(z_{d}-v_{d}))^{\top}$, where $\sigma(s)=\max\{0,s\}$ and $\bm{z}=(z_{1},\ldots,z_{d})^{\top}\in\mathbb{R}^{d}$. ReLU activation function enjoys both theoretical and computational advantages. The projection property $\sigma\circ\sigma=\sigma$ can facilitate the proof of consistency, while ReLU activation function can help avoid vanishing gradient problem. The ReLU feedforward neural network $\bm{f}(\bm{z},W,v)$ is given by

$$\bm{f}(\bm{z},W,v)=W_{L}\sigma_{\bm{v}_{L}}\ldots W_{1}\sigma_{\bm{v}_{1}}W_{0}\bm{z},\quad\bm{z}\in\mathbb{R}^{p_{0}},$$

(2)

where $\{(W_{0},\ldots,W_{L}):W_{l}\in\mathbb{R}^{p_{l+1}\times p_{l}},0\leq l\leq L\}$ is the collection of weight matrices, $\{(\bm{v}_{1},\ldots,\bm{v}_{L}):\bm{v}_{l}\in\mathbb{R}^{p_{l}},1\leq l\leq L\}$ is the collection of so-called biases (in the neural network literature), and $\sigma_{\bm{v}_{l}}(\cdot)$, $1\leq l\leq L$, are the ReLU shifted activation functions. Here, $L$ measures the number of hidden layers, i.e., the length of the network, while $p_{j}$ is the number of units in each layer, i.e., the depth of the network. When using a ReLU feedforward neural network to estimate the regression problem (1), we need to have $p_{0}=d$ and $p_{L+1}=1$; and the parameters that need to be estimated are the weight matrices $(W_{j})_{j=0,\ldots,L}$ and the biases $(\bm{v}_{j})_{j=1,\ldots,L}$.

By definition, a ReLU feedforward neural network can be written as a composition of simple nonlinear functions; that is,

$$\bm{f}(\bm{z},W,v)=\bm{g}_{L}\circ\ldots\circ\bm{g}_{1}\circ\bm{g}_{0}(\bm{z}),\quad\bm{z}\in\mathbb{R}^{p_{0}},$$

where $\bm{g}_{i}(\bm{z})=W_{i}\sigma_{\bm{v}_{i}}(\bm{z})$, $\bm{z}\in\mathbb{R}^{p_{i}}$, $i=1,\ldots,L$, and $\bm{g}_{0}(\bm{z})=W_{0}\bm{z},\bm{z}\in\mathbb{R}^{p_{0}}$. Unlike traditional function approximation theory where a complex function is considered as an infinite sum of simpler functions (such as Tayler series, Fourier Series, Chebyshev approximation, etc.), deep neural networks approximate a complex function via compositions, i.e., approximating the function by compositing simpler functions (Lu et al.,, 2020; Farrell et al.,, 2021; Yarotsky,, 2017). Thus, a composite function can be well approximated by a feedforward neural network. That is why we assume the true mean function $f_{0}$ has a compositional structure.

In practice, the length and depth of the networks can be extremely large, thereby easily causing overfitting. To overcome this problem, a common practice in deep learning is to randomly set some neurons to zero, which is called dropout. Therefore, it is natural to assume the network space is sparse and all the parameters are bounded by one, where the latter can be achieved by dividing all the weights by the maximum weight (Bauer and Kohler,, 2019; Schmidt-Hieber,, 2020; Kohler and Langer,, 2021). As such, we consider the following sparse neural network class with bounded weights

$$\mathcal{F}(L,\bm{p},\tau,F)=\left\{f\textrm{ is of form }\eqref{fnn}:\max_{j=0,\ldots,L}\|W_{j}\|_{\infty}+|\bm{v}_{j}|_{\infty}\leq 1,\sum_{j=0}^{L}(\|W_{j}\|_{0}+|\bm{v}_{j}|_{0})\leq\tau,\|f\|_{\infty}\leq F\right\},$$

(3)

where $\bm{p}=(p_{0},\ldots,p_{L+1})$ with $p_{0}=d$ and $p_{L+1}=1$, and $\bm{v}_{0}$ is a vector of zeros. This class of neural networks is also adopted in Schmidt-Hieber, (2020), Liu et al., (2020), and Li et al., (2021).

Suppose that the process $y(\cdot)$ is observed at a finite number of spatial locations $\{\bm{s}_{1},\ldots,\bm{s}_{n}\}$ in $\mathcal{S}$. The desired DNN estimator of $f_{0}$ in Model (1) is a sparse neural network in $\mathcal{F}(L,\bm{p},\tau,F)$ with the smallest empirical risk; that is,

$$\widehat{f}_{\textrm{global}}(\widehat{W},\widehat{v})=\mathop{\mathrm{argmin}}_{f\in\mathcal{F}(L,\bm{p},\tau,F)}n^{-1}\sum_{i=1}^{n}(y(\bm{s}_{i})-f(\bm{x}(\bm{s}_{i}))^{2}.$$

(4)

For simplicity, we sometimes write $\widehat{f}_{\textrm{global}}$ for $\widehat{f}_{\textrm{global}}(\widehat{W},\widehat{v})$ if no confusion arises.

Because (4) does not have an exact solution, we use a stochastic gradient descent (SGD)-based algorithm to optimize (4). In contrast to a gradient descent algorithm which requires a full dataset to estimate gradients in each iteration, SGD or mini-batch gradient descent only needs an access to a subset of observations during each update, which is capable of training relatively complex models for large datasets and computationally feasible with streaming data. Albeit successful applications in machine learning and deep learning, SGD still suffers from some potential problems. For example, the rate of convergence to the minima is slow; the performance is very sensitive to tuning parameters. To circumvent these problems, various methods have been proposed, such as RMSprop and Adam (Kingma and Ba,, 2014). In this paper, we use Adam optimizer to solve (4).

During the training process, there are many hyper-parameters to tune in our approach: the number of layers $L$, the number of neurons in each layer $\bm{p}$, the sparse parameter $\tau$, and the learning rate. These hyper-parameters play an important role in the learning process. However, it is challenging to determine the values of hyper-parameters without domain knowledge. In particular, it is challenging to control the sparse parameter $\tau$ directly in the training process. Thus, we add an $\ell_{1}$-regularization penalty to control the number of inactive neurons in the network. The idea of adding a sparse regularization to hidden layers in deep learning is very common; see, for instance, Scardapane et al., (2017) and Lemhadri et al., (2021). In this paper, we use a $5$-fold cross-validation to select tuning parameters.

Recall that the minimizer of (4), $\widehat{f}_{\textrm{global}}$, is practically unattainable and we use an SGD-based algorithm to minimize the objective function (4), which may converge to a local minimum. The actual estimator obtained by minimizing (4) is denoted by $\widehat{f}_{\textrm{local}}\in\mathcal{F}(L,\bm{p},\tau,F)$. We define the difference between the expected empirical risks of $\widehat{f}_{\textrm{global}}$ and $\widehat{f}_{\textrm{local}}$ as

	$\displaystyle\Delta_{n}(\widehat{f}_{\textrm{local}})$	$\displaystyle\doteq\bm{E}_{f_{0}}\left[\frac{1}{n}\sum_{i=1}^{n}(y(\bm{s}_{i})-\widehat{f}_{\textrm{local}}(\bm{x}(\bm{s}_{i}))^{2}-\inf_{\tilde{f}\in\mathcal{F}(L,\bm{p},\tau,F)}\frac{1}{n}\sum_{i=1}^{n}(y(\bm{s}_{i})-\tilde{f}(\bm{x}(\bm{s}_{i}))^{2}\right]$
		$\displaystyle=\bm{E}_{f_{0}}\left[\frac{1}{n}\sum_{i=1}^{n}(y(\bm{s}_{i})-\widehat{f}_{\textrm{local}}(\bm{x}(\bm{s}_{i}))^{2}-\frac{1}{n}\sum_{i=1}^{n}(y(\bm{s}_{i})-\widehat{f}_{\textrm{global}}(\bm{x}(\bm{s}_{i}))^{2}\right],$		(5)

where $\bm{E}_{f_{0}}$ stands for the expectation with respect to the true regression function $f_{0}$. For any $\widehat{f}\in\mathcal{F}(L,\bm{p},\tau,F)$, we consider the following estimation error:

$$R_{n}(\widehat{f},f_{0})\doteq\bm{E}_{f_{0}}\left[\frac{1}{n}\sum_{i=1}^{n}\big{(}\widehat{f}(\bm{x}(\bm{s}_{i}))-f_{0}(\bm{x}(\bm{s}_{i}))\big{)}^{2}\right].$$

(6)

The oracle-type theorem below gives an upper bound for the estimation error.

The convergence rate in Theorem 1 is determined by three components. The first component $\inf_{\tilde{f}\in\mathcal{F}(L,\bm{p},\tau,F)}\|\tilde{f}-f_{0}\|^{2}_{\infty}$ measures the distance between the neural network class $\mathcal{F}(L,\bm{p},\tau,F)$ and $f_{0}$, i.e., the approximation error. The second term $\zeta_{n,\varepsilon,\delta}$ pertains to the estimation error, and $\Delta_{n}(\widehat{f})$ is owing to the difference between $\widehat{f}$ and the oracle neural network estimator $\widehat{f}_{\textrm{global}}$. It is worth noting that the upper bound in Theorem 1 does not depend on the network architecture parameter $\bm{p}$, i.e., the width of the network, in that the network is sparse and its “actual” width is controlled by the sparsity parameter $\tau$. To see this, after removing all the inactive neurons, it is straightforward to show that $\mathcal{F}(L,\bm{p},\tau,F)=\mathcal{F}(L,(p_{0},p_{1}\wedge\tau,\ldots,p_{L}\wedge\tau,p_{L+1}),\tau,F)$ (Schmidt-Hieber,, 2020).

Next, we turn to the consistency of the DNN estimator $\widehat{f}_{\textrm{local}}$ for $f_{0}\in\mathcal{CS}(L_{*},\bm{r},\bm{\tilde{r}},\bm{\beta},\bm{a},\bm{b},\bm{C})$. In nonparametric regression, the estimation convergence rate is heavily affected by the smoothness of the function. Consider the class of composite functions $\mathcal{CS}(L_{*},\bm{r},\bm{\tilde{r}},\bm{\beta},\bm{a},\bm{b},\bm{C})$. Let $\beta_{i}^{*}=\beta_{i}\prod_{s=i+1}^{L_{*}}(\beta_{s}\wedge 1)$ for $i=0,\ldots,L_{*}$ and $i^{*}=\mathop{\mathrm{argmin}}_{0\leq i\leq L_{*}}\beta_{i}^{*}/\tilde{r}_{i}$, with the convention $\prod_{s=L_{*}+1}^{L_{*}}(\beta_{s}\wedge 1)=1$. Then $\beta^{*}=\beta_{i^{*}}^{*}$ and $r^{*}=\tilde{r}_{i^{*}}$ are known as the intrinsic smoothness and intrinsic dimension of $f\in\mathcal{CS}(L_{*},\bm{r},\bm{\tilde{r}},\bm{\beta},\bm{a},\bm{b},\bm{C})$. Similar definitions could be found in literature Kurisu, (2019); Nakada and Imaizumi, (2020); Cloninger and Klock, (2021); Kohler et al., (2022); Wang et al., (2023). These quantities play an important role in controlling the convergence rate of the estimator. To better understand $\beta_{i}^{*}$ and $i^{*}$, think about the composite function from the $i$th to the last layer, i.e., $\bm{h}_{i}(\bm{z})=\bm{g}_{L_{*}}\circ\ldots\circ\bm{g}_{i+1}\circ\bm{g}_{i}(\bm{z}):[a_{i},b_{i}]^{r_{i}}\to\mathbb{R}$; then $\beta_{i}^{*}$ can be viewed as the smoothness of $\bm{h}_{i}$ and $i^{*}$ is the layer of the least smoothness after rescaled by the respective number of “active” variables $\tilde{r}_{i}$, $i=0,\ldots,L_{*}$. The following theorem establishes the consistency of $\widehat{f}_{\textrm{local}}$ as an estimator of $f_{0}$ and its convergence rate in the presence of spatial dependence.

The consistency of $\widehat{f}_{\textrm{local}}$ can be achieved by, for instance, letting $L\asymp\log(n),N\asymp n^{\frac{r^{*}}{2\beta^{*}+r^{*}}},\operatorname{tr}(\Gamma_{n}^{2})=o(n^{\frac{4\beta^{*}+r^{*}}{2\beta^{*}+r^{*}}}(\log n)^{-3})$, and $\Delta_{n}(\widehat{f}_{\textrm{local}})=o(1)$, as a result of which $\varsigma_{n}\asymp n^{-\frac{2\beta^{*}}{2\beta^{*}+r^{*}}}(\log n)^{3}+\Delta_{n}(\widehat{f}_{\textrm{local}})=o(1)$. As expected, the rate of convergence is affected by the intrinsic smoothness and intrinsic dimension of $\mathcal{CS}(L_{*},\bm{r},\bm{\tilde{r}},\bm{\beta},\bm{a},\bm{b},\bm{C})$, the architecture of the neural network $\mathcal{F}(L,\bm{p},\tau,F)$, and the magnitude of the spatial dependence.

We also compare the proposed DNN estimator with various estimators in the literature. The first estimator of $f_{0}$ is based on the Gaussian process-based spatially varying coefficient model (GP-SVC) which is given by

$$y(\bm{s})=\beta_{1}(\bm{s})x_{1}(\bm{s})+\ldots,+\beta_{p}(\bm{s})x_{p}(\bm{s})+\epsilon,\,\,\,\epsilon\sim\mathcal{N}(0,\tau^{2}),\,\,\,\bm{s}\in\mathcal{S},$$

and the spatially varying coefficient $\beta_{j}(\cdot)$ is the sum of a fixed effect and a random effect. That is, $\beta_{j}(\bm{s})=\mu_{j}+\eta_{j}(\bm{s})$, where $\mu_{j}$ is a non-random fixed effect and $\eta_{j}(\cdot)$ is a zero-mean Gaussian process with an isotropic covariance function $c(\cdot;\bm{\theta}_{j})$. In this work, we use the popular Matérn covariance function defined as

$$c\left(r;\rho,\nu,\sigma^{2}\right)=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{r}{\rho}\right)^{\nu}K_{\nu}\left(\sqrt{2\nu}\frac{r}{\rho}\right),$$

where $\rho>0$ is the range parameter, $\nu>0$ is the smoothness parameter, and $K_{\nu}(\cdot)$ is the modified Bessel function of second kind with order $\nu$.

The second estimator is the Nadaraya-Watson (N-W) kernel estimator for spatially dependent data discussed in Robinson, (2011), which considers the following spatial regression model

$$y(\bm{s}_{i})=f_{0}(\bm{x}(\bm{s}_{i}))+\sigma(\bm{x}(\bm{s}_{i}))V_{i},\,\,\,V_{i}=\sum_{j=1}^{\infty}a_{ij}\epsilon_{j},\,\,\,i=1,\ldots,n,$$

where $f_{0}(\bm{x}):[0,1]^{d}\rightarrow\mathbb{R}$ and $\sigma(\bm{x}):[0,1]^{d}\rightarrow[0,\infty)$ are the mean and variance functions, respectively, $\epsilon_{j}$ are independent random variables with zero mean and unit variance, and $\sum_{j=1}^{\infty}a_{ij}^{2}=1$. Robinson, (2011) introduces the following Nadaraya-Watson kernel estimator for $f_{0}$:

$$\hat{f}(\bm{x})=\frac{\hat{\nu}(\bm{x})}{\widehat{g}(\bm{x})},$$

where

$$\widehat{g}(\bm{x})=\frac{1}{nh_{n}^{d}}\sum_{i=1}^{n}K_{i}(\bm{x}),\quad\hat{\nu}(\bm{x})=\frac{1}{nh_{h}^{d}}\sum_{i=1}^{n}y_{i}K_{i}(\bm{x}),$$

with

$$K_{i}(\bm{x})=K\left(\frac{\bm{x}-\bm{x}(\bm{s}_{i})}{h_{n}}\right),$$

and $h_{n}$ is a scalar, positive bandwidth sequence satisfying $h_{n}\rightarrow 0$ as $n\rightarrow\infty$.

The third estimator of $f_{0}$ is based on the generalized additive model (GAM) mentioned in Example 1. That is, we assume that

$$f_{0}(\bm{x}(\bm{s}))=\Psi\left(\sum_{j=1}^{d}g_{j}(x_{j}(\bm{s}))\right),$$

where $g_{j}(\cdot):[0,1]\rightarrow\mathbb{R}$ and $\Psi(\cdot):\mathbb{R}\rightarrow\mathbb{R}$ are some smooth functions. In this model, spatial dependence is not assumed. The fourth estimator we consider is the linear regression (LR) model which assumes a linear relationship between the response variable and the predictor variables, without incorporating spatial dependence explicitly.

To evaluate the performance of each method, we generate additional $m=n/10$ observations at new locations, treated as a test set. Similar to Chu et al., (2014), we adopt mean squared estimation error (MSEE) and mean squared prediction error (MSPE) to evaluate the estimation and prediction performance, where MSEE and MSPE are defined as

$$\textrm{MSEE}=m^{-1}\sum_{i=1}^{m}(\widehat{f}(\bm{x}(s_{i}))-f_{0}(\bm{x}(s_{i})))^{2},\quad\textrm{MSPE}=m^{-1}\sum_{i=1}^{m}(\widehat{f}(\bm{x}(s_{i}))-y(s_{i}))^{2},$$

and $\widehat{f}(\bm{x}(s_{i}))$ is an estimator of $f_{0}(\bm{x}(s_{i}))$. The mean and standard deviation of MSEE and MSPE over the 100 independent replicates are summarized in Tables 1 – 4.

Tables 1 and 2 pertain to Simulation Design 1, for fixed and expanding domains, respectively. For each combination of the sample size $n$ and the spatial dependence $\rho$, we highlight the estimator in boldface that yields the smallest MSEE and MSPE. Overall, GAM, GP-SVC, N-W, and LR methods perform similar to each other. The proposed DNN estimator produces a smaller estimation error and prediction error than the others in all cases except when $n=200$ and $\rho=0.1$ in the fixed-domain case, GAM yields the smallest MSPE of 1.26. But the MSPE produced by DNN is close. Despite that spatial dependence has an adverse impact on the performance, when $n$ increases (and $D$ increases for the expanding-domain case), both estimation error and prediction error decrease as expected.

We depict in Figure 2 the estimated mean functions $\widehat{f}(\bm{x}(s))$ via our method with $n=100$ and $\rho=0.5$ from the 100 replications along with the $95\%$ pointwise confidence intervals for both fixed and expanding domains. Here, the $95\%$ pointwise intervals are defined as

$$\left(2^{-1}(\widehat{f}_{(2)}(\bm{x}(s_{i}))+\widehat{f}_{(3)}(\bm{x}(s_{i}))),2^{-1}(\widehat{f}_{(97)}(\bm{x}(s_{i}))+\widehat{f}_{(98)}(\bm{x}(s_{i})))\right),~{}i=1,2,\ldots,n,$$

where $\widehat{f}_{(k)}(\bm{x}(s_{i}))$ is the $k$th smallest value of $\{\widehat{f}_{[j]}(\bm{x}(s_{i})):j=1,\ldots,100\}$, and $\widehat{f}_{[j]}(\bm{x}(s_{i}))$ is the estimator of $f_{0}(\bm{x}(s_{i}))$ from the $j$th replicate.

Tables 3 and 4 report the results for Simulation Design 2. For both fixed and expanding domains, our method performs the best among the five methods and N-W comes next. This is mainly because LR, GAM and GP-SVC treat $f_{0}$ to be linear and cannot handle complex interactions and nonlinear structures in $f_{0}$.