Introduction
With the release of the R package varycoef on CRAN (see) we enable you to analyze your spatial data and in a simple, yet versatile way. The underlying idea are spatially varying coefficients, short SVC, which have been studied over the last decades. These models are based on the linear regression model and are therefore very easy to interpret. Due to their design, SVC models offer a high flexibility. These two properties allow for better understanding of the data, gaining new insights and, ultimately, better predictive performance.
In this blog post, we will show what SVC models are and how to define them. Afterwards, we give a short and illustrative example on how to apply these tools in varycoef using the well known data set meuse from the package sp.
Disclaimer This analysis is meant to be a show case of what is possible using the package varycoef, not to write the most rigorous statistical analysis of a data set.
But first, let us start with the set up and let me show you where to find help:
# install from CRAN
# install.packages("varycoef")
# attach package
library(varycoef)
# general package help file
help("varycoef")
# vignette for detailed information on how to use varycoef
vignette("example", package = "varycoef")## Warning: vignette 'example' not foundThere is also a version of varycoef on Github from which you can get the latest devel version.
Preliminaries
Before we start, we want to give some prerequisites in order to understand the analysis below. Beside a classical linear regression model, we require the knowledge of geostatistics and tools thereof, like:
- what is the advantage of geostatistics?
- two-dimensional Gaussian processes, i.e. Gaussian random fields
- covariance functions and how the range and variance parameter influence them
- maximum likelihood estimation
What are Spatially Varying Coefficient Models?
SVC models are a generalization of the classical linear regression model. In matrix notation with a response vector \(\mathbf y\), a data matrix \(\mathbf X\), and a vector \(\boldsymbol \varepsilon\) containing the errors, a linear regression model is defined as
\[\mathbf y = \mathbf X \boldsymbol \beta + \boldsymbol \varepsilon.\]
In geostatistics, the error term usually is called the nugget. SVC models are introduced through a random effect \(\boldsymbol \eta(s)\) that is depending on the location \(s\). The underlying model takes the form
\[\mathbf y = \mathbf X \boldsymbol \beta + \mathbf W \boldsymbol \eta (s) + \boldsymbol \varepsilon\] where \(\mathbf W\) parses the covariates to the random effect \(\boldsymbol \eta(s)\). Specifically, each SVC is defined by a Gaussian process, more specifically a Gaussian random field (GRF), with mean zero and some stationary covariance function. This function usually depends on some marginal variance (sill) and a range. An example of 3 SVC, i.e. 3 GRFs, and the nugget \(\boldsymbol \varepsilon\) on a regular grid is given in the figure below.
# setting seed
set.seed(123)
# number of SVC
p <- 3
# sqrt of total number of observations
m <- 20
# covariance parameters
(pars <- data.frame(var = c(0.1, 0.2, 0.3),
scale = c(0.3, 0.1, 0.2)))## var scale
## 1 0.1 0.3
## 2 0.2 0.1
## 3 0.3 0.2nugget.var <- 0.05
# function to sample SVCs
sp.SVC <- fullSVC_reggrid(m = m, p = p,
cov_pars = pars,
nugget = nugget.var)
library(sp)
# visualization of sampled SVC
spplot(sp.SVC, colorkey = TRUE)
Just by looking at the figure, one can see that
- SVC_1, SVC_2 and SVC_3 have some spatial clustering…
- …, whereas the nugget does not.
- The largest variance can be observed in SVC_2 and SVC_3.
- The largest spatial clusters are in SVC_1.
- They all vary around 0, which is expected as they are all zero mean GRF.
Goal of SVC Modelling
Given the response \(\mathbf y\), the data matrices \(\mathbf X\) and \(\mathbf W\), we want to estimate the coefficients \(\boldsymbol \beta\) as we would with a linear regression model, and further, try to understand the properties of \(\boldsymbol \eta(s)\). Another important aspect is that we are able to predict \(\boldsymbol \eta(s')\) for some new locations \(s'\). This is what we are about to do in the following sections.
