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 \in \mathbb{R}\) defined on the real line around the unit circle to make its domain adhere to the interval \([0,2\pi]\). The result is a circular – or wrapped – random variable \(X \in \lbrack 0,2\pi)\). Then, the unwrapped variable \(Y\) can be expressed as \[ Y = X + 2 \pi K \] where \(K \in \mathbb{Z}\) measures the number of “turns” \(Y\) has be wrapped around the unit circle.

Let \(\mathbf{s}\) denote a spatial index variable within a spatial domain \(S \subset \mathbb{R}^2\). 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(\mathbf{s}_i) = \beta_0 + \mathbf{z}(\mathbf{s}_i)' \boldsymbol{\beta} +\gamma(\mathbf{s}_i, \mathbf{z}^\kappa(\mathbf{s}_i)) + 2 \pi K(\mathbf{s}_i) + \varepsilon_i ,\] where \(i=1, \ldots,n\), \(\gamma(\mathbf{s}_i)\) is a zero-mean Gaussian random field (GRF) and \(\varepsilon_i \sim N(0, \sigma^2_\varepsilon)\) are the i.i.d. error terms, with the following novelties:

1. The mean depends on \(\mathbf{z}(\mathbf{s}_i)' \boldsymbol{\beta}\), i.e. it is non-stationary,

2. The covariance depends on \(\mathbf{z}^\kappa(\mathbf{s}_i)\), i.e. it is non-stationary,

3. We include small-scale noise \(\varepsilon_i\)s.

We focus on 1. and 2.

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:

\[\begin{align*} X(\mathbf{s}) &= \eta(\mathbf{s}) + 2 \pi K(\mathbf{s}) + \varepsilon, \\ &= \beta_0 + \mathbf{z}(\mathbf{s})' \boldsymbol{\beta} + \gamma(\mathbf{s}, \mathbf{z}^\kappa(\mathbf{s})) + 2 \pi K(\mathbf{s}) + \varepsilon,\\ \beta_0 &\sim N(0, 10) ,\mbox{ } \boldsymbol{\beta} \vert \xi^2 \sim N(0, \xi^2 \mathbf{I}), \mbox{ }\boldsymbol{\gamma} \vert \boldsymbol{\theta} \sim N(\mathbf{0}, \mathbf{Q}^{-1}(\boldsymbol{\theta})) \\ \theta_0^\tau &\sim U(-10, 3) ,\mbox{ } \theta_0^\kappa \sim U(1, 6) ,\mbox{ } \boldsymbol{\theta}^\kappa_z \vert \zeta^2 \sim N(\mathbf{0}, \zeta^2 \mathbf{I}) \\ \xi^2 &\sim PC\left(\pi, \alpha \right) ,\mbox{ } \zeta^2 \sim PC\left(3, \alpha_\kappa \right) \\ \mathbf{\varepsilon} \vert \sigma^2_\varepsilon &\sim N(\mathbf{0}, \sigma^2_\varepsilon \mathbf{I}) ,\mbox{ } \sigma^2_\varepsilon \sim IG (0.001, 0.001)\\ p(K(\mathbf{s})) &= \begin{cases} \frac{1}{3},& \text{if } K(\mathbf{s}) \in \lbrace -1, 0, 1 \rbrace\\ 0, & \text{otherwise} \end{cases} \end{align*}\]

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 \(\mathbf{z}(\mathbf{s})\) are maximum wind speed, average air humidity, and temperature at 10 meters height. For \(\mathbf{z^\kappa}(\mathbf{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.

About the authors

• Isa Marques
• Thomas Kneib
• Nadja Klein