Area-weighted index standardization

class: center, middle, inverse, title-slide

.title[
# Area-weighted index standardization
]
.subtitle[
## IMR sdmTMB workshop
]
.author[
### 
]
.date[
### May 23–25 2023
]

---

# Why model-based index standardization?

* Design-based methods (stratified means) are straightforward but assume perfect survey implementation

* Model-based approaches can account for (some) aspects of imperfect implementation (ref WKUSER2)

* Model-based approaches can leverage spatial correlation to improve precision

* Model-based approaches can integrate additional data sources, multiple gear types, accounting for gear calibration, etc.

---

# Design- vs. model-based estimators

* Design-based: fish are homogenously distributed within strata; noise comes from sampling

* Model-based: fish are distributed according to a mathematical model; noise from uncertainty in fish distribution and sampling

---

# Calculating an area-weighted population index using sdmTMB

1. Fit a sensible model: .blue[`sdmTMB()`]

2. Predict on a grid covering the survey domain: .blue[`predict(..., newdata = ...)`]

3. Sum up the density multiplied by cell area to calculate total abundance or biomass: .blue[`get_index()`]

---

# Fit a sensible model

```r
pcod$year_f <- as.factor(pcod$year)
mesh <- make_mesh(pcod, xy_cols = c("X", "Y"), cutoff = 10)
fit <- sdmTMB(
* density ~ 0 + year_f,
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "on",
  time = "year_f",
* spatiotemporal = "iid",
  silent = FALSE # show progress!
)
```

---

# Predict over the survey domain

prediction grid needs to be much smaller than the estimated spatial range
.small[

```r
qcs_grid1 <- purrr::map_dfr(seq_len(length(unique(pcod$year))), ~ qcs_grid)
qcs_grid1$year = rep(unique(pcod$year), each=nrow(qcs_grid))
qcs_grid1$year_f = as.factor(qcs_grid1$year)
pred <- predict(
  fit, newdata = qcs_grid1, return_tmb_object = TRUE
)
select(pred$data, year, X, Y, est) |> 
  head()
#>   year   X    Y      est
#> 1 2003 456 5636 1.639821
#> 2 2003 458 5636 1.682348
#> 3 2003 460 5636 1.724875
#> 4 2003 462 5636 1.767403
#> 5 2003 464 5636 1.889935
#> 6 2003 466 5636 2.361575
```
]

---

# Predict over the survey domain

---

# Sum up the density multiplied by cell area

(and calculate standard errors)

```r
index <- get_index(pred, area = 4, bias_correct = TRUE)
head(index)
#>   year_f       est       lwr       upr  log_est        se
#> 1   2003  936175.6  653699.1 1340715.9 13.74956 0.1832462
#> 2   2004 1832130.4 1359020.1 2469942.7 14.42099 0.1524087
#> 3   2005 1757227.0 1224189.4 2522360.5 14.37925 0.1844208
#> 4   2007  452112.1  328785.5  621698.2 13.02169 0.1625155
#> 5   2009  722981.6  518712.7 1007691.5 13.49114 0.1694080
#> 6   2011 1357885.0 1028850.9 1792146.6 14.12144 0.1415770
```

---

# The resulting standardized index

.xsmall[

```r
ggplot(index, aes(as.numeric(year_f), est)) + geom_line() +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4) +
  xlab('Year') + ylab('Biomass estimate (kg)')
```

<img src="07-index-standardization_files/figure-html/pcod-st-index-plot-1.png" width="600px" style="display: block; margin: auto;" />
]

---

# Remember...

There is *one* design-based index.

A good model-based estimated depends on the modeller constructing a reasonable model!

There are *many* possible model-based indexes. Thankfully, results are often (but not always!) qualitatively similar across configurations.

.tiny[
Commander, C.J.C, L.A.K. Barnett, E.J. Ward, S.C. Anderson, T.E. Essington. 2022. The shadow model: how and why small choices in spatially explicit species distribution models affect predictions. PeerJ. 10: e12783. <https://doi.org/10.7717/peerj.12783>
]