Introduction

The Fisheries Integrated Management Systems aims to offer a system of modular, maintainable, and extensible models for next generation stock assessments. The choice of software platform used to develop FIMS will heavily influence the success of the initiative. While a thorough review of FIMS requirements is currently underway, a concurrent exploration of candidate software platforms will provide insight that will guide final decisions once software requirements have been established.

This project aims to produce a comprehensive review of available statistical computing platforms and summarize their suitability with respect to FIMS. The cross-comparison study will focus on state-space and explicit spatio-temporal modeling. Models will be developed and implemented across the following software platforms: TMB, C++, Julia, Greta tensorflow, and Stan. Inference pathways for both frequentist and Bayesian inference will be outlined for each platform when applicable. Under consideration will be performance, speed, the usability of the software with respect to both the user and developer, and software features that meet FIMS requirements under consideration.

Models

The state-space comparison looks at a logistic growth model:

\[\begin{align} \eta_{t} &= u_{t-1} + ru_{t-1}(1-u_{t-1}/K) \\ u_{t} &\sim LN(log(\eta_{t}),\sigma_{proc}) \\ y_{t} &\sim LN(log(u_{t}), \sigma_{obs}) \end{align}\],

where \(r\) represents the growth rate and \(K\) represents the density-dependent carrying capacity.

Software Platforms

TMB/R

Fast random effects, reliance on small community of developers

Base code written in C++, models implemented in R
CppAD used for forward/reverse autodifferentiation
Derivatives out to third order calculated for Laplace approximation
LA used to integrate out random effects
Sparse matrix algorithms used to automatically detect sparsity in Hessian and calculate sparse cholesky decompostion
Can fit models with large number of parameters and random effects
Memory bottleneck when translating from R to C++ (MakeADFun) and returning parameter/derived standard errors (sdreport)
?Modularity
Do not need to be expert C++ programmer to write TMB models, although contributors need more advanced skills to write efficient code or contribute to TMB base
Well tested AD capabilities
Small developers community

Julia Fast, Code written and executed in single language, still in infancy (may be bugs in base code, small community of developers)

New dynamic programming language with Just-In-Time (JIT) Compiler
Fast although first time to plot can be slow
Avoids the ‘two-language problem’ – all code is written and executed with Julia
Modular language – development occurs within packages
Multiple packages useful to numerical optimization: Distributions.jl, SPDEPrecisionMatrices.jl, GaussianMarkovRandomFields.jl, Optim.jl, SparsityDetection.jl
Autodifferentiation packages: ForwardDiff.jl, ReverseDiff.jl, Zygote.jl
Bayesian package: Turing.jl
Higher order differentiation to use Laplace approximation to integrate out parameters is still under development

C++/Rcpp Examples: r4MAS, CASAL, GADGET, MultifanCL

Stan HMC sampler, Laplace approximation using the adjoint method, variational Bayes

TMBstan
Writing code in Stan versus interacting with C++

Greta Machine learning approach to Bayesian inference

Relies on tensorflow to computationally automate Baye’s Rule in its application to generative DAGs
greta is a CRAN package that interfaces between R and tensorflow
tensorflow is a python package that implements machine learning

Project Deliverables

Table detailing whether or not each platform meets list of FIMS requirements
Speed test results
Outline inference pathways

Results

Platform Comparisons

State-Space Gompertz Results

Parameter estimates for state-space gompertz model (n=100)
	true	tmb	tmbstan	stanP0	stanP1
\(\theta_{1}\)	2.0	2.151	2.266	2.253	2.191
\(\theta_{2}\)	0.8	0.786	0.775	0.776	0.782
\(\sigma_{proc}\)	0.1	0.193	0.201	0.201	0.065
\(\sigma_{obs}\)	0.5	0.705	0.708	0.708	0.745

Standard error estimates for state-space gompertz model (n=100)
	tmb	tmbstan	stanP0	stanP1
\(\theta_{1}\)	2.151	0.012	0.370	0.315
\(\theta_{2}\)	0.786	0.001	0.038	0.032
\(\sigma_{proc}\)	0.079	0.030	0.368	1.194
\(\sigma_{obs}\)	0.062	0.003	0.090	0.092

Run times for state-space gompertz model (n=100)
	time..min.
tmb	0.001
tmbstan	0.972
stanP0	1.386
stanP1	2.283

State-Space logistic growth results

Parameter estimates for state-space logistic model (n=100)
	true	tmb	tmbstan	stanP0	stanP1
r	0.200	0.2124133	0.2093303	0.2072028	0.2081566
K	100.000	104.9323592	105.2542708	105.9495336	104.9503777
\(\sigma_{proc}\)	0.010	0.0817817	0.0917259	0.0905077	0.0976050
\(\sigma_{obs}\)	0.001	0.0406975	0.0128433	0.0262478	0.0031561

Standard error estimates for state-space logistic model (n=100)
	tmb	tmbstan	stanP0	stanP1
r	0.027	0.012	0.029	0.031
K	4.704	0.002	6.009	5.837
\(\sigma_{proc}\)	0.012	0.010	0.011	0.008
\(\sigma_{obs}\)	0.015	0.249	0.016	0.005

Run times for state-space logistic model (n=100)
	time..min.
tmb	0.002
tmbstan	0.501
stanP0	1.071
stanP1	0.393

Benchmark results for logistic growth model

## Warning: Removed 14 row(s) containing missing values (geom_path).

Inference Pathways

TMB

Frequentist: -> fn() and gr() passed to nlminb in R -> quasi-Newton optimization routine (Fox, Hall, and Schryer, 1978; Fox, 1997)
- fn(): inner optimization step finds optimum random effects values; Laplace approximation calculates the marginal likelihood after integrating out random effects
- gr(): Gradient functions are determined using forward/reverse automatic differentiation from cppAD and tinyAD
Bayesian -> fn() and gr() passed to Stan -> HMC algorithm (using tmbStan)

C++

Frequentist
Bayesian -> fn() and gr() passed to Stan -> HMC algorithm (using tmbStan)
- -> fn() and gr() passed to Stan -> HMC algorithm (using tmbStan)
- ? -> fn() passed to Stan -> Laplace/HMC ?
- ? -> pass to python and tensorflow ?
- Framework for acquiring gradients?

Performance results (eg. stability of source code) * Software platform pros and cons * Ease of use summary * Summary of input/output procedures

FIMS software platform exploratory analysis

Andrea Havron

7/20/2021