Time-varying coefficients with sdmTMB

class: center, middle, inverse, title-slide

.title[
# Time-varying coefficients with sdmTMB
]
.subtitle[
## DFO TESA sdmTMB workshop
]
.author[
### 
]
.date[
### January 16–19 2023
]

---

.small[
# Why might we want time-varying effects?

* Time-varying slopes: 
  * To allow for evolving responses to covariates (e.g., species moving deeper over time)
  
  * Example use: [English et al. (2021) Fish and Fisheries](https://doi.org/10.1111/faf.12613)  
    Modelled groundfish density with depth; didn't want to constrain fish if they were moving deeper when water was warmer
  
* Time-varying intercepts:
  * To allow variable means across time with constraints
  * To have a model to interpolate or forecast over time
]

---

# Time-varying intercepts

Several ways in sdmTMB:

* factors: `as.factor(year)` (independent)
* random effects: ` + (1 | year)` (drawn from normal distribution)
* smooth: ` + s(year)`
* as random walk (shown next)

---

#  Random walk covariates in sdmTMB

Random walk:

$$
`\begin{aligned}
x_t &= x_{t-1} + \epsilon_t\\
\epsilon &\sim \mathrm{Normal(0, \sigma^2)}
\end{aligned}`
$$

Defined by `time_varying` argument

Takes a *one-sided* formula, e.g. `~ 1` or `~ 0 + depth`

Note: initial coefficient is unconstrained, i.e. **do not place the same covariate in 
the `formula` argument** (this includes the intercept)

---

# Time-varying intercept

Note: a `0` or `-1` in formula for suppressing global intercept

Otherwise, both the main effects and time-varying effects would have the same parameter and this can't be estimated.

.small[

```r
mesh <- make_mesh(pcod, xy_cols = c("X", "Y"), cutoff = 10)
fit <- sdmTMB(
  density ~ 0 + s(depth, k = 5), 
* time_varying = ~ 1,
  data = pcod, mesh = mesh,
  time = "year",  
  family = tweedie(link = "log")
)
```
]

---

# Getting coefficients

Return with

.small[

```r
print(fit)
```
]

.small[

```r
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: density ~ 0 + s(depth, k = 5)
#> Time column: "year"
#> ...
*#> Time-varying parameters:
*#>                  coef.est coef.se
*#> (Intercept)-2003     1.96    0.29
*#> (Intercept)-2004     2.31    0.27
*#> (Intercept)-2005     2.06    0.27
*#> ...
*#> (Intercept)-2015     2.07    0.27
*#> (Intercept)-2017     1.55    0.29
#> ...
```
]

---

# Getting coefficients

Or by digging into `fit$sd_report`

(Not yet in `tidy.sdmTMB()`.)

```r
library(TMB)
est <- as.list(fit$sd_report, "Est")
est_se <- as.list(fit$sd_report, "Std. Error")
cbind(est$b_rw_t, est_se$b_rw_t)
#>           [,1]      [,2]
#>  [1,] 1.958251 0.2889057
#>  [2,] 2.313664 0.2724793
#>  [3,] 2.064983 0.2707960
#>  [4,] 1.232808 0.3014059
#>  [5,] 1.511154 0.2758021
#>  [6,] 1.932893 0.2710203
#>  [7,] 2.159133 0.2672042
#>  [8,] 2.073629 0.2725238
#>  [9,] 1.549745 0.2933550
```

---

# Other approaches to modeling time-varying intercepts

.small[

```r
density ~ s(depth) + 0 + as.factor(year)
```
]
.small[

```r
density ~ s(depth) + (1 | year)
```
]
.small[

```r
density ~ s(depth) + s(year)
```
]
---

# These approaches are similar but subtly different

---

# Time-varying coefficients

Time-varying (random walk) effect of depth

Intercept in this model NOT time-varying

```r
fit_tv <- sdmTMB(
  density ~ 1, 
* time_varying = ~ 0 + depth_scaled + depth_scaled2,
  data = pcod, mesh = mesh,
  time = "year",
  family = tweedie(link = "log"),
  spatial = "on",
  spatiotemporal = "iid",
  silent = FALSE
)
```

---

# Time-varying coefficients

Time-varying (random walk) effect of depth

---

# Time-varying coefficient notes

* `time_varying` is a formula for coefficients that follow a random walk over time

* Make sure a coefficient isn't in both `formula` and `time_varying`, this includes the intercept

* The `time_varying` formula cannot have smoothers `s()` in it! Instead:
  * Polynomials: `time_varying = ~ x + I(x^2)`
  * `formula = s(depth, by = factor_year)` (independent smooths) 
  * `formula = s(depth, year)` (2D smooth)