Spatial Regression Analysis Continued¶

These are the libraries we will be using

And this is some code that will become handy at some point during the exposition. The files will be uploaded in the course web page.

Some initial adjustments to define the plots size, and we are ready to go...

iid and beyond¶

Let $y_i = \beta_0 + \beta_1 x_i + f_i$, with $f_i \sim \mathcal{N}(0, 1)$.

Say $\beta_0 = 2$ and $\beta_1 = 1/3$. For a random realization of $f_i$'s we would have something like this:

Now we are going to remove $\beta_0$ and $\beta_1 * x$ and focus only on what is going on with $f$.

Not very interesting since all the $f_i$ are iid.

We can think of the vector $f$ as being generated by a Multivariate Normal, with mean zero and covariance matrix: $I \sigma^2$.

The covariance matrix would look like this:

Let's assume that $f = f(x)$. This assumption doesn't change anything so far.

But what would happen if $f_i$ were not iid., and their pairwise correlations dependend on the values $x_i, x_j$?

Exponentiated Quadrtic Kernel¶

Also know as RBF kernel. For any pair of points $x_i, x_j$, the covariance is given by the function:

$k(x_i, x_j) = \sigma^2 exp \big(- \frac{(x_i - x_j)^2}{2 \ell ^2} \big)$¶

The plot above is using $\sigma = 1$ and $\ell = 1$. What happens with different values of $\sigma$?

Now different values of $\ell$.

In our initial example, we had these randomly generated errors:

Think what would happen if instead of a covariance given by $I \sigma^2$, we had a covariance given by the Squared Exponential kernel. How do random realizations of $f$ would look like?

Periodic kernel¶

$k(x_i, x_j) = \sigma^2 exp \big( \frac{ - 2 \sin( \pi |x_i - x_j|)^2}{\ell^2} \big)$¶

Let's look at some random realizatios of $f$.

Linear kernel¶

$k(x_i, x_j) = b + \sigma^2 (x_i - c) (x_j - c)$¶

You can find out more about kernels here: http://www.cs.toronto.edu/~duvenaud/cookbook/

Kernel functions in 2D¶

This is a class on spatial statistics, so we need to consider at least 2 dimensions.

The kernels we saw before are defined in the same way. The only changes we need are the following:

• Instead of differences $x_i - x_j$ we use euclidean distances $\sqrt{(x_{i1} - x_{j1})^2 + (x_{i2} - x_{j2})^2}$
• Instead of products $x_i x_j$ we use inner products $x_{i1} x_{j1} + x_{i2} x_{j2}$

Consider a grid of points in 2D.

A random realization of $f(x)$ under the iid assumption would look like:

Squared Exponential kernel:

Periodic kernel:

Depending on the kernel used (and its parameters) $f$ is able to represent different functions across the space.

Gaussian Process Regression¶

So far we have only spoken about random realizations of multivariate normals given a kernel, but we haven't actually linked it to an inference method.

Definition¶

Let $f(x)$ be a function such that any finite collection of its realizations $\{f_1, \ldots, f_n\}$ follows a multivariate Gaussian distribution, with mean $m$ and covariance function $K$. Then $f$ is known as a Gaussian process.

Say we have some a set of observations $(x_i, y_i)$ and we are interested in knowing which functions $f(x)$ are able to represent our observations. We use the following model:

$y_i = f_i(x_i) + \epsilon_i$, with $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$.

$f_i$ are never observed. They are latent variables of our model that we need to learn. For every observation $(x_i, y_i)$, we need to infer the corresponding $f_i$.

We do this using a Bayesian approach as follows:

$p(f|y,x) = \frac{p(y|f,x)p(f|x)}{p(y|x)}$¶

As prior $p(f|x)$ we use a multivariate normal with covariance kernel K. The type of kernel depends also on our knowledge or assumptions of the relation between $x$ and $y$. Let the mean of our prior be zero, as this simplifies the explication.

Then we have that $p(f|x) = N(0, K)$.

Now if $f$ is given, we have that $p(y|f,x) = N(f, I\sigma^2)$.

And the posterior of $f$ becomes:

$p(f|y,x) = N \big( K (K + I\sigma^2)^{-1} y, K - K (K + I\sigma^2)^{-1} K \big)$¶

The proof of this formulation goes beyond the scope of todays discussion, so we will skip it. However you can find out more about the technical details of Gaussian processes here: http://www.gaussianprocess.org/gpml/

The solution of this inference problem is analytically tractable. And this is not only for $f(x)$, but for predictions $y_\ast$ at new input locations $x_\ast$ as well. The predictive distribution for $y_\ast$ is given by:

$p(y_\ast|y,x, x_\ast) = N \big( K_{\ast x} (K_{xx} + I\sigma^2)^{-1} y, K_{\ast\ast} - K_{\ast x} (K_{xx} + I\sigma^2)^{-1} K_{x\ast} \big) + I\sigma^2$¶

where $K_{xx} = k(x,x)$ and $K_{\ast x} = k(x_\ast, x)$.

Now let's see what this means in terms of the inference process. We will us an exponentiated quadratic kernel for this example.

Example: Earthquakes around the world¶

Things to improve?

• How to avoid negative values?
• How to select the kernel parameters (hyperparameters)?

Laplace Approximation¶

We mentioned before that $p(f|y, x)$ has a Gaussian form. This is the result of it being the product of two Gaussian distributions: $p(y|f,x)$ and $p(f|x)$.

If we assume that $y$ is not Gaussian (e.g. is positive or discrete), the posterior distribution will no longer be Gaussian. Hence to know of

$p(f|y,x) = \frac{p(y|f,x)p(f|x)}{p(y|x)}$¶

will require calculating the denominator

$p(y|x) = \int p(y|f,x) p(f|x) df$¶

which might render intractable.

The Laplace Approximation estimates the posterior distribution with a Gaussian distribution that resembles it the closest.

https://en.wikipedia.org/wiki/Laplace%27s_method

http://www.sumsar.net/blog/2013/11/easy-laplace-approximation/

Hyperparameters¶

We can define priors over them:

• MCMC https://chi-feng.github.io/mcmc-demo/
• Maximum likelihood type II

The big $n$ problem¶

Although we have a closed analytical expression for computing predictions, one disadvantage of these methods is the computational cost of the algebraic operations where the covariance function is involved.

Both the predictive mean and the predictive covariance, require the computation of the factor:

$(K + I\sigma^2)^{-1}$¶

This is known to have a computational cost of order $O(n^3)$. Which means a cubic increase in the number of operations needed, related to the number of data points modeled.

In particular, for applications in epidemiology, where we it is needed to make inference on large regions at, but with enough detail at local level, this approach can be infeasible.

SPDE approach¶

In previous sessions, the topic of Gaussian Markov Random Fields (GMRF) has been discussed. GMRF are commonly applied for inference on areal data. The neighboring structure leads to sparse precision matrices, and thus allows computational benefits.

Why?

• Elements in the precision matrix are non-zero only for neighbour areas and diagonal elements.
• Sparse matrices are easer to invert

A continuous field can be represented using a GMRF indexed across space. The countinuous process $f(x)$ is replaced by a discrete field using a combination of weights $w$ and discrete basis functions $\psi(x)$:

$f(x) \rightarrow \sum_j w_j \psi_j(x)$¶

$w$ is found by solving the stocastic partial differential equation (SPDE):

$(\kappa^2 - \Delta)^{\alpha/2}\tau f(x) = W(x)$¶

where $\Delta$ is the Laplacian, $\alpha$ controls the smoothness, $\kappa$ is a scale parameter and $\tau$ is a variance rescaling parameter.

The solution to this equation is the a spatial process with Matern covariance:

$k(x_i, x_j) = \frac{\sigma^2}{\Gamma(\lambda) 2^{\lambda - 1}} (\kappa (||x_i - x_j||)^\lambda K_\lambda (\kappa||x_i - x_j||)$¶

The core ideas of this approach are:

• Construct a discretized spatial field.
• Find a GMRF with local neighborhoud and precision matrix $Q$ that represents the Gaussian field.
• Do the computations using the GMRF representation taking advantage of the sparse matrices.

R-INLA¶

We will dedicate the rest of the session to show some how to use the INLA package in R to implement geostatistical models. We will just cover some basic steps. A detailed tutorial of R-INLA can be found here: http://www.math.ntnu.no/inla/r-inla.org/tutorials/spde/html/

INLA stands for Integrated Nested Laplace Approximation. The implementation of geostatistical models is done by using the SPDE approximation.

The SPDE representation also makes it possible to use the Integrated Nested Laplace Approximation (INLA) to handle the non-linear transformation on the latent variable. This algorithm computes the posterior marginals of the latent field through a step-wise implementation of a simplified Laplace approximation with a cost of O(n^2 log n).

You can also check these publications:

http://www.statslab.cam.ac.uk/~rjs57/RSS/0708/Rue08.pdf

http://www.maths.ed.ac.uk/~flindgre/rinla/isbaspde.pdf

Example: Burkina Faso¶

We will use an example of Malaria prevalence in Burkina Faso. We have data from 109 villages sampled.

At each village, the number of positves can be modeled as:

$y \sim Binomial(p, n)$¶

where p is the rate of positives and $n$ is the number of examined people.

Say that $p$ varies across space. We can define a Gaussian process $f$, that is linked to $p$ through the transformation:

$p = \frac{1}{1 + e^{-f}}$¶

Notice that while $f \in (-\infty, \infty)$, $p \in (0, 1)$.

First we need to define a mesh. A "good" mesh:

• has few very sharp angles
• has similar sized triangles
• extends outside the observation area to reduce edge effects