pbs-assess / tmbpop

(Experimental) R package for implementing TMB age-structured POP model
1 stars 0 forks source link

Existing code into R package format #1

Open andrew-edwards opened 6 years ago

seananderson commented 6 years ago

Got a first draft of this running without changing any input/output formats yet.

# install.packages("remotes")
remotes::install_github("pbs-assess/tmb-pop")
library(tmbpop)

Documentation:

?sim
?build

See example model fitting:

?fit

Example with built-in pop_example data:

sim_dat <- sim(
  survey1 = pop_example$S1t,
  survey2 = pop_example$S2t,
  survey3 = pop_example$S3t,
  paa.catch.female = pop_example$patC1,
  paa.catch.male = pop_example$patC2,
  n.trips.paa.catch = pop_example$ntC,
  paa.survey1.female = pop_example$patS11,
  paa.survey1.male = pop_example$patS12,
  n.trips.paa.survey1 = pop_example$ntS1,
  catch = pop_example$Ct,
  paa.mature = pop_example$ma,
  weight.female = pop_example$wa1,
  weight.male = pop_example$wa2,
  misc.fixed.param = pop_example$MiscFixedParam,
  theta.ini = NULL, # if NULL, uses default values
  lkhd.paa = "binomial",
  var.paa.add = TRUE,
  enable.priors = TRUE
)
model <- fit(sim_dat)

Doesn't always converge that well yet.

MCMC with NUTS/Stan (that also doesn't converge very well yet):

library(tmbstan)
# options(mc.cores = parallel::detectCores()) # for parallel processing
pop_mcmc <- tmbstan(
  model$obj,
  chains = 1, # using only 1 chain and...
  iter = 600, # only 600 iterations for a quick example
  init = list("last.par.best"),
  control = list(adapt_delta = 0.9, max_treedepth = 20L) # as needed, see ?stan
)

pars <- c("logR0", "logM1", "logh", "logmuC", "deltaC", "logqS1") # a selection
bayesplot::mcmc_trace(as.array(pop_mcmc), pars = pars)
bayesplot::mcmc_hist(as.array(pop_mcmc), pars = pars)
bayesplot::mcmc_pairs(as.array(pop_mcmc), pars = pars)

Continuous integration testing with Travis: https://travis-ci.org/pbs-assess/tmbpop

Could use a better package name, function names, data format, and argument names still.