---
title: "GPs for GLMs<br/>+ Spatial Data"
subtitle: "Lecture 16"
author: "Dr. Colin Rundel"
footer: "Sta 344/644 - Fall 2023"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute:
  echo: true
  warning: false
  collapse: true
---

```{r setup}
#| include: false
library(tidyverse)
library(patchwork)
library(dukestm)

library(sf)
library(stringr)

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

options(width=60)

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

set.seed(20221028)
```


# GPs and GLMs

## Logistic Regression
 
A typical logistic regression problem uses the following model,

$$
\begin{aligned}
y_i &\sim \text{Bern}(p_i)\\
\text{logit}(p_i) 
  &= \boldsymbol{X}\,\boldsymbol{\beta} \\
  &= \beta_0 + \beta_1 \, x_{i1} + \cdots + \beta_k \, x_{ik}
\end{aligned}
$$

. . .

there is no reason that the linear equation above can't contain thing like random effects or GPs

$$
\begin{aligned}
y_i &\sim \text{Bern}(p_i) \\
\text{logit}(p_i) 
  &= \eta_i = \boldsymbol{X}\,\boldsymbol{\beta} + w(\boldsymbol{x}) \\
w(\boldsymbol{x}) &\sim N(0,\Sigma)
\end{aligned}
$$

## A toy example

```{r echo=FALSE}
n=200

eta = function(x) 2.5*sin(2.1*pi*(x-0.25)) 

d = tibble(x=runif(n)) %>%
  mutate(
    eta = eta(x),
    p = 1/(1+exp(-eta)),
    y = rbinom(n(), size=1, prob = p)
  )

(ggplot(d, aes(x=x, y=eta)) + geom_line() +
 ggplot(d, aes(x=x, y=p)) + geom_line()) / 
(ggplot(d, aes(x=x, y=y)) + geom_jitter(height = 0.05))
```

## A standard GLM

```{r}
(g = glm(y~x, family="binomial", data=d))
```

##

```{r echo=FALSE}
newdata = data.frame(
  x=seq(0,1,length.out=101)
)

p1 = g |> broom::augment(newdata=newdata, type.predict = "link")
p2 = g |> broom::augment(newdata=newdata, type.predict = "response")


(
  ggplot(d, aes(x=x, y=eta)) + 
    geom_line() + 
    geom_line(data=p1, aes(y=.fitted), color='red') +
 ggplot(d, aes(x=x, y=p)) + 
    geom_line() + 
    geom_line(data=p2, aes(y=.fitted), color='red')
) / (
  ggplot(d, aes(x=x, y=y)) + 
    geom_jitter(height = 0.05) +
    geom_line(data=p2, aes(y=.fitted), color='red')
)
```


## A quadratic GLM

```{r}
(g2 = glm(y~poly(x,2), family="binomial", data=d))
```

##

```{r echo=FALSE}
newdata = data.frame(
  x=seq(0,1,length.out=101)
)

p1 = g2 |> broom::augment(newdata=newdata, type.predict = "link")
p2 = g2 |> broom::augment(newdata=newdata, type.predict = "response")


(
  ggplot(d, aes(x=x, y=eta)) + 
    geom_line() + 
    geom_line(data=p1, aes(y=.fitted), color='red') +
 ggplot(d, aes(x=x, y=p)) + 
    geom_line() + 
    geom_line(data=p2, aes(y=.fitted), color='red')
) / (
  ggplot(d, aes(x=x, y=y)) + 
    geom_jitter(height = 0.05) +
    geom_line(data=p2, aes(y=.fitted), color='red')
)
```


## Model fitting

::: {.small}
```{r eval=FALSE}
m = gpglm(
  chains = 2,
  y~1, family="binomial",
  data = d, coords = c("x"),
  cov_model = "exponential",
  n_batch = 200,
  batch_len = 100,
  starting=list(
    "beta"=0, "phi"=3/0.5, "sigma.sq"=5, "w"=0
  ),
  priors=list(
    "beta.Normal"=list(0,0.01),
    "phi.unif"=c(3/0.75, 3/0.01),
    "sigma.sq.ig"=c(2, 1)
  ),
  tuning=list(
    "beta"=0.5, "phi"=0.5, "sigma.sq"=0.5, "w"=0.5
  ),
  verbose=TRUE
)
```
:::

```{r echo=FALSE}
file = "Lec16_gpglm_logistic_model.rds"
  
if (!file.exists(file)) {
  m = gpglm(
    y~1, family="binomial",
    data = d, coords = c("x"),
    cov_model = "exponential",
    n_batch = 200,
    batch_len = 100,
    starting=list(
      "beta"=0, "phi"=3/0.1, "sigma.sq"=3, "w"=0
    ),
    priors=list(
      "beta.Normal"=list(0,1),
      "phi.unif"=c(3/0.75, 3/0.01),
      "sigma.sq.ig"=c(2, 1)
    ),
    tuning=list(
      "beta"=0.5, "phi"=0.5, "sigma.sq"=0.5, "w"=0.5
    ),
    verbose=TRUE
  )
  
  newdata = data.frame(
    x=seq(0,1,length.out=101)
  )
  
  p = predict(m, newdata=newdata, coords="x", thin=10)
  
  saveRDS(
    list(
      m = m,
      newdata = newdata,
      p = p
    ),
    file=file
  )
  
  
} else {
  res = readRDS(file)
  m = res$m
  p = res$p
}
```

## Model diagnostics

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


## Model predictions

::: {.medium}
```{r eval=FALSE}
newdata = data.frame(
  x=seq(0,1,length.out=101)
)

(p = predict(m, newdata=newdata, coords="x", thin=10))
```

```{r echo=FALSE}
p
```
:::

## Predicted y

::: {.small}
```{r}
p |>
  tidybayes::gather_draws(y[i]) |>
  ggplot2::ggplot(ggplot2::aes(x=i/max(i),y=.value)) +
    tidybayes::stat_lineribbon(alpha=0.25) +
    geom_line(data=d |> arrange(x), aes(x=x, y=p), color='red')
```
:::


## Predicted w

::: {.small}
```{r}
p |>
  tidybayes::gather_draws(w[i]) |>
  ggplot2::ggplot(ggplot2::aes(x=i/max(i),y=.value)) +
    tidybayes::stat_lineribbon(alpha=0.25) +
    geom_line(data=d |> arrange(x), aes(x=x, y=eta), color='red')
```
:::



# Count data

## Polio cases

`Polio` from the `glarma` package.

> This data set gives the monthly number of cases of poliomyelitis in the U.S. for the years 1970–1983 as reported by the Center for Disease Control.


```{r echo=FALSE, fig.height=3.5}
data(Polio, package="glarma")

polio = tibble(
  year = seq(1970, 1984-1/12, by=1/12),
  cases = Polio$Cases
)

ggplot(polio, aes(x=year, y=cases)) +
  geom_line() +
  geom_point()
```
 
## Polio Model

#### Model:
$$
\begin{aligned}
y_i &\sim \text{Pois}(\lambda_i) \\
\text{log}(\lambda_i) 
  &= \beta_0 + w(\boldsymbol{t}) \\
\\
w(\boldsymbol{t}) &\sim N(0,\Sigma) \\
\{\boldsymbol{\Sigma}\}_{ij} &= \sigma^2 \exp(-| l \, d_{ij}|)
\end{aligned}
$$


#### Priors:

$$
\begin{aligned}
\beta_0 &\sim N(0,1) \\
\phi &\sim \text{Unif}\left(\frac{3}{6}, \frac{3}{1/12}\right)\\
\sigma^2 &\sim \text{Inv-Gamma}(2,1)
\end{aligned}
$$


## Model fitting

::: {.small}
```{r eval=FALSE}
m = gpglm(
  cases~1, family="poisson",
  data = polio, coords = c("year"),
  cov_model = "exponential",
  starting=list(
    "beta"=0, "phi"=3/2, "sigma.sq"=1, "w"=0
  ),
  tuning=list(
    "beta"=0.5, "phi"=0.5, "sigma.sq"=0.5, "w"=0.5
  ),
  priors=list(
    "beta.Normal"=list(0,1),
    "phi.unif"=c(3/6, 3/(1/12)),
    "sigma.sq.ig"=c(2, 1)
  ),
  n_batch = 100,
  batch_len = 100,
  verbose = FALSE
)
```
:::

```{r echo=FALSE}
file = "Lec16_gpglm_poisson_model.rds"

if (!file.exists(file)) {
  m = gpglm(
  cases~1, family="poisson",
  data = polio, coords = c("year"),
    cov_model = "exponential",
    starting=list(
      "beta"=0, "phi"=3/2, "sigma.sq"=1, "w"=0
    ),
    tuning=list(
      "beta"=0.5, "phi"=0.5, "sigma.sq"=0.5, "w"=0.5
    ),
    priors=list(
      "beta.Normal"=list(0,1),
      "phi.unif"=c(3/6, 3/(1/12)),
      "sigma.sq.ig"=c(2, 1)
    ),
    n_batch = 100,
    batch_len = 100,
    chains = 2,
    verbose = TRUE
  )
  
  newdata = data.frame(
    year = seq(1970, 1984, by=0.1) |> jitter()
  )

  p = predict(m, newdata=newdata, coords="year", thin=5)
  
  saveRDS(
    list(
      m=m, p=p, newdata=newdata
    ), 
    file=file
  )
} else {
  res = readRDS(file)
  m = res$m
  p = res$p
}
```



## Model diagnostics

```{r}
plot(m, thin=5)
```


## Model fit

::: {.medium}
```{r eval=FALSE}
newdata = data.frame(
  year = seq(1970, 1984, by=0.1) |> jitter()
)
(p = predict(m, newdata=newdata, coords="year", thin=5))
```

```{r echo=FALSE}
p
```
:::

## Predicted y

::: {.small}
```{r}
p |>
  tidybayes::gather_draws(y[i]) |>
  ggplot2::ggplot(ggplot2::aes(x=i/max(i),y=.value)) +
    tidybayes::stat_lineribbon(alpha=0.25) +
    geom_line(data=polio, aes(x=(year-min(year))/(max(year)-min(year)), y=cases), color='red', alpha=0.5)
```
:::


## Predicted w

::: {.small}
```{r}
p |>
  tidybayes::gather_draws(w[i]) |>
  ggplot2::ggplot(ggplot2::aes(x=i/max(i),y=.value)) +
    tidybayes::stat_lineribbon(alpha=0.25)
```
:::



# Spatial data in R


## Packages for geospatial data in R

::: {.small}
R has a rich package ecosystem for read/writing, manipulating, and analyzing geospatial data. Some core packages:

* `sp` - core classes for handling spatial data, additional utility functions - **Deprecated**

* `rgdal` - R interface to `gdal` (Geospatial Data Abstraction Library) for reading and writing spatial data - **Deprecated**

* `rgeos` - R interface to `geos` (Geometry Engine Open Source) library for querying and manipulating spatial data. Reading and writing WKT. - **Deprecated**

* `raster` - classes and tools for handling spatial raster data.

* `sf` - Combines the functionality of `sp`, `rgdal`, and `rgeos` into a single package based on tidy simple features.

* `stars` - Reading, manipulating, writing and plotting spatiotemporal arrays (rasters)

* `terra` - Methods for spatial data analysis with vector (points, lines, polygons) and raster (grid) data. Replaces `raster`. 

See more - [Spatial task view](http://cran.r-project.org/web/views/Spatial.html)
:::

## Installing `sf`

This is the hardest part of using the `sf` package, difficulty comes from its dependence on several external libraries (`geos`, `gdal`, `proj`, and `udunits2`).

* *Windows* - installing from source works when Rtools is installed (system requirements are downloaded from rwinlib)

* *MacOS* - install dependencies via homebrew: `gdal2`, `geos`, `proj`, `udunits2`.

* *Linux* - Install development packages for GDAL (>= 2.0.0), GEOS (>= 3.3.0), Proj.4 (>= 4.8.0), and udunits2 from your package manager of choice.


More specific details are included in the [README on github](https://github.com/r-spatial/sf).



## Simple Features

```{r, echo=FALSE}
par(mar=c(1,1,2,1), mfrow=c(2,4))

## Single Geometries

pt = st_point(c(30, 10))
ls = st_linestring(matrix(c(30, 10, 10, 30, 40, 40), byrow=TRUE, ncol=2))
poly = st_polygon(list(matrix(c(30, 10, 40, 40, 20, 40, 10, 20, 30, 10), ncol=2, byrow=TRUE)))
poly_hole = st_polygon(
              list(
                matrix(c(35, 10, 45, 45, 15, 40, 10, 20, 35, 10), ncol=2, byrow=TRUE),
                matrix(c(20, 30, 35, 35, 30, 20, 20, 30), ncol=2, byrow=TRUE)
              )
            )

plot(pt, axes=FALSE, main="Point", pch=16); box()
plot(ls, axes=FALSE, main="Linestring");    box()
plot(poly, axes=FALSE, col="lightgrey", main="Polygon");  box()
plot(poly_hole, axes=FALSE, col="lightgrey", main="Polygon w/ Hole(s)");  box()


## Multi Geometries

pts = st_multipoint(matrix(c(10, 40, 40, 30, 20, 20, 30, 10), ncol=2, byrow=TRUE))
lss = st_multilinestring(list( 
        matrix(c(10, 10, 20, 20, 10, 40), ncol=2, byrow=TRUE),
        matrix(c(40, 40, 30, 30, 40, 20, 30, 10), ncol=2, byrow=TRUE)
      ))
        
polys = st_multipolygon(list(
          list(matrix(c(30, 20, 45, 40, 10, 40, 30, 20), ncol=2, byrow=TRUE)),
          list(matrix(c(15, 5, 40, 10, 10, 20, 5, 10, 15, 5), ncol=2, byrow=TRUE))
        ))
polys_hole = st_multipolygon(list(
                list(matrix(c(40, 40, 20, 45, 45, 30, 40, 40), ncol=2, byrow=TRUE)),
                list(
                  matrix(c(20, 35, 10, 30, 10, 10, 30, 5, 45, 20, 20, 35), ncol=2, byrow=TRUE),
                  matrix(c(30, 20, 20, 15, 20, 25, 30, 20), ncol=2, byrow=TRUE)
                )
              ))

plot(pts, axes=FALSE, main="Multipoint", pch=16); box()
plot(lss, axes=FALSE, main="Multilinestring");    box()
plot(polys, axes=FALSE, col="lightgrey", main="Multipolygon");  box()
plot(polys_hole, axes=FALSE, col="lightgrey", main="Multipolygon w/ Hole(s)");  box()
```


## Reading, writing, and converting

- `sf`
    * `st_read()` / `st_write()` - Shapefile, GeoJSON, KML, ...
    * `read_sf()` / `write_sf()` - Same, supports tibbles ...
    * `st_as_sfc()` / `st_as_wkt()` - sf <-> WKT
    * `st_as_sfc()` / `st_as_binary()` - sf <-> WKB
    * `st_as_sfc()` / `as(x, "Spatial")` - sf <-> sp

::: {.aside}
See [sf vignette #2 - Reading, Writing and Converting Simple Features](https://cran.r-project.org/web/packages/sf/vignettes/sf2.html).
:::

# Geospatial data in<br/>the real world

## Projections

```{r projs, echo=FALSE, message=FALSE}
lat_lines  = map(seq(9.999, 89.999, length.out = 9), ~ cbind(seq(-179.999, -9.999, length.out=100), .))
long_lines = map(seq(-179.999, -9.999, length.out = 17), ~ cbind(., seq(9.999, 89.999, length.out=100))) 

lat_long = c(lat_lines, long_lines) %>% 
  st_multilinestring() %>%
  st_sfc() %>%
  st_set_crs("+proj=longlat +datum=WGS84 +no_defs")


world = st_as_sf(s2::s2_data_tbl_countries) |>
  st_set_crs("+proj=longlat +datum=WGS84 +no_defs")

NAm = world %>% filter(continent %in% "North America")
NAm_google = st_transform(NAm, "+init=epsg:3857")


par(mar=c(3,2,2,1),mfrow=c(2,3))
plot(lat_long, col=adjustcolor("grey",alpha.f = 0.5), axes=TRUE, main="Lat/Long (epsg:4326)")
plot(st_geometry(NAm), col="black", add=TRUE)

plot(st_transform(lat_long, "+init=epsg:3857"), col=adjustcolor("grey",alpha.f = 0.5), axes=TRUE, main="Google / Web Mercator (epsg:3857)", ylim=c(0, 2e7))
plot(st_transform(NAm, "+init=epsg:3857") %>% st_geometry(), col="black", add=TRUE)

lcc = "+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
plot(st_transform(lat_long, lcc), col=adjustcolor("grey",alpha.f = 0.5), axes=TRUE, main="Lambert Conformal Conic:")
plot(st_transform(NAm, lcc) %>% st_geometry(), col="black", add=TRUE)

aea = "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
plot(st_transform(lat_long, aea), col=adjustcolor("grey",alpha.f = 0.5), axes=TRUE, main="Alberts Equal Area")
plot(st_transform(NAm, aea) %>% st_geometry(), col="black", add=TRUE)

robinson = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
plot(st_transform(lat_long, robinson), col=adjustcolor("grey",alpha.f = 0.5), axes=TRUE, main="Robinson")
plot(st_transform(NAm, robinson) %>% st_geometry(), col="black", add=TRUE)

mollweide = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
plot(st_transform(lat_long, mollweide), col=adjustcolor("grey",alpha.f = 0.5), axes=TRUE, main="Mollweide")
plot(st_transform(NAm, mollweide) %>% st_geometry(), col="black", add=TRUE)
```

## Distance on a Sphere

```{r echo=FALSE}
states = st_read("data/us/states-unfiltered.shp", quiet = TRUE, stringsAsFactors = FALSE) %>%
  filter(!STATE %in% c("Alaska", "Hawaii", "Puerto Rico", "U.S. Virgin Islands"))

plot(st_geometry(states))

durham = c(-78.8986, 35.9940)
la = c(-118.243, 34.0522)
locs = rbind(durham, la)

gc = geosphere::gcIntermediate(durham, la, n=50, addStartEnd=TRUE)

points(locs, col='blue', pch=16, cex=2)
lines(locs, col='red', lty=2, lwd=2)
lines(gc, col='red', lty=3, lwd=2)
```

## Dateline

How long is the flight between the Western most and the Eastern most points in the US?

. . .

```{r echo=FALSE, message=FALSE, fig.height=5, fig.width=8}
par(mar=c(3,3,1,1))

ak = st_read("data/ak/states.shp", quiet = TRUE, stringsAsFactors = FALSE)
ak_geom = st_geometry(ak)

west_hem = st_polygon(list(matrix(c(-180,90, -180,-90, 0,-90, 0,90, -180,90), ncol=2,byrow=TRUE))) %>% 
  st_sfc() %>%
  st_set_crs("+proj=longlat +datum=WGS84")

east_hem = st_polygon(list(matrix(c(180,90, 180,-90, 0,-90, 0,90, 180,90), ncol=2,byrow=TRUE))) %>% 
  st_sfc() %>%
  st_set_crs("+proj=longlat +datum=WGS84")

ak_west = st_intersection(ak_geom, west_hem)
ak_east = st_intersection(ak_geom, east_hem)

ak_east_shift = (ak_east - c(360,0)) %>% st_set_crs("+proj=longlat +datum=WGS84")

ak_shift = st_union(ak_east_shift, ak_west)

plot(ak_shift, axes=TRUE, col="black", border=NA, xlim=c(-190, -130))
points(c(-360+179.776,-179.146), c(51.952,51.273),col='red')
abline(v=-180,col='blue',lty=2)
```

##

```{r echo=FALSE}
plot(ak_shift, axes=TRUE, col="black", border=NA, xlim=c(-190, -170), ylim=c(50, 55))
points(c(-360+179.776,-179.146), c(51.952,51.273),col='red')
abline(v=-180,col='blue',lty=2)
```


##

```{r}
path = geosphere::gcIntermediate(
  c(179.776, 51.952), c(-179.146, 51.273), 
  n=50, addStartEnd=TRUE
)
```

```{r echo=FALSE, fig.align="center", fig.width=8, fig.height=4}
par(mar=c(0,0,0,0))
plot(st_geometry(world), col="black", ylim=c(-90,90), axes=TRUE)
lines(path,col='red',lwd=2,lty=3)
```



# Using `sf`

## Example data

::: {.small}
```{r}
nc  = read_sf("data/gis/nc_counties/", quiet=TRUE)
air = read_sf("data/gis/airports/", quiet=TRUE)
hwy = read_sf("data/gis/us_interstates/", quiet=TRUE)
```

```{r}
nc
```
:::

##

::: {.small}
```{r}
air
```
:::

##

::: {.small}
```{r}
hwy
```
:::

## `sf` structure

::: {.small}
```{r}
str(nc)
```
:::

## `sf` classes

:::{.small}
```{r}
class(nc)

class(nc$geometry)

class(nc$geometry[[1]])
```
:::


## Projections

::: {.small}
```{r}
st_crs(nc)
```
:::

::: {.aside}
For EPSG details see [here](https://en.wikipedia.org/wiki/EPSG_Geodetic_Parameter_Dataset)
:::

##

::: {.small}
```{r}
st_crs(hwy)
```
:::

## UTM Zones

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


## Lat/Long

```{r}
nc_ll = nc
air_ll = air
hwy_ll = st_transform(hwy, st_crs(nc))
```

```{r echo=FALSE, fig.align="center", fig.height=3.5, fig.width=8}
par(mar=c(3,3,3,0.1), mfrow=c(1,3))
plot(st_geometry(nc_ll),  axes=TRUE, main="nc")
plot(st_geometry(air_ll), axes=TRUE, pch=16, col="blue", main="air")
plot(st_geometry(hwy_ll), axes=TRUE, col="red", main="hwy")
```


## UTM

```{r}
nc_utm = st_transform(nc, st_crs(hwy))
air_utm = st_transform(air, st_crs(hwy))
hwy_utm = hwy
```


```{r echo=FALSE, fig.align="center", fig.height=3.5, fig.width=8}
par(mar=c(3,3,3,0.1), mfrow=c(1,3))
plot(st_geometry(nc_utm),  axes=TRUE, main="nc")
plot(st_geometry(air_utm), axes=TRUE, pch=16, col="blue", main="air")
plot(st_geometry(hwy_utm), axes=TRUE, col="red", main="hwy")
```

## Comparison

::: {.panel-tabset .small}
### Lat/Long

```{r}
par(mar=c(3,5,0.1,0.1), las=1)
plot(st_geometry(nc_ll),  axes=TRUE)
plot(st_geometry(hwy_ll), col="red", add=TRUE)
plot(st_geometry(air_ll), pch=16, col="blue", add=TRUE)
```

### UTM

```{r}
par(mar=c(3,5,0.1,0.1), las=1)
plot(st_geometry(nc_utm),  axes=TRUE)
plot(st_geometry(hwy_utm), col="red", add=TRUE)
plot(st_geometry(air_utm), pch=16, col="blue", add=TRUE)
```
:::