---
title: "Gaussian Processes models for point-referenced spatial data, part 2"
subtitle: "Spatial Statistics - BIOSTAT 696/896"
author: "Michele Peruzzi"
institute: "University of Michigan"
execute:
  cache: true
format:
  revealjs:
    auto-stretch: true
    center: true
    html-math-method: katex
    transition: none
    slide-number: true
    callout-appearance: minimal
    code-line-numbers: true
    fontsize: 16pt
    output-file: l05_gaussianprocesses2.html
    theme: [simple, custom.css] #, web/pulse_colors.css] # alternative themes (subset): default, night, dark
    embed-resources: false
    echo: true
    fig-dpi: 150
    # incremental: true  # Remove comment if you like incremental bullet points
    footer: "[Home page](index.html)"
knitr:
  opts_chunk:
    out.width: "60%"
editor: 
  markdown: 
    wrap: 72
editor_options: 
  chunk_output_type: console
---

```{r echo=FALSE, cache=FALSE}
library(sf)
library(rnaturalearth)
library(tidyverse)
library(gridExtra)
library(magrittr)
library(scico)
```

## Likelihood inference: MLE estimation of covariance parameters

-  Let's consider the Greenness data
-  Small preprocessing: remove sample mean

```{r echo=FALSE}
df <- read.csv("data/michigan_greenness.csv")
df %<>% mutate(Greenness = Greenness-mean(Greenness))

mainmap <- ne_states(country = c("united states of america", "canada"), returnclass = "sf")
michiganplus <- mainmap %>% dplyr::filter(name %in% c("Illinois", "Michigan", "Wisconsin", "Ontario"))

  (data_plot <- ggplot() +
    geom_sf(data = michiganplus, fill = "white", color = "black") +
      coord_sf(ylim = c(42, 46), xlim = c(-87, -82)) +
      geom_raster(data = df, aes(x=Longitude, y=Latitude, fill=Greenness)) +
      scale_fill_scico(palette="batlowK") + 
      theme_minimal())
```

## Likelihood inference: MLE estimation of covariance parameters

-   Call $\cal T$ our set of observed locations. The data vector at $\cal T$ is $Y_{\cal T}$, or simply $Y$
-   Let's use an exponential covariance model for a mean-zero GP
-   Because we are directly assuming that the data are a realization of a GP, this is called a GP model of the *response* or response model
$$Y(\cdot) \sim GP(0, C_{\theta}(\cdot)) $$
For our $n$-dimensional vector of observations, we have
$$Y \sim N(0, C_{\theta}) $$ where the $(i,j)$ element of the $C_{\theta}$ matrix is $C_{\theta}(i,j) = \sigma^2 \exp \{ -\phi \|s_i - s_j \| \}$ where $s_i, s_j \in \cal T$.


## Likelihood inference: MLE estimation of covariance parameters

-  Recall from preliminaries that because of your MVN assumption, we can write

$$\log p(Y \mid \theta) = -\frac{n}{2} \log(2\pi) - \sum_{i=1}^n\log(L_{ii}) - \frac{1}{2}Y^T L^{-T} L^{-1} Y$$

- [The matrix $L$ is... ?]{.orange}
- Our goal is to find the MLE of $\theta$

$$\hat{\theta}_{MLE} = \arg\max_{\theta} \log p(Y \mid \theta)$$


## Likelihood inference: MLE estimation of covariance parameters

-  Recall from preliminaries that because of your MVN assumption, we can write

$$\log p(Y \mid \theta) = -\frac{n}{2} \log(2\pi) - \sum_{i=1}^n\log(L_{ii}) - \frac{1}{2}Y^T L^{-T} L^{-1} Y$$

- The matrix $L$ is the Cholesky decomposition of $C_{\theta}$, i.e. $C_{\theta} = L L^\top$
- $L$ is where $\theta$ appears in the likelihood function
- Our goal is to find the MLE of $\theta$

$$\hat{\theta}_{MLE} = \arg\max_{\theta} \log p(Y \mid \theta)$$

## Likelihood inference: MLE estimation of covariance parameters

-  Recall from preliminaries that because of your MVN assumption, we can write

$$\log p(Y \mid \theta) = -\frac{n}{2} \log(2\pi) - \sum_{i=1}^n\log(L_{ii}) - \frac{1}{2}Y^T L^{-T} L^{-1} Y$$

- The matrix $L$ is the Cholesky decomposition of $C_{\theta}$, i.e. $C_{\theta} = L L^\top$
- $L$ is where $\theta$ appears in the likelihood function
- Our goal is to find the MLE of $\theta$

$$\hat{\theta}_{MLE} = \arg\min_{\theta} \left( \sum_{i=1}^n\log(L_{ii}) + \frac{1}{2}Y^T L^{-T} L^{-1} Y \right)$$


## Likelihood inference: MLE estimation of covariance parameters

-  Let's do this in R in a simple way

```{r}
set.seed(696)

coords <- df %>% dplyr::select(Longitude, Latitude) %>% as.matrix()
Y <- df$Greenness %>% as.vector()

Covariance_make <- function(theta){
  sigmasq <- exp(theta[1])
  phi <- exp(theta[2])
  Cmatrix <- sigmasq * exp(-phi * as.matrix(dist(coords)) )
  return(Cmatrix)
}

mvn_eval <- function(theta){
  Ctheta <- Covariance_make(theta)
  Ltheta <- t(chol(Ctheta))
  Linv_Y <- solve(Ltheta, Y) # = L^-1 Y
  
  neg_logdens <- sum(log(diag(Ltheta))) + 0.5 * crossprod(Linv_Y)
  return(neg_logdens)
}

(save_time <- system.time({ results <- optim(c(0,0), mvn_eval) }))
```

## Likelihood inference: MLE estimation of covariance parameters

```{r echo=FALSE}
results
hatssq <- exp(results$par[1])
hatphi <- exp(results$par[2])
```

-  `convergence=0` means the optimization algorithm found a (local) minimum
-  Because of the way we parametrized, $\hat{\sigma}^2=`r hatssq`$ and $\hat{\phi} =`r hatphi`$
-  Let's try different starting values:

```{r}
results2 <- optim(c(-1,1), mvn_eval)
```

```{r echo=FALSE}
results2
hatssq <- exp(results2$par[1])
hatphi <- exp(results2$par[2])
```

-  $\hat{\sigma}^2=`r hatssq`$ and $\hat{\phi} =`r hatphi`$

## Key take-aways from the previous slides

### Scalability of likelihood-based inference of GP models

-  With $n= `r nrow(df)`$ observations the optimization routine takes about 40 seconds for `r results$counts["function"]` evaluations of the likelihood function
-  That is $\approx 0.5$s/iteration: *what if we had many more observations?*
-  Presumably, optimizing a more complex covariance function will require more likelihood evaluations

### Covariance estimation (aka kernel learning)

-  We changed the starting values and we got slightly different optima
-  The likelihood values are very similar
-  In general, GP models have relatively *flat* likelihoods
-  Many different values of elements of $\theta$ lead to very similar $p(Y | \theta)$
-  This will be one justification for using a *Bayesian model* in which we include prior information on the parameters


## Making predictions based on a GP response model

-  We have a point-estimate $\hat{\theta}$ for our covariance parameters
-  In other words, we have found what values of $\theta$ in our covariance model best fit our data
-  We do not have any other unknowns in our model
-  [How do we make predictions?]{.orange}


## Making predictions based on a GP response model

Back to our theory of Gaussian Processes:

-  Our GP assumption says that *any finite set of spatially-indexed random variables are jointly Gaussian*
-  We have $\cal T$ as our observed locations, $Y_{\cal T}$ the realization of a GP at those locations
-  Suppose we now have $\cal U$ as the set of locations where we want to make predictions, i.e. $Y_{\cal U}$ is unobserved
-  $\cal T$ and $\cal U$ are both *finite* sets of locations, hence

$$\begin{bmatrix}Y_{\cal T} \\ Y_{\cal U}\end{bmatrix} \sim \text{????}$$

## Making predictions based on a GP response model

Back to our theory of Gaussian Processes:

-  Our GP assumption says that *any finite set of spatially-indexed random variables are jointly Gaussian*
-  With locations $\cal L$, we are saying $Y_{\cal L} \sim N(0, C_{\cal L})$ 
-  Key: *for any finite set* $\cal L$

-  In our applied problem we have $\cal T$ as our observed locations, $Y_{\cal T}$ the realization of a GP at those locations
-  Suppose we now have $\cal U$ as the set of locations where we want to make predictions, i.e. $Y_{\cal U}$ is unobserved
-  $\cal T$ and $\cal U$ are both *finite* sets of locations, hence we can build $\cal L = \cal T \cup \cal U$, and then:

$$Y_{\cal L} = \begin{bmatrix}Y_{\cal T} \\ Y_{\cal U}\end{bmatrix} \sim N\left( 0, \begin{bmatrix} C_{\cal T} & C_{\cal T, \cal U} \\ C_{\cal T, \cal U}^ \top & C_{\cal U} \end{bmatrix} \right) = N(0, C_{\cal L})$$

- [We should now be able to make predictions. How do we continue?]{.orange}


## Making predictions based on a GP response model

$$Y_{\cal L} = \begin{bmatrix}Y_{\cal T} \\ Y_{\cal U}\end{bmatrix} \sim N\left( 0, \begin{bmatrix} C_{\cal T} & C_{\cal T, \cal U} \\ C_{\cal T, \cal U}^ \top & C_{\cal U} \end{bmatrix} \right) = N(0, C_{\cal L})$$

We use the joint density above to find the conditional density below:

$$Y_{\cal U} | Y_{\cal T} \sim N\left( m_{\cal U|\cal T}, C_{\cal U|\cal T} \right)$$
where
\begin{align*}
m_{\cal U|\cal T} &= C_{\cal U, \cal T} C_{\cal T}^{-1} Y_{\cal T}  \\
C_{\cal U|\cal T} &= C_{\cal U} - C_{\cal U, \cal T} C_{\cal T}^{-1} C_{\cal T, \cal U}
\end{align*}
The conditional mean is our point prediction at new locations; the conditional variance is giving us uncertainty at the new locations
\begin{align*}
m_{\text{new}|\text{obs}} &= C_{\text{new}, \text{obs}} C_{\text{obs}}^{-1} Y_{\text{obs}}  \\
C_{\text{new}|\text{obs}} &= C_{\text{new}} - C_{\text{new}, \text{obs}} C_{\text{obs}}^{-1} C_{\text{obs}, \text{new}}
\end{align*}


## Making predictions based on a GP response model

```{r echo=FALSE, cache=FALSE}

mycppcode <- "
//[[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

//[[Rcpp::export]]
arma::mat Cf(const arma::mat& cx, const arma::mat& cy, const arma::vec& param){
  double sigmasq = param(0);
  double phi = param(1);
  
  arma::mat result = arma::zeros(cx.n_rows, cy.n_rows);
  for(unsigned int i=0; i<cx.n_rows; i++){
    for(unsigned int j=0; j<cy.n_rows; j++){
      arma::rowvec h = cx.row(i) - cy.row(j);
      double distance = arma::conv_to<double>::from( sqrt(h * arma::trans(h)) );
      result(i, j) = sigmasq * exp(- phi * distance);
    }
  }
  return result;
}
"
Rcpp::sourceCpp(code=mycppcode)

pred_points <- expand.grid(list(Longitude=seq(-86.6, -82.5, length.out=50),
                           Latitude=seq(41.8, 46, length.out=50))) 
outcoords <- as.matrix(pred_points)

(data_plot <- ggplot() +
  geom_sf(data = michiganplus, fill = "white", color = "black") +
    coord_sf(ylim = c(42, 46), xlim = c(-87, -82)) +
    geom_raster(data = df, aes(x=Longitude, y=Latitude, fill=Greenness), alpha=.6) +
    geom_point(data=pred_points, aes(x=Longitude, y=Latitude), color="magenta", size=.4) +
    scale_fill_scico(palette="batlowK") + 
    theme_minimal())

m_hat <- Cf(outcoords, coords, c(hatssq, hatphi)) %*% solve(Cf(coords, coords, c(hatssq, hatphi)), df$Greenness)
```

## Making predictions based on a GP response model

```{r echo=FALSE}
df_preds <- data.frame(pred_points, Predictions=m_hat)

(pred_plot <- ggplot() +
  geom_sf(data = michiganplus, fill = "white", color = "black") +
    coord_sf(ylim = c(42, 46), xlim = c(-87, -82)) +
    geom_raster(data = df_preds, aes(x=Longitude, y=Latitude, fill=Predictions)) +
    scale_fill_scico(palette="batlowK") + 
    theme_minimal())
```


## Making predictions based on a GP response model

```{r echo=FALSE}

pred_points <- expand.grid(list(Longitude=seq(-86.6, -82.5, length.out=100),
                           Latitude=seq(41.8, 46, length.out=100))) 
outcoords <- as.matrix(pred_points)

(data_plot <- ggplot() +
  geom_sf(data = michiganplus, fill = "white", color = "black") +
    coord_sf(ylim = c(42, 46), xlim = c(-87, -82)) +
    geom_raster(data = df, aes(x=Longitude, y=Latitude, fill=Greenness), alpha=.6) +
    geom_point(data=pred_points, aes(x=Longitude, y=Latitude), color="magenta", size=.2) +
    scale_fill_scico(palette="batlowK") + 
    theme_minimal())

m_hat <- Cf(outcoords, coords, c(hatssq, hatphi)) %*% solve(Cf(coords, coords, c(hatssq, hatphi)), df$Greenness)
```

## Making predictions based on a GP response model

```{r echo=FALSE}
df_preds <- data.frame(pred_points, Predictions=m_hat)

(pred_plot <- ggplot() +
  geom_sf(data = michiganplus, fill = "white", color = "black") +
    coord_sf(ylim = c(42, 46), xlim = c(-87, -82)) +
    geom_raster(data = df_preds, aes(x=Longitude, y=Latitude, fill=Predictions)) +
    scale_fill_scico(palette="batlowK") + 
    theme_minimal())
```


## Making predictions based on a GP response model

-  Uncertainty of predictions: we can use 95% credible bands point-wise. At 1 new location $s^ *$, we have

$$y(s^*) | Y_{\cal T} \sim N\left( m_{\text{pred}}, C_{\text{pred}} \right)$$
where
\begin{align*}
m_{\text{pred}} &= C(s^*, {\cal T}) C_{\cal T}^{-1} Y_{\cal T}  \\
C_{\text{pred}} &= \sigma^2 - C(s^*, {\cal T}) C_{\cal T}^{-1} C({\cal T}, s^*)
\end{align*}

-  Therefore $U_{95} = \Phi^{-1}(0.975, m_{\text{pred}}, C_{\text{pred}}) - \Phi^{-1}(0.025, m_{\text{pred}}, C_{\text{pred}})$
-  $\Phi^{-1}$ is the inverse CDF of the (univariate) Gaussian. In R, `qnorm`.
-  Because $\cal U$ can be large, applying the above formula for each element of $\cal U$ saves on memory relative to extracting the diagonal elements from $C_{\cal U | \cal T}$, but that would also work

## Making predictions based on a GP response model

```{r echo=FALSE}
Cband_hat <- Cf(outcoords, outcoords, c(hatssq, hatphi)) - Cf(outcoords, coords, c(hatssq, hatphi)) %*%
  solve(Cf(coords, coords, c(hatssq, hatphi)), Cf(coords, outcoords, c(hatssq, hatphi)))

Cband_hat_diag <- diag(Cband_hat)
UQsize <- qnorm(0.975, m_hat, sqrt(Cband_hat_diag)) - qnorm(0.025, m_hat, sqrt(Cband_hat_diag))

df_uq <- data.frame(pred_points, UBands=UQsize)

(pred_plot <- ggplot() +
  geom_sf(data = michiganplus, fill = "white", color = "black") +
    coord_sf(ylim = c(42, 46), xlim = c(-87, -82)) +
    geom_raster(data = df_uq, aes(x=Longitude, y=Latitude, fill=UBands)) +
    scale_fill_scico(palette="batlowK") + 
    theme_minimal())
```


## Making predictions based on a GP response model

-  As we move away from the observed data
      - The predictive mean gets closer to 0 (which is the prior mean we have assumed for the GP)
      - Uncertainty bands become wider 
-  This matches the intuition that spatial dependence is determined by distance


## More on point predictions and uncertainty bands

-  We considered the Greenness dataset which is a gridded dataset in 2 dimensions
-  Let's simplify and consider a 1-dimensional spatial domain and a small dataset
-  [Explain what the code below is doing]{.orange}
-  [What covariance function are we using]{.orange}

```{r}
set.seed(696)
n <- 10
x <- runif(n)
sigmasq <- 1
phi <- 12
tausq <- 0.0001
C <- sigmasq * exp(-phi * as.matrix(dist(x))^2) + tausq * diag(n)

y <- t(chol(C)) %*% rnorm(n)

df <- data.frame(x=x, y=y)
df %>% ggplot(aes(x,y)) +
  geom_point() +
  theme_minimal()
```

## More on point predictions 

-  Let's make predictions on a fine grid 
-  The GP interpolates!

```{r echo=FALSE}
make_gp_1d <- function(n, xout, sigmasq, phi, tausq, uq=FALSE){
  set.seed(696)
  #n <- 10
  x <- runif(n)
  #sigmasq <- 1
  #phi <- 12
  #tausq <- 0.0001
  C <- sigmasq * exp(-phi * as.matrix(dist(x))^2) + tausq * diag(n)
  
  y <- t(chol(C)) %*% rnorm(n)
  
  #n_out <- nout
  #x_out <- seq(-0.2, 1, length.out=n_out)
  x_out <- xout 
  n_out <- length(x_out)
  
  x_all <- c(x, x_out)
  n_all <- n + n_out
  C_all <- sigmasq * exp(-phi * as.matrix(dist(x_all))^2) + tausq * diag(n_all)
  C_io <- C_all[1:n, (n+1):n_all]
  y_pred <- t(C_io) %*% solve(C, y)
  
  C_oo <- C_all[(n+1):n_all, (n+1):n_all]
  C_pred <- C_oo - t(C_io) %*% solve(C, C_io)
  uupp <- qnorm(0.975, y_pred, sqrt(diag(C_pred))) 
  ulow <- qnorm(0.025, y_pred, sqrt(diag(C_pred)))
  
  df <- data.frame(x=x, y=y)
  df_out <- data.frame(x=x_out, y=y_pred, uq_low=ulow, uq_upp=uupp)
  
  visual <- ggplot() +
    geom_point(data=df, aes(x=x, y=y))
  
  if(uq) {visual <- visual +
    geom_ribbon(data=df_out, aes(x=x, ymin = uq_low, ymax=uq_upp), fill="darkblue", alpha=0.2) 
  }

  visual <- visual +
  geom_line(data=df_out, aes(x=x, y=y_pred), color="orange") +
  theme_minimal()

  
  return(list(data_in = df, data_out=df_out, visual=visual))
}

gp1 <- make_gp_1d(10, seq(-0.1, 1, length.out=200), 1, 12, 0.0001)
gp1$visual
```




## More on point predictions 

-  Far from data? Prediction = prior mean of the GP (zero, usually)

```{r echo=FALSE}
gp2 <- make_gp_1d(10, seq(-3, 3, length.out=200), 1, 12, 0.0001)
gp2$visual
```


## More on point predictions 

```{r echo=FALSE}
tausq <- 0.05
```

-  Notice what happens when we increase `tausq` (now at `r tausq`)

```{r echo=FALSE}
gp3 <- make_gp_1d(10, seq(-0.1, 1, length.out=200), 1, 12, tausq)
gp3$visual
```


## More on uncertainty bands

-  Uncertainty increases as we move away from the observations

```{r echo=FALSE}
gp1u <- make_gp_1d(10, seq(-0.1, 1, length.out=200), 1, 12, 0.0001, TRUE)
gp1u$visual
```


## More on uncertainty bands

-  Far from data? Uncertainty band = prior variance

```{r echo=FALSE}
gp2u <- make_gp_1d(10, seq(-3, 3, length.out=200), 1, 12, 0.0001, TRUE)
gp2u$visual
```


## More on point predictions and uncertainty bands

```{r echo=FALSE}
tausq <- 0.05
```

-  Notice what happens when we increase `tausq` (now at `r tausq`)

```{r echo=FALSE}
gp3u <- make_gp_1d(10, seq(-0.1, 1, length.out=200), 1, 12, tausq, TRUE)
gp3u$visual
```


## Summary

-  The same things we see in 1 dimension occur in 2 or more dimensions!
-  As we move away from the observed data
      - The predictive mean gets closer to the prior mean we have assumed for the GP
      - Uncertainty bands become wider, up to being equal to the prior variance
-  This matches the intuition that spatial dependence is determined by distance
-  When we are far from data, we do not have information to make predictions
-  We thus revert to our prior information (i.e. our GP assumption)
-  Because the GP interpolates at the observed locations, it is a flexible model for any dataset
-  The *Gaussian* assumption is **difficult to criticize when we only have 1 replicate of the data** because for any dataset with 1 replicate, the GP will "connect all the points" via interpolation
-  In the presence of noise, other assumptions can be criticized, such as the choice of the covariance function and its parameters


