---
title: "Lecture 22" 
subtitle: "Spatio-temporal Models"
author: "Colin Rundel"
date: "11/28/2018"
fontsize: 11pt
output: 
  beamer_presentation:
    theme: metropolis
    highlight: pygments
    fig_caption: false
    latex_engine: xelatex
    keep_tex: true
    includes:
      in_header: ../settings.tex
---

```{r setup, include=FALSE}
library(raster)
library(dplyr)
library(ggplot2)
library(patchwork)
library(sf)
library(rstan)

set.seed(20180405)

knitr::opts_chunk$set(
  collapse = TRUE,
  fig.width=7,
  fig.height=4.5,
  out.width="\\textwidth",
  fig.align="center",
  echo=TRUE,
  warning=FALSE
)

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

rstan::rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

source("../util.R")
```


# Spatial Models with AR time dependence

## Example - Weather station data

\footnotesize

Based on Andrew Finley and Sudipto Banerjee's notes from [National Ecological Observatory Network (NEON) Applied Bayesian Regression Workshop, March 7 - 8, 2013](http://blue.for.msu.edu/NEON/SC/) [Module 6](http://blue.for.msu.edu/NEON/SC/exercises/exercise-6/initial-exploration-spDynLM.pdf)

`NETemp.dat` - Monthly temperature data (Celsius) recorded across the Northeastern US starting in January 2000.

\scriptoutput

```{r echo=FALSE}
data("NETemp.dat", package = "spBayes")
ne_temp = NETemp.dat %>%
  filter(UTMX > 5.5e6, UTMY > 3e6) %>%
  mutate(x=UTMX/1000, y=UTMY/1000) %>%
  select(-UTMX, -UTMY) %>%
  select(x, y, 1:25) %>%
  tbl_df()

names(ne_temp) = stringr::str_replace(names(ne_temp), "y\\.", "t_")

ne_temp
```

##

```{r echo=FALSE}
ne_temp %>%
  select(x, y, t_1, t_4, t_7, t_10) %>%
  tidyr::gather(month, temp, -x, -y) %>%
  ggplot() +
    geom_point(aes(x=x, y=y, color=temp)) +
    facet_wrap(~forcats::as_factor(month))
```

##

```{r echo=FALSE}
ne_temp %>%
  mutate(station = as.factor(as.character(1:n()))) %>%
  slice(c(1, 10, 20, 30)) %>%
  tidyr::gather(month, temp, -x, -y, -station, -elev) %>%
  mutate(month = stringr::str_replace(month,"t_","") %>% as.integer()) %>%
  ggplot() +
    geom_point(aes(x=month, y=temp, color=station)) +
    geom_line(aes(x=month, y=temp, color=station, group=station))
```


## Dynamic Linear / State Space Models (time)

$$ 
\begin{aligned}
{{y}_t} &= \underset{1 \times p}{\symbf{F}'_t} ~ \underset{p \times 1}{\symbf{\theta}_t} + {{v}_t} 
&\qquad\qquad\text{observation equation}\\
\underset{p\times 1}{\symbf{\theta}_t} &= \underset{p \times p}{\symbf{G}_t} ~ \underset{p \times 1}{\symbf{\theta}_{t-1}}+ \underset{p \times 1}{\symbf{\omega}_t}
&\qquad\qquad\text{evolution equation}\\ 
\end{aligned}
$$

$$ 
\begin{aligned}
\symbf{v}_t &\sim \mathcal{N}(0,\symbf{V}_t) \\
\symbf{\omega}_t &\sim \mathcal{N}(0,\symbf{W}_t) \\
\end{aligned}
$$

## DLM vs ARMA {.t}

ARMA / ARIMA are a special case of a dynamic linear model, for example an $AR(p)$ can be written as
$$ F_t' = (1, 0, \ldots, 0) $$
$$
G_t = \begin{pmatrix}
\phi_1 & \phi_2 & \cdots & \phi_{p-1} & \phi_p \\
1      & 0      & \cdots & 0          & 0      \\
0      & 1      & \cdots & 0          & 0      \\
\vdots & \vdots & \ddots & \vdots     & 0      \\
0      & 0      & \cdots & 1          & 0      \\
\end{pmatrix}
$$
$$ 
\begin{aligned}
\omega_t &= (\omega_1, 0, \ldots, 0), 
\quad &\omega_1 \sim \mathcal{N}(0,\,\sigma^2)
\end{aligned}
$$

. . .

$$
\begin{aligned}
y_t &= \theta_t + v_t,
  \quad &v_t \sim \mathcal{N}(0,\, \sigma^2_v) \\
\theta_t &= \sum_{i=1}^p \phi_i\, \theta_{t-i} + \omega_1, 
  \quad &\omega_1 \sim \mathcal{N}(0,\, \sigma^2_\omega) \\
\end{aligned}
$$


## Dynamic spatio-temporal model

\vspace{2mm}

The observed temperature at time $t$ and location $s$ is given by $y_t(s)$ where,
\footnotesize
$$
\begin{aligned}
y_t(\symbf{s}) & = \symbf{x}_t(\symbf{s})\symbf{\beta}_t + u_t(\symbf{s}) + \epsilon_t(\symbf{s}) \\
\epsilon_t(\symbf{s}) &\stackrel{ind.}\sim \mathcal{N}(0,\tau_{t}^2) \\
\\
\symbf{\beta}_t & = \symbf{\beta}_{t-1} + \symbf{\eta}_t \\
\symbf{\eta}_t &\stackrel{i.i.d.}\sim \mathcal{N}(0,\symbf{\Sigma}_{\eta}) \\
\\
u_t(\symbf{s}) &= u_{t-1}(\symbf{s}) + w_t(\symbf{s}) \\
w_t(\symbf{s}) &\stackrel{ind.}{\sim} \mathcal{N}\left(\symbf{0}, \Sigma_t(\phi_t, \sigma^2_t)\right)
\end{aligned}
$$

\vspace{3mm}

. . .

\normalsize
Additional assumptions for $t=0$, 
\footnotesize
$$
\symbf{\beta}_{0} \sim \mathcal{N}(\symbf{\mu}_0, \symbf{\Sigma}_0)
$$
$$
u_{0}(\symbf{s}) = 0
$$

## Variograms by time

```{r echo=FALSE, message=FALSE, fig.height=5}
coords = ne_temp %>% 
  select(x, y) %>% 
  as.matrix()

par(mfrow=c(2,2), mar=c(2,2,2,2))

lm(t_1~elev, data=ne_temp)$residuals %>%
  geoR::variog(data=., coords=coords, uvec=seq(0,200, length.out=30), messages=FALSE) %>% 
  plot(main="Jan 2000", ylim=c(0,5))

lm(t_4~elev, data=ne_temp)$residuals %>%
  geoR::variog(data=., coords=coords, uvec=seq(0,200, length.out=30), messages=FALSE) %>% 
  plot(main="Apr 2000", ylim=c(0,5))

lm(t_7~elev, data=ne_temp)$residuals %>%
  geoR::variog(data=., coords=coords, uvec=seq(0,200, length.out=30), messages=FALSE) %>% 
  plot(main="Jul 2000", ylim=c(0,5))

lm(t_10~elev, data=ne_temp)$residuals %>%
  geoR::variog(data=., coords=coords, uvec=seq(0,200, length.out=30), messages=FALSE) %>% 
  plot(main="Oct 2000", ylim=c(0,5))
```

## Data and Model Parameters

**Data**:
\scriptoutput
```{r}
max_d = coords %>% dist() %>% max()
n_t = 24
n_s = nrow(ne_temp)
```

**Parameters**:
\scriptoutput
```{r}
n_beta = 2
starting = list(
  beta = rep(0, n_t * n_beta), phi = rep(3/(max_d/4), n_t),
  sigma.sq = rep(1, n_t), tau.sq = rep(1, n_t), 
  sigma.eta = diag(0.01, n_beta)
)
tuning = list(phi = rep(1, n_t))
priors = list(
  beta.0.Norm = list(rep(0, n_beta), diag(1000, n_beta)), 
  phi.Unif = list(rep(3/(0.9 * max_d), n_t), rep(3/(0.05 * max_d), n_t)), 
  sigma.sq.IG = list(rep(2, n_t), rep(2, n_t)), 
  tau.sq.IG = list(rep(2, n_t), rep(2, n_t)),
  sigma.eta.IW = list(2, diag(0.001, n_beta))
)
```

## Fitting with `spDynLM` from `spBayes`

\scriptoutput
```{r eval=FALSE}
n_samples = 10000
models = lapply(paste0("t_",1:24, "~elev"), as.formula)

m = spBayes::spDynLM(
  models, data = ne_temp, coords = coords, get.fitted = TRUE,
  starting = starting, tuning = tuning, priors = priors,
  cov.model = "exponential", n.samples = n_samples, n.report = 1000)

m = clean_spdynlm(m, n_samples/2+1, n_samples, (n_samples/2)/1000)
save(m, file="dynlm.Rdata")

##  ----------------------------------------
##  	General model description
##  ----------------------------------------
##  Model fit with 34 observations in 24 time steps.
##  
##  Number of missing observations 0.
##  
##  Number of covariates 2 (including intercept if specified).
##  
##  Using the exponential spatial correlation model.
##  
##  Number of MCMC samples 10000.
##
##  ...
```

```{r echo=FALSE}
load(file="dynlm.Rdata")
```

## Posterior Inference - $\beta$s

\vspace{4mm}

```{r echo=FALSE, fig.height=4}
betas = m$p.beta.samples %>% 
  post_summary() %>%
  mutate(
    month = stringr::str_extract(param,"[0-9]+") %>% as.integer(),
    param = stringr::str_extract(param,"\\(.*\\)|elev")
  )

ggplot(betas, aes(x=month, y=post_mean, color=param)) +
  geom_point(size=2) +
  geom_linerange(aes(ymin = post_lower, ymax = post_upper)) +
  facet_wrap(~param, scale="free_y")
```

\vvfill


[Lapse Rate](https://en.wikipedia.org/wiki/Lapse_rate) $\approx -9.8~^\circ C/km$.


## Posterior Inference - $\theta$

\vspace{4mm}

```{r echo=FALSE, fig.height=4}
theta = m$p.theta.samples %>% 
  post_summary() %>%
  mutate(
    month = stringr::str_extract(param,"[0-9]+") %>% as.integer(),
    param = stringr::str_extract(param,"sigma\\.sq|tau\\.sq|phi")
  )
eff_range = theta %>% 
  filter(param == "phi") %>%
  mutate(
    param = "Eff. Range",
    post_mean = 3/post_mean,
    post_med = NA,
    tmp = post_lower,
    post_lower = 3/post_upper,
    post_upper = 3/tmp,
    tmp = NULL
  )

theta = rbind(theta, eff_range)

ggplot(theta, aes(x=month, y=post_mean, color=param)) +
  geom_point(size=2) +
  geom_linerange(aes(ymin = post_lower, ymax = post_upper)) +
  facet_wrap(~param, scale="free_y", nrow=3)
```


## Posterior Inference - Observed vs. Predicted

\vspace{4mm}

```{r echo=FALSE, fig.height=5, fig.width=7, fig.align="center",out.width="\\textwidth"}
y_hat = m$p.y.samples %>% 
  t() %>% coda::as.mcmc() %>%
  post_summary() %>%
  bind_cols(
    tidyr::gather(ne_temp, month, temp, t_1:t_24),
    .
  )

ggplot(y_hat, aes(x=temp, y=post_mean)) +
  geom_point(size=1, alpha=0.5) +
  ylab("posterior mean") +
  geom_abline(slope = 1, intercept=0, color="black", alpha=0.2, size=2) +
  ylim(-20,25) + xlim(-20,25)
```

## Prediction

`spPredict` does not support `spDynLM` objects but it will impute missing values.

```{r}
r = raster(xmn=5750, xmx=6300, ymn=3000, ymx=3550, nrow=20, ncol=20)

pred = xyFromCell(r, 1:length(r)) %>% 
  as.data.frame() %>%
  mutate(type="pred") %>%
  bind_rows(
    ne_temp %>% mutate(type = "obs"),
    .
  )
```  

##

```{r echo=FALSE}
coords_pred = pred %>% select(x, y) %>% as.matrix()
n_p = nrow(pred)

n_beta = 1
starting = list(
  beta = rep(0, n_t * n_beta), phi = rep(3/(max_d/4), n_t),
  sigma.sq = rep(1, n_t), tau.sq = rep(1, n_t), 
  sigma.eta = diag(0.01, n_beta)
)

tuning = list(phi = rep(1, n_t))

priors = list(
  beta.0.Norm = list(rep(0, n_beta), diag(1000, n_beta)), 
  phi.Unif = list(rep(3/(0.9 * max_d), n_t), rep(3/(0.05 * max_d), n_t)), 
  sigma.sq.IG = list(rep(2, n_t), rep(2, n_t)), 
  tau.sq.IG = list(rep(2, n_t), rep(2, n_t)),
  sigma.eta.IW = list(2, diag(0.001, n_beta))
)
```

```{r eval=FALSE}
models_pred = lapply(paste0("t_",1:n_t, "~1"), as.formula)

n_samples = 5000
m_pred = spBayes::spDynLM(
  models_pred, data = pred, coords = coords_pred, get.fitted = TRUE,
  starting = starting, tuning = tuning, priors = priors,
  cov.model = "exponential", n.samples = n_samples, n.report = 1000)

m_pred = clean_spdynlm(m_pred, n_samples/2+1, n_samples, thin = 5)
save(m_pred, file="dynlm_pred.Rdata")

## ----------------------------------------
## 	General model description
## ----------------------------------------
## Model fit with 434 observations in 24 time steps.
## 
## Number of missing observations 9600.
## 
## Number of covariates 1 (including intercept if specified).
## 
## Using the exponential spatial correlation model.
## 
## Number of MCMC samples 5000.
```

```{r echo=FALSE}
load("dynlm_pred.Rdata")
```

## 

```{r echo=FALSE}
y_hat_pred = m_pred$p.y.samples %>% 
  t() %>% coda::as.mcmc() %>%
  post_summary() %>%
  bind_cols(
    tidyr::gather(pred, month, temp, t_1:t_24),
    .
  )

ggplot(y_hat_pred) +
  geom_point(aes(y=post_mean, x=temp)) +
  geom_abline(slope = 1, intercept=0, color="black", alpha=0.2, size=2) +
  ylim(-15,20) + xlim(-15,20)
```

##

```{r echo=FALSE, fig.width=6, fig.align="center", out.width="0.8\\textwidth"}
r_t1 = r
r_t1[] = y_hat_pred %>% filter(type == "pred", month=="t_1") %>% .$post_mean

r_t4 = r
r_t4[] = y_hat_pred %>% filter(type == "pred", month=="t_4") %>% .$post_mean

r_t7 = r
r_t7[] = y_hat_pred %>% filter(type == "pred", month=="t_7") %>% .$post_mean

r_t10 = r
r_t10[] = y_hat_pred %>% filter(type == "pred", month=="t_10") %>% .$post_mean


par(mfrow=c(2,2), mar=c(1,1,3,4))

s = raster::stack(r_t1, r_t4, r_t7, r_t10)
names(s) = c("Jan 2000","Apr 2000","Jul 2000","Oct 2000")

rasterVis::levelplot(s, contour=TRUE)
```


## Out-of-sample validation

```{r echo=FALSE}
oos_full = mutate(ne_temp, type = "train", station=1:n())
oos_full$type[c(1,7,19,24,33)] = "test"

oos = oos_full %>%
  tidyr::gather(month, temp, t_1:t_24) %>%
  mutate(temp = ifelse(type=="train", temp, NA)) %>%
  tidyr::spread(month, temp) %>%
  arrange(station)
  
oos
```  

```{r eval=FALSE, echo=FALSE}
n_samples = 10000
m_oos = spBayes::spDynLM(
  lapply(paste0("t_",1:n_t, "~1"), as.formula),
  data = oos, coords = coords, get.fitted = TRUE,
  starting = starting, tuning = tuning, priors = priors,
  cov.model = "exponential", n.samples = n_samples, n.report = 1000)

m_oos = clean_spdynlm(m_oos, n_samples/2+1, n_samples, 10)
save(m_oos, file="dynlm_oos.Rdata")
```

```{r echo=FALSE}
load("dynlm_oos.Rdata")
```

## 

```{r echo=FALSE}
m_oos$p.y.samples %>% 
  t() %>% coda::as.mcmc() %>%
  post_summary() %>%
  bind_cols(
    tidyr::gather(oos_full, month, temp, t_1:t_24),
    .
  ) %>%
ggplot() +
  geom_point(aes(y=post_mean, x=temp, color=type)) +
  geom_abline(slope = 1, intercept=0, color="black", alpha=0.2, size=2) +
  ylim(-15,20) + xlim(-15,20)
```


# Spatio-temporal models for continuous time

## Additive Models {.t}

In general, spatiotemporal models will have a form like the following,

$$
\begin{aligned}
y(\symbf{s},{t}) 
  &= \underset{\text{mean structure}}{\mu(\symbf{s},{t})} + \underset{\text{error structure}}{{e}(\symbf{s},{t})} \\
  &= \underset{\text{Regression}}{\symbf{x}(\symbf{s},{t}) \, \symbf{\beta}(\symbf{s},{t})} + \underset{\text{Spatiotemporal RE}}{{w}(\symbf{s},{t})} + \underset{\text{White Noise}}{\epsilon(\symbf{s},{t})}
\end{aligned} 
$$

. . .

\vspace{5mm}

The simplest possible spatiotemporal model is one were assume there is no dependence between observations in space and time,

$$
w(\symbf{s},t) = \alpha(t) + \omega(\symbf{s})
$$

these are straight forward to fit and interpret but are quite limiting (no shared information between space and time).


## Spatiotemporal Covariance

Lets assume that we want to define our spatiotemporal random effect to be a single stationary Gaussian Process (in 3 dimensions$^\star$),
$$ 
\symbf{w}(\symbf{s},\symbf{t}) \sim \mathcal{N}\big(\symbf{0}, \symbf{\Sigma}(\symbf{s},\symbf{t})\big) 
$$
where our covariance function depends on both $\lVert s-s'\rVert$ and $\lvert t-t'\rvert$,
$$
\text{cov}(\symbf{w}(\symbf{s},\symbf{t}), \symbf{w}(\symbf{s}',\symbf{t}')) = c(\lVert s-s'\rVert, \lvert t-t'\rvert)
$$

. . .

* Note that the resulting covariance matrix $\Sigma$ will be of size $n_s \cdot n_t \times n_s \cdot n_t$.

    + Even for modest problems this gets very large (past the point of direct computability).

    + If $n_t = 52$ and $n_s = 100$ we have to work with a $5200 \times 5200$ covariance matrix


## Separable Models

One solution is to use a seperable form, where the covariance is the product of a valid 2d spatial and a valid 1d temporal covariance / correlation function,
$$
\text{cov}(\symbf{w}(\symbf{s},\symbf{t}), \symbf{w}(\symbf{s}',\symbf{t}')) = \sigma^2 \, \rho_1(\lVert \symbf{s}-\symbf{s}'\rVert;\symbf{\theta}) \, \rho_2(\lvert \symbf{t}-\symbf{t}' \rvert; \symbf{\phi})
$$
\pause
If we define our observations as follows (stacking time locations within spatial locations)
\footnotesize
$$
\symbf{w}(\symbf{s},\symbf{t}) = \big(
  w(\symbf{s}_1,t_1)     ,\, \cdots ,\, w(\symbf{s}_1,t_{n_t}) ,\,
  \cdots ,\,
  w(\symbf{s}_{n_s},t_1) ,\, \cdots ,\, w(\symbf{s}_{n_s},t_{n_t}) \big)^t
$$
\normalsize
then the covariance can be written as
\footnotesize
$$
\underset{n_s n_t \,\times\, n_s n_t}{\symbf{\Sigma}_w(\sigma^2, \theta, \phi)} = \sigma^2 \, \underset{n_s \,\times\, n_s}{\symbf{H}_s(\theta)} \otimes \underset{n_t \,\times\, n_t}{\symbf{H}_t(\phi)}
$$
\normalsize
where $\symbf{H}_s(\theta)$ and $\symbf{H}_t(\theta)$ are correlation matrices defined by
\footnotesize
$$
\begin{aligned}
\{\symbf{H}_s(\theta)\}_{ij} &= \rho_1(\lVert \symbf{s}_i - \symbf{s}_j \rVert; \theta) \\
\{\symbf{H}_t(\phi)\}_{ij} &= \rho_2(\lvert t_i - t_j \rvert; \phi) \\
\end{aligned}
$$


## Kronecker Product

Definition:
$$
\begin{aligned}
\underset{[m \times n]}{\symbf{A}} \otimes \underset{[p \times q]}{\symbf{B}} = \underset{[m \cdot p \times  n \cdot q]}{\begin{pmatrix}
a_{11} \symbf{B} & \cdots & a_{1n} \symbf{B} \\
\vdots        & \ddots & \vdots        \\
a_{m1} \symbf{B} & \cdots & a_{mn} \symbf{B} \\
\end{pmatrix}}
\end{aligned}
$$

\vspace{4mm}

. . .

Properties:
$$
\begin{aligned}
\symbf{A} \otimes \symbf{B}       &\ne \symbf{B} \otimes \symbf{A}  \qquad\text{(usually)} \\ 
(\symbf{A} \otimes \symbf{B})^t   &= \symbf{A}^t \otimes \symbf{B}^t \\ \\
\det(\symbf{A} \otimes \symbf{B}) &= \det(\symbf{B} \otimes \symbf{A}) \\ 
&=\det(\symbf{A})^{\text{rank}(\symbf{B})} \det(\symbf{B})^{\text{rank}(\symbf{A})} \\ \\
(\symbf{A} \otimes \symbf{B})^{-1} &= \symbf{A}^{-1} \otimes \symbf{B}^{-1}
\end{aligned}
$$



## Kronecker Product and MVN Likelihoods {.t}

If we have a spatiotemporal random effect with a separable form,
$$
\symbf{w}(\symbf{s},\symbf{t}) \sim \mathcal{N}(\symbf{0},\, \symbf{\Sigma}_w)
$$
$$
\symbf{\Sigma}_w = \sigma^2 \, \symbf{H}_s \otimes \symbf{H}_t
$$

then the likelihood for $\symbf{w}$ is given by

$$
-\frac{n}{2}\log 2\pi - \frac{1}{2} \log |\Sigma_w| - \frac{1}{2} \symbf{w}^t \symbf{\Sigma_w}^{-1} \symbf{w}
$$
$$
= -\frac{n}{2}\log 2\pi - \frac{1}{2} \log \left[ (\sigma^2)^{n_t \cdot n_s} |H_s|^{n_t} |H_t|^{n_s}\right] - \frac{1}{2\sigma^2} \symbf{w}^t (\symbf{H}_s^{-1} \otimes \symbf{H}_t^{-1}) \symbf{w}
$$

## Non-seperable Models {.t}

* Additive and separable models are still somewhat limiting

\vspace{2mm}

* Cannot treat spatiotemporal covariances as 3d observations

\vspace{2mm}

* Possible alternatives:

    * Specialized spatiotemporal covariance functions, i.e.
$$ 
\gamma(\symbf{s},\symbf{s}', t,t') 
= \sigma^2 (\lvert t - t'\rvert+1)^{-1} \exp\big(-\lVert\symbf{s}-\symbf{s}'\rVert (\lvert t-t' \rvert + 1)^{-\beta/2}\big)
$$
  
    * Mixtures of separable covariances, i.e.
$$
w(\symbf{s},t) = w_1(\symbf{s},t) + w_2(\symbf{s},t)
$$
