---
title: "Spatial GLM + Point Reference Spatial Data"
subtitle: "Lecture 21"
author: "Dr. Colin Rundel"
footer: "Sta 344 - Fall 2022"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
execute:
  echo: true
  warning: false
  collapse: true
---

```{r setup, include=FALSE}
library(sta344)
library(sf)
library(raster)

library(tidyverse)
library(patchwork)

knitr::opts_chunk$set(
  fig.align = "center"
)

options(width=50)

set.seed(20221109)

ggplot2::theme_set(ggplot2::theme_bw())
```


# Spatial GLM Models

## Scottish Lip Cancer Data

```{r echo=FALSE}
lip_cancer = readRDS("data/lip_cancer.rds")

lip_cancer |>
  pivot_longer(Observed:Expected, names_to = "var", values_to = "value") |>
  ggplot() +
    geom_sf(aes(fill=value), color=NA) +
    facet_wrap(~forcats::as_factor(var))
```

##

```{r echo=FALSE}
( ggplot(lip_cancer) +
    geom_sf(aes(fill=Observed/Expected), color=NA) + 
    labs(title="Obs/Exp",fill="")
) +
( ggplot(lip_cancer) +
    geom_sf(aes(fill=pcaff), color=NA) +
    labs(title="% Agg Fish Forest",fill="")
)
```

## Neighborhood / weight matrix

::: {.small}
```{r}
A = (st_distance(lip_cancer) |> unclass()) < 1e-6
listw = spdep::mat2listw(A)
```
:::

```{r echo=FALSE}
plot(st_geometry(lip_cancer), asp=1) 
plot(listw, coords = st_coordinates(st_centroid(lip_cancer)), add=TRUE, col="blue", pch=16)
```


## Moran's I

:::: {.columns .small}
::: {.column width='50%'}
```{r}
spdep::moran.test(lip_cancer$Observed, listw)
```
:::

::: {.column width='50%'}
```{r}
spdep::moran.test(lip_cancer$Observed / lip_cancer$Expected, listw)
```
:::
::::


## GLM

::: {.small}
```{r}
l = glm(Observed ~ offset(log(Expected)) + pcaff, 
        family="poisson", data=lip_cancer)
summary(l)
```
:::


## GLM Fit

```{r echo=FALSE} 
lip_cancer = lip_cancer %>%
  mutate(
    glm_pred = predict(l, type="response"),
    glm_resid = Observed - glm_pred
  )

( ggplot(lip_cancer) +
    geom_sf(aes(fill=Observed), color=NA) + 
    labs(title="Observed Cases",fill="")) +
( ggplot(lip_cancer) +
    geom_sf(aes(fill=glm_pred), color=NA) + 
    labs(title="GLM Predicted Cases",fill=""))
```

## GLM Fit 

```{r echo=FALSE}
ggplot(lip_cancer) +
  geom_abline(intercept=0, slope=1, color="grey") +
  geom_point(aes(x=Observed, y=glm_pred)) 
```

## GLM Residuals


```{r echo=FALSE} 
( ggplot(lip_cancer) +
    geom_sf(aes(fill=glm_pred), color=NA) + 
    labs(title="GLM Predicted Cases",fill="")) +
( ggplot(lip_cancer) +
    geom_sf(aes(fill=glm_resid), color=NA) +
    labs(title="GLM Residuals",fill=""))
```

## Model Results

::: {.small}
```{r}
#RMSE
yardstick::rmse_vec(lip_cancer$Observed, lip_cancer$glm_pred)
```
:::

. . .

::: {.small}
```{r}
#Moran's I
spdep::moran.test(lip_cancer$glm_resid, listw)
```
:::

## A hierachical model for lip cancer

We have observed counts of lip cancer for 56 districts in Scotland. Let $y_i$ represent the number of lip cancer for district $i$.

$$\begin{aligned}
y_i &\sim \text{Poisson}(\lambda_i) \\
\\
\log(\lambda_i) &= \log(E_i) + x_i \beta + \omega_i \\
\\
\boldsymbol{\omega} &\sim \mathcal{N}(\boldsymbol{0},~\sigma^2(\boldsymbol{D}-\phi\,\boldsymbol{A})^{-1})
\end{aligned}$$

where $E_i$ is the expected counts for each region (and serves as an offset).


## Data prep & CAR model

::: {.small}
```{r}
X = model.matrix(~scale(lip_cancer$pcaff))
offset = lip_cancer$Expected
y = lip_cancer$Observed
```

```{r}
car_model = "
data {
  int<lower=0> N;
  int<lower=0> p;
  int<lower=0> y[N];
  matrix[N,N] A;
  matrix[N,p] X;
  vector[N] offset;
}
transformed data {
  vector[N] nb = A * rep_vector(1, N);
  matrix[N,N] D = diag_matrix(nb);
}
parameters {
  vector[N] w_s;
  vector[p] beta;
  real<lower=0> sigma2;
  real<lower=0,upper=1> phi;
}
transformed parameters {
  vector[N] eta = log(offset) + X * beta + w_s; 
}
model {
  matrix[N,N] Sigma_inv = (D - phi * A) / sigma2;
  w_s ~ multi_normal_prec(rep_vector(0,N), Sigma_inv);

  beta ~ normal(0,10);
  sigma2 ~ cauchy(0,5);
  
  y ~ poisson_log(eta);
}
"
```
:::

## CAR Fitting

```{r include=FALSE, cache=TRUE}
car = rstan::stan_model(model_code = car_model)

car_m = rstan::sampling(
  car, iter=5000, cores=4,
  data = list(N=nrow(X), A=A, X=X, p=ncol(X), offset=offset, y=y)
)
```

```{r eval=FALSE}
car = rstan::stan_model(model_code = car_model)

car_m = rstan::sampling(
  car, iter=5000, cores=4,
  data = list(N=nrow(X), A=A, X=X, p=ncol(X), offset=offset, y=y)
)
```

##

```{r echo=FALSE}
car_m %>%
  tidybayes::gather_draws(sigma2, phi, beta[i]) %>%
  mutate(
    .variable = ifelse(is.na(i), .variable, paste0(.variable,"_",i)),
    .chain = as.factor(.chain)
  ) %>%
  ggplot(aes(x=.iteration, y=.value, color=.chain)) +
    geom_line() +
    facet_wrap(~.variable, scales = "free_y")
```

##

```{r echo=FALSE}
car_m %>%
  tidybayes::gather_draws(sigma2, phi, beta[i]) %>%
  mutate(
    .variable = ifelse(is.na(i), .variable, paste0(.variable,"_",i)),
    .chain = as.factor(.chain)
  ) %>%
  ggplot(aes(x=.value, fill=.chain)) +
    geom_density(alpha=0.2) +
    facet_wrap(~.variable, scales = "free")
```

## CAR Predictions ($\hat\lambda$)

```{r echo=FALSE}
car_pred = tidybayes::gather_draws(car_m, eta[i], w_s[i]) %>%
  group_by(.variable, i) %>%
  summarize(
    mean = mean(.value),
    med = median(.value),
    q025 = quantile(.value, 0.025),
    q975 = quantile(.value, 0.975),
    .groups = "drop"
  )


car_lip_cancer = lip_cancer %>% 
  mutate(i = seq_len(n())) %>%
  left_join(
    car_pred %>% filter(.variable == "eta"),
    by = "i"
  ) %>%
  mutate(y_pred = exp(mean))

( ggplot(car_lip_cancer) +
    geom_sf(aes(fill=Observed), color=NA) + 
    labs(title="Observed Cases",fill="")) +
( ggplot(car_lip_cancer) +
    geom_sf(aes(fill=y_pred), color=NA) + 
    labs(title="Predicted Cases",fill=""))
```

## CAR Predictions

```{r echo=FALSE}
ggplot(car_lip_cancer, aes(x=Observed, y=y_pred)) +
  geom_abline(intercept=0, slope=1, color="grey") +
  geom_point() +
  geom_errorbar(aes(ymin=exp(q025), ymax=exp(q975)), alpha=0.25) +
  labs(y="ʎ")
```

## CAR Residuals

```{r echo=FALSE}
( ggplot(car_lip_cancer) +
    geom_sf(aes(fill=y_pred), color=NA) + 
    labs(title="Predicted Cases",fill="")) +
( ggplot(car_lip_cancer) +
    geom_sf(aes(fill=Observed-y_pred), color=NA) +
    labs(title="Residuals",fill=""))
```


## CAR Results

::: {.small}
```{r}
#RMSE
yardstick::rmse_vec(car_lip_cancer$Observed, car_lip_cancer$y_pred)
```
:::

. . .

::: {.small}
```{r}
#Moran's I
spdep::moran.test(car_lip_cancer$Observed - car_lip_cancer$y_pred, listw)
```
:::


## Latent spatial process

```{r echo=FALSE}
lip_cancer %>% 
  mutate(i = seq_len(n())) %>%
  left_join(
    car_pred %>% filter(.variable == "w_s"),
    by = "i"
  ) %>% ggplot() +
    geom_sf(aes(fill=mean), color=NA) +
    labs(title="w(s)")
```


## Intrinsic Autoregressive Model (IAR)

```{r iar, include=FALSE, cache=TRUE}
rownames(A) = lip_cancer$District

iar_m = brms::brm(
  Observed~offset(Expected)+pcaff+car(A, gr=District, type="icar"), 
  data=lip_cancer, data2=list(A=A), 
  family = poisson, cores=4, iter=5000
)
```

```{r eval=FALSE, ref.label="iar"}
```

##

```{r echo=FALSE}
plot(iar_m)
```



## Predictions

```{r echo=FALSE}
iar_pred = sta344::epred_draws_fix(iar_m, newdata=lip_cancer) %>%
  filter(.chain == 1) %>%
  as_tibble() %>%
  group_by(District) %>%
  summarize(
    Observed = Observed[1],
    y_pred = mean(.epred),
    q025 = quantile(.epred, 0.025),
    q975 = quantile(.epred, 0.975),
    resid =  Observed - y_pred,
    geometry = geometry[1],
    .groups = "drop"
  ) %>%
  st_as_sf()

( ggplot(iar_pred) +
    geom_sf(aes(fill=Observed), color=NA) + 
    labs(title="Observed Cases",fill="")) +
( ggplot(iar_pred) +
    geom_sf(aes(fill=y_pred), color=NA) + 
    labs(title="Predicted Cases",fill=""))
```

## Observed vs predicted

```{r echo=FALSE}
ggplot(iar_pred, aes(x=Observed, y=y_pred)) +
  geom_abline(intercept=0, slope=1, color="grey") +
  geom_point() +
  geom_errorbar(aes(ymin=q025, ymax=q975), alpha=0.25) +
  labs(y="ʎ")
```

## Residuals

```{r echo=FALSE}
( ggplot(iar_pred) +
    geom_sf(aes(fill=y_pred), color=NA) + 
    labs(title="Predicted Cases",fill="")) +
( ggplot(iar_pred) +
    geom_sf(aes(fill=resid), color=NA) +
    labs(title="Residuals",fill=""))
```

## IAR Results

::: {.small}
```{r}
#RMSE
yardstick::rmse_vec(iar_pred$Observed, iar_pred$y_pred)

#Moran's I
spdep::moran.test(iar_pred$resid, listw)
```
:::


# Point Referenced Data

## Example - PM2.5 from CSN

The Chemical Speciation Network are a series of air quality monitors run by EPA (221 locations in 2007). We'll look at a subset of the data from Nov 11th, 2007 (n=191) for just PM2.5. 

```{r include=FALSE}
states = st_read("data/us/") %>%
  filter(!STATE %in% c("Hawaii", "Alaska", "Puerto Rico", "U.S. Virgin Islands")) %>%
  select(state = STATE) %>%
  filter(!duplicated(state))

usa = states |>
  summarize()
```

```{r echo=FALSE, fig.height=3.5}
csn = readRDS("data/epa/csn.rds")

csn = csn %>% 
  select(site:date, pm25) %>% 
  na.omit() %>% 
  filter(date == lubridate::ymd("2007-11-14")) %>% 
  filter(longitude > -140) %>% # Remove Hawaii Location
  as_tibble()

csn_plot = ggplot() +
  geom_sf(data=states, color="grey", fill=NA) +
  geom_point(data=csn, aes(x=longitude, y=latitude, color=pm25), alpha=1, size=3) +
  scale_color_viridis_c() 

csn_plot
```




##

```{r}
csn
```

## Aside - Splines

![](imgs/Lec21/spline1.png){fig-align="center" width="80%"}

##

![](imgs/Lec21/spline2.png){fig-align="center" width="100%"}


## Splines in 1d - Smoothing Splines

These are a mathematical analogue to the drafting splines represented using a penalized regression model.

. . .

We want to find a function $f(x)$ that best fits our observed data $\boldsymbol{y} = y_1, \ldots, y_n$ while being *smooth*.

$$ \underset{f(x)}{\arg\min} ~ \sum_{i=1}^n\left(y_i - f(x_i)\right)^2 + \lambda \int_{-\infty}^\infty f''(x)^2 ~ dx $$

. . .

Interestingly, this minimization problem has an exact solution which is given by a mixture of weighted natural cubic splines (cubic splines that are linear in the tails) with knots at the observed data locations ($x$s).


## Splines in 2d - Thin Plate Splines

Now imagine we have observed data of the form $(x_i, y_i, z_i)$ where we wish to predict $z_i$ given $x_i$ and $y_i$ for all $i$. We can extend the smoothing spline model in two dimensions,

$$ \underset{f(x,y)}{\arg\min} ~~ \sum_{i=1}^n (z_i-f(x_i,y_i))^2 + \lambda  \int_{-\infty}^\infty  \int_{-\infty}^\infty \left(\frac{\partial^2 f}{\partial x^2} + 2 \frac{\partial^2 f}{\partial x \, \partial y} + \frac{\partial^2 f}{\partial y^2} \right) dx\, dy$$

. . .

The solution to this equation has a natural representation using a weighted sum of *radial basis functions* with knots at the observed data locations ($\boldsymbol{x_i}$)
 
$$ f(\boldsymbol{x}) = \sum_{i=1}^n w_i ~ d(\boldsymbol{x}, \boldsymbol{x_i})^2 \log d(\boldsymbol{x}, \boldsymbol{x_i}).$$

## Prediction locations

::: {.small}
```{r}
r_usa = stars::st_rasterize(
  usa, 
  stars::st_as_stars(st_bbox(usa), 
    nx = 100, ny = 50, values=NA_real_)
)
plot(r_usa)
```
:::

## Fitting a TPS

::: {.small}
```{r}
coords = select(csn, long=longitude, lat=latitude) |>
  as.matrix()
(tps = fields::Tps(x = coords, Y=csn$pm25, lon.lat=TRUE))
```
:::


## Predictions

::: {.small}
```{r echo=FALSE}
pred = r_usa |>
  as.data.frame() |>
  filter(!is.na(ID)) %>%
  select(longitude = x, latitude = y) %>%
  as_tibble()

tps_pred = pred %>%
  mutate(
    .,
    pred = predict(tps, as.matrix(.[,1:2]))
  ) %>%
  stars::st_as_stars()

plot(tps_pred)
```
:::



# Gaussin Process Models / Kriging

## Variogram

::: {.small}
```{r}
coords = csn %>% select(latitude, longitude) %>% as.matrix()
d = fields::rdist(coords)

geoR::variog(
  coords = coords, data = csn$pm25, messages = FALSE, 
  uvec = seq(0, max(d)/2, length.out=50)
) %>% 
  plot()
```
:::


##

::: {.small}
```{r}
geoR::variog(
  coords = coords, data = csn$pm25, messages = FALSE,
  uvec = seq(0, max(d)/4, length.out=50)
) %>% plot()
```
:::


## Isotropy / Anisotropy

::: {.small}
```{r message=FALSE}
v4 = geoR::variog4(
  coords = coords, data = csn$pm25, messages = FALSE,
  uvec = seq(0, max(d)/4, length.out = 50)
)
plot(v4)
```
:::

## GP Spatial Model

If we assume that our data is *stationary* and *isotropic* then we can use a Gaussian Process model to fit the data. We will assume an exponential covariance structure. 

$$ \boldsymbol{y} \sim \mathcal{N}(\boldsymbol{\mu},~\Sigma) $$
$$ \{\Sigma\}_{ij} = \sigma^2 \exp(- l \, \lVert s_i - s_j\lVert) + \sigma^2_n \, 1_{i=j} $$

. . .

we can also view this as a spatial random effects model where

$$ 
y(\boldsymbol{s}) = \mu(\boldsymbol{s}) + w(\boldsymbol{s}) + \epsilon(\boldsymbol{s}) \\
w(\boldsymbol{s}) \sim \mathcal{N}(0,\Sigma') \\
\epsilon(s_i) \sim \mathcal{N}(0,\sigma^2_n) \\
\{\Sigma'\}_{ij} = \sigma^2 \exp(- r \, \lVert s_i - s_j\lVert) 
$$

## Fitting with `gplm()` (`spBayes`)

```{r eval=FALSE}
max_range = max(dist(csn[,c("longitude", "latitude")])) / 4

m = gplm(
  pm25~1, data = csn, coords=c("longitude", "latitude"),
  cov_model = "exponential",
  starting = list(phi = 3/3, sigma.sq = 33, tau.sq = 17),
  tuning = list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1),
  priors = list(
    phi.Unif = c(3/max_range, 3/(0.5)), 
    sigma.sq.IG = c(2, 2), 
    tau.sq.IG = c(2, 2)
  ),
  thin=10,
  verbose=TRUE
)
```

```{r include=FALSE}
m = readRDS("Lec21_csn_gp.rds")
```



##

```{r}
m
```


## Parameter values

```{r}
plot(m)
```


## Predictions

```{r eval=FALSE}
(p = predict(m, newdata=pred, coords=c("longitude", "latitude")))
```

```{r echo=FALSE}
(p = readRDS("Lec21_csn_gp_pred.rds"))
```


##

```{r echo=FALSE}
gp_pred = tidybayes::gather_draws(p, y[i]) %>%
  filter(.chain == 1) %>%
  group_by(.chain, i) %>%
  summarize(
    mean = mean(.value),
    med = median(.value),
    sd = sd(.value),
    .groups = "drop"
  ) %>%
  left_join(
    pred %>% mutate(i = seq_len(n())),
    ., 
    by = "i"
  )

gp_pred %>% 
  stars::st_as_stars() %>% 
  select(mean) %>% 
  plot()
```

##

```{r echo=FALSE}
gp_pred %>% 
  stars::st_as_stars() %>% 
  select(sd) %>% 
  plot()
```