Geostatistical methods for spatial and spatio-temporal data are in great demand from fields such as earth and climate sciences, epidemiology and agriculture. A comprehensive overview can be found in cressie2011statistics. Spatially stationary processes are commonly used as models in geostatistical applications, but often the assumption of stationarity and isotropic covariance functions are difficult to hold in real applications; see for example, guttorp1994space, le2006statistical, damian2001bayesian. A recent review of methods that allow nonhomogenous covariance models is schmidt2020flexible.
We propose a semiparametric method of nonstationary spatial covariance function that expands upon the work of sampson1992nonparametric. Data observed in a set
2 Spatial Deformation Model (SDM)
Let be spatial locations with for all in a geostatistical domain and with for all Moreover, let such that for all and are bijective, differentiable functions, for
Let be a Gaussian random field with spatial covariance function For any pair of spatial sites we have with given by
where and is a stationary, isotropic correlation function with parameters We remark that this model may not be identifiable, since is invariant for shifts and rigid rotations in Moreover, scaling by a constant will produce the same if is also scaled by
Nevertheless the model can be used for Kriging and spatial interpolation without issues.
For example, consider and given by the swirl transformation shown in Figure 1. We sampled from a spatio-temporal process on a regular grid of points, with and separable covariance function
where and A realization of a random field on a fine mesh grid is also shown in Figure 1, obtained with the RandomFields package (schlather2015analysis).
3 Spatial Deformation Estimation
then is the semivariogram of the random field . In particular, for a stationary, isotropic
where is a positive non-decreasing function of . Note valid functions are conditionally negative definite (cressie1993statistics).
If there are ways to obtain a sample covariance matrix, such as temporal replicates, then we can construct a sample variogram as
If the variogram is isotropic on a set of artificial coordinates , then it must be have form
where Therefore, the variogram entries can be seen as a dispersion metric, with
3.1 Non-metric multidimensional scaling
sampson1992nonparametric first step is to consider a non-metric multidimensional scaling approach (hereafter, nMDS; see kruskal1964multidimensional; mardia1979multivariate; cox2000multidimensional). From a set of dispersions we seek a monotone transformation such that
and an artificial set of coordinates with interpoint Euclidean distances that minimizes the stress criterion
Here we need to impose constraints on such that the function is conditionally positive definite. kruskal1964multidimensional only finds an isotonic transformation, so sampson1992nonparametric adapted the optimization algorithm as follows:
Set k = 1. Find an initial configuration using, for example, Kruskal’s isotonic nMDS.
Fit a variogram model to the data using the coordinates obtaining
Obtain via MDS performed on
Set and return to step 2 until convergence.
3.2 Spline smoothing
The second step in sampson1992nonparametric approach is to find an approximation for such that This can be done if we approximate the functions with, for example, thin-plate splines (wahba1980some). Thin-plate splines are the minimizers of the variational problem
The solution for a fixed is
where and are estimating by plug-in.
3.3 Tensor product of B-splines
We propose a different approach for the spline smoothing by making use of low rank approximation (finite approximation) for the deformation functions. Specifically, let and assume we have a collection of points be a sequence of points in . Define a map such that
where the functions , form a partition of unity.
Particularly, suppose that are the B-splines basis functions and there is a integer number such that a map can be represented by
where and Let for some index set such that then the functions are well approximated by
where are fixed positive integers and are B-spline basis functions (see, e.g., ramsay2005functional).
To guarantee that the functions do not fold onto themselves, we must guarantee that is locally invertible and differentiable. In fact, a diffeomorphism is desirable (see, e.g., perrin1999modelling). A necessary condition is that everywhere in the domain,
but in practice such constraint is difficult to implement, requiring the evaluation of the Jacobian for every pair We will show in Section 4 that such condition can be translated into a condtion onto the coefficients for Note however that it is enough to ensure since a change in signs would imply a discontinuous Jacobian.
3.4 Simultaneous estimation of covariance function and deformation
The covariance parameters
as well as the mean vectorare estimated by maximizing the profile log-likelihood. Without loss of generality set Thus our proposed B-spline approach can be estimated as
Given solve the optimization problem
where corresponds to the optimization target seen as a function of only, with parameters fixed.
Given obtain such that
where are the row vector of basis functions and is the Kronecker product, and for some choice of such that if
Repeat steps 1 and 2 until convergence.
We remark that the matrices are sparse when using B-splines, therefore solutions to step 2 are scalable. To evaluate the inequality consider the following: write
or simply where is a vector of B-spline basis functions evaluated at , similarly is a vector of B-spline basis evaluated at , and is a matrix of spline coefficients for the tensor product approximation of the -th function . Similarly,
is the partial derivative of with respect to and if we write then This allows us to see that
so the determinant of the Jacobian as a function of is an inner product of
, weighted by the skew-Symmetric matrixEnsuring the inequality for all values of remains a difficult task, but we can chose a set of basis for which does not depend on
4 Constrained spatial deformation estimation
We will use an approach similar to musse2001topology. Consider B-splines of degree 1 on . Assume that there are equally spaced inner knots, where . Since the knots are equally spaced, they can be written as where In this case, the B-spline bases are given by
Consider where and is between and (where and ). Then there are only 4 bases that evaluate to non-zero values, indexed by , , and , so
and therefore, looking only at the non-zero pairs , we have
note is therefore proportional to
The above equations describe a plane in with coefficients depending on . Now, since and , where (and similarly for ), we have four restrictions to consider:
When and ,
When and ,
When and ,
When and ,
This collection of constraints, for and together imply in a non-folding deformation map, and can be enforced with constrained optimization routines. We have employed svanberg2002class constrained optimization algorithm, avaliable in the nloptr package (johnson2020nloptr). The code is available as an R package in https://github.com/guiludwig/bsplinedef.
5 Simulation study
To evaluate the performance of the algorithm, we conducted a simulation study based on the swirl function shown in Figure 1. Each sample has spatial points on a regular grid in the geographical domain We simulated from a Gaussian random field with mean function and covariance function with parameters and First, consider a single realization of the spatial random field. The estimated deformation maps are shown in Figure 2. We have obtained the constrained B-spline deformations (hereafter, bdef) with and basis functions. The estimated deformation function does not fold, even though the true deformation function is difficult to be recovered. The case when performs better than or
indicating that there are features in the deformation map that need a large number of degrees of freedom to be estimated. On the other hand, the functions usingsampson1992nonparametric (SG) method have smoothing parameters and set to match the bdef approach. They start showing folding at and cannot recover deformation maps that require a number of degrees of freedom larger than Note that the SG maps were stretched or shrunk to fit the plot area, but the bdef maps did not require this step.
The comparison of estimated covariance matrices allow us to overlook the identifiability issues with rotations, shifts and scaling of the estimated maps. In Figure 3 we show a scatterplot of the upper-diagonal entries of the estimated covariance matrices for the data, versus the true covariance matrices. We remark that the bdef method shows no apparent bias and becomes more accurate as the number of degrees of freedom increase. On the other hand, SG has good performance at a small number of degrees of freedom (more smoothing), but tends to overestimate matrix entries as the number of degrees of freedom for the thin-plate spline increases.
6 Case study: Rainfall data in southeastern Brazil
The dataset we use to illustrate our method comes from meteorological surveys conducted by INMET – Instituto Nacional de Meteorologia, Brazil. The measurements are made at every 15 minutes, and daily accumulated values are made available. Following rozante2010combining, we grouped the rainfall data in periods of 10 days.
Since we seek temporally stationary data, we decided to restrain the data collection to 2018-01-01 to 2018-03-30 (rainfall season), for a total of 9 time periods (period 1: 2018-01-01 to 2018-01-10, …, period 9: 2018-03-22 to 2018-03-31).
We selected the 50 stations of the southeastern region of Brazil that had complete observations available during the period. The stations are shown in the left panel of Figure 4. Data can be obtained at http://www.inmet.gov.br/portal/index.php?r=bdmep/bdmep.
In the right panel of Figure 4, we show the estimated deformation map, using B-spline basis functions. The estimated deformation map reveals topographical features shown in the left panel of the same Figure. Weather stations located in northwestern flat lands of Minas Gerais are treated as closer to each other than stations near the rough coast of Rio de Janeiro and São Paulo states, where elevation changes are abrupt.
In Figure 5 we perform conditional simulation of the Gaussian random field (Kriging) for the 10-day periods starting in 2018–01–01 and 2018–01–11, showing the resulting prediction maps.
Ronaldo Dias was supported by the FAPESP foundation, grant 2018/04654-9. Guilherme Ludwig was supported by FAPESP grants 2018/04654-9 and 2019/03517-0, as well as FAEPEX grant 2253/19.
Appendix A Derivation of the non-folding constraint
We can derive the constraint presented in Section 4 as follows: Observe that
Further simplifying leads to