Evaluate or simulate the tinyAM model

Evaluates the population and observation equations for a parameter list. Most users should use fit_tam() to fit a model and sim_tam() to generate simulated data. This lower-level function is useful for checking model equations or simulating from chosen parameter values.

Usage

nll_fun(par, dat, simulate = FALSE)

Arguments

par: Parameter list with the structure produced by make_par().
dat: Data and settings returned by prepare_tam().
simulate: Logical; generate process states and observations instead of returning the likelihood? Defaults to FALSE.

Value

With simulate = FALSE, a scalar joint negative log-likelihood. When used inside RTMB::MakeADFun(), reported population and observation quantities are available through the resulting object's report() method; selected log population summaries also receive uncertainty via sdreport(). With TRUE, a list containing log_r, log_obs, log_f when its process is active, and applicable log_n0, log_n, and log_m values. Missing observations remain NA; predictions and population states are calculated for the full model grid.

Details

See tinyAM-model for the complete state equations, process distributions, likelihood, initial conditions, and derived quantities.

Latent-state convention: log_r, log_n0, log_n, log_f, and log_m represent absolute latent quantities on the log scale. Full model surfaces are log_N, log_F, and log_M. Process deviations are separate quantities: recruitment uses its formula's conditional predictions (eta_R), N uses cohort predictions (eta_log_N), and F/M use their mean log surfaces (eta_log_f, eta_log_m). In particular, log_f = log_mu_F + eta_log_f and log_m = log_mu_M + eta_log_m on their represented years and ages/blocks. The mean is not added again when constructing mortality from the latent state.

With simulate = FALSE, the result is the joint negative log-likelihood. It includes recruitment, active N/F/M processes, random initial abundance, and Gaussian densities for observed log observations. Random effects are integrated out by fit_tam(), not by this function itself. Normal densities are evaluated on log observations, so predictions on the natural scale are conditional medians rather than arithmetic means. Catchability uses the link chosen by index_settings$q_link: the log link exponentiates the formula predictor, while the logit link applies the inverse logit and restricts q to between zero and one. Both estimation and simulation use the same predictor, including any mono() increments.

With simulate = TRUE, mortality states are drawn first. Recruitment and N are constructed chronologically, using the same simulated parent SSB for stock-recruit relationships. This is followed by predictions and observation draws with the SD for each matching row. Derived quantities therefore use the same realization as the returned states. If biological inputs extend above the modeled plus age, hidden age abundances are reconstructed within each year before SSB is calculated. Effective terminal W and P preserve biomass and mature biomass; see tinyAM-model. No extra parameters or process penalties are introduced. N_plus, W, and P are then included in report(). Initial older-age abundance, log_r0, and additional early recruitment boundary states retain their supplied values during process simulation. See recruitment_formulas for fixed versus integrated initialization. RW processes retain their starting state because it has no process density. For N, the first log_n row is retained and its starting residual is computed against the newly simulated cohort prediction. Subsequent residuals follow the temporal walk. No new F states are drawn for projection years: projected F scales terminal historical F by a supplied multiplier or a multiplier solved from projected catch biomass. Both paths use the same chronological population. Active M states already extend into projected years under their process density.

Examples

dat <- prepare_tam(
  cod_obs,
  years = 1983:2024,
  ages = 2:14,
  N_settings = list(process = "iid", init = "exp"),
  F_settings = list(process = "rw", mu_form = NULL),
  M_settings = list(process = "off", mu_supplied = ~ I(0.3)),
  catch_settings = list(sd_form = ~ 1),
  index_settings = list(sd_form = ~ 1, q_form = ~ q_block)
)
par <- make_par(dat)
set.seed(1)
simulated <- nll_fun(par, dat, simulate = TRUE)
head(exp(simulated$log_r))
#>       1984       1985       1986       1987       1988       1989 
#> 0.28583061 0.10415661 0.22080651 0.05967620 0.10113673 0.05931925