---
title: "Point referenced data (pt. 2)"
subtitle: "Lecture 22"
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=70)

set.seed(20221109)

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


# Loa Loa Example

## Loa Loa

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


## Data

::: {.small}
```{r}
loaloa = PrevMap::loaloa %>%
  as_tibble() %>% 
  setNames(., tolower(names(.))) %>%
  rename(elev=elevation)

loaloa
```
:::

## Spatial Distribution

```{r echo=FALSE, out.width="0.9\\textwidth"}
w_africa = read_sf("data/AfricanCountries/", quiet = TRUE) %>%
  filter(COUNTRY %in% c("Cameroon", "Nigeria", "Central African Republic", "Gabon", "Congo", "Equatorial Guinea", "Chad"))

ggplot() +
  geom_sf(data=w_africa) +
  geom_point(data=loaloa, aes(x=longitude, y=latitude, size=no_exam, color=no_inf/no_exam), alpha=0.8) +
  xlim(7.5,16) + ylim(2,12) +
  scale_size(range=c(0.1,5))
```


## Normalized Difference Vegetation Index (NDVI)

![](imgs/Lec22/ndvi_cameroon.png){fig-align="center" width="50%"}


## Paper / Data summary

Original paper - Diggle, et. al. (2007). *Spatial modelling and prediction of Loa loa risk: decision making under uncertainty*. Annals of Tropical Medicine and Parasitology, 101, 499-509.

* `no_exam` and `no_inf` - Collected between 1991 and 2001 by NGOs (original paper mentions 168 villages and 21,938 observations)

* `elev` - USGS gtopo30 (1km resolution)

* `mean9901` to `stdev9901` - aggregated data from 1999 to 2001 from the Flemish Institute for Technological Research (1 km resolution)


## Diggle's Model

$$ 
\begin{aligned}
\log \left( \frac{p(s)}{1-p(s)} \right) = \alpha &+ f_1(\text{elev}(s)) \\
&+ f_2(\text{MAX.NDVI}(s)) \\
&+ f_3(\text{SD.NDVI}(s)) + w(s) 
\end{aligned}
$$

where

$$ w(s) \sim \mathcal{N}(0, \Sigma) $$
$$ \{\Sigma\}_{ij} = \sigma^2 \, \exp(-d \,\phi) $$


## EDA

```{r echo=FALSE}
logit = function(p) log(p/(1-p))
loaloa = loaloa %>%
  mutate(prop = no_inf/no_exam) %>%
  mutate(logit_prop = logit(prop)) %>%
  mutate(logit_prop = ifelse(is.finite(logit_prop), logit_prop, NA))

( ggplot(loaloa, aes(x=elev, y=logit_prop)) + geom_point() + 
  geom_vline(xintercept=c(1000,1300), color="red")) /
( ggplot(loaloa, aes(x=max9901, y=logit_prop)) + geom_point() + 
    geom_vline(xintercept=c(0.8), color="red")) / 
( ggplot(loaloa, aes(x=stdev9901, y=logit_prop)) + geom_point())
```

## Diggle's EDA

![](imgs/Lec22/diggle_eda.png){fig-align="center" width="33%"}

## Feature engineering

::: {.small}
```{r}
loaloa = loaloa %>% 
  mutate(
    elev_f = cut(elev, breaks=c(0,1000,1300,2000), dig.lab=5),
    max_f  = cut(max9901, breaks=c(0,0.8,1))
  )
loaloa %>% select(elev, elev_f, max9901, max_f)
```
:::

## Model Matrix

::: {.small}
```{r}
model.matrix(
  ~ elev:elev_f - 1, 
  data = loaloa
) %>%
  as_tibble()
```
:::


## OOS Validation

```{r message=FALSE}
set.seed(12345)
loaloa_test = loaloa %>% slice_sample(prop=0.20)
loaloa = anti_join(loaloa, loaloa_test, quiet=TRUE)
```

```{r echo=FALSE}
( ggplot() +
    geom_sf(data=w_africa) +
    geom_point(data=loaloa, aes(x=longitude, y=latitude, color=no_inf/no_exam), alpha=0.8) +
    xlim(7.5,16) + ylim(2,7.2) +
    labs(title="Training")
) +
( ggplot() +
    geom_sf(data=w_africa) +
    geom_point(data=loaloa_test, aes(x=longitude, y=latitude, color=no_inf/no_exam), alpha=0.8) +
    xlim(7.5,16) + ylim(2,7.2) +
    labs(title="Testing")
)
```

## Model

::: {.small}
```{r}
g = glm(no_inf/no_exam ~ elev:elev_f + max9901:max_f + stdev9901, 
        data=loaloa, family=binomial, weights=loaloa$no_exam)
summary(g)
```
:::

## Predictions - Training

```{r echo=FALSE}
loaloa = loaloa %>% 
  mutate(glm_pred = predict(g, type="response"))

loaloa_test = loaloa_test %>% 
  mutate(glm_pred = predict(g, newdata=loaloa_test, type="response"))
```


```{r echo=FALSE}
( ggplot() +
    geom_sf(data=w_africa) +
    geom_point(data=loaloa, aes(x=longitude, y=latitude, color=no_inf/no_exam), alpha=0.8) +
    xlim(7.5,16) + ylim(2,7.2) +
    labs(title="Data")
) +
( ggplot() +
    geom_sf(data=w_africa) +
    geom_point(data=loaloa, aes(x=longitude, y=latitude, color=glm_pred), alpha=0.8) +
    xlim(7.5,16) + ylim(2,7.2) +
    labs(title="GLM Prediction")
)
```

## Predictions - Testing

```{r echo=FALSE}
( ggplot() +
    geom_sf(data=w_africa) +
    geom_point(data=loaloa_test, aes(x=longitude, y=latitude, color=no_inf/no_exam), alpha=0.8) +
    xlim(7.5,16) + ylim(2,7.2) +
    labs(title="Data")
) +
( ggplot() +
    geom_sf(data=w_africa) +
    geom_point(data=loaloa_test, aes(x=longitude, y=latitude, color=glm_pred), alpha=0.8) +
    xlim(7.5,16) + ylim(2,7.2) +
    labs(title="GLM Prediction")
)
```

## Fit - Training

```{r echo=FALSE}
ggplot(loaloa, aes(x=no_inf/no_exam, y=glm_pred)) + 
  geom_point() +
  geom_abline(slope = 1, intercept = 0)
```

## Fit - Testing

```{r echo=FALSE, out.width="\\textwidth"}
ggplot(loaloa_test, aes(x=no_inf/no_exam, y=glm_pred)) + 
  geom_point() +
  geom_abline(slope = 1, intercept = 0)
```

## Fit - RMSE

```{r}
# Training
yardstick::rmse_vec(loaloa$no_inf/loaloa$no_exam, loaloa$glm_pred)

# Testing
yardstick::rmse_vec(loaloa_test$no_inf/loaloa_test$no_exam, loaloa_test$glm_pred)
```

## Spatial Structure?

```{r message=FALSE}
geoR::variog(coords = cbind(loaloa$longitude, loaloa$latitude), 
       data = loaloa$prop - loaloa$glm_pred,
       uvec = seq(0, 4, length.out = 50)) %>% plot()
```

## `gpglm` model

::: {.small}
```{r eval=FALSE}
ll_gp = gpglm(
  no_inf ~ scale(elev):elev_f + scale(max9901):max_f + scale(stdev9901), 
  data = loaloa, family="binomial", weights=loaloa$no_exam,
  coords = c("longitude", "latitude"),
  cov_model="exponential",
  starting = list(
    beta=rep(0,7), 
    phi=3, sigma.sq=1, w=0
  ),
  priors = list(
    beta.Normal=list(rep(0,7), rep(10,7)),
    phi.unif=c(3/4, 3/0.25), sigma.sq.ig=c(2, 2)
  ),
  tuning = list(
    "beta"=rep(0.1, 7),
    "phi"=0.6, "sigma.sq"=0.3, "w"=0.1
  ),
  n_batch = 400,
  batch_len = 50,
  verbose = TRUE,
  n_report = 10,
  chains=4
)

saveRDS(ll_gp, file="Lec22_loaloa.rds")
```

```{r eval=FALSE, include=FALSE}
ll_gp = gpglm(
  no_inf ~ scale(elev):elev_f + scale(max9901):max_f + scale(stdev9901), 
  data = loaloa, family="binomial", weights=loaloa$no_exam,
  coords = c("longitude", "latitude"),
  cov_model="exponential",
  starting = list(
    beta=rep(0,7), 
    phi=3, sigma.sq=1, w=0
  ),
  priors = list(
    beta.Normal=list(rep(0,7), rep(10,7)),
    phi.unif=c(3/4, 3/0.25), sigma.sq.ig=c(2, 2)
  ),
  tuning = list(
    "beta"=rep(0.1, 7),
    "phi"=0.6, "sigma.sq"=0.3, "w"=0.1
  ),
  n_batch = 400,
  batch_len = 50,
  verbose = TRUE,
  n_report = 10,
  chains=4
)


saveRDS(ll_gp, file="Lec22_loaloa.rds")
```
:::

```{r echo=FALSE}
ll_gp = readRDS("Lec22_loaloa.rds")
```

##

::: {.small}
```{r}
ll_gp
```
:::

## Diagnostics

```{r}
plot(ll_gp, vars=1:3)
```

##

```{r}
plot(ll_gp, vars=4:6)
```

##

```{r}
plot(ll_gp, vars=7:9)
```

## Prediction (training)

::: {.small}
```{r message=FALSE}
ll_gp_pred = predict(
  ll_gp, 
  newdata=loaloa, 
  coords = c("longitude", "latitude"),
  thin = 25,
  verbose=FALSE
)

ll_gp_pred_y = tidybayes::gather_draws(ll_gp_pred, y[i]) %>%
  group_by(.chain, i) %>%
  summarize(
    post_p = mean(.value),
    q025 = quantile(.value, 0.025),
    q975 = quantile(.value, 0.975)
  )
```
:::

##

```{r echo=FALSE}
ll_gp_pred_y %>%
  inner_join(
    loaloa %>% mutate(i = seq_len(n())),
    by = "i"
  ) %>%
  ggplot(aes(x=no_inf/no_exam, y=post_p, color=as.factor(.chain))) + 
    geom_point() +
    geom_abline(slope = 1, intercept = 0) +
    geom_errorbar(aes(ymin=q025, ymax=q975), alpha=0.1)+
    guides(color="none")
```



## Prediction - Testing

::: {.small}
```{r message=FALSE}
ll_gp_test_pred = predict(
  ll_gp, 
  newdata=loaloa_test, 
  coords = c("longitude", "latitude"),
  thin = 25,
  verbose=FALSE
)

ll_gp_test_pred_y = tidybayes::gather_draws(ll_gp_test_pred, y[i]) %>%
  group_by(.chain, i) %>%
  summarize(
    post_p = mean(.value),
    q025 = quantile(.value, 0.025),
    q975 = quantile(.value, 0.975)
  )
```
:::

##

```{r echo=FALSE}
ll_gp_test_pred_y %>%
  inner_join(
    loaloa %>% mutate(i = seq_len(n())),
    by = "i"
  ) %>%
  ggplot(aes(x=no_inf/no_exam, y=post_p, color=as.factor(.chain))) + 
    geom_point() +
    geom_abline(slope = 1, intercept = 0) +
    geom_errorbar(aes(ymin=q025, ymax=q975), alpha=0.1) +
    guides(color="none")
```

## Diggle's Predictive Surface

![](imgs/Lec22/diggle_fig2.png){fig-align="center"}


## Exceedance Probability

```{r echo=FALSE}
village_i = c(13, 78, 6, 151)

village_pred = predict(
  ll_gp, 
  newdata=loaloa[village_i,], 
  coords = c("longitude", "latitude"),
  thin = 25,
  verbose=FALSE
) %>%
  tidybayes::gather_draws(y[i]) %>%
  mutate(village = paste0("Village ", i))

ggplot(village_pred, aes(x=.value, fill=village)) +
    geom_density() +
    facet_wrap(~village) +
    guides(fill="none")
```

## Exceedance Probability Predictive Surface

![](imgs/Lec22/diggle_fig4.png){fig-align="center"}


# Spatial Assignment of Migratory Birds

## Background

Using intrinsic markers (genetic and isotopic signals) for the purpose of inferring migratory connectivity.

* Existing methods are too coarse for most applications

* Large amounts of data are available ( \textgreater 150,000 feather samples from \textgreater 500 species)

* Genetic assignment methods are based on Wasser, et al. (2004)

* Isotopic assignment methods are based on Wunder, et al. (2005)



## Data - DNA microsatellites and $\delta {}^2{H}$


| Hermit Thrush<br/>(*Catharus guttatus*) | Wilson's Warbler<br/>(*Wilsonia pusilla*) |
|:------------------------------------|:--------------------------------------|
| 138 individuals                     | 163 individuals                       |
| 14 locations                        | 8 locations                           |
| 6 loci                              | 9 loci                                |
| 9-27 alleles / locus                | 15-31 alleles / locus                 |

:::: {.columns}
::: {.column width='50%'}
![](imgs/Lec22/hermit_thrush.jpeg){fig-align="center" width="50%"}
:::

::: {.column width='50%'}
![](imgs/Lec22/wilsons_warbler.jpeg){fig-align="center" width="50%"}
:::
::::



## Sampling Locations

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

## Allele Frequency Model

For the allele $i$, from locus $l$, at location $k$

$$
\begin{aligned}
\boldsymbol{y}_{\cdot l k}|\boldsymbol{\Theta} &\sim \mathcal{N}\left(\textstyle\sum_i y_{ilk},\: \boldsymbol{f}_{\cdot l k}\right) \\
\\
f_{ilk} &= \frac{\exp(\Theta_{ilk})}{\sum_i \exp(\Theta_{ilk})} \\
\\
\boldsymbol{\Theta}_{il}|\boldsymbol{\alpha},\boldsymbol{\mu} &\sim \mathcal{N}( \boldsymbol{\mu}_{il},\, \boldsymbol{\Sigma_{}}) \\
\end{aligned}
$$

$$
\left\{\Sigma\right\}_{ij} = \sigma^2 \, \exp \Big(-(\{d\}_{ij}\, r)^{\psi} \Big) + \sigma^2_n \, {1}_{i=j}
$$

## Predictions by Allele (Locus 3)

![](imgs/Lec22/allele3/Med-Al3-1.png){fig-align="center" width="25%"}
![](imgs/Lec22/allele3/Med-Al3-2.png){fig-align="center" width="25%"}
![](imgs/Lec22/allele3/Med-Al3-3.png){fig-align="center" width="25%"}
![](imgs/Lec22/allele3/Med-Al3-4.png){fig-align="center" width="25%"}
![](imgs/Lec22/allele3/Med-Al3-5.png){fig-align="center" width="25%"}
![](imgs/Lec22/allele3/Med-Al3-6.png){fig-align="center" width="25%"}


## Genetic Assignment Model

Assignment model assuming Hardy-Weinberg equilibrium and allowing for genotyping ($\delta$) and single amplification ($\gamma$) errors.

$$
\begin{aligned}
P(S_G|\boldsymbol{f},k) &= \prod_l P(i_l, j_l | \boldsymbol{f},k) \\
\\
P(i_l, j_l | \boldsymbol{f},k) &= 
\begin{cases}
  \gamma P(i_l|\boldsymbol{f},k) + (1-\gamma)P(i_l|\boldsymbol{\tilde f},k)^2 & \text{if $i=j$} \\
  (1-\gamma) P(i_l|\boldsymbol{f},k) P(j_l|\boldsymbol{f},k)                  & \text{if $i \ne j$}
\end{cases} \\
\\
P(i_l|\boldsymbol{f},k) &= (1-\delta) f_{lik} + \delta / m_l
\end{aligned}
$$


## Combined Model

<br/>

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



## Model Assessment

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



## Migratory Connectivity

<br/>

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