sel_comp: sel_comp

In nikolaifish/bamExtras: Functions and Data in Support of Stock Assessments with the Beaufort Assessment Model (BAM)

sel_comp

Description

Estimate selectivity from age or length compositions using optim()

Usage

sel_comp(
  par_init = c(a1 = 0.5, b1 = 3, a2 = 0.5, b2 = 6),
  par_lo = par_init * 0.001,
  par_up = par_init * 1000,
  fn_sel = "dbllgs",
  optim_method = "Nelder-Mead",
  optim_control = list(reltol = sqrt(.Machine$double.eps), maxit = 1000),
  age_pop,
  age_agec,
  P_agec = NULL,
  P_lenc = NULL,
  lenprob,
  F_full = 0.2,
  M_age,
  N0 = 1,
  penalty_out_of_bounds = 100,
  dbeta_args = list(shape1 = 1e-05, shape2 = 1e-05),
  value_type = "SSQ",
  output_type = "value",
  pass_dat_par_trace = TRUE,
  penalty_scale = 1,
  fit_to_zeros = FALSE,
  plot_format = list(col = list(init = "green2", obs = "purple", pr = "black"), lty =
    list(init = 2, obs = 0, pr = 1), lwd = list(init = 2, obs = 2, pr = 2), pch =
    list(init = NA, obs = 1, pr = NA), type = list(init = "l", obs = "p", pr = "l"))
)

Arguments

par_init	named vector of initial guesses for parameters to estimate
par_lo	lower bounds for par
par_up	upper bounds for par
fn_sel	name of selectivity function which can be called with two arguments, "x" and "par", where "x" is a vector of ages and "par" is a named numeric vector where the names correspond the names of the parameters in fn_sel.
optim_method	method argument to pass to optim
optim_control	control argument to pass to optim
age_pop	vector of ages in population (often includes 0)
age_agec	vector of ages in age comp data (often does not include 0)
P_agec	vector of proportion-at-age in age comp data
P_lenc	vector of proportion-at-age in length comp data
lenprob	age-length conversion matrix. Only used if P_lenc is specified
F_full	fixed value of full F. numeric vector of length 1.
M_age	vector of natural mortality at age, corresponding to age_agec
N0	initial population size
penalty_out_of_bounds	large value to add as penalty to objective function value (NLL or SSQ) if parameter estimates go outside of bounds
dbeta_args	list of arguments to pass to dbeta() function for defining penalty function. The shape parameters should be equal for symmetrical U-shaped penalty profile, which approach infinity as the parameter nears the bounds. Smaller values of the shape parameters make the U-shape increasingly flat-bottomed while larger values result in more curvature.
value_type	"SSQ" to return sum of squared deviations (default). "NLL" to return negative log-likelihood.
output_type	"value" (default) returns objective function value (NLL or SSQ). Use "value" with optim(). "predict" to return predicted age comp values.
pass_dat_par_trace	logical. Should dat_par_trace be passed to the global environment? Allows the user to follow the optimization of par.
penalty_scale	value to multiply by penalties added to objective function value. Defaults to 1 which leaves penalties as is. A value of zero effectively removes all penalties.
fit_to_zeros	logical. Should observed comps with values of zero be included in the objective function?
plot_format	list of plot formatting defaults

Details

As of 2023-04-20, this function seems to work pretty well based on testing with red snapper data (see example). However fitting to age comps works better than length comps. Overall, fitting seems to be sensitive to things like dbeta_args, penalty_out_of_bounds, and the variability of length at age in lenprob. So, at this point I'm going to use it with caution and continue to improve it.

Author(s)

Nikolai Klibansky

Examples

## Not run: 
rdat <- rdat_RedSnapper
agec <- tail(rdat$comp.mats$acomp.sCT.ob,1)
lenc <- tail(rdat$comp.mats$lcomp.cHL.ob,1)

# Estimate selectivity
sela1 <- sel_comp(age_pop=rdat$a.series$age,
                    age_agec=as.numeric(dimnames(agec)[[2]]),
                    P_agec=as.numeric(agec),
                    M_age=rdat$a.series$M
)
# Run the code again with the estimates from sela1 for par_init
sela2 <- sel_comp(par_init = sela1$fit$par,
                    age_pop=rdat$a.series$age,
                    age_agec=as.numeric(dimnames(agec)[[2]]),
                    P_agec=as.numeric(agec),
                    M_age=rdat$a.series$M
)

# Build length-at-age probability matrix
lenprob <- agelenprob(
  age=rdat$a.series$age,
  len=rdat$a.series$length,
  binmid_lenc=as.numeric(colnames(rdat$comp.mats$lcomp.cHL.ob)),
  cv_lenc=rdat$parm.cons$len.cv.val.L[8],
  plot=TRUE)

sell1 <- sel_comp(
  age_pop=rdat$a.series$age,
  age_agec=as.numeric(dimnames(agec)[[2]]),
  P_lenc=as.numeric(lenc),
  M_age=rdat$a.series$M,
  lenprob=lenprob
)

sell2 <- sel_comp(par_init = sell1$fit$par,
                    age_pop=rdat$a.series$age,
                    age_agec=as.numeric(dimnames(agec)[[2]]),
                    P_lenc=as.numeric(lenc),
                    M_age=rdat$a.series$M,
                    lenprob=lenprob
)

## End(Not run)

nikolaifish/bamExtras documentation built on July 21, 2023, 8:26 a.m.