Bayesian Regression with Gaussian Processes

Spatial Statistics - BIOSTAT 696/896

Michele Peruzzi

University of Michigan

From response model to regression with latent GP

  • We have considered (ML) inference of the following model

Y(\cdot) \sim GP(0, C_{\theta}(\cdot, \cdot))

  • Optimizing becomes more difficult as we increase the dimension of \theta
  • We may want to model a non-zero mean based on covariate matrix X
  • We may want to make the model more complex or flexible

Bayesian linear regression with latent GP

  • Linear regression model with spatial random effects
  • Spatial random effects = unobserved error term with spatial correlation
  • Equivalent to having correlated measurement error
  • Separate purely spatial error from purely noise error terms

At location s, we model

Y(s) = X(s)^\top\beta + w(s) + \varepsilon(s)

  • X(s) is a p\times 1 vector of covariates at location s
  • \beta is a constant vector of regression coefficients
  • w(s) is the realization at s of a GP
  • \varepsilon(s) is independent noise

Model DAG

Model DAG

Bayesian linear regression with latent GP

At location s, we model

Y(s) = X(s)^\top\beta + w(s) + \varepsilon(s)

  • X(s) is a p\times 1 vector of covariates at location s

  • \beta is a constant vector of regression coefficients

  • w(s) is the realization at s of a GP: we assume w(\cdot) \sim GP(0, \tau^2 C_{\theta}(\cdot, \cdot) ) where we have included \tau^2 for conjugacy

  • Example: C_{\theta}(s, s') = s^2 \exp\{ -\phi \| s-s' \| \} and the process variance can be computed as

    Var(w) = Cov(w(s), w(s)) = \tau^2 s^2

  • Notice: if we denote Var(w) = \sigma^2 then \sigma^2 = \tau^2 s^2 and we interpret s^2 as the signal-to-noise ratio

\text{SNR:} \qquad s^2 = \frac{\sigma^2}{\tau^2}

  • \varepsilon(s) is iid noise \varepsilon(s) \sim N(0, \tau^2)
  • We do not need a nugget term in C_{\theta}(\cdot, \cdot) as we are including iid noise in our model

Bayesian linear regression with latent GP

Take \cal T as the set of observed locations, then the model can be written as

Y = X\beta + w + \varepsilon where Y, X, w and \varepsilon are vectors and matrices with appropriate dimensions.

  • We call this a latent GP model because we are explicitly introducing latent (unobserved) variables (i.e., each element of w)
  • Let’s assign conjugate priors to most unknown parameters:
\begin{align*} \tau^2 &\sim Inv.G(a, b) \\ \beta | \tau^2 &\sim N(0, \tau^2 V_{\beta}) \end{align*}
  • GP assumption implies w \sim N(0, \tau^2 C_{\theta})
  • The prior on \theta is not conjugate: write generically \theta \sim p(\theta) (we will specify later)

Bayesian linear regression with latent GP

Priors:

\begin{align*} \theta &\sim p(\theta) \\ \tau^2 &\sim Inv.G(a, b) \\ \beta | \tau^2 &\sim N(0, \tau^2 V_{\beta}) \\ \text{at}\ \cal T \ \text{:}\qquad w &\sim N(0, \tau^2 C_{\theta}) \qquad \text{because} \qquad w(\cdot) &\sim GP(0, \tau^2 C_{\theta}(\cdot)) \end{align*}

Model/Likelihood:

Y | \beta, w, \tau^2 \sim N(X \beta + w, \tau^2 I_n)

Posterior we want to find:

p(\beta, \theta, \tau^2 | Y) = \text{???}

Bayesian linear regression with latent GP

Posterior we want to find:

p(\beta, \theta, \tau^2 | Y) = \text{???}

  • w are “nuisance” parameters. They are “in the way” when estimating \theta a posteriori
  • p(w|Y) is the posterior of the GP, which we will need to make predictions at new locations

p(w | Y) = \int p(w, \beta, \theta, \tau^2 | Y)d\beta, d\theta, d\tau^2 = \int p(w | Y,\beta, \theta, \tau^2)p(\beta, \theta, \tau^2 | Y)d\beta, d\theta, d\tau^2

  • Therefore, once again we first need p(\beta, \theta, \tau^2 | Y)
  • Notice: p(w | Y, \beta, \theta, \tau^2) is the conditional posterior for fixed values of \beta, \theta, \tau^2
  • We have a roadmap, now we will work towards finding p(\beta, \theta, \tau^2 | Y)

Finding p(\beta, \theta, \tau^2 | Y)

  • We marginalize over w (aka integrate out w) and get the marginal likelihood

p(Y | \beta, \theta, \tau^2) =\int p(Y, w| \beta, \theta, \tau^2) dw = \int p(Y | \beta, w, \theta, \tau^2) p(w | \theta) dw

  • With Gaussian things, this is easy! If U \sim N(m_U, S_U) independent of V \sim N(m_V, S_V) and Z = U + V + \varepsilon where \varepsilon \sim N(0, a I_n) then Z|U,V \sim N(U + V, aI_n), Z|U \sim N(U + m_V, S_V + aI_n), Z \sim N(m_U + m_V, S_U + S_V + aI_n).

  • In our case, the Gaussisan things are \beta, w, and Y, therefore:

Y | \beta, \theta, \tau^2 \sim N(X \beta, \tau^2 C_{\theta} + \tau^2 I_n) = N(X \beta, \tau^2 (C_{\theta} + I_n))

  • Notice: the marginal model is just the likelihood of the usual linear regression model, with correlated errors! (See e.g. slides on preliminaries, page 13, set V_y = C_{\theta} + I_n)

Finding p(\beta, \theta, \tau^2 | Y)

Having the Marginal Likelihood, our target posterior can be written as

p(\beta, \theta, \tau^2 | Y) \propto p(Y | \beta, \theta, \tau^2) p(\beta, \theta, \tau^2)

This is not a conjugate model because of \theta, however:

p(\beta, \theta, \tau^2 | Y) \propto p(\beta, \tau^2 | Y, \theta) p(\theta | Y)

  • p(\beta, \tau^2 | Y, \theta) is the posterior of the marginal model when we condition (or fix) a value of \theta
  • Because we chose a Normal-Inverse Gamma prior on (\tau^2, \beta), we know that p(\beta, \tau^2 | Y, \theta) is also a Normal-Inverse Gamma, with updated parameters.

Finding p(\beta, \theta, \tau^2 | Y)

For a fixed and known \theta value, set V_y = C_{\theta} + I_n.

The Bayesian hierarchical model for our latent GP spatial regression is:

\begin{align*} \tau^2 &\sim Inv.G(a, b) \\ \beta | \tau^2 &\sim N(0, \tau^2 V_{\beta}) \\ Y | \beta, \tau^2 &\sim N(X \beta, \tau^2 V_y) \end{align*}

Then the posterior is (see prelim slides here):

p(\tau^2, \beta | Y) = Inv.G(\tau^2; a^*,b^*) N(\beta; m_{\beta}^*, \tau^2 V_{\beta}^*)

where we just need to compute the updated parameters

\begin{align*} V_{\beta}^* &= (V_{\beta}^{-1} + X^T V_y^{-1} X)^{-1} \\ m_{\beta}^* &= V_{\beta}^* X^T V_y^{-1} y \\ a^* = a + n/2 \quad\text{and}\quad b^* &= b + \frac{1}{2}(y^TV_y^{-1} y - m_{\beta}^* V_{\beta}^{* -1} m_{\beta}^*) \end{align*}

Finding p(\beta, \theta, \tau^2 | Y)

For a fixed and known \theta value,

the Bayesian hierarchical model for our latent GP spatial regression is:

\begin{align*} \tau^2 &\sim Inv.G(a, b) \\ \beta | \tau^2 &\sim N(0, \tau^2 V_{\beta}) \\ Y | \beta, \tau^2 &\sim N(X \beta, \tau^2 \textcolor{teal}{(C_{\theta} + I_n)}) \end{align*}

Then the posterior is (see prelim slides here):

p(\tau^2, \beta | Y) = Inv.G(\tau^2; a^*,b^*) N(\beta; m_{\beta}^*, \tau^2 V_{\beta}^*)

where we just need to compute the updated parameters

\begin{align*} V_{\beta}^* &= (V_{\beta}^{-1} + X^T \textcolor{teal}{(C_{\theta} + I_n)}^{-1} X)^{-1} \\ m_{\beta}^* &= V_{\beta}^* X^T \textcolor{teal}{(C_{\theta} + I_n)}^{-1} y) \\ a^* = a + n/2 \quad\text{and}\quad b^* &= b + \frac{1}{2}(y^T\textcolor{teal}{(C_{\theta} + I_n)}^{-1} y - m_{\beta}^* V_{\beta}^{* -1} m_{\beta}^*) \end{align*}

Finding p(\beta, \theta, \tau^2 | Y)

Recap so far:

  • We want to compute p(\beta, \theta, \tau^2 | Y)
  • We have found that if we fixed \theta, then the “remaining” posterior is conjugate
  • We know everything about p(\beta, \tau^2 | Y, \theta)!
  • For example, we can sample from it, i.e. we can build

\{ \tau^2_{(t)}, \beta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \tau^2_{(t)}, \beta_{(t)} \sim p(\beta, \tau^2 | Y, \theta)

  • Unfortunately for us, p(\beta, \theta, \tau^2 | Y) \neq p(\beta, \tau^2 | Y, \theta)
  • The posterior we really want is

p(\beta, \theta, \tau^2 | Y) \propto p(\beta, \tau^2 | Y, \theta) p(\theta | Y)

Finding p(\beta, \theta, \tau^2 | Y)

  • Our task is to find p(\beta, \theta, \tau^2 | Y)
  • Notice that we can write it “the other way”, i.e. p(\beta, \theta, \tau^2 | Y) \propto p(\theta | Y, \beta, \tau^2) p(\beta, \tau^2 | Y)
  • In the term p(\theta | Y, \beta, \tau^2), we are looking at the posterior of \theta after fixing a value for the pair (\beta, \tau^2)
  • The term p(\beta, \tau^2 | Y) is not the N-Inv.G we had before. Notice this does not condition on \theta

Make a big leap of faith (for now):

What if we could extract samples from p(\theta | Y, \beta, \tau^2)?

i.e., suppose we have a procedure that allows us to build

\{ \theta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \theta_{(t)} \sim p(\theta | Y, \beta, \tau^2)

Finding p(\beta, \theta, \tau^2 | Y)

Putting things together:

We know how to get

\{ \tau^2_{(t)}, \beta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \tau^2_{(t)}, \beta_{(t)} \sim p(\beta, \tau^2 | Y, \theta)

We “know” how to get

\{ \theta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \theta_{(t)} \sim p(\theta | Y, \beta, \tau^2)

Therefore, we can build a Gibbs sampler!

Gibbs sampling

  • Suppose we want to extract samples from the joint distribution p(A, B)
  • We know how to extract samples from p(A | B) and p(B | A)
  • We then build a Markov chain \Theta_{t} = (A_t, B_t) on discrete time
  • Fix some initial value \Theta_0 by setting A_0 = a, B_0 = b
  • At time t we sample a new \Theta_{t} by updating its components: A_t \sim p(A | B = b_{t-1}) \qquad B_t \sim p(B | A = a_t)
  • The stationary distribution for this Markov chain is p(A, B)
  • In other words, we are alternating sampling from each conditional distribution to obtain (with large t) a (correlated) sample from the joint distribution

Gibbs sampling: toy example

  • Suppose we want to extract samples from the joint distribution N\left(\begin{bmatrix} x_1 \\ x_2\end{bmatrix}; \begin{bmatrix} 0 \\ 0\end{bmatrix}, \begin{bmatrix}1 & \rho \\ \rho & 1 \end{bmatrix}\right)
  • We know how to extract samples from p(X_1 | X_2=x_2) and p(X_2 | X_1=x_1):

\begin{align*} p(X_1 | X_2=x_2) &= N(x_1; \rho x_2, 1-\rho^2) \\ p(X_2 | X_1=x_1) &= N(x_2; \rho x_1, 1-\rho^2) \end{align*}

set.seed(696)

rho <- 0.6

gibbs_step <- function(x){
  x[1] <- rnorm(1, rho * x[2], sqrt( 1-rho^2 ) )
  x[2] <- rnorm(1, rho * x[1], sqrt( 1-rho^2 ) )
  return(x)
}

num_samples <- 10000
x_all_samples <- matrix(NA, ncol=2, nrow=num_samples)
colnames(x_all_samples) <- c("X1", "X2")

x <- c(0,0) # initial value

# run the Gibbs sampler
for(i in 1:num_samples){
  x <- gibbs_step(x)
  x_all_samples[i,] <- x
}

Gibbs sampling: toy example

Finding p(\beta, \tau^2, \theta | Y)

We know how to get

\{ \tau^2_{(t)}, \beta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \tau^2_{(t)}, \beta_{(t)} \sim p(\beta, \tau^2 | Y, \theta)

We “know” how to get

\{ \theta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \theta_{(t)} \sim p(\theta | Y, \beta, \tau^2)

Therefore, we build a Gibbs sampler and get

\{ \tau^2_{(t)}, \beta_{(t)}, \theta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \tau^2_{(t)}, \beta_{(t)}, \theta_{(t)} \sim p(\beta, \tau^2, \theta | Y)

We learn p(\beta, \tau^2, \theta | Y) by summarizing samples from it.

Marginal posteriors:

Need samples from p(\beta | Y), p(\tau^2 | Y), p(\theta | Y)?

Just pick the element we need, e.g. extract \{ \beta_{(t)} \}_{t=1}^T from \{ \tau^2_{(t)}, \beta_{(t)}, \theta_{(t)} \}_{t=1}^T for a sample from p(\beta | Y)

With samples from p(\beta, \tau^2, \theta | Y)

One important thing we still need is p(w | Y), the posterior for the spatial process

\begin{align*} p(w | Y) &= \int p(\beta, \tau^2, \theta, w | Y) d\beta d\tau^2 d\theta \\ &= \int p(w | Y, \beta, \tau^2, \theta) p(\beta, \tau^2, \theta | Y) d\beta d\tau^2 d\theta \end{align*}

  • This means that for each sample of \tau^2_{(t)}, \beta_{(t)}, \theta_{(t)}, all we need to do is independently sample from p(w | Y, \beta, \tau^2, \theta) to obtain a sample from the joint posterior.
  • p(w | Y, \beta, \tau^2, \theta) is our next target
  • It will be simple!

Finding p(w | Y, \beta, \tau^2, \theta)

  • Before we start, take a moment to notice what we are trying to do
  • We want a posterior (we condition on Y) but we take everything else as known

Prior: w | \theta \sim N(0, C_{\theta})

Model: Y = X\beta + w + \varepsilon,\quad \varepsilon \sim N(0, \tau^2 I_n)

Rewrite the model: Y^* = I_n w + \varepsilon \qquad \text{where}\quad Y^* = Y-X\beta

(full conditional) Posterior

How do we continue?

Finding p(w | Y, \beta, \tau^2, \theta)

  • Before we start, take a moment to notice what we are trying to do
  • We want a posterior (we condition on Y) but we take everything else as known

Prior: w | \theta \sim N(0, \tau^2 C_{\theta})

Model: Y = X\beta + w + \varepsilon,\quad \varepsilon \sim N(0, \tau^2 I_n)

Rewrite the model: Y^* = I_n w + \varepsilon \qquad \text{where}\quad Y^* = Y-X\beta

(full conditional) Posterior

This is just the usual linear regression. Here, w plays the role of \beta. It is a Gaussian model with Gaussian prior, so the posterior is Gaussian (conjugate):

\begin{align*} &w | Y, \beta, \tau^2, \theta \sim N(m_w, \tau^2 \Sigma_w)\\ &\Sigma_w^{-1} = C_{\theta}^{-1} + I_n\\ &m_w = \Sigma_w (Y - X\beta) \end{align*}

Predictions

  • The last thing we are left to figure out is how to make predictions at new locations
  • We know how to do that for a GP response model (How?), what about now?
  • We are looking for Y(s^*) where s^* was not originally in our dataset

Y(s^*) \sim N( X(s^*)^\top \beta + w(s^*), \tau^2)

  • First, note that p(w(s^*) | w, Y) = p(w(s^*) | w)
  • p(w(s^*) | w) is the GP prediction as we saw it here
  • With samples \tau^2_{(t)}, \beta_{(t)}, \theta_{(t)}, w_{(t)}, we just need to compute p(w(s^*) | w) and then extract from Y(s^*) \sim N( X(s^*)^\top \beta + w(s^*), \tau^2)

Why? Because:

\begin{align*} p(Y^*, w^*, w, \beta, \tau^2, \theta | Y) &= p(Y^*, w^* | Y, w, \beta, \tau^2, \theta) p(w, \beta, \tau^2, \theta | Y) \\ &= p(Y^* | w^*, Y, w, \beta, \tau^2, \theta) p(w^* | Y, w, \beta, \tau^2, \theta) p(w, \beta, \tau^2, \theta | Y) \\ &= p(Y^* | w^*, \beta, \tau^2) p(w^* | w, \theta)p(w, \beta, \tau^2, \theta | Y) \end{align*}

Back to the DAG

Back to the DAG

We’re done!

  • We have outlined all steps for estimating parameters and making predictions
  • The basic univariate linear spatial model can be extended in many ways
  • Make it a latent-Gaussian model of non-Gaussian data (spatial GLMM)
  • The model we have considered is also referred to as a spatially-varying intercept model
  • Assume Z(s) \in \Re is a covariate not already in X(s). Then,

Y(s) = X(s) \beta + Z(s) w(s) + \varepsilon(s) is a spatially-varying coefficient (SVC) model. The regression coefficient that explains how Z(s) is associated to Y(s) is a spatial process (basically, we can draw a map of it) rather than a constant

But we are not done!!!

We have not explained how we are getting samples

\{ \theta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \theta_{(t)} \sim p(\theta | Y, \beta, \tau^2)

For that, we need to outline the Metropolis-Hastings algorithm.

Metropolis-Hastings algorithm

  • Goal: obtain samples from a target density p(x)
  • Idea: make a Markov chain that has p(x) as its stationary distribution

Iteration t of the MH algorithm

  • Given the current value x_{(t-1)}, propose a new value x^* by sampling it from density q(\cdot | x_{(t-1)})
  • Compute (by evaluating the densities) \alpha = \frac{p(x^*)}{p(x_{(t-1)})} \frac{q(x_{(t-1)}| x^*)}{q(x^* | x_{(t-1)})}
  • Accept and therefore set x_{(t)} = x^* with probability \min\{1, \alpha\}
  • For example, sample U \sim U(0,1) and accept if U < \alpha, reject otherwise
  • We have built a Markov chain, and it can be shown to converge to p(\cdot)
  • Output: (correlated) samples \{ x_{(t)} \}_{t=1}^T

In our case

  • Our target is p(\theta | Y, \beta, \tau^2) \propto p(Y | \theta, \beta, \tau^2) p(\theta)
  • Because we are taking ratios of this target, the normalizing constant cancels out
  • This is very convenient: MH works even if we can only evaluate the target density up to a proportionality constant

Metropolis-Hastings algorithm

  • Goal: obtain samples from p(\theta | Y, \beta, \tau^2) \propto p(Y | \theta, \beta, \tau^2) p(\theta)
  • Idea: make a Markov chain that has p(Y | \theta, \beta, \tau^2) p(\theta) as its stationary distribution

Iteration t of the MH algorithm

  • Given the current value \theta_{(t-1)}, propose a new value \theta^* by sampling it from density q(\cdot | \theta_{(t-1)})

  • Compute (by evaluating the densities) \alpha = \frac{p(Y | \theta^*, \beta, \tau^2) p(\theta^*)}{p(Y | \theta_{(t-1)}, \beta, \tau^2) p(\theta_{(t-1)})} \frac{q(\theta_{(t-1)}| \theta^*)}{q(\theta^* | \theta_{(t-1)})}

  • Accept and therefore set \theta_{(t)} = \theta^* with probability \min\{1, \alpha\}

  • For example, sample U \sim U(0,1) and accept if U < \alpha, reject otherwise

  • Output: (correlated) samples \{ \theta_{(t)} \}_{t=1}^T from p(\theta | Y, \beta, \tau^2)

Metropolis-Hastings algorithm & Gibbs samplers

  • A Gibbs sampler assumes we can sample from all the full conditional distributions
  • We would then alternate sampling from each of them to obtain a correlated sample from the joint posterior
  • If even one of the full conditionals is not available, we can replace that step with a Metropolis-Hastings step
  • The resulting Metropolis-within-Gibbs algorithm is one of the most important algorithms in Bayesian statistics because of its extremely wide applicability
  • The Gibbs sampler itself is a special case of Metropolis-Hastings
  • Many generalizations of Metropolis-Hastings (adaptive MH, Hamiltonian MC, …)

Back to finding p(\beta, \tau^2, \theta | Y)

We use a conditionally conjugate step to get

\{ \tau^2_{(t)}, \beta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \tau^2_{(t)}, \beta_{(t)} \sim p(\beta, \tau^2 | Y, \theta)

We use Metropolis-Hastings to get

\{ \theta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \theta_{(t)} \sim p(\theta | Y, \beta, \tau^2)

Therefore, we build a Metropolis-within-Gibbs sampler and get

\{ \tau^2_{(t)}, \beta_{(t)}, \theta_{(t)} \}_{t=1}^T \qquad \text{where} \qquad \tau^2_{(t)}, \beta_{(t)}, \theta_{(t)} \sim p(\beta, \tau^2, \theta | Y)

We learn p(\beta, \tau^2, \theta | Y) by summarizing samples from it.

Alternative sampling algorithms

  • In the previous discussion, we have marginalized out w
  • This is convenient because w is high dimensional
  • The “remaining” joint posterior is lower dimensional
  • Lower dimension = easier to explore and learn
  • We could also integrate out \beta (R package spBayes does this)

Non-marginalized sampler

  • What if we do not integrate out w?

  • We are then looking at the joint posterior p(w, \beta, \tau^2, \theta | Y)

  • Here, a Gibbs sampler (with a Metropolis step) works this way:

    • Sample from p(\tau^2 | Y, w, \beta, \theta)
    • Sample from p(w | Y, \beta, \tau^2, \theta)
    • Sample from p(\beta | Y, w, \tau^2, \theta)
    • Update \theta with a Metropolis step targeting p(\theta | Y, w, \tau^2, \beta)

Non-marginalized sampler

  • What if we do not integrate out w?

  • We are then looking at the joint posterior p(w, \beta, \tau^2, \theta | Y)

  • Here, a Gibbs sampler (with a Metropolis step) works this way:

    • Sample from p(\tau^2 | Y, w, \beta, \theta) = p(\tau^2 | Y, w, \beta) why?
    • Sample from p(w | Y, \beta, \tau^2, \theta)
    • Sample from p(\beta | Y, w, \tau^2, \theta) = p(\beta | Y, w, \tau^2) why?
    • Update \theta with a Metropolis step targeting p(\theta | Y, w, \tau^2, \beta) = p(\theta | w) why?

Non-marginalized sampler

  • What if we do not integrate out w?

  • We are then looking at the joint posterior p(w, \beta, \tau^2, \theta | Y)

  • Here, a Gibbs sampler (with a Metropolis step) works this way:

    • Sample from p(\tau^2 | Y, w, \beta)
    • Sample from p(w | Y, \beta, \tau^2, \theta)
    • Sample from p(\beta | Y, w, \tau^2)
    • Update \theta with a Metropolis step targeting p(\theta | w) \propto p(w | \theta) p(\theta) what is p(w | \theta)?

Non-marginalized sampler, step 1

  • Sample from p(\tau^2 | Y, w, \beta)
  • We condition on everything else, so Y^* = Y-X\beta - w can be computed
  • Treat Y^* as the dataset
  • Then the prior-model pair is

\begin{align*} \tau^2 &\sim Inv.G(a_{\tau}, b_{\tau})\\ Y^* &\sim N(0, \tau^2 I_n) \end{align*}

  • How do we find the (full conditional) posterior in this case?

Non-marginalized sampler, step 1

  • Sample from p(\tau^2 | Y, w, \beta)
  • We condition on everything else, so Y^* = Y-X\beta - w can be computed
  • Treat Y^* as the dataset
  • Then the prior-model pair is

\begin{align*} \tau^2 &\sim Inv.G(a_{\tau}, b_{\tau})\\ Y^* &\sim N(0, \tau^2 I_n) \end{align*}

  • In the usual linear regression model, drop everything related to \beta
  • Inverse Gamma full conditional

Non-marginalized sampler, step 2

  • Sample from p(w | Y, \beta, \tau^2, \theta)
  • Treat Y^* = Y - X\beta as the dataset
  • Then the prior-model pair is

\begin{align*} w &\sim N(0, C_{\theta})\\ Y^* &\sim N(I_n w, \tau^2 I_n) \end{align*}

  • How do we find the (full conditional) posterior in this case?

Non-marginalized sampler, step 2

  • Sample from p(w | Y, \beta, \tau^2, \theta)
  • Treat Y^* = Y - X\beta as the dataset
  • Then the prior-model pair is

\begin{align*} w &\sim N(0, C_{\theta})\\ Y^* &\sim N(I_n w, \tau^2 I_n) \end{align*}

  • This is a regression model with known error variance
  • Gaussian prior, Gaussian model, therefore Gaussian posterior
  • p(w | Y, \beta, \tau^2, \theta) = N(w; \mu_w, \Sigma_w) where

\begin{align*} \Sigma_w^{-1} &= C_{\theta}^{-1} + 1/\tau^2\\ \mu_w &= \Sigma_w (Y - X\beta)/\tau^2 \end{align*}

Non-marginalized sampler, step 3

  • Sample from p(\beta | Y, w, \tau^2)
  • Treat Y^* = Y - w as the dataset
  • Then the prior-model pair is

\begin{align*} \beta &\sim N(0, V)\\ Y^* &\sim N(X\beta, \tau^2 I_n) \end{align*}

  • How do we find the (full conditional) posterior in this case?

Non-marginalized sampler, step 3

  • Sample from p(\beta | Y, w, \tau^2)
  • Treat Y^* = Y - w as the dataset
  • Then the prior-model pair is

\begin{align*} \beta &\sim N(0, V)\\ Y^* &\sim N(X\beta, \tau^2 I_n) \end{align*}

  • This is a regression model with known error variance
  • Gaussian prior, Gaussian model, therefore Gaussian posterior
  • p(\beta | Y, w, \tau^2, \theta) = N(\beta; \mu_{\beta}, \Sigma_{\beta}) where

\begin{align*} \Sigma_{\beta}^{-1} &= V^{-1} X^\top X/\tau^2\\ \mu_{\beta} &= \Sigma_{\beta} X^\top (Y - w)/\tau^2 \end{align*}

Non-marginalized sampler, step 4

  • Update p(\theta | w) \propto p(w | \theta) p(\theta)
  • p(w | \theta) = N(w; 0, C_{\theta})

Alternative methods

  • Because Y = X\beta + w + \varepsilon can be written as Y = \begin{bmatrix} X & I_n\end{bmatrix} \begin{bmatrix} \beta \\ w \end{bmatrix} + \varepsilon

  • We can call X^* = \begin{bmatrix} X & I_n\end{bmatrix} and \beta^* = \begin{bmatrix} \beta \\ w \end{bmatrix} and then update as a block via

p(\beta^* | Y, \tau^2, \theta) = N(\beta^*; \mu_{\beta^*}, \Sigma_{\beta^*})

  • \Sigma_{\beta^*}^{-1} = V_{\beta^*}^{-1} + X^{*^\top} X^*/\tau^2
  • \mu_{\beta^*} = \Sigma_{\beta^*} X^{*^\top} Y/\tau^2
  • Once again, this follows from the linear regression posterior we saw in the preliminaries

Using package spBayes for spatial linear regression with latent GPs

Using package spBayes for spatial linear regression with latent GPs

  • First, specify how we want to set things up.
  • The model we are implementing is

Y = X\beta + w + \varepsilon, \qquad \text{where}\quad \varepsilon \sim N(0, \tau^2 I_n)

Priors:

  • \beta \sim N(0, 10^3 I_p) where p is the number of columns of X
  • w(\cdot) \sim GP(0, C_{\sigma^2, \phi}(\cdot)) and we let C_{\sigma^2, \phi}(s, s') = \sigma^2 \exp\{ -\phi \| s-s' \| \}
  • \tau^2 \sim Inv.G(2, 0.1) and \sigma^2 \sim Inv.G(2,2)
  • \phi \sim U(3, 30). It is common to let \phi have a uniform prior where the range of possible values is set to give “reasonable” spatial dependence
library(spBayes)

n.samples <- 2000

starting <- list("phi"=3/0.5, "sigma.sq"=50, "tau.sq"=1)
tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1)

priors <- list("beta.Norm"=list(rep(0,ncol(X)), diag(1000,ncol(X))),
                 "phi.Unif"=c(3/1, 3/0.1), "sigma.sq.IG"=c(2, 2),
                 "tau.sq.IG"=c(2, 0.1))

cov.model <- "exponential"

Using package spBayes for spatial linear regression with latent GPs

We fit the GP regression model in 2 steps

system.time({
  model_fit <- spLM(Y~X-1, coords=coords, starting=starting,
            tuning=tuning, priors=priors, cov.model=cov.model,
            n.samples=n.samples, verbose=FALSE, n.report=500) 
  })
   user  system elapsed 
  7.602   0.010   7.629 
system.time({
  model_fit <- model_fit %>% spRecover(start=round(n.samples/3), verbose=FALSE)
})
   user  system elapsed 
 35.201   0.913  25.636

Using package spBayes for spatial linear regression with latent GPs

  • We chose to run MCMC for 2000 iterations
  • If we are satisfied with convergence, we drop some initial samples (burn-in period) because they are likely not sampled at convergence
  • Recall with MCMC we get correlated samples from the joint posterior after reaching convergence

Using package spBayes for spatial linear regression with latent GPs

  • spBayes calls theta the vector of (\sigma^2, \tau^2, \phi)
  • In this example, the MCMC chains look poor
  • All chains mix poorly, or in other terms they show high autocorrelation
  • Inference based on these chains should not be trusted
  • We need to improve MCMC tuning and possibly run the Markov chain longer

Using package spBayes for spatial linear regression with latent GPs

  • At each step, a new proposed value for \phi is generated
  • tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1) sets the standard deviation of the proposal for phi to 0.1
  • In other words, this determines the jump size
  • If the jump is too big, the proposal will probably get rejected
  • In our case: 16.4% acceptance rate
  • We target an acceptance rate of about 24% so we should lower tuning$phi

Using package spBayes for spatial linear regression with latent GPs

  • Let’s take the previous MCMC run (even though we did not like it too much)
  • We can now predict at a grid of locations
predictions <- spPredict(model_fit,pred.coords = coords_out, pred.covars = X_out)
----------------------------------------
    General model description
----------------------------------------
Model fit with 1000 observations.

Prediction at 900 locations.

Number of covariates 3 (including intercept if specified).

Using the exponential spatial correlation model.

-------------------------------------------------
        Sampling
-------------------------------------------------
Sampled: 100 of 2000, 4.95%
Sampled: 200 of 2000, 9.95%
Sampled: 300 of 2000, 14.95%
Sampled: 400 of 2000, 19.95%
Sampled: 500 of 2000, 24.95%
Sampled: 600 of 2000, 29.95%
Sampled: 700 of 2000, 34.95%
Sampled: 800 of 2000, 39.95%
Sampled: 900 of 2000, 44.95%
Sampled: 1000 of 2000, 49.95%
Sampled: 1100 of 2000, 54.95%
Sampled: 1200 of 2000, 59.95%
Sampled: 1300 of 2000, 64.95%
Sampled: 1400 of 2000, 69.95%
Sampled: 1500 of 2000, 74.95%
Sampled: 1600 of 2000, 79.95%
Sampled: 1700 of 2000, 84.95%
Sampled: 1800 of 2000, 89.95%
Sampled: 1900 of 2000, 94.95%
Sampled: 2000 of 2000, 99.95%

Using package spBayes for spatial linear regression with latent GPs

  • Calculate residuals from the model fit
w_post <- model_fit$p.w.recover.samples
b_post <- model_fit$p.beta.recover.samples
XB_post <- X %*% t(b_post)
XBW_mean <- (XB_post + w_post) %>% apply(1, mean)
model_residual <- Y - XBW_mean

df <- df %>% mutate(residual = model_residual)
residual_map <- ggplot(df, aes(xcoord, ycoord, color=residual)) +
  geom_point() + 
  scale_color_viridis_c() +
  theme_minimal() +
  ggtitle("Residual")

sv <- df %>% with(geoR::variog(data = residual, coords = cbind(xcoord, ycoord), messages=FALSE))

sv_df <- data.frame(dists = sv$u, variogram = sv$v, npairs = sv$n, sd = sv$sd)
sv_plot <- ggplot(sv_df, aes(x=dists, y=variogram)) + geom_point(size=2, shape=8) +
  theme_minimal()

grid.arrange(residual_map, sv_plot, nrow=1)

Why does spBayes need two steps

Priors:

  • \tau^2 \sim Inv.G(a_{\tau}, b_{\tau})
  • \beta \sim N(0, \tau^2 V)
  • s^2 \sim Inv.G(a_{s}, b_s)
  • \phi \sim U(l_{\phi}, u_{\phi})
  • w(\cdot) \sim GP(0, C_{s, \phi}(\cdot)) where C_{s, \phi}(s, s') = \sigma^2 \exp\{ -\phi \| s-s' \| \}

Model:

Y(s) = X(s)^\top \beta + w(s) + \varepsilon(s),\ \text{where}\ \varepsilon(s) \sim N(0, \tau^2) \quad \text{(iid)} In vector form: Y = X\beta + w + \varepsilon, where \varepsilon \sim N(0, \tau^2 I_n).

Marginalized model:

  • In addition to integrating out w as we have seen before, we can also do that with \beta
  • In spBayes, this is the marginal model that is fit via MCMC

Y | \tau^2, s^2, \phi \sim N(0, \tau^2 (X V X^\top + I_n) + C_{\sigma^2, \phi})

Why does spBayes need two steps

Marginalized model:

  • In addition to integrating out w as we have seen before, we can also do that with \beta
  • In spBayes, this is the marginal model that is fit via MCMC

Y | \tau^2, s^2, \phi \sim N(0, \tau^2 (X V X^\top + I_n) + C_{\sigma^2, \phi})

  • It is possible to update \tau^2 via a Gibbs step or jointly with \phi and \sigma^2 as a Metropolis update
  • The Metropolis-Hastings algorithm targets the posterior p(\tau^2, \sigma^2, \phi | Y) \propto p(Y | \tau^2, \sigma^2, \phi) p(\tau^2, \sigma^2, \phi)
  • Once we have samples from p(\tau^2, \sigma^2, \phi | Y) we can sample from p(\beta | Y, \tau^2, \sigma^2, \phi)
  • Once we have samples from p(\tau^2, \sigma^2, \phi | Y) we can sample from p(w | Y, \tau^2, \sigma^2, \phi)
  • Both of these steps are performed in spRecover

Moving forward

  • spBayes::spLM fits a GP regression model
  • Computations take a considerable amount of time for n=1000
  • Datasets with n=10000 or larger are quite common
  • Let’s see what happens with spLM at increasing sample sizes

Moving forward

  • Things do not look promising
  • Next, we will look into why this is the performance profile we get
  • We will then consider some solutions to the scalability problem

Home page