Non-stationary wrapped Gaussian spatial response model
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Background
Circular data, i.e., data defined on the unit circle, can be found in many areas of science. The unique nature of these data means that conventional methods for non-circular data are not valid for these. At the same time, advances in geographical information and global positioning systems have generated large amounts of spatial data and, consequently, have increased the need for spatial models that accurately describe it. An example of circular, spatial data are the average wind directions in Germany from April 1 to 10, 2019, plotted in space (left) and on a circular histogram (right):
It was Jona-Lasinio, Gelfand, and Jona-Lasinio (2012) and Wang and Gelfand (2014) that first brought the wrapped and projected Gaussian distributions to space and space-time settings. Nonetheless, today,circular data models still need to catch up on the latest advances in spatial statistics. In particular, models developed so far assume stationarity in the mean and the covariance of the response modeled; i.e., these are constant in space. Our work in Marques, Kneib, and Klein (2022) aims to bridge this gap and allow for non-stationarity in the mean and covariance of distributions for circular responses defined on spatial domains.
The response distribution
To construct distributions for circular data, we can rely on the construction principle of wrapping the distribution of a random variable Y∈R defined on the real line around the unit circle to make its domain adhere to the interval [0,2π].
The result is a circular – or wrapped – random variable X∈[0,2π).
Then, the unwrapped variable Y can be expressed as Y=X+2πK
Let s denote a spatial index variable within a spatial domain S⊂R2.
Jona-Lasinio, Gelfand, and Jona-Lasinio (2012) first showed that this wrapping approach can be extended to the context of spatial data.
We extend their framework and consider
X(si)=β0+z(si)′β+γ(si,zκ(si))+2πK(si)+εi,
The mean depends on z(si)′β, i.e. it is non-stationary,
The covariance depends on zκ(si), i.e. it is non-stationary,
We include small-scale noise εis.
We focus on 1. and 2.
1. Non-stationarity in the mean
Linear covariates included in the mean of a circular response might induce a circular likelihood for the regression coefficients that has infinitely many maxima. The typical solution is to use a link function for fixed effects with co-domain defined in R and of length equal to the circular variable’s period of 2π, e.g. the inverse tangent link function.
We argue that:
The winding number K(s) can govern the behavior of the combination of GRF and fixed effects such that no link function is needed,
We use prior information that shrinks towards a stationary mean and enforces the length of the fixed effects to be approximately 2π.
Concretely, β0∼N(0,10) represents a diffuse prior for the intercept and we let β|ξ2∼N(0,ξ2I). One can show that if the base model is such that ξ2→0 and the prior is constructed according to the PC-prior principles (Simpson et al. 2017), then ξ2 has a Weibull prior with shape 12 and scale λ, i.e., ξ2∼Weibull(12,λ) (Klein and Kneib 2016).
For α∈(0,1), λ is chosen such that
P(maxs∈S∣z(s)′β∣≤π)≥1–α,
2. Non-stationarity in the covariance
The stochastic partial differential equation (SPDE) approach (see Lindgren, Rue, and Lindström 2011) allows us to model non-stationarity in the covariance of the GRF in a computationally efficient way. Consider the parameters of the Matérn covariance with smoothness 1 and κ and τ controlling the range and the marginal variance of the GRF, respectively.
We let
log(τ)=θτ0 and log(κ(s))=θκ0+zκ(s)′θκz,
P(maxs∈S∣zκ(s)′θκz∣≤c)≥1–ακ,
More concretely, we have θτ0∼U(−10,3), θκ0∼U(1,6) and θκz|ζ2∼N(0,ζ2I), with ζ2∼Weibull(12,λκ).
The full-hierarchical model
The SPDE is typically used in combination with integrated nested Laplace approximations (INLA). However, INLA does not perform well when combined with shrinkage priors for non-stationarity in the covariance as defined here. Thus, we use Markov Chain Monte Carlo with full-hierarchical model:
X(s)=η(s)+2πK(s)+ε,=β0+z(s)′β+γ(s,zκ(s))+2πK(s)+ε,β0∼N(0,10), β|ξ2∼N(0,ξ2I), γ|θ∼N(0,Q−1(θ))θτ0∼U(−10,3), θκ0∼U(1,6), θκz|ζ2∼N(0,ζ2I)ξ2∼PC(π,α), ζ2∼PC(3,ακ)ε|σ2ε∼N(0,σ2εI), σ2ε∼IG(0.001,0.001)p(K(s))={13,if K(s)∈{−1,0,1}0,otherwise
Outcome
One can use our model to estimate wind directions in the dataset above and test the model’s performance. We randomly select a holdout set consisting of 20% of the locations and use the remaining 80% as training data. The covariates considered for z(s) are maximum wind speed, average air humidity, and temperature at 10 meters height. For zκ(s) we consider altitude and latitude.
The figure below shows the posterior predictive mean (left) and standard deviation (right) for each test location for a fully stationary, i.e., stationary mean and covariance and fully non-stationary models. The black arrows represent the true wind directions (left) and the corresponding standard deviation of zero (right). The figure shows that the fully non-stationary model reaches smaller bias than the fully stationary model throughout the whole domain. In the north where this bias is close to zero, our model also reaches lower uncertainty.
References
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.