Mean-dependent nonstationary spatial models

TL;DR

This paper introduces a mean-dependent nonstationary spatial model that simplifies covariance estimation, improves prediction in certain cases, and includes a test for nonstationarity, demonstrated on Puerto Rico precipitation data.

Contribution

The paper proposes a novel mean-dependent nonstationary spatial model that reduces complexity and enhances predictive performance compared to traditional models.

Findings

Model improves predictions over stationary models in specific scenarios.

Proposed approximation maintains predictive accuracy with lower computational cost.

The nonstationarity test reliably detects nonstationary spatial data.

Abstract

Nonstationarity is a major challenge in analyzing spatial data. For example, daily precipitation measurements may have increased variability and decreased spatial smoothness in areas with high mean rainfall. Common nonstationary covariance models introduce parameters specific to each location, giving a highly-parameterized model which is difficult to fit. We develop a nonstationary spatial model that uses the mean to determine the covariance in a region, resulting in a far simpler, albeit more specialized, model. We explore inferential and predictive properties of the model under various simulated data situations. We show that this model in certain circumstances improves predictions compared to a standard stationary spatial model. We further propose a computationally efficient approximation that has comparable predictive accuracy. We also develop a test for nonstationary data and show it reliably identifies nonstationarity. We apply these methods to daily precipitation in Puerto Rico.

Figures

Tables

Table 1. Table 1 : Estimation Accuracy for Simulated Data: Mean squared error (with standard errors) for covariance parameter estimates, averaged over all simulations. Values are totaled over the 6 elements in 𝜼 𝜼 \eta and multiplied by 1000 for clarity. See Appendix A for MSEs for individual elements.

Training Data			Models
True Covariance	Sites	Days	Stationary	One-step	Full MLE
Stationary	50	5	218.3 (1.0)	1241.6 (36.1)	1330.2 (34.5)
		10	198.1 (0.6)	311.9 (4.6)	312.4 (4.5)
	100	5	96.5 (0.4)	113.2 (0.5)	114.0 (0.5)
		10	87.0 (0.3)	95.6 (0.3)	95.2 (0.3)
	200	5	56.0 (0.3)	65.6 (0.3)	65.8 (0.3)
		10	52.5 (0.2)	56.3 (0.2)	56.2 (0.2)
Nonstationary	50	5	1628.9 (4.6)	683.7 (8.8)	1767.6 (14.6)
		10	1304.2 (2.5)	502.6 (1.2)	1735.7 (10.9)
	100	5	746.2 (2.1)	147.4 (0.7)	182.8 (0.9)
		10	541.7 (1.0)	99.0 (0.4)	108.9 (0.5)
	200	5	1133.2 (2.0)	78.4 (0.4)	69.6 (0.5)
		10	868.2 (0.9)	56.1 (0.2)	34.7 (0.2)

Table 2. Table 2 : Mean computation time (in seconds) for simulated data.

Training Data			Models
True Covariance	Sites	Days	Stationary	One-step	Full MLE
Stationary	50	5	8.9	18.7	81.3
		10	18.3	38.8	231.8
	100	5	27.8	61.2	172.4
		10	54.3	112.5	443.4
	200	5	132.4	264.3	686.3
		10	241.6	444.9	1442.0
Nonstationary	50	5	8.0	15.6	219.0
		10	15.7	27.5	1186.7
	100	5	25.3	71.2	531.6
		10	53.0	146.9	2919.3
	200	5	109.6	328.3	2128.6
		10	242.0	665.0	11589.0

Table 3. Table 3 : Evaluation of Prediction Distributions for Simulated Data: Percent improvement in mean squared errors (top) and differences in deviance (bottom) between the nonstationary and stationary models. Positive percentages indicate that the nonstationary model was better than the stationary model. Values in parentheses are standard errors.

Training Data			Models
True Covariance	Sites	Days	One-step	Full MLE
Stationary	50	5	-0.55 (0.12)	-0.82 (0.15)
		10	-0.22 (0.04)	-0.31 (0.07)
	100	5	-0.26 (0.05)	-0.27 (0.06)
		10	-0.14 (0.02)	-0.16 (0.03)
	200	5	-0.16 (0.04)	-0.17 (0.04)
		10	-0.07 (0.02)	-0.07 (0.02)
Nonstationary	50	5	1.04 (0.20)	1.24 (0.37)
		10	0.98 (0.13)	1.74 (0.24)
	100	5	1.12 (0.26)	2.44 (0.02)
		10	1.32 (0.17)	2.89 (0.23)
	200	5	1.54 (0.42)	2.28 (0.44)
		10	2.10 (0.24)	3.05 (0.27)

Table 4. Table 4 : Covariance parameter estimates and log-likelihoods for Puerto Rico precipitation data. The nugget and sill are fitted under stationary, mean-dependent nonstationarity (MDNS), or log-mean-dependent nonstationarity (L-MDNS) conditions. Range is fitted under stationary or mean-dependent nonstationarity (MDNS) conditions.

Model				Nugget		Sill		Range
$p$	Nugget/Sill	Range	Log-Likelihood	$a_1$	$b_1$	$a_2$	$b_2$	$a_3$	$b_3$
4	Stationary	Stationary	-1575.7	-0.91	0	-2.22	0	2.70	0
	MDNS	Stationary	-1436.1	-1.72	1.08	-5.24	-9.27	54.78	0
	MDNS	MDNS	-1431.9	-1.80	1.15	-4.26	-7.34	48.20	11.18
	LMDNS	Stationary	-1435.2	-2.56	1.95	-1.80	0.72	1.88	0
	LMDNS	MDNS	-1421.2	-2.60	1.99	-4.62	8.94	2.62	-12.15
7	Stationary	Stationary	-1580.1	-1.48	0	-1.52	0	-2.35	0
	MDNS	Stationary	-1464.1	-7.55	-11.09	-1.70	0.92	-1.80	0
	MDNS	MDNS	-1405.9	-17.79	-24.30	-1.36	0.65	5.93	-9.27
	LMDNS	Stationary	-1492.8	-2.36	1.45	-0.93	1.11	-1.64	0
	LMDNS	MDNS	-1452.1	-2.27	1.37	3.03	4.29	10.10	-11.27

Table 5. Table 5 : Five-fold cross-validation results for Puerto Rico precipitation data. For predictive accuracy, we calculate the score, mean-squared error, and coverage of the 95% prediction interval using the prediction distribution. The quantiles of the standard errors for the test sites illustrate how the nonstationary models compare to the stationary models in terms of precision.

Model			Prediction Accuracy			Standard Errors
$p$	Nugget/Sill	Range	Score	MSE	Coverage	$Q_{.05}$	$Q_{.50}$	$Q_{.95}$
4	Stationary	Stationary	-2.555	0.518	93.5%	0.69	0.72	0.72
	MDNS	Stationary	-0.335	0.505	96.0%	0.38	0.59	1.23
	MDNS	MDNS	-0.431	0.508	95.7%	0.39	0.58	1.27
	LMDNS	Stationary	-0.424	0.513	95.8%	0.35	0.62	1.12
	LMDNS	MDNS	-0.060	0.505	96.3%	0.33	0.62	1.13
7	Stationary	Stationary	-1.305	0.451	94.8%	0.64	0.67	0.69
	MDNS	Stationary	+0.388	0.450	96.2%	0.39	0.57	1.14
	MDNS	MDNS	+1.715	0.441	95.3%	0.22	0.61	1.02
	LMDNS	Stationary	-0.264	0.453	96.0%	0.37	0.62	1.00
	LMDNS	MDNS	+0.546	0.444	95.9%	0.34	0.62	1.00

Training Data			Models
True Covariance	Sites	Days	Stationary	One-step	Full MLE
Stationary	50	5	654.2 (2.8)	644.6 (2.8)	645.2 (2.9)
		10	496.1 (2.3)	516.3 (2.4)	499.2 (2.4)
	100	5	111.0 (0.7)	155.5 (1.0)	154.3 (1.0)
		10	67.9 (0.5)	99.1 (0.6)	98.9 (0.6)
	200	5	57.9 (0.4)	72.2 (0.5)	72.2 (0.5)
		10	33.6 (0.2)	37.7 (0.3)	38.1 (0.3)
Nonstationary	50	5	1412.4 (8.3)	1278.2 (1.6)	1198.3 (1.9)
		10	1313.1 (1.0)	1294.8 (1.3)	1189.6 (1.7)
	100	5	1454.8 (9.6)	210.0 (1.5)	199.3 (1.2)
		10	916.7 (5.5)	111.3 (0.7)	103.3 (0.7)
	200	5	313.3 (2.2)	92.1 (0.7)	70.2 (0.6)
		10	167.9 (1.5)	47.1 (0.3)	29.1 (0.2)

Training Data			Models
True Covariance	Sites	Days	Stationary	One-step	Full MLE
Stationary	50	5	1024.5 (7.9)	1040.9 (8.0)	1077.9 (8.3)
		10	702.1 (4.9)	734.7 (5.5)	780.5 (5.9)
	100	5	699.5 (3.8)	705.7 (3.8)	701.9 (3.8)
		10	590.7 (3.0)	568.0 (3.2)	560.1 (3.1)
	200	5	416.5 (2.5)	453.6 (2.8)	450.9 (2.8)
		10	421.9 (2.0)	441.5 (2.1)	434.4 (2.1)
Nonstationary	50	5	56032.7 (226.5)	2503.7 (15.0)	1003.8 (8.3)
		10	39556.9 (96.0)	1841.9 (10.2)	643.5 (4.4)
	100	5	26599.4 (116.9)	1156.4 (14.1)	829.4 (4.8)
		10	16148.4 (54.2)	526.3 (4.9)	754.9 (4.4)
	200	5	44648.2 (122.4)	746.8 (7.8)	772.3 (4.4)
		10	30649.4 (54.4)	319.9 (4.0)	493.0 (3.0)

Training Data			Models
True Covariance	Sites	Days	Stationary	One-step	Full MLE
Stationary	50	5	11421.5 (57.5)	12167.9 (65.3)	12099.2 (67.2)
		10	10688.2 (35.1)	10788.3 (35.3)	10582.3 (35.5)
	100	5	4981.9 (24.7)	5290.1 (271)	5182.4 (27.1)
		10	4561.4 (18.8)	4702.6 (19.3)	4627.8 (19.4)
	200	5	2887.9 (14.2)	3108.9 (151)	3055.3 (15.0)
		10	2694.4 (9.3)	2756.5 (10.2)	2724.6 (10.3)
Nonstationary	50	5	35189.6 (165.6)	20653.4 (92.7)	27312.4 (138.8)
		10	32284.0 (116.1)	20253.4 (63.2)	20338.2 (69.2)
	100	5	11614.9 (50.7)	5166.4 (32.1)	6438.8 (37.8)
		10	10337.6 (36.4)	3954.5 (20.0)	3939.0 (22.9)
	200	5	17930.9 (39.4)	3291.7 (18.8)	2589.1 (20.3)
		10	16175.1 (29.4)	2617.6 (11.4)	1261.9 (8.9)

Equations

$${\mu }_t(\mathbf{s})=\sum _{j=1}^J{Z}_j(\mathbf{s}){\beta }_{t,j},$$

$$\mathbb{C}\text{ov}[Y_{t}(\mathbf{s}_{i}),Y_{t}(\mathbf{s}_{j})]=\mbox{\boldmath$\Sigma$}_{t,ij}=\tau_{t}^{2}(\mathbf{s}_{i})\cdot\mathbb{I}(i=j)+\sigma_{t}(\mathbf{s}_{i})\;\sigma_{t}(\mathbf{s}_{j})\;\mathbf{R}_{t}^{NS}(\mathbf{s}_{i},\mathbf{s}_{j}),$$

$${\mathbf{R}}_t^{NS}(i,j)=|{\mathbf{C}}_{t,i}{|}^{1/4}|{\mathbf{C}}_{t,j}{|}^{1/4}{\left|\frac{{\mathbf{C}}_{t,i}+{\mathbf{C}}_{t,j}}{2}\right|}^{-1/2}{\mathbf{R}}^S(\sqrt{Q_{t,ij}})$$

$$Q_{t,ij}=({\mathbf{s}}_i-{\mathbf{s}}_j{)}^T{\left(\frac{{\mathbf{C}}_{t,i}+{\mathbf{C}}_{t,j}}{2}\right)}^{-1}({\mathbf{s}}_i-{\mathbf{s}}_j)$$

$${\mathbf{R}}_t^{NS}(i,j)={\left\{\frac{4{\rho }_t({\mathbf{s}}_i){\rho }_t({\mathbf{s}}_j)}{[{\rho }_t({\mathbf{s}}_i)+{\rho }_t({\mathbf{s}}_j){]}^2}\right\}}^{p/4}{\mathbf{R}}^S\left\{h_{ij}/\sqrt{\left[{\rho }_t({\mathbf{s}}_i)+{\rho }_t({\mathbf{s}}_j)\right]/2}\right\}$$

$${\mathbf{R}}_t^{NS}(i,j)={\mathbf{R}}^S\left(h_{ij}/\sqrt{{\rho }_t(\mathbf{s})}\right),$$

$$\tau_{t}^{2}(\mathbf{s})=g_{1}[\mu_{t}(\mathbf{s});\mbox{\boldmath$\eta$}];\;\sigma_{t}(\mathbf{s})=g_{2}[\mu_{t}(\mathbf{s});\mbox{\boldmath$\eta$}];\;\rho_{t}(\mathbf{s})=g_{3}[\mu_{t}(\mathbf{s});\mbox{\boldmath$\eta$}]$$

$$L(\mbox{\boldmath$\eta$};\mbox{\boldmath$\beta$}_{1},\dots,\mbox{\boldmath$\beta$}_{m})=\prod_{t=1}^{m}\phi[\mathbf{Y}_{t}|\mbox{\boldmath$\beta$}_{t},\mbox{\boldmath$\Sigma$}(\mbox{\boldmath$\beta$}_{t},\mbox{\boldmath$\eta$})]\cdot\phi(\mbox{\boldmath$\beta$}_{t}|\mbox{\boldmath$\beta$}_{0},\mbox{\boldmath$\Omega$})$$

$$\hat{\mbox{\boldmath$\beta$}}_{t}(\hat{\mbox{\boldmath$\eta$}}_{0})=\left(\mathbf{Z}^{T}\mbox{\boldmath$\Sigma$}(\hat{\mbox{\boldmath$\eta$}}_{0})^{-1}\mathbf{Z}+\mbox{\boldmath$\Omega$}^{-1}\right)^{-1}\mathbf{Z}^{T}\mbox{\boldmath$\Sigma$}(\hat{\mbox{\boldmath$\eta$}}_{0})^{-1}\mathbf{Y}_{t}.$$

$$H_0:$$

$$H_A:$$

$${\chi }^2=\sum _{t=1}^m$$

$$\begin{pmatrix}\mathbf{Y}_{0}\\ \mathbf{Y}_{1}\end{pmatrix}\sim\mbox{Normal}\left[\begin{pmatrix}\mbox{\boldmath$\mu$}_{0}\\ \mbox{\boldmath$\mu$}_{1}\end{pmatrix},\begin{pmatrix}\mbox{\boldmath$\Sigma$}_{00}&\mbox{\boldmath$\Sigma$}_{01}\\ \mbox{\boldmath$\Sigma$}_{10}&\mbox{\boldmath$\Sigma$}_{11}\end{pmatrix}\right]$$

$$(\mathbf{Y}_{0}|\mathbf{Y}_{1},\hat{\mbox{\boldmath$\beta$}},\hat{\mbox{\boldmath$\eta$}})\sim\mbox{Normal}\left(\mbox{\boldmath$\mu$}_{0}+\mbox{\boldmath$\Sigma$}_{01}\mbox{\boldmath$\Sigma$}_{11}^{-1}(\mathbf{Y}_{1}-\mbox{\boldmath$\mu$}_{1}),\mbox{\boldmath$\Sigma$}_{00}-\mbox{\boldmath$\Sigma$}_{01}\mbox{\boldmath$\Sigma$}_{11}^{-1}\mbox{\boldmath$\Sigma$}_{10}\right),$$

$$\text{\% Improvement in MSE}=100\cdot \frac{MSE(\text{Stationary})-MSE(\text{Nonstationary})}{MSE(\text{Stationary})}$$

Full text

Key words: spatial statistics, nonstationary, continuous regimes, environmental statistics, precipitation, hypothesis testing

1 Introduction

Modeling daily precipitation poses unique challenges to typical spatio-temporal models. Sites experiencing higher amounts of precipitation often exhibit larger variability than sites experiencing low amounts of precipitation, as seen in Figure 1. Furthermore, sites in regions experiencing similar weather patterns will be more correlated than sites experiencing different weather patterns, so the spatial covariances will change over the spatial domain. Given that the most widely used long-term daily precipitation time series do not account for this, and as the regions that experience high precipitation changes daily, we need a flexible model that accounts for changes in the daily mean and spatial covariance structures. The mean structure for a given day can be modeled using specified covariates, but we need a model that allows the spatial covariances between sites to shift both from day to day and from region to region.

Unlike stationary models, where the covariance structure is consistent across the spatial domain, nonstationary models allow for different levels of variability and spatial smoothness within the same domain. Nonstationary models have been developed through multiple approaches, which are summarized in Sampson [2010]. Typical techniques for nonstationary models are spatial deformation methods [Anderes and Stein, 2008; Sampson and Guttorp, 1992; Schmidt and O’Hagan, 2003], basis function models [Holland et al., 1999; Nychka and Saltzman, 1998; Nychka et al., 2002], kernel based methods [Fuentes, 2001; Fuentes and Smith, 2001; Higdon, 1998; Reich et al., 2013], and spectral methods [Fuentes, 2002; Guinness and Fuentes, 2015].

Regime-switching models are another class of the nonstationary spatio-temporal covariance models. Regime-switching models use a latent process, called the spatial regime, to determine the mean and spatial covariance of each site in the domain [Paciorek and Schervish, 2006]. Discrete regime models are common in hidden Markov models [Baum and Petrie, 1966; Barbu and Limnios, 2006; Leroux, 1992; Lindgren, 1978], econometrics [Gray, 1996; Zhou and Yin, 2003], wind forecasting [Gneiting et al., 2006; Pinson and Madsen, 2012; Neto et al., 2014], and wind time series [Shamshad et al., 2005; Kazor and Hering, 2015]. In these cases, the underlying physics of the process switches between specific modes, such as bull versus bear markets or directional shifts of major wind currents. However, discrete regimes are not always flexible enough to account for geographic and temporal shifts in the size and shape of regimes.

Instead of discrete regimes, continuous regime processes allow for more complex patterns in the nonstationary spatial model. Several authors have already extended discrete spatial regimes into continuously-indexed spatial regimes. Kleiber and Nychka [2012] extended the method for multivariate spatial processes, although without considering detailed modeling of the latent regime process. Modeling the regime process using covariates have been used extensively [Calder, 2008; Reich et al., 2011; Schmidt et al., 2011; Risser and Calder, 2015]. These covariate-driven models create continuously-varying regimes that directly inform the spatial covariance structure. In particular, temporal covariates can flexibly account for each day’s weather, whether it is clear skies, light rain, or even hurricane-level storms.

This paper builds on continuous regime-switching models but deviates from other approaches by indexing the regimes using the mean process. That is, we allow a separate spatial covariance for generally dry regions versus generally wet regions, but assume that all dry spatiotemporal regions have the same spatial covariance. While the proposed method cannot capture all forms of nonstationarity, it is useful in the specific case where mean/covariance relationship are thought to be critical. The data and the resulting spatial model is then used to create a predictive distribution across the entire spatial domain with appropriately varying standard errors. We discuss computation of both stationary and nonstationary models within this framework and develop a test statistic for significant departures from stationarity. We illustrate the model’s predictive capabilities using daily precipitation data from Puerto Rico and confirm that as the mean precipitation in an area on a given day increases, the variability increases and spatial smoothness decreases.

2 Nonstationary spatial model from conditionally stationary regimes

Let $Y_{t}(\mathbf{s})$ be the response at spatial location $\mathbf{s}$ and time $t$. The data are observed at locations $\mathbf{s}_{1},\dots,\mathbf{s}_{n}$ and days $t_{1},\dots,t_{m}$. The response is decomposed as $Y_{t}(\mathbf{s})=\mu_{t}(\mathbf{s})+e_{t}(\mathbf{s})$, where $\mu_{t}(\cdot)$ is the mean process and $e_{t}(\cdot)$ is, conditionally on $\mu_{t}(\cdot)$, a zero-mean Gaussian error process. The mean function is modeled as a linear combination of $J$ known covariates $Z_{1}(\mathbf{s}),\dots,Z_{J}(\mathbf{s})$:

[TABLE]

where $\beta_{t,j}$ is the time-varying effect of covariate $Z_{j}$ on day $t$. Potential covariates can be geographic variables such as elevation, spatiotemporal variables such as temperature, or polynomial functions of $\mathbf{s}$ and $t$ to model large-scale trends. The coefficient vectors \mbox{\boldmath\beta}_{t}=[\beta_{t,1},\dots,\beta_{t,J}]^{T}\overset{iid}{\sim}\mbox{Normal}(\mbox{\boldmath\beta}_{0},\mbox{\boldmath\Omega}) are modeled as random effects with mean \mbox{\boldmath\beta}_{0} and covariance $\Omega$ since each day has a different weather pattern. Conditioned on the mean, the errors $\mathbf{e}_{t}=[e_{t}(\mathbf{s}_{1}),\dots,e_{t}(\mathbf{s}_{n})]^{T}$ are assumed to be independent \mbox{Normal}(\vec{0},\mbox{\boldmath\Sigma}_{t}) with a spatial covariance matrix \mbox{\boldmath\Sigma}_{t}, that will depend on $\mu_{t}$ as described below. Note we are assuming conditional normality but not marginal normality, except for trivial cases where the covariance is not dependent on the mean. Furthermore we do not assume temporal covariances between the \mbox{\boldmath\beta}_{t} covariates to focus on the spatial model, though inclusion would likely help estimation and forecasting the mean process.

The error process can have different variances and correlations across multiple spatial regimes. While there are many approaches to index the regimes (typically based on large-scale geographic features), we index the regimes by the mean process. Using the mean to determine the regime allows the covariances to vary smoothly over space and time in a simple but effective way. Let \mbox{\boldmath\Sigma}_{t} be the spatial covariance matrix on day $t$, such that

[TABLE]

where $\tau_{t}^{2}(\mathbf{s}_{i})$ is the spatially-varying nugget effect, $\sigma^{2}_{t}(\mathbf{s}_{i})$ is the spatially-varying partial sill, and $\mathbf{R}_{t}^{NS}(\mathbf{s}_{i},\mathbf{s}_{j})$ is a nonstationary spatial correlation function.

For the nonstationary correlation matrix $\mathbf{R}_{t}^{NS}$, we use the method developed in Paciorek and Schervish [2006] to construct a positive-definite, nonstationary correlation function using local kernel matrices. Define $\mathbf{C}_{t,i}$ to be a positive-definite local kernel matrix for the correlation between $\mathbf{s}_{i}$ and nearby locations. For locations $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$, their nonstationary correlation is determined by their local kernel matrices:

[TABLE]

where $\mathbf{R}^{S}(\cdot)$ is a stationary correlation function of distance, and

[TABLE]

is the nonstationary extension of the Mahalanobis distance. We assume that the local kernel matrix is isotropic, so that $\mathbf{C}_{t,i}=\rho_{t}(\mathbf{s}_{i})\mathbf{I}_{2}$ and

[TABLE]

where $h_{ij}=||\mathbf{s}_{i}-\mathbf{s}_{j}||$. Note that if the regimes are equivalent, $\mu_{t}(\mathbf{s}_{i})=\mu_{t}(\mathbf{s}_{j})$, then

[TABLE]

implying that $\rho_{t}(\cdot)$ determines the range for the correlation within a regime. The Matèrn class of functions is a typical choice for $\mathbf{R}^{S}(r)$, of which we use the exponential function $\mathbf{R}^{s}(r)=\exp\{-r\}$ as a simple special case.

For spatially-varying nugget, partial sill, and range, we elect to use mean-dependent link functions

[TABLE]

that depend upon a set of unknown parameters $\eta$. Other methods of obtaining nonstationarity would be to estimate parameters through the covariates [Risser and Calder, 2015], but we present here a flexible ableit highly specialized model where the mean is the random spatial process varying across days to determine the degree of nonstationarity. While any nonnegative link function would suffice, we elect to use an exponential link function. The link function does not need to be monotonic, and much more sophisticated link functions could be used. For the exposition and simulation study, we assume the link is linear in the mean, such that g_{i}[\mu_{t}(\mathbf{s});\mbox{\boldmath\eta}]=\exp[a_{i}+b_{i}\mu_{t}(\mathbf{s})]. For the data analysis, we also explore the link that is linear in the log-mean g_{i}[\mu_{t}(\mathbf{s});\mbox{\boldmath\eta}]=\exp[a_{i}+b_{i}\log\mu_{t}(\mathbf{s})], which leads to a better fit.

In both cases, there are six parameters \mbox{\boldmath\eta}=(a_{1},b_{1},a_{2},b_{2},a_{3},b_{3}) to be estimated. The $a_{i}$ parameters affect the covariance parameters where either at the baseline mean equal to zero or (if $b_{i}=0$) under stationary assumptions. The $b_{i}$ parameters indicate how the mean affects the local covariance function, and because the mean is the only spatially varying component of the covariance function, the $b_{i}$ parameters control the degree of nonstationarity in the model. The $(a_{1},b_{1})$ parameters affect the nugget of the covariance function, the $(a_{2},b_{2})$ parameters affect the sill, and the $(a_{3},b_{3})$ parameters affect the range of correlated sites. However, the entire set of parameters should be considered when interpreting how the local covariance function changes with respect to the mean.

3 Computation

In vector form, \mathbf{Y}_{t}=\mbox{\boldmath\mu}_{t}+\mathbf{e}_{t}=\mathbf{Z}\mbox{\boldmath\beta}_{t}+\mathbf{e}_{t}, where the design matrix $\mathbf{Z}$ contains the covariates $\mathbf{Z}_{ij}=Z_{j}(\mathbf{s}_{i})$. In our representaion, the covariances change over time only via \mbox{\boldmath\beta}_{t}, so we denote the covariance matrix \mathbb{C}\text{ov}(\mathbf{e}_{t}|\mbox{\boldmath\beta}_{t})=\mbox{\boldmath\Sigma}(\mbox{\boldmath\beta}_{t},\mbox{\boldmath\eta}). The penalized likelihood function is then

[TABLE]

where $\phi$ is the multivariate normal density function. While construction of this likelihood function is well-reasoned, we acknowledge that we cannot guarantee the likelihood is not multimodal as parameters are only weakly identifiable, particularly when the nonstationary spatial covariances are relatively small. As a result, we carefully restrict parameter estimation in this study through profile maximum likelihood estimation methods, as described below. Bayesian methods could be explored to address uncertainty in the parameter estimations, but maximum likelihood is chosen because of the computational costs.

3.1 The Stationary Model

When $b_{1}=b_{2}=b_{3}=0$, the covariance function in Equation (3) loses dependence on the mean and becomes stationary. For the stationary model, we have \mbox{\boldmath\eta}\equiv\mbox{\boldmath\eta}_{0}=(a_{1},0,a_{2},0,a_{3},0) and \mathbb{C}\text{ov}(\mathbf{Y}_{t}|\mbox{\boldmath\beta}_{t})\equiv\mbox{\boldmath\Sigma}(\mbox{\boldmath\eta}_{0}) for all $t$. Furthermore, conditioned on estimated covariance parameters \hat{\mbox{\boldmath\eta}}_{0}, the maximum penalized likelihood estimate of \mbox{\boldmath\beta}_{t} is

[TABLE]

To estimate $a_{1}$, $a_{2}$, and $a_{3}$, we numerically maximize (8) with $b_{1}$, $b_{2}$, and $b_{3}$ set to 0 and iterate between updating \mbox{\boldmath\eta}_{0} and \mbox{\boldmath\beta}_{t} until the likelihood converges to its maximum.

3.2 The Nonstationary Model

For the full nonstationary model without the restriction that $b_{1}=b_{2}=b_{3}=0$, the optimization is challenging because \mbox{\boldmath\beta}_{t} appears in the covariance and a closed-form solution like (9) is no longer available. However, the estimator in (9) provides a reasonable approximation and propose a computationally-efficient one-step process. The one-step approximation fixes \mbox{\boldmath\beta}_{t} equal to the estimated \hat{\mbox{\boldmath\beta}}_{t} from the stationary model in Section 3.1 and optimizes (8) with respect to \mbox{\boldmath\eta}=(a_{1},b_{1},a_{2},b_{2},a_{3},b_{3}). We choose not to iterate between estimating \mbox{\boldmath\beta}_{t} and $\eta$, making this method a one-step approximation to the profile maximum likelihood estimate.

We found that multiple iterations of refitting the mean given the previously estimated covariance matrix did not guarantee convergence in the estimates, tended to overfit the mean process, and negatively affected predictive accuracy. However, the ultimate goal is to compare the stationary and nonstationary models to test for significant differences. Both the stationary and the one-step approximation have the same mean process but different spatial covariance models that can be compared to determine if there exists significant nonstationarity. Conversely, multiple iterations would lead to different mean processes and covariances, so the increased degrees of freedom make testing for nonstationarity more complicated.

To obtain the full maximum likelihood estimates for the nonstationary model, we also numerically optimize the likelihood in (8) over the full (\mbox{\boldmath\eta};\mbox{\boldmath\beta}_{1},\dots,\mbox{\boldmath\beta}_{m}) parameter set. Numerically maximizing over all of the \mbox{\boldmath\beta}_{t} coefficients is time-consuming, especially as the number of days increases, so we use the one-step approximation estimates as starting values to achieve convergence within fewer steps. In practice, this full maximization would be performed after deeming the one-step approximation is deemed significantly different from the stationary model.

4 Evaluation of models

Once models are fit to the data, we evaluate the nonstationary model to determine if it significantly improves upon the stationary model. First, we derive a hypothesis test to determine if the mean-dependent terms of the nonstationary model are significant. Then, we obtain the prediction distribution that uses the fitted model to estimate the response over the spatial domain.

4.1 Test for Nonstationarity

Under the restriction $b_{1}=b_{2}=b_{3}=0$, the spatial covariance no longer depends on the mean, leading to a stationary covariance model. We can hence test the hypotheses

[TABLE]

using a likelihood ratio test. Let \hat{\mbox{\boldmath\eta}}_{0} and \hat{\mbox{\boldmath\beta}}_{t(0)} be the parameter estimates under stationarity, and \hat{\mbox{\boldmath\eta}}_{1} and \hat{\mbox{\boldmath\beta}}_{t(1)} be the parameter estimates under nonstationarity. The Wilks’ test statistic \chi^{2}=2\cdot[\log L(\hat{\mbox{\boldmath\eta}}_{1};\hat{\mbox{\boldmath\beta}}_{1(1)},\dots,\hat{\mbox{\boldmath\beta}}_{m(1)})-\log L(\hat{\mbox{\boldmath\eta}}_{0};\hat{\mbox{\boldmath\beta}}_{1(0)},\dots,\hat{\mbox{\boldmath\beta}}_{m(0)})] is proportional to the difference in the log-likelihood:

[TABLE]

Under $H_{0}$, the Wilks’ test statistic [Wilks, 1938] should follow a Chi-squared distribution with degrees of freedom equal to the difference in dimensionality of the parameter space between the full and reduced model. Since stationarity occurs by restricting three parameters, the difference in dimensionality would be three, so we would reject the hypothesis that the data are stationary at the $\alpha$ significance level if the test statistic exceeds the $\alpha$-level critical value for the Chi-squared distribution with $df=3$ degrees of freedom. This method could then be generalized to determine the significance of a wide variety of nested models, depending on the defined link functions in (7) and which parameters are restricted.

In theory, the test statistic converges in distribution to the Chi-squared with $df=3$ degrees of freedom as the number of data points approaches infinity. In practice, the convergence rate may be too slow to assume the test statistic follows the Chi-squared. Furthermore, the inherent connection between the covariance and the mean indicate that the computation method from Section LABEL:label:sec:computation could affect the degrees of freedom of the test statistic. Under the one-step approximation, \mbox{\boldmath\Sigma}_{t(0)} and \mbox{\boldmath\Sigma}_{t(1)} would be the same, but under the full MLE computation, they would be different. Hence, the test statistic would most closely follow the Chi-squared distribution with $df=3$ degrees of freedom when using the one-step approximation method to estimate the covariance parameters. In Section 5.2, we explore the convergence of the Type I error and power for this test statistic using both computation methods and increasing sample sizes.

4.2 Prediction Distribution

Let $\mathbf{Y}_{0}$ be the unknown values of the responses at a set of prediction sites, and $\mathbf{Y}_{1}$ be the known values of the responses at a set of training sites for a given day. Following the general format of spatial kriging, the training and prediction sites are stacked into a joint multivariate normal distribution:

[TABLE]

where \mbox{\boldmath\mu}_{0}=\mathbf{Z}_{0}\mbox{\boldmath\beta} and \mbox{\boldmath\Sigma}_{00}=\mbox{\boldmath\Sigma}(\mbox{\boldmath\mu}_{0},\mbox{\boldmath\eta}) are the marginal mean and covariance matrix for the prediction sites, \mbox{\boldmath\mu}_{1}=\mathbf{Z}_{1}\mbox{\boldmath\beta} and \mbox{\boldmath\Sigma}_{11}=\mbox{\boldmath\Sigma}(\mbox{\boldmath\mu}_{1},\mbox{\boldmath\eta}) are the marginal mean and covariance for the training sites, and \mbox{\boldmath\Sigma}_{01} and \mbox{\boldmath\Sigma}_{10} are the cross-covariance matrices between the prediction and training sites.

After the training data is used to estimate the mean function coefficients \hat{\mbox{\boldmath\beta}} and covariance parameters \hat{\mbox{\boldmath\eta}}, we calculate the predicted responses $\hat{\mathbf{Y}}_{0}$ and their standard errors from the mean and variances in the conditional distribution of the prediction responses given the training responses:

[TABLE]

In (12), the prediction distribution is evaluated on a given day $t$, yet the notation has been suppressed for brevity. The prediction mean in (12) is our best prediction for the responses, with standard errors calculated from the square root of the diagonals of the covariance matrix. Note that training responses could be forecasted values using data from previous days at fixed sites, but the covariance matrix will need to incorporate temporal or autoregressive terms.

5 Simulation study

In this simulation study, we generate spatial data from the model, estimate parameters, and evaluate the predictive distribution. The purpose of the study is to compare the computational algorithms and understand the frequentist properties of the methods, with a type of stationary data as a baseline comparison. We show that we can estimate the covariance parameters accurately, detect nonstationarity, and predict more accurately than the stationary models.

5.1 Simulated data generation

The simulated data are split into a training set to fit the model and a testing set to evaluate predictions. The training data contains $n$ spatial locations with observations on $m$ days. The testing set has 100 randomly-selected locations for each day. All training and testing locations are uniformly distributed over a rectangular domain with longitudes in the range $(-67.3,-65.7)$ and latitudes in the range $(17.9,18.5)$, which is approximately the geographic area of Puerto Rico used in Section 6 (see Figure 5). We use $J=3$ covariates for the mean, using an intercept and spatial coordinates $\mathbf{s}=(s_{1},s_{2})^{T}$. When generating data, the coefficients \mbox{\boldmath\beta}_{t}=(\beta_{t1},\beta_{t2},\beta_{t3}) are fixed such that the intercept is constant, $\beta_{t1}=1$, and the spatial effects increase over time, $\beta_{t2}=\beta_{t3}=2\frac{t-1}{m-1}$.

Given the \mbox{\boldmath\beta}_{t} values for each day, the responses are generated using the covariance function in (2). The mean-independent covariance parameters remain fixed at $a_{1}=-1.0$, $a_{2}=0.5$, and $a_{3}=4.0$. For nonstationary data, we set the mean-dependent parameters to be $b_{1}=0.1$, $b_{2}=0.5$, and $b_{3}=-0.5$; for stationary data, these parameters are set to zero. Figure 6 illustrates how the simulations generate stationarity and nonstationarity under the chosen parameters. We consider both stationary and nonstationary data sets with $n=$ 50, 100 or 200 locations and $m=$ 5 or 10 days. For each set of experimental conditions, we create 200 replicates of the data simulations.

The parameters for the stationary and nonstationary models are estimated using the methods described in Section 3. Since there are only $J=3$ parameters and the true \mbox{\boldmath\beta}_{t} values are fixed for each simulation, we set \mbox{\boldmath\beta}_{0}=\vec{0} and \mbox{\boldmath\Omega}=e^{10}\mathbf{I}_{3} to effectively remove the random effects from the penalized likelihood function.

5.2 Type I error and power rate for nonstationarity test

We explore how effectively the test described in Section 4.1 determines a data set is nonstationary. Data sets are generated for a fixed $m=5$ days, with a varying number of spatial locations $n$. For each data set, we fit the data generated from the stationary model (i.e. $b_{1}=b_{2}=b_{3}=0$) and use the Chi-squared test statistic in (10) to test for non-stationarity at the 5% significance level. The number of days is fixed at $m=5$, and the number of spatial locations $n$ is allowed to vary.

To evaluate Type I error in Figure 7(a), data is generated from a stationary model, with $b_{1}=b_{2}=b_{3}=0$. The one-step approximation converges to the expected 5% error rate quicker than the full MLE. The full MLE generally has a larger Type I error rate than the expected 5%, which indicates that the test will favor using the nonstationary model even when the data are stationary.

To evaluate power in Figure 7(b), data is generated from a range of nonstationarity, with $b_{1}=0.1*c$, $b_{2}=0.5*c$, and $b_{3}=-0.5*c$ and $c$ varied from 0.1 to 1.0. The power increases as the number of spatial locations increases from 50 to 100, and the one-step approximation typically has lower power than the full MLE.

5.3 Parameter estimation accuracy

Table 1 contains the mean squared error (MSE) for the parameter estimates summed over all of the covariance parameters; the MSE for each individual $\eta$ parameter appears in Appendix A. As expected, for the stationary data, the stationary model is the most accurate, and for the nonstationary data, the nonstationary models are most accurate. The full MLE is the most accurate for all nonstationary cases, except for the smallest sample size. Table 2 contains the average times needed to fit the model. Between these tables we can see that the one-step approximation is faster than the full MLE while still being competitively accurate.

5.4 Prediction distribution accuracy

To measure prediction accuracy, we use prediction mean squared error, averaged over 100 prediction locations. The MSE is computed as the average squared difference between the testing data at the prediction locations (generated concurrently with the training data, as described in Section 5.1) and the predicted mean in (12) obtained from each model for each simulation. We calculate the percent improvement in MSE between each nonstationary model and the stationary model

[TABLE]

averaged over the simulations. In addition to MSE for point prediction, we evaluate the deviance of the test data from the negative log-likelihood of the prediction distribution and calculate the differences between each nonstationary model and the stationary model, averaged over all of the simulations. For both measures, positive differences indicate that the nonstationary model yielded better predictions of the test data.

Table 3 contains the percent improvement in MSE and log-likelihood for the simulated data. We observe that the nonstationary model fitted using full MLE is the best predictive model, but the one-step approximation is competitive. For the stationary data, the predictions from the nonstationary models are only slightly less accurate than the simpler stationary model (most likely due to overfitting to the training data), but the accuracy of the nonstationary data becomes comparable as the sample size increases. For the nonstationary data, the maximum likelihood estimates for the model have the most accurate predictive means and best overall fitting scores, but the one-step approximation is comparable.

6 Data application: Precipitation data in Puerto Rico

We apply the mean-dependent nonstationary model to daily rain gauge measurements (in millimeters) from the island of Puerto Rico. The dataset comes from the Global Historical Climatology Network-Daily database [Menne et al., 2012]. We use the square-root of the measurements as the response throughout the data analysis. We consider two time periods: one during the dry season (January 2013) and one during the rainy season (May 2013). Figure 4 depicts the daily mean and maximum values for the precipitation measurements over $n=47$ reporting weather stations. The rainy season has consistently higher means than the dry season. These two time periods were selected to illustrate the seasonal differences in Puerto Rico’s climate, so the resulting best models may be different for other months or years.

Figure 2 shows significant differences in empirical correlograms, so we will fit the mean-dependent nonstationary model separately for each month. For the covariance parameters in (7), we perform an exploratory analysis to determine the proper link functions. For each day in 2013, we estimated a simple stationary model with exponential covariance function to obtain daily estimates for the nugget, sill, and range. These estimates are plotted against the daily averages in Figure 3. The log-transformations on nugget and sill indicate that exponential link functions should be used, $g_{i}(\mu)=\exp(a_{i}+b_{i}\mu)$, which we will refer to as the mean-dependent nonstationary (MDNS) fitted model in the analysis. However, log-transformations on the average lead to better linear relationships with the nugget and sill, so we will also consider $g_{i}(\mu)=\exp[a_{i}+b_{i}\log(1+\mu)]$, which we will refer to as the log-mean-dependent nonstationary (L-MDNS) fitted model. For the range, the log-transformation indicates that MDNS fit would work, but the relationship to the mean may be too weak. As a result, we will also consider fitting the model where $b_{3}=0$ is fixed, which would make the range stationary.

To fit the regime process, we use the longitude $s_{1}$, latitude $s_{2}$, and elevation $e(\mathbf{s})$ (see Figure 5). The mean is then regressed against either $p=4$ linear predictors $Z(\mathbf{s})=[1,s_{1},s_{2},e(\mathbf{s})]^{T}$ or $p=7$ quadratic predictors $Z(\mathbf{s})=[1,s_{1},s_{2},e(\mathbf{s}),s_{1}^{2},s_{1}s_{2},s_{2}^{2}]^{T}$. For the random-effects distribution, we compute the estimated $\beta$ for each day and set \mbox{\boldmath\beta}_{0} and $\Omega$ to the sample mean and covariance of the estimates. We did not consider temporal effects mainly because there were no clear temporal trends, outside of the seasonal differences observed in Figure 4.

We used five-fold cross-validation to evaluate the methods. For each fold, a random selection of locations are removed as testing data. The model is fit using the remaining sites as the training data. The fitted model is then used to construct the predictive distribution at the test sites. We evaluate the accuracy of the predictive distributions using three measures: the average log-likelihood score of the prediction distribution evaluates overall fit of the distribution, the mean squared error MSE measures the difference between the observed values and prediction mean, and the percent coverage of the observed values using 95% probability prediction intervals. We also calculate the 5%, 50%, and 95% quantiles of the standard errors to illustrate how the nonstationary models compare to the stationary models in terms of precision.

Negative values can be a potential issue with the raw precipitation data analysis. Our model may predict negative values, which is not physically possible. This issue could be remedied using a zero-truncated normal distribution [Stein, 1992; Militino and Ugarte, 1999], but this distribution is too computationally inefficient to use on such a large dataset. For practical purposes, we simply threshold negative predictions to zeroes.

6.1 January 2013 data

From the parameter estimates in Table 4(a) for the January 2013 data, we observe that the best fitting model is the full MDNS model with $p=7$ predictors. The highly negative values for $a_{1}$ and $b_{1}$ indicate that very little nugget effect is present. The nonstationary models are significantly better fits than the stationary model ($\chi^{2}\geq 139.6$ for all $p=4$ fits and $\chi^{2}\geq 87.3$ for all $p=7$ fits). The positive $b_{2}$ values and negative $b_{3}$ value indicate that, as the mean increases, the covariability increases and the range of the spatial correlation decreases, as depicted in Figure 8.

The cross-validation results in Table 5(a) indicate that the predictions from the full MDNS model with $p=7$ predictors are the most accurate in terms of low MSE and high prediction score. Also, the small quantiles for the standard error indicate high precision in the predictions. The coverages indicate that the 95% prediction intervals of the fully stationary model fail contain all of the removed values.

Figure 11 depicts the predicted values and standard errors on 26 January 2013, using MDNS and $p=7$ spatial predictors. The raw data appears in Figure 10. The maps predict high precipitation in the northeastern half of the island and practically zero precipitation in the southwestern half, with a possible outlier in the northwest corner. The gaps in the maps indicate regions where the predicted precipitation was below zero, which are treated as zero values for practical purposes.

6.2 May 2013 data

From the parameter estimates in Table 4(b) for the May 2013 data, we see that the L-MDNS model with $p=4$ predictors is the best fitting model. The nonstationary models are significantly better fits than the stationary model ($\chi^{2}\geq 57.2$ for all $p=4$ fits and $\chi^{2}\geq 76.6$ for all $p=7$ fits). For the range, the difference in log-likelihood between stationary and nonstationary fits are not very significant. The positive values for $b_{1}$, $b_{2}$, and $b_{3}$ indicate that the variability and range increases with the mean, as illustrated in Figure 9. The highly negative $a_{2}$ and $b_{2}$ values for the MDNS model with $p=7$ predictors indicate insignificant spatial covariance in the fitted model, and the highly negative $a_{3}$ value and positive $b_{3}$ for the L-MDNS model with $p=7$ predictors indicate that the spatial covariance is significant in areas with high mean precipitation.

From the cross-validation results in Table 5(b), we see that the predictions from full L-MDNS model with $p=7$ predictors are more accurate in terms of low MSE and high prediction score. However, there are insignificant differences in predictions between the stationary and nonstationary range fits. The coverages indicate that the 95% prediction intervals of the fully stationary model fail contain all of the removed values.

Figure 12 depicts the predicted values and standard errors on 9 May 2013, using L-MDNS with $p=4$ spatial predictors. The raw data appears in Figure 10. From the observed values, there was higher precipitation along the northern coastline and around the eastern Pico El Yunque mountain and lower precipitation in the southeastern coastal plains and central mountains. The predictions depict these features quite well, while the high standard errors indicate that most of the uncertainty lies in the transitory areas between the northern coast and the central mountains.

7 Conclusion

From the simulation study, we see that mean-dependent nonstationary covariance models significantly improve predictions if the data are nonstationary. Testing for nonstationarity appears to be best suited for the single step least-squares method, which is faster to estimate and converges quicker to the asymptotic distribution, and if the test is significant, the full MLE method should be used to compute the most accurate predictions. Using the full MLE for the nonstationarity test is more likely to favor the nonstationary model over the stationary model, but this is not a significant practical problem with the method due to the nonstationary model’s overall better predictive fit. Ideally, the full MLE should be used for the test, but the computation cost would have to be drastically reduced through more complex computing methods, such as Stein et al. [2004] or Vecchia [1988].

From the Puerto Rico precipitation data analysis, the nonstationary model significantly improves the accuracy and precision of the prediction distribution over the stationary model. For the dry season month of January 2013, the fully mean-dependent model with quadratic spatial predictors created the best prediction distribution. However, assuming local normality while the dry season data has a preponderance of zeroes may result in an ill-behaved likelihood and unrealistic predictions. Future work could include zero-censored normal distributions to address this issue, though this model was designed to be flexible enough to reasonably model both wet and dry seasons without having to use the computation costs of calculating zero-censored multivariate distributions. For the wet season month of May 2013, the fully log-mean-dependent model with quadratic predictors creates the best prediction distributions, but no significant spatial covariance is presence in the fitted model. The log-mean-dependent model with linear predictors has better internal verifiability in terms of log-likelihood, and the spatial covariance is more clearly present in the fitted model.

The results agree with the notion that precipitation measurements are more volatile during heavy rain events. Accounting for this volatility improves prediction and uncertainty quantification, particularly in regards to the possibility of extreme values in areas with high mean values. Improving prediction distributions will be a significant benefit in weather forecasting, reanalysis models, and for developing long-term time series that blend remotely sensed (satellite or radar) and land-based records, particularly for regions with highly varied climates like Puerto Rico. Further research would incorporate temporal or autoregressive effects to predict the mean process for the future days from previous climate data, then utilize the mean- or log-mean-dependent nonstationary spatial model to fully predict the possible range of precipitation values.

8 Disclaimer

Use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.

Appendix A Covariance parameter MSEs from Simulation Study