---
title: "GLMs +
Point referenced data"
subtitle: "Lecture 22"
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(dukestm)
library(sf)
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
{fig-align="center" width="100%"}
## Data
::: {.small}
```{r}
( loaloa = PrevMap::loaloa |>
as_tibble() |>
rename_with(tolower) |>
rename(elev=elevation)
)
```
:::
## Spatial Distribution
```{r}
#| echo: false
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
{fig-align="center" width="45%"}
::: {.aside}
See [here](https://www.usgs.gov/land-resources/nli/landsat/landsat-normalized-difference-vegetation-index) for more details on NDVI.
:::
## 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`, `max9901`, `min9901`, `stdev9901` - aggregated NDVI data from 1999 to 2001 from the Flemish Institute for Technological Research (1 km resolution)
## Diggle's Model
$$
\begin{aligned}
y(s) &\sim \text{Binom}(p(s), n(s)) \\
\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) \\
\\
w(s) &\sim N(0, \Sigma) \\
\{\Sigma\}_{ij} &= \sigma^2 \, \exp(-d \,\phi)
\end{aligned}
$$
## 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
{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,
data = loaloa
) |>
as_tibble()
```
:::
## OOS Validation
::: {.small}
```{r message=FALSE}
set.seed(12345)
loaloa_test = loaloa |> slice_sample(prop=0.20)
loaloa = anti_join(loaloa, loaloa_test, by="row")
```
:::
```{r echo=FALSE}
bind_rows(
loaloa |> mutate(type = "Training"),
loaloa_test |> mutate(type = "Testing")
) |>
mutate(
type = as_factor(type)
) |>
ggplot() +
geom_sf(data=w_africa) +
geom_point(aes(x=longitude, y=latitude, color=no_inf/no_exam), alpha=0.8) +
xlim(7.5,16) + ylim(2,7.2) +
facet_wrap(~type)
```
## 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
#### Training
```{r}
yardstick::rmse_vec(loaloa$no_inf/loaloa$no_exam, loaloa$glm_pred)
```
#### Testing
```{r}
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
)
```
```{r}
#| eval: false
#| include: false
saveRDS(ll_gp, file="Lec21_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="Lec21_loaloa.rds")
```
:::
```{r echo=FALSE}
ll_gp = readRDS("Lec21_loaloa.rds")
```
##
::: {.small}
```{r}
ll_gp
```
:::
## Diagnostics
```{r}
vars = colnames(ll_gp$mcmc)
plot(ll_gp, vars=vars[1:4])
```
##
```{r}
plot(ll_gp, vars=vars[5: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") +
coord_fixed()
```
## 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") +
coord_fixed()
```
## Diggle's Predictive Surface
{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
{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
(*Catharus guttatus*) | Wilson's Warbler
(*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%'}
{fig-align="center" width="50%"}
:::
::: {.column width='50%'}
{fig-align="center" width="50%"}
:::
::::
## Sampling Locations
{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} | f_{\cdot l k} &\sim \text{Multinom}\left(\boldsymbol{f}_{\cdot l k}, n_{l k}\right) \\
\\
f_{ilk} &= \frac{\exp(\Theta_{ilk})}{\sum_i \exp(\Theta_{ilk})} \\
\\
\boldsymbol{\Theta}_{il}|\boldsymbol{\alpha},\boldsymbol{\mu} &\sim 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)
{fig-align="center" width="25%"}
{fig-align="center" width="25%"}
{fig-align="center" width="25%"}
{fig-align="center" width="25%"}
{fig-align="center" width="25%"}
{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
{fig-align="center" width="100%"}
## Model Assessment
{fig-align="center" width="100%"}
## Migratory Connectivity
{fig-align="center" width="100%"}