Example: Spatial modeling of Pacific cod

class: center, middle, inverse, title-slide

.title[
# Example: Spatial modeling of Pacific cod
]
.subtitle[
## DFO TESA sdmTMB workshop
]
.author[
### 
]
.date[
### January 16–19 2023
]

---

# 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="03-pcod-spatial_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"
)
```

---

# Checking convergence

* Learning curve of TMB can be steep  
* Models can be very complex (hard to diagnose issues)
* `sanity()` function tries to help

.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 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)
#> # A tibble: 3 × 3
#>   term                 estimate std.error
#>   <chr>                   <dbl>     <dbl>
#> 1 (Intercept)            -0.450     0.431
#> 2 poly(log(depth), 2)1  -74.4       8.87 
#> 3 poly(log(depth), 2)2 -107.        8.76

# Variance-related terms:
tidy(fit, "ran_pars", conf.int = TRUE)
#> # 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)
#> # 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="03-pcod-spatial_files/figure-html/unnamed-chunk-9-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="03-pcod-spatial_files/figure-html/unnamed-chunk-11-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="03-pcod-spatial_files/figure-html/unnamed-chunk-12-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="03-pcod-spatial_files/figure-html/unnamed-chunk-13-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="03-pcod-spatial_files/figure-html/pcod-sp-depth-pred-plot-1.png" width="500px" style="display: block; margin: auto;" />
]

---

# Plotting depth effect with visreg

.small[

```r
visreg::visreg(fit, xvar = "depth", gg = TRUE)
```

<img src="03-pcod-spatial_files/figure-html/pcod-sp-depth-visreg-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Plotting depth effect with visreg

.small[

```r
visreg::visreg(fit, xvar = "depth", scale = "response")
```

<img src="03-pcod-spatial_files/figure-html/pcod-sp-depth-visreg-resp-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Plotting depth effect with ggeffects

.small[

```r
g <- ggeffects::ggeffect(fit, "depth [10:500 by=2]")
plot(g)
```

<img src="03-pcod-spatial_files/figure-html/pcod-sp-depth-ggeffects-1.png" width="700px" style="display: block; margin: auto;" />
]