SS_profile: Run a likelihood profile in Stock Synthesis.

In r4ss: R Code for Stock Synthesis

Run a likelihood profile in Stock Synthesis.

Description

Iteratively changes the control file using SS_changepars.

Usage

SS_profile(
  dir = "C:/myfiles/mymodels/myrun/",
  masterctlfile = "control.ss_new",
  newctlfile = "control_modified.ss",
  linenum = NULL,
  string = NULL,
  profilevec = NULL,
  usepar = FALSE,
  globalpar = FALSE,
  parfile = "ss.par",
  parlinenum = NULL,
  parstring = NULL,
  dircopy = TRUE,
  exe.delete = FALSE,
  model = "ss",
  extras = "-nox",
  systemcmd = FALSE,
  saveoutput = TRUE,
  overwrite = TRUE,
  whichruns = NULL,
  version = "3.30",
  prior_check = TRUE,
  read_like = TRUE,
  verbose = TRUE
)

Arguments

dir	Directory where input files and executable are located.
masterctlfile	Source control file. Default = "control.ss_new"
newctlfile	Destination for new control files (must match entry in starter file). Default = "control_modified.ss".
linenum	Line number of parameter to be changed. Can be used instead of string or left as NULL. Can be a vector if you are profiling multiple parameters at the same time.
string	String partially matching name of parameter to be changed. Can be used instead of linenum or left as NULL. Can be a vector if you are profiling multiple parameters at the same time.
profilevec	Vector of values to profile over. If you are profileing over multiple parameters at the same time this should be a data.frame or matrix with a column for each parameter.
usepar	Use PAR file from previous profile step for starting values?
globalpar	Use global par file ("parfile_original_backup.sso", which is automatically copied from original parfile) for all runs instead of the par file from each successive run
parfile	Name of par file to use (for 3.30 models, this needs to remain 'ss.par'). When globalpar=TRUE, the backup copy of this is used for all runs.
parlinenum	Line number in par file to change (if usepar = TRUE). Can be a vector if you are profiling multiple parameters at the same time.
parstring	String in par file preceding line number to change as an alternative to parlinenum (only needed if usepar = TRUE). Can be a vector if you are profiling multiple parameters at the same time.
dircopy	Copy directories for each run? NOT IMPLEMENTED YET.
exe.delete	Delete exe files in each directory? NOT IMPLEMENTED YET.
model	Name of executable. Default = "ss".
extras	Additional commands to use when running SS. Default = "-nox" will reduce the amount of command-line output.
systemcmd	Should R call SS using "system" function instead of "shell". This may be required when running R in Emacs. Default = FALSE.
saveoutput	Copy output .SSO files to unique names. Default = TRUE.
overwrite	Overwrite any existing .SSO files. Default = TRUE. If FALSE, then some runs may be skipped.
whichruns	Optional vector of run indices to do. This can be used to re-run a subset of the cases in situations where the function was interrupted or some runs fail to converge. Must be a subset of 1:n, where n is the length of profilevec.
version	SS version number. Currently "3.24" or "3.30" are supported, either as character or numeric values (noting that numeric 3.30 = 3.3). version = NULL is no longer the default or an allowed entry. The default is version = "3.30".
prior_check	Check to make sure the starter file is set to include the prior likelihood contribution in the total likelihood. Default = TRUE.
read_like	Read the table of likelihoods from each model as it finishes. Default = TRUE. Changing to FALSE should allow the function to play through even if something is wrong with reading the table.
verbose	Controls amount of info output to command line. Default = TRUE.

Note

The starting values used in this profile are not ideal and some models may not converge. Care should be taken in using an automated tool like this, and some models are likely to require rerunning with alternate starting values.

Also, someday this function will be improved to work directly with the plotting function SSplotProfile(), but they don't yet work well together. Thus, even if SS_profile() is used, the output should be read using SSgetoutput() or by multiple calls to SS_output() before sending to SSplotProfile().

Author(s)

Ian Taylor

See Also

SSplotProfile(), SSgetoutput(), SS_changepars(), SS_parlines()

Examples

## Not run: 
# note: don't run this in your main directory
# make a copy in case something goes wrong
mydir <- "C:/ss/Simple - Copy"

# the following commands related to starter.ss could be done by hand
# read starter file
starter <- SS_readstarter(file.path(mydir, "starter.ss"))
# change control file name in the starter file
starter[["ctlfile"]] <- "control_modified.ss"
# make sure the prior likelihood is calculated
# for non-estimated quantities
starter[["prior_like"]] <- 1
# write modified starter file
SS_writestarter(starter, dir = mydir, overwrite = TRUE)

# vector of values to profile over
h.vec <- seq(0.3, 0.9, .1)
Nprofile <- length(h.vec)

# run SS_profile command
profile <- SS_profile(
  dir = mydir, # directory
  # "NatM" is a subset of one of the
  # parameter labels in control.ss_new
  model = "ss",
  masterctlfile = "control.ss_new",
  newctlfile = "control_modified.ss",
  string = "steep",
  profilevec = h.vec
)


# read the output files (with names like Report1.sso, Report2.sso, etc.)
profilemodels <- SSgetoutput(dirvec = mydir, keyvec = 1:Nprofile)
# summarize output
profilesummary <- SSsummarize(profilemodels)

# OPTIONAL COMMANDS TO ADD MODEL WITH PROFILE PARAMETER ESTIMATED
MLEmodel <- SS_output("C:/ss/SSv3.24l_Dec5/Simple")
profilemodels[["MLE"]] <- MLEmodel
profilesummary <- SSsummarize(profilemodels)
# END OPTIONAL COMMANDS

# plot profile using summary created above
SSplotProfile(profilesummary, # summary object
  profile.string = "steep", # substring of profile parameter
  profile.label = "Stock-recruit steepness (h)"
) # axis label

# make timeseries plots comparing models in profile
SSplotComparisons(profilesummary, legendlabels = paste("h =", h.vec))


###########################################################################

# example two-dimensional profile
# (e.g. over 2 of the parameters in the low-fecundity stock-recruit function)
base_dir <- "c:/mymodel"

dir_profile_SR <- file.path(base_dir, "Profiles/Zfrac_and_Beta")

# make a grid of values in both dimensions Zfrac and Beta
# vector of values to profile over
Zfrac_vec <- seq(from = 0.2, to = 0.6, by = 0.1)
Beta_vec <- c(0.5, 0.75, 1.0, 1.5, 2.0)
par_table <- expand.grid(Zfrac = Zfrac_vec, Beta = Beta_vec)
nrow(par_table)
## [1] 25
head(par_table)
##   Zfrac Beta
## 1   0.2 0.50
## 2   0.3 0.50
## 3   0.4 0.50
## 4   0.5 0.50
## 5   0.6 0.50
## 6   0.2 0.75

# run SS_profile command
# requires modified version of SS_profile available via
# remotes::install_github("r4ss/r4ss@profile_issue_224")
profile <- SS_profile(
  dir = dir_profile_SR, # directory
  masterctlfile = "control.ss_new",
  newctlfile = "control_modified.ss",
  string = c("Zfrac", "Beta"),
  profilevec = par_table,
  extras = "-nohess"
)

# get model output
profilemodels <- SSgetoutput(
  dirvec = dir_profile_SR,
  keyvec = 1:nrow(par_table), getcovar = FALSE
)
n <- length(profilemodels)
profilesummary <- SSsummarize(profilemodels)

# add total likelihood (row 1) to table created above
par_table[["like"]] <- as.numeric(profilesummary[["likelihoods"]][1, 1:n])

# reshape data frame into a matrix for use with contour
like_matrix <- reshape2::acast(par_table, Zfrac ~ Beta, value.var = "like")

# make contour plot
contour(
  x = as.numeric(rownames(like_matrix)),
  y = as.numeric(colnames(like_matrix)),
  z = like_matrix
)

## End(Not run)

r4ss documentation built on May 28, 2022, 1:11 a.m.