Skip to contents

Comparison with the ADMB model

Source: vignettes/sbt_vs_admb.Rmd
sbt_vs_admb.Rmd

Introduction

This vignette compares outputs from the ADMB model to the TMB model compiled within the sbt package. This is done using a fixed parameter run and then by optimising. In the fixed parameter run the parameter values from a single cell of the ADMB model are used as fixed value inputs to sbt and outputs from sbt are generated.

Fixed parameter run

Comparison table

The table below includes transformed parameters (i.e., in natural space), derived quantities, priors, penalties, likelihoods, and objective function values for the ADMB model, the TMB model, and the difference and percent difference between the two models:

ADMB TMB Difference Percent difference
Parameters
B0 1.08358e+07 1.083578e+07 18.0862526 0.0001669
psi 1.50000e+00 1.500000e+00 0.0000000 0.0000000
sigmaR 6.00000e-01 6.000000e-01 0.0000000 0.0000000
h 5.50000e-01 5.500000e-01 0.0000000 0.0000000
q HSP 1.00000e+00 1.000000e+00 0.0000000 0.0000000
m0 4.00000e-01 4.000000e-01 0.0000000 0.0000000
m4 1.67050e-01 1.670507e-01 -0.0000007 0.0004073
m10 6.50000e-02 6.500000e-02 0.0000000 0.0000000
m30 4.57410e-01 4.574100e-01 0.0000000 0.0000039
Derived quantities
R0 8.74439e+06 8.744386e+06 3.9612256 0.0000453
alpha 1.09929e+07 1.099294e+07 -42.4487450 0.0003861
beta 2.78634e+06 2.786344e+06 -3.9206779 0.0001407
tau_ac2 6.44716e-01 6.447156e-01 0.0000004 0.0000590
Priors & penalties
kludge 0.00000e+00 0.000000e+00 0.0000000 0.0000000
sel.change 7.06247e+01 7.062467e+01 0.0000305 0.0000431
sel.smooth 3.42249e+01 3.422494e+01 -0.0000408 0.0001192
rec -3.03390e+01 -3.033900e+01 0.0000045 0.0000147
m0 0.00000e+00 0.000000e+00 0.0000000 0.0000000
m10 0.00000e+00 0.000000e+00 0.0000000 0.0000000
steep 0.00000e+00 0.000000e+00 0.0000000 0.0000000
omega 0.00000e+00 0.000000e+00 0.0000000 0.0000000
expl 0.00000e+00 0.000000e+00 0.0000000 0.0000000
sel.init 2.00000e-07 2.000000e-07 0.0000000 0.0002443
hstar 7.63360e+00 7.633597e+00 0.0000030 0.0000390
Likelihoods & objective function
LL1 1.94607e+02 1.946072e+02 -0.0002198 0.0001130
LL2 3.09202e+01 3.092031e+01 -0.0001054 0.0003409
LL3 3.48981e+01 3.489813e+01 -0.0000286 0.0000819
LL4 3.82627e+01 3.826276e+01 -0.0000591 0.0001546
Indo 9.30514e+01 9.305138e+01 0.0000202 0.0000217
Aus 4.54265e+01 4.542650e+01 -0.0000019 0.0000041
CPUE -6.46131e+01 -6.461295e+01 -0.0001494 0.0002312
Tags 1.76553e+02 1.765530e+02 0.0000081 0.0000046
Aerial 3.98139e+00 3.981441e+00 -0.0000505 0.0012676
Troll -1.20015e+01 -1.200157e+01 0.0000656 0.0005466
POP 1.75465e+03 1.754654e+03 -0.0037102 0.0002114
HSP 2.19884e+03 2.199098e+03 -0.2582337 0.0117441
GT 1.45904e+03 1.459035e+03 0.0046028 0.0003155
ObjF 6.03576e+03 6.036018e+03 -0.2577550 0.0042705

Comparison figures

The following figures compare outputs from the ADMB and TMB models:

Average weight (kg) at age.

Average weight (kg) at age.

Natural mortality (M) at age.

Natural mortality (M) at age.

Initial numbers at age in the population.

Initial numbers at age in the population.

Recruitment each year.

Recruitment each year.

Total number of inidividuals each year.

Total number of inidividuals each year.

Spawning biomass (tonnes).

Spawning biomass (tonnes).

Catch per unit effort (CPUE).

Catch per unit effort (CPUE).

Phi at age for a subset of years.

Phi at age for a subset of years.

Selectivity at age for a subset of years.

Selectivity at age for a subset of years.

LF observations for a subset of years.

LF observations for a subset of years. 

#> [1] 3.526193e-11
AF observations for a subset of years.

AF observations for a subset of years. 

#> [1] -6.342165e-07  6.605903e-07
AF predictions for a subset of years.

AF predictions for a subset of years.

Optimising

Now the model is optimised using nlminb to see if TMB produces the same result. First, make the initial values a little different:

Params <- list(par_log_B0 = 17,
               par_log_psi = log(data_par1$psi),
               par_log_m0 = log(data_par1$m0), 
               par_log_m4 = log(0.1),
               par_log_m10 = log(data_par1$m10), 
               par_log_m30 = log(0.5),
               par_log_h = log(data_par1$steep), 
               par_log_sigma_r = log(lr$sigma.r), 
               par_rdev_y = rnorm(n = data_par1$Reps, sd = 0.1),
               par_sels_init_i = data_par1$par_sels_init_i, 
               par_sels_change_i = data_par1$par_sels_change_i,
               par_log_cpue_q = -0.02,
               par_log_cpue_sigma = log(data_par1$sigma_cpue),
               par_log_cpue_omega = log(data_par1$cpue_omega),
               par_log_aerial_tau = log(data_par1$tau_aerial),
               par_log_aerial_sel = data_par1$ln_sel_aerial,
               par_log_troll_tau = log(0.4),
               par_log_hsp_q = data_par1$lnqhsp, 
               par_logit_hstar_i = qlogis(exp(data_par1$par_log_hstar_i)),
               par_log_tag_H_factor = log(data_par1$tag_H_factor)
)
Map <- list()
Map[["par_log_psi"]] <- factor(NA)
Map[["par_log_m0"]] <- factor(NA)
Map[["par_log_m10"]] <- factor(NA)
Map[["par_log_h"]] <- factor(NA)
Map[["par_log_sigma_r"]] <- factor(NA)
Map[["par_log_cpue_sigma"]] <- factor(NA)
Map[["par_log_cpue_omega"]] <- factor(NA)
Map[["par_log_aerial_tau"]] <- factor(NA)
Map[["par_log_aerial_sel"]] <- factor(rep(NA, 2))
Map[["par_log_hsp_q"]] <- factor(NA)
Map[["par_log_tag_H_factor"]] <- factor(NA)
obj <- MakeADFun(data = Data, parameters = Params, map = Map, random = c(), 
                 hessian = TRUE, inner.control = list(maxit = 1000), 
                 DLL = "sbt_v1")
bnd <- get_bounds(obj = obj)
opt <- nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, 
              lower = bnd$lb, upper = bnd$ub,
              control = list(eval.max = 1000,  # deaults to 200
                             iter.max = 1000)) # deaults to 150
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation

opt$message

# obj$par[1:10]
# obj$report()$spawning_biomass_y
# opt$par[1:10]
# obj$par <- opt$par

The nlminb function finds the same optimum as ADMB:

Spawning biomass (tonnes).

Spawning biomass (tonnes).