Forecasting with sdmTMB

# Forecasting with sdmTMB
## NOAA PSAW Seminar Series
### 
### March 9, 2022

---

# Types of prediction we might be interested in

* Time  
  * e.g., interpolating over missed survey year
  * e.g., forecasting future year
  
* Extrapolating in space
  * e.g., unsampled area (MPA, area beyond existing domain)

---

# Forecasting in time

---

# Predicting to missing/future years

* Need a model for time: e.g., can't predict with years as factors

* Options: 
  * AR(1) or random walk random fields
  * Time-varying intercept
  * Smoother on year (`s(year)`)
  * Ignore time (fixed)
  * Some combination of these

---

# AR(1) spatiotemporal field

```r
# missing and forecasted years:
extra_years <- c(
  2006, 2008, 2010, 2012, 2014, 2016,
  2018:2025
)
fit_ar1 <- sdmTMB(
  density ~ depth_scaled + depth_scaled2,
  time = "year",
* extra_time = extra_years,
* spatiotemporal = "AR1",
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "off",
  silent = FALSE
)
```
]

---

# Random walk spatiotemporal field

```r
fit_rw <- sdmTMB(
  density ~ depth_scaled + depth_scaled2,
  time = "year",
* extra_time = extra_years,
* spatiotemporal = "RW",
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "off",
  silent = FALSE
)
```
]

---

# Random walk intercept + AR(1) fields

```r
fit_rw_ar1 <- sdmTMB(
  density ~ 0 + depth_scaled + depth_scaled2,
  time = "year",
* time_varying = ~1,
* extra_time = extra_years,
* spatiotemporal = "AR1",
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "off",
  silent = FALSE
)
```
]

---

# Smoother on year + AR(1) fields

```r
fit_sm <- sdmTMB(
* density ~ s(year) + depth_scaled + depth_scaled2,
  time = "year",
* extra_time = extra_years,
* spatiotemporal = "AR1",
  data = pcod,
  mesh = mesh,
  family = tweedie(link = "log"),
  spatial = "off",
  silent = FALSE
)
```
]

---

# Comparing predicted density at a point in space

---

# AR(1) spatiotemporal fields evolve towards mean zero

---

# Random walk fields do not evolve towards mean

---

# Spatiotemporal field uncertainty grows without data

(Here AR(1) fields; random walk similar.)

---

# Forecasting and interpolating in time summary

* Use `extra_time` argument to fill in or forecast  
  .xsmall[(a more flexible interface for forecasting may be forthcoming)]
* Need a model for time
* AR(1) field processes revert towards mean
* Random walk field processes do not revert towards mean
* Smoothers should be used with caution when forecasting  
  .xsmall[(they continue whatever the basis functions were doing)]

---

# Spatial extrapolation

---

# Previously explored BEI tree dataset

---

# Mesh more uncertain in unsampled areas

---

# Re-fitting the model

Intercept only, 1 time slice

```r
fit <- sdmTMB(
  n ~ 1,
  data = dat,
  mesh = mesh,
  family = truncated_nbinom2(link = "log"),
)
```

---

# Predict to grid

* Slower, but include `se_fit = TRUE`
* Could also use `predict(..., nsim = 500)` and summarize matrix

```r
# makes all combinations of x and y:
newdf <- expand.grid(
  x = seq(min(dat$x), max(dat$x), 5),
  y = seq(min(dat$y), max(dat$y), 5)
)
p <- predict(fit,
  newdata = newdf, se_fit = TRUE
)
```

---
 
# Visualizing uncertainty

* Why is uncertainty higher at vertices?  
  * Fewer neighbors, e.g. [this tutorial](https://ourcodingclub.github.io/tutorials/spatial-modelling-inla/)
  
<img src="10-forecasting_files/figure-html/vis-vert-1.png" width="700px" style="display: block; margin: auto;" />

---
 
# Predicting outside the survey domain

* Predict into border area