Example: Spatiotemporal modeling of Pacific cod

class: center, middle, inverse, title-slide

.title[
# Example: Spatiotemporal modeling of Pacific cod
]
.subtitle[
## IMR sdmTMB workshop
]
.author[
### 
]
.date[
### May 23–25 2023
]

---

# Extending to a spatiotemporal model

Reminder, the Pacific Cod dataset:

.small[

```r
pcod %>% 
  select(year, X, Y, depth, density, present) %>% 
  head()
#> # A tibble: 6 × 6
#>    year     X     Y depth density present
#>   <int> <dbl> <dbl> <dbl>   <dbl>   <dbl>
#> 1  2003  446. 5793.   201   113.        1
#> 2  2003  446. 5800.   212    41.7       1
#> 3  2003  449. 5802.   220     0         0
#> 4  2003  437. 5802.   197    15.7       1
#> 5  2003  421. 5771.   256     0         0
#> 6  2003  418. 5772.   293     0         0

mesh <- make_mesh(pcod, xy_cols = c("X", "Y"), cutoff = 10)

pcod$fyear <- as.factor(pcod$year)# year factor
```
]

---

# Our previous spatial model fit

```r
fit <- sdmTMB(
  present ~ poly(log(depth), 2),
  data = pcod,
  mesh = mesh,
  family = binomial(link = "logit"),
  spatial = "on"
)
```

---

# Switching to density

```r
fit <- sdmTMB(
* density ~ poly(log(depth), 2),
  data = pcod,
  mesh = mesh,
* family = tweedie(link = "log"),
  spatial = "on"
)
```

---

# Switching to a smoother + annual mean

```r
fit <- sdmTMB(
* density ~ s(depth) + 0 + fyear,
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "on"
)
```

---

# Adding spatiotemporal fields

```r
fit <- sdmTMB(
  density ~ s(depth, k = 5) + 0 + fyear,
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "on",
* time = "year",
* spatiotemporal = "iid",
  silent = FALSE # show progress!
)
```

---

# Inspecting model convergence

.small[

```r
sanity(fit)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
```
]

---

# Inspecting the model fit

.small[

```r
fit
```

```r
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: density ~ s(depth) + 0 + as.factor(year)
#> Time column: "year"
#> Mesh: mesh
#> Data: pcod
#> Family: tweedie(link = 'log')
#>                     coef.est coef.se
*#> as.factor(year)2003     1.79    0.29
*#> as.factor(year)2004     2.36    0.27
#> ...
#> 
#> Dispersion parameter: 10.78
#> Tweedie p: 1.49
*#> Matern range: 12.85
#> Spatial SD: 1.81
*#> Spatiotemporal SD: 1.78
#> ML criterion at convergence: 6246.433
```
]

---

# Model residuals

.xsmall[
Warning: these residuals are fast but might look off even if the model is fine. Also see MCMC residuals. See the ['Residual checking' vignette](https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html) and the example in the exercises.
]

```r
set.seed(1)
rq_res <- residuals(fit) # randomized quantile residuals
qqnorm(rq_res);qqline(rq_res)
```

---

# Model residuals in space

.small[

```r
pcod$resids <- residuals(fit)
filter(pcod, year %in% c(2015, 2017)) %>% 
* ggplot(aes(X, Y, colour = resids)) +
  geom_point() +
  facet_wrap(~year) +
  scale_colour_gradient2() +
  coord_fixed()
```

<img src="05-pcod-spatiotemporal_files/figure-html/resid2-1.png" width="500px" style="display: block; margin: auto;" />
]

---

# Predicting on the survey grid

```r
data("qcs_grid")
all_grid <- replicate_df(qcs_grid, 
            "year", unique(pcod$year))
all_grid$fyear <- as.factor(all_grid$year)
p <- predict(fit, newdata = all_grid)
```

* `est`: Overall estimate in link space (log here)
* `est_non_rf`: Estimate of non-random-field components
* `est_rf`: Estimate of random-field components
* `omega_s`: Spatial random field
* `epsilon_st`: **Spatiotemporal random field**

---

# Plotting overall predictions

.xsmall[

```r
ggplot(p, aes(X, Y, fill = exp(est))) +
  geom_raster() +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~year) +
  coord_fixed()
```

<img src="05-pcod-spatiotemporal_files/figure-html/pcod-st-plot-est-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Plotting overall predictions (truncated)

.xsmall[

```r
max_est <- quantile(p$est, 0.999)
p <- mutate(p, est_trim = if_else(est > max_est, max_est, est))
```

<img src="05-pcod-spatiotemporal_files/figure-html/pcod-st-plot-est3-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Plotting main effect contributions

.xsmall[

```r
ggplot(p, aes(X, Y, fill = exp(est_non_rf))) +
  geom_raster() +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~year) +
  coord_fixed()
```

<img src="05-pcod-spatiotemporal_files/figure-html/pcod-st-plot-non-rf-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Plotting spatial random effects

.xsmall[

```r
ggplot(p, aes(X, Y, fill = omega_s)) +
  geom_raster() +
  scale_fill_gradient2() +
  facet_wrap(~year) +
  coord_fixed()
```

<img src="05-pcod-spatiotemporal_files/figure-html/pcod-st-plot-omega-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Plotting spatiotemporal random effects

.xsmall[

```r
ggplot(p, aes(X, Y, fill = epsilon_st)) +
  geom_raster() +
  scale_fill_gradient2() +
  facet_wrap(~year) +
  coord_fixed()
```

<img src="05-pcod-spatiotemporal_files/figure-html/pcod-st-plot-eps-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Visualizing spatial uncertainty

By sampling from the joint precision matrix

.small[

```r
*psims <- predict(fit, newdata = all_grid, nsim = 50)
dim(psims)
#> [1] 65826    50
psims[1:3, 1:4]
#>           [,1]      [,2]      [,3]      [,4]
#> 2003 -3.327548 -2.105384 -2.718930 -3.555146
#> 2003  2.056039  2.606377  1.703947  1.019264
#> 2003  3.102490  3.511679  2.837596  2.080360

*p$sd <- apply(psims, 1, sd)
```
]

---

# Visualizing spatial uncertainty

.small[

```r
filter(p, year %in% c(2013, 2015, 2017)) %>% 
  ggplot(aes(X, Y, fill = sd)) +
  geom_raster() +
  facet_wrap(~year) +
  scale_fill_viridis_c(trans = "log10") +
  coord_fixed()
```

<img src="05-pcod-spatiotemporal_files/figure-html/pcod-st-sims-plot-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Considering anisotropy

* Default is isotropic: correlation decays in all directions at same rate
* Anisotropic: directionally dependent spatial correlation

```r
fit_aniso <- update(fit, anisotropy = TRUE)
```

```r
AIC(fit, fit_aniso)
#>           df      AIC
#> fit       16 12524.87
#> fit_aniso 18 12528.46
```

.xsmall[
* Not favoured here; often important on narrow continental shelves
]

---

# Plotting anisotropy

.small[
Plot shows range in all directions from zero in the center

```r
plot_anisotropy(fit_aniso)
```

<img src="05-pcod-spatiotemporal_files/figure-html/pcod-aniso-plot-1.png" width="700px" style="display: block; margin: auto;" />
]