Example: Spatial modeling of Pacific cod

class: center, middle, inverse, title-slide

# Example: Spatial modeling of Pacific cod
## NOAA PSAW Seminar Series
### March 2, 2022

---

# The Pacific Cod dataset

DFO trawl survey data from Queen Charlotte Sound off BC

.small[

```r
library(dplyr)
library(ggplot2)
library(sdmTMB)
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
```
]

---

# Building a mesh

.small[

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

<img src="04-pcod-example_files/figure-html/unnamed-chunk-2-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Adding spatial random fields

```r
fit <- sdmTMB(
  present ~ 1, # intercept only
  data = pcod,
  mesh = mesh,
  family = binomial(link = "logit"),
* spatial = "on"
)
```

---

# Should we include depth?

* Predictions here in *logit* space

---

# Add depth as a quadratic effect

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

---

# Inspecting the model output

.small[

```r
fit
*#> Spatial model fit by ML ['sdmTMB']
#> Formula: present ~ poly(log(depth), 2)
#> Mesh: mesh
#> Data: pcod
#> Family: binomial(link = 'logit')
#>                      coef.est coef.se
#> (Intercept)             -0.45    0.43
#> poly(log(depth), 2)1   -74.39    8.87
#> poly(log(depth), 2)2  -107.16    8.76
#> 
#> Matern range: 43.54
#> Spatial SD: 1.65
#> ML criterion at convergence: 1042.157
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
```
]

---

# Inspecting the model output

.small[

```r
fit
#> Spatial model fit by ML ['sdmTMB']
*#> Formula: present ~ poly(log(depth), 2)
*#> Mesh: mesh
*#> Data: pcod
*#> Family: binomial(link = 'logit')
#>                      coef.est coef.se
#> (Intercept)             -0.45    0.43
#> poly(log(depth), 2)1   -74.39    8.87
#> poly(log(depth), 2)2  -107.16    8.76
#> 
#> Matern range: 43.54
#> Spatial SD: 1.65
#> ML criterion at convergence: 1042.157
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
```
]

---

# Inspecting the model output

.small[

```r
fit
#> Spatial model fit by ML ['sdmTMB']
#> Formula: present ~ poly(log(depth), 2)
#> Mesh: mesh
#> Data: pcod
#> Family: binomial(link = 'logit')
*#>                      coef.est coef.se
*#> (Intercept)             -0.45    0.43
*#> poly(log(depth), 2)1   -74.39    8.87
*#> poly(log(depth), 2)2  -107.16    8.76
#> 
#> Matern range: 43.54
#> Spatial SD: 1.65
#> ML criterion at convergence: 1042.157
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
```
]
---

# Inspecting the model output

.small[

```r
fit
#> Spatial model fit by ML ['sdmTMB']
#> Formula: present ~ poly(log(depth), 2)
#> Mesh: mesh
#> Data: pcod
#> Family: binomial(link = 'logit')
#>                      coef.est coef.se
#> (Intercept)             -0.45    0.43
#> poly(log(depth), 2)1   -74.39    8.87
#> poly(log(depth), 2)2  -107.16    8.76
#> 
*#> Matern range: 43.54
*#> Spatial SD: 1.65
*#> ML criterion at convergence: 1042.157
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
```
]

---

# Extracting parameters in a data frame

```r
# Main effects:
tidy(fit)
#>                   term     estimate std.error
#> 1          (Intercept)   -0.4501517 0.4308226
#> 2 poly(log(depth), 2)1  -74.3932665 8.8714966
#> 3 poly(log(depth), 2)2 -107.1584439 8.7597670

# Variance-related terms:
tidy(fit, effects = "ran_pars")
#>      term  estimate std.error
#> 1   range 43.536586        NA
#> 3 sigma_O  1.650969        NA
```

---

# Extracting parameters in a data frame

.small[

```r
tidy(fit, "ran_pars", conf.int = TRUE) %>% 
  as_tibble()
#> # A tibble: 2 × 5
#>   term    estimate std.error conf.low conf.high
#>   <chr>      <dbl> <lgl>        <dbl>     <dbl>
#> 1 range      43.5  NA           28.4      66.8 
#> 2 sigma_O     1.65 NA            1.27      2.14
```
]

---

# Making predictions

.small[

```r
*p <- predict(fit)
select(p, X, Y, est:omega_s) %>% 
  as_tibble()
#> # A tibble: 2,143 × 6
#>        X     Y     est est_non_rf  est_rf omega_s
#>    <dbl> <dbl>   <dbl>      <dbl>   <dbl>   <dbl>
#>  1  446. 5793.  0.0701      0.241 -0.171  -0.171 
#>  2  446. 5800. -0.184      -0.169 -0.0158 -0.0158
#>  3  449. 5802. -0.429      -0.480  0.0510  0.0510
#>  4  437. 5802.  1.21        0.383  0.826   0.826 
#>  5  421. 5771. -2.34       -1.98  -0.352  -0.352 
#>  6  418. 5772. -3.49       -3.64   0.147   0.147 
#>  7  408. 5771. -7.77       -9.01   1.25    1.25  
#>  8  408. 5767. -6.76       -7.96   1.21    1.21  
#>  9  414. 5761. -3.09       -3.27   0.184   0.184 
#> 10  400. 5749.  0.467      -2.60   3.07    3.07  
#> # … with 2,133 more rows
```
]

---

# Making predictions

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

---

# Making predictions on new data

.small[
* `qcs_grid`: a 2x2 km grid extending over the full survey domain
]
.small[

```r
nd <- filter(qcs_grid, year == 2017) # pick any year
ggplot(nd, aes(X, Y, fill = depth)) + 
  geom_raster() +
  coord_fixed()
```

<img src="04-pcod-example_files/figure-html/unnamed-chunk-10-1.png" width="500px" style="display: block; margin: auto;" />
]

---

# Making predictions on new data

```r
head(select(nd, X, Y, depth))
#>     X    Y    depth
#> 1 456 5636 347.0834
#> 2 458 5636 223.3348
#> 3 460 5636 203.7408
#> 4 462 5636 183.2987
#> 5 464 5636 182.9998
#> 6 466 5636 186.3892

*p <- predict(fit, newdata = nd)
```

---

# Plotting predictions on new data

.small[

```r
ggplot(p, aes(X, Y, fill = plogis(est))) +
  geom_raster() +
  scale_fill_viridis_c() +
  coord_fixed()
```

<img src="04-pcod-example_files/figure-html/unnamed-chunk-12-1.png" width="600px" style="display: block; margin: auto;" />
]

---

# Plotting main effect contributions

.small[

```r
ggplot(p, aes(X, Y, fill = plogis(est_non_rf))) +
  geom_raster() +
  scale_fill_viridis_c() +
  coord_fixed()
```

<img src="04-pcod-example_files/figure-html/unnamed-chunk-13-1.png" width="600px" style="display: block; margin: auto;" />
]

---

# Plotting spatial random effects

.small[

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

<img src="04-pcod-example_files/figure-html/unnamed-chunk-14-1.png" width="600px" style="display: block; margin: auto;" />
]

---

# Plotting depth effect

.small[

```r
nd <- data.frame(depth = seq(30, 500, length.out = 300))
p <- predict(
  fit, 
  newdata = nd, 
* se_fit = TRUE,
* re_form = ~ 0
)
p$lwr <- plogis(p$est - 1.96 * p$est_se)
p$upr <- plogis(p$est + 1.96 * p$est_se)
```

<img src="04-pcod-example_files/figure-html/pcod-sp-depth-pred-plot-1.png" width="500px" style="display: block; margin: auto;" />
]

---

# 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
```
]

---

# 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"
)
```

---

# Adding spatiotemporal fields

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

---

# 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 might look off even if the model is fine. Also see `dharma_residuals()` or MCMC residuals. See the ['Residual checking' vignette](https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html).
]

```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="04-pcod-example_files/figure-html/resid2-1.png" width="500px" style="display: block; margin: auto;" />
]

---

# Predicting on the survey grid

```r
p <- predict(fit, newdata = qcs_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="04-pcod-example_files/figure-html/pcod-st-plot-est-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="04-pcod-example_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="04-pcod-example_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="04-pcod-example_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 = qcs_grid, nsim = 500)
dim(psims)
#> [1] 65826   500
psims[1:3, 1:4]
#>           [,1]      [,2]      [,3]       [,4]
#> 2003 -4.281085 -2.569879 -3.618530 -3.6595682
#> 2003  1.475317  1.390593  1.308786  0.9383632
#> 2003  2.659430  2.241994  2.265373  1.9791928

*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="04-pcod-example_files/figure-html/pcod-st-sims-plot-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Calculating an area-weighted population index

.xsmall[

```r
*p2 <- predict(fit, newdata = qcs_grid, area = 4, return_tmb_object = TRUE)
*index <- get_index(p2, bias_correct = FALSE)
ggplot(index, aes(year, est)) + geom_line() +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4) +
  xlab('Year') + ylab('Biomass estimate (kg)')
```

<img src="04-pcod-example_files/figure-html/pcod-st-index1-1.png" width="600px" style="display: block; margin: auto;" />
]

---

# Calculating the center of gravity

.small[

```r
*cog <- get_cog(p2, bias_correct = FALSE, format = "wide")
ggplot(cog, aes(est_x, est_y, colour = year)) + 
  geom_linerange(aes(xmin = lwr_x, xmax = upr_x)) +
  geom_linerange(aes(ymin = lwr_y, ymax = upr_y)) +
  coord_equal()
```

<img src="04-pcod-example_files/figure-html/pcod-st-cog-1.png" width="600px" style="display: block; margin: auto;" />
]