Related packages, troubleshooting, tips

class: center, middle, inverse, title-slide

.title[
# Related packages, troubleshooting, tips
]
.subtitle[
## DFO TESA sdmTMB workshop
]
.author[
### 
]
.date[
### January 16–19 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

.center[
<img src="images/table_comparison.png" width="550px" height = "500px"/>
]

---

# 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))

---

class: center, middle, inverse

# Troubleshooting

---

# Example of non-converging model

.small[

```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:

.small[

```r
fit
#> Spatial model fit by ML ['sdmTMB']
#> Formula: present ~ depth * year
#> Mesh: mesh
#> Data: pcod
#> Family: binomial(link = 'logit')
#>  
#>             coef.est coef.se
#> (Intercept)    41.49   24.32
#> depth           0.02    0.00
*#> year           -0.02    0.01
#> depth:year      0.00     NaN
#> 
#> Matern range: 37.06
#> Spatial SD: 2.20
#> ML criterion at convergence: 1116.725
#> 
#> See ?tidy.sdmTMB to extract these values as a data frame.
```
]

---

# Sanity

* `NaN` checks alone can't diagnose all problems

---

# Inspecting output

We could also look directly at output from TMB `TMB::sdreport()`:

.small[

```r
fit$sd_report
#> sdreport(.) result
#>               Estimate   Std. Error
#> b_j       4.148809e+01 24.316054458
#> b_j       2.349202e-02  0.002301295
#> b_j      -1.860347e-02  0.012087876
*#> b_j      -2.357971e-05          NaN
#> ln_tau_O  5.202919e-01  0.153541429
#> ln_kappa -2.572838e+00  0.088994367
*#> Warning:
*#> Hessian of fixed effects was not positive definite.
*#> Maximum gradient component: 8.067217
```
]

---

# Looking at gradients

Log likelihood gradient with respect to fixed effects:

.small[

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

* 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)

---

# Interpretation

* 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"
)
```

---
 
# 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, plot the predictions, plot the random fields, ...]

* 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]

* Consider if observation error variance >> process variance   
  .xsmall[Known to be problematic for state-space models]