Related packages, troubleshooting, tips

.title[
# Related packages, troubleshooting, tips
]
.subtitle[
## IMR sdmTMB workshop
]
.author[
### 
]
.date[
### May 23–25 2023
]

---

# Getting help

* sdmTMB GitHub repository:  
  <https://github.com/pbs-assess/sdmTMB>
  
* sdmTMB documentation:  
  <https://pbs-assess.github.io/sdmTMB/index.html>

*  New features to suggest? Bugs?  
  <https://github.com/pbs-assess/sdmTMB/issues>
  
*  Resources and reference material: 
  <https://github.com/pbs-assess/sdmTMB/wiki/resources>

---

# Related modelling software

---

# sdmTMB limitations

sdmTMB fits a variety of models with univariate responses

sdmTMB: 
* does not fit multivariate models ([see VAST](https://github.com/James-Thorson-NOAA/VAST))
* does not conduct spatial (dynamic) factor analysis ([see VAST](https://github.com/James-Thorson-NOAA/VAST))

---

# Troubleshooting

---

# Example of non-converging model

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

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

```r
#> Warning message:
#> The model may not have converged: 
#> non-positive-definite Hessian matrix. 
```
]
---

# Who is Hessian and why are they not "positive definite"?

* A Hessian matrix is a matrix of second derivatives of a function (here the log likelihood surface)

* The inverse of the negative Hessian is the parameter covariance matrix

* Warning means the curvature of the log-likelihood surface is inconsistent with the model having found the best fit

* Overparameterized? Variance estimated near zero?

* See [vignette("troubleshooting", "glmmTMB")](https://cran.r-project.org/web/packages/glmmTMB/vignettes/troubleshooting.html)

---

# Inspecting output

One or more standard errors being `NaN` is a problem:

```r
fit
#> Spatial model fit by ML ['sdmTMB']
#> Formula: present ~ depth * year
#> Mesh: mesh (isotropic covariance)
#> Data: pcod
#> Family: binomial(link = 'logit')
#>  
#>             coef.est coef.se
#> (Intercept)    41.66   24.32
#> depth           0.02    0.00
#> year           -0.02    0.01
*#> depth:year      0.00     NaN
#> 
#> Matérn range: 37.06
#> Spatial SD: 2.20
#> ML criterion at convergence: 1116.725
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
#> 
#> **Possible issues detected! Check output of sanity().**
```
]

---

# Inspecting output

We could also look directly at the TMB `sdreport`:

```r
fit$sd_report
#> sdreport(.) result
#>               Estimate   Std. Error
#> b_j       4.166376e+01 24.316746431
#> b_j       2.349096e-02  0.002301342
#> b_j      -1.869076e-02  0.012088203
*#> b_j      -2.357976e-05          NaN
#> ln_tau_O  5.202341e-01  0.153542815
#> ln_kappa -2.572805e+00  0.088995040
*#> Warning:
*#> Hessian of fixed effects was not positive definite.
*#> Maximum gradient component: 131.1179
```
]

---

# Looking at gradients

Log likelihood gradient with respect to fixed effects:

```r
max(fit$gradients)
#> [1] 0.0006444363
```
]

* Gradient becomes smaller as likelihood surface becomes flatter at the location in parameter space

* Flatter surface = closer to point of maximum likelihood

* So, small values are consistent with convergence (no hard rule, perhaps at least < 0.001 but depends on scale of log likelihood)

---

# Sanity function - an useful wrapper

* the `sanity()` function checks all of the above mentioned inspection points

```r
sanity(fit)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✖ Non-positive-definite Hessian matrix: model may not have converged
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✖ `b_j` standard error is NA
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
```
]

---

# Interpretation of the above model

* Linear interaction term seems tricky to estimate

* Solutions:
  * make sure the model is identifiable
  * drop linear interaction (maybe not important)
  * try re-scaling predictors
  * try a different mesh setup
  * simplify some other part of the model
  * try a related fixed-effect structure, e.g., non-linear smooths
  
---

# Smooth effects by year

For example, separate smooths by year:

```r
pcod$fyear <- as.factor(pcod$year)
fit <- sdmTMB(
* present ~ s(depth, by = fyear),
  data = pcod,
  mesh = mesh,
  family = binomial(link = "logit"),
  spatial = "on"
)
```

---

# Tips with mesh

Unfortunately, there is no "fit-for-all" mesh

It all depends on data (distance between survey points, etc) but check that the estimated `Matern range > mesh size (cutfoff)`

```r
#Spatial model fit by ML ['sdmTMB']
#[...]
*#Matern range: 37.06

# data 
quantile(dist(fit$data[,c("X","Y")]), c(0.01,0.05,0.1))
#>        1%        5%       10% 
#>  8.901301 21.907257 32.712026
# mesh
quantile(dist(fit$spde$mesh$loc[,1:2]), c(0.01,0.05,0.1))
#>       1%       5%      10% 
#> 14.07544 30.25935 42.70809
```

---

# Residuals diagnostics

---

# Randomized quantile residuals - default

```r
set.seed(1)
rq_res <- residuals(fit) 
qqnorm(rq_res)
qqline(rq_res)
```

<img src="04-troubleshooting_files/figure-html/resid1-1.png" width="500px" style="display: block; margin: auto;" />
]

.xsmall[
Warning: these residuals are fast but might look off even if the model is fine. Also see MCMC residuals. See the ['Residual checking' vignette](https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html).
]

---

# Model residuals in space

```r
pcod$resids <- residuals(fit)
filter(pcod, year %in% c(2011, 2013, 2015, 2017)) %>% 
* ggplot(aes(X, Y, colour = resids)) +
  geom_point() +
  facet_wrap(~year) +
  scale_colour_gradient2() +
  coord_fixed()
```

<img src="04-troubleshooting_files/figure-html/resid2-1.png" width="700px" style="display: block; margin: auto;" />
]

---

# Model residuals in space (continue)

```r
pred_fixed <- fit$family$linkinv(predict(fit)$est)
s_fit = DHARMa::createDHARMa(
  simulatedResponse = simulate(fit, nsim=500) ,
  observedResponse = pcod$present,
  fittedPredictedResponse = pred_fixed
)
DHARMa::testSpatialAutocorrelation(s_fit, x= pcod$X, y=pcod$Y, plot=FALSE)
#> 
#> 	DHARMa Moran's I test for distance-based autocorrelation
#> 
#> data:  s_fit
#> observed = -0.00186286, expected = -0.00046685, sd = 0.00228332, p-value = 0.5409
#> alternative hypothesis: Distance-based autocorrelation
```
]

---
  
# Other warning messages

* `Matrix` package, compiler warnings
  * generally OK, but see glmmTMB page:  
    <https://glmmtmb.github.io/glmmTMB/>

* TMB users forum  
  <https://groups.google.com/g/tmb-users/>

---

# Words of wisdom

**Start simple!**

These are complex models:

* May take a while to fit (`silent = FALSE`)

* Easy to create configurations that won't converge

* Can also be hard to wrap your head around

---

# Words of wisdom

Options for starting simple:

* Start intercept only

* Start with a spatial not spatiotemporal model

* Start without random fields:  
  `spatial = 'off', spatiotemporal = 'off'`
  
--

* Start with coarse meshes, maybe increase later  
  .xsmall[Mesh resolution greatly impacts speed and finer is not always better; can affect convergence.]

---

# More words of wisdom

* Make lots of plots  
  .xsmall[Plot the raw data, predictions, random fields, residuals, ...]

* Don't get too caught up in model configuration  
  .xsmall[Many are only subtly different and can produce similar predictions]

* Consider the scale of the response and predictors
  .xsmall[Keep parameters not too big or small; standardizing variables?]

* Consider if observation error variance (phi) >> process variance (sigma_O and/or sigma_E)  
  .xsmall[Known to be problematic for state-space models]
  
---

# And a few more

* Make sure the estimated spatial range > mesh size (cut-off value)

* large SE in the TMB `sdreport` might suggest overparametrization