sd report gradients vs optim

[Samuel Johnson]

Hi friends,

I'm getting some strange gradients in my sd report as compared to the gradients getting returned by optim().

For example, when optimising, I get a final MGC of ~ 0.09 (you can argue that this is too high, but it's small enough for illustrative purposes).

After optimising, I call sdreport(obj) to calculate standard errors and check if the Hessian (covariance matrix) is PD.

The final gradients after calculating standard errors are sometimes much higher, for some parameters they are around -3, and for two around -20.

I am concerned that the delta method procedure used for calculating standard errors is ending in some strange place. But the Hessian comes out PD and all the standard errors look fine, so I'm not sure what to conclude, except there's something funky going on.

Since I had a covariance matrix to sample initial parameters, I also ran HMC via tmbstan. The resulting posteriors all have reasonable convergence metrics and the marginal densities all look good.

Sam

SessionInfo:

R version 4.1.2 (2021-11-01)
Platform: aarch64-apple-darwin21.1.0 (64-bit)
Running under: macOS Monterey 12.4

Matrix products: default
BLAS:   /opt/homebrew/Cellar/openblas/0.3.20/lib/libopenblasp-r0.3.20.dylib
LAPACK: /opt/homebrew/Cellar/r/4.1.2/lib/R/lib/libRlapack.dylib

locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base    

other attached packages:
 [1] knitr_1.36           corrplot_0.90        devtools_2.4.2      
 [4] usethis_2.1.3        MASS_7.3-54          beepr_1.3          
 [7] scales_1.1.1         wesanderson_0.3.6    Rcpp_1.0.8          
[10] RColorBrewer_1.1-2   here_1.0.1           tmbstan_1.0.2      
[13] rstan_2.21.2         StanHeaders_2.21.0-7 TMB_1.9.0          
[16] reshape2_1.4.4       forcats_0.5.1        stringr_1.4.0      
[19] dplyr_1.0.7          purrr_0.3.4          readr_2.0.2        
[22] tidyr_1.1.4          tibble_3.1.5         ggplot2_3.3.5      
[25] tidyverse_1.3.1      kableExtra_1.3.4     tictoc_1.0.1        
[28] bookdown_0.24        coda_0.19-4        

loaded via a namespace (and not attached):
 [1] matrixStats_0.61.0 fs_1.5.0           lubridate_1.8.0    webshot_0.5.2    
 [5] httr_1.4.2         rprojroot_2.0.2    tools_4.1.2        backports_1.3.0  
 [9] utf8_1.2.2         R6_2.5.1           DBI_1.1.1          colorspace_2.0-2  
[13] withr_2.4.2        tidyselect_1.1.1   gridExtra_2.3      prettyunits_1.1.1
[17] processx_3.5.2     curl_4.3.2         compiler_4.1.2     cli_3.1.0        
[21] rvest_1.0.2        xml2_1.3.2         desc_1.4.0         callr_3.7.0      
[25] systemfonts_1.0.3  digest_0.6.28      rmarkdown_2.11     svglite_2.0.0    
[29] pkgconfig_2.0.3    htmltools_0.5.2    sessioninfo_1.2.1  highr_0.9        
[33] dbplyr_2.1.1       fastmap_1.1.0      rlang_0.4.12       readxl_1.3.1      
[37] rstudioapi_0.13    farver_2.1.0       jquerylib_0.1.4    generics_0.1.1    
[41] jsonlite_1.7.2     inline_0.3.19      magrittr_2.0.1     loo_2.4.1        
[45] Matrix_1.3-4       munsell_0.5.0      fansi_0.5.0        lifecycle_1.0.1  
[49] yaml_2.2.1         stringi_1.7.5      pkgbuild_1.2.0     plyr_1.8.6        
[53] grid_4.1.2         crayon_1.4.2       lattice_0.20-45    haven_2.4.3      
[57] hms_1.1.1          ps_1.6.0           pillar_1.6.4       pkgload_1.2.3    
[61] codetools_0.2-18   stats4_4.1.2       reprex_2.0.1       glue_1.4.2        
[65] evaluate_0.14      V8_3.5.0           remotes_2.4.1      RcppParallel_5.1.4
[69] modelr_0.1.8       vctrs_0.3.8        tzdb_0.2.0         testthat_3.1.0    
[73] cellranger_1.1.0   gtable_0.3.0       assertthat_0.2.1   cachem_1.0.6      
[77] xfun_0.28          broom_0.7.10       viridisLite_0.4.0  audio_0.1-8      
[81] memoise_2.0.0      ellipsis_0.3.2   

[Cole Monnahan - NOAA Federal]

Hi Sam, I think the first thing I'd try is to see if you can use some Newton steps to get those gradients to 0 and see if the problem persists (e.g., in fit_tmb). Do those parameters with large gradients post uncertainty calcs have a much smaller scale relative to the SE? If so could rescaling it help? I have never looked at the gradients like that so don't have a sense of how well they behave in practice, but it does seem concerning if calculation of the Hessian involves something with those gradients. -Cole