After fitting multiple WHAM (or ASAP) models, compare_wham_models produces plots
and a table of AIC and Mohn's rho to aid model comparison.
Usage
compare_wham_models(
mods,
do.table = TRUE,
do.plot = TRUE,
fdir = getwd(),
compare.opts = NULL,
table.opts = NULL,
plot.opts = NULL,
fname = NULL,
sort = NULL,
calc.rho = NULL,
calc.aic = NULL,
do.print = NULL
)Arguments
- mods
(named) list of fit WHAM/ASAP models. To read in ASAP model output, use
read_asap3_fit. If no names are given, 'm1', 'm2', ... will be used.- do.table
T/F, produce table of AIC and/or Mohn's rho? Default = TRUE.
- do.plot
T/F, produce plots? Default = TRUE.
- fdir
character, path to directory to save table and/or plots. Default =
getwd().- compare.opts
list of options to generate comparison results:
$stockinteger, which stock to include in results. Default = 1.
$regioninteger, which region to include in results. Default = 1.
- table.opts
list of options for AIC/rho table:
$fnamecharacter, filename to save CSV results table (.csv will be appended). Default =
'model_comparison'.$sortT/F, sort by AIC? Default = TRUE.
$calc.rhoT/F, calculate Mohn's rho? Retrospective analysis must have been run for all modes. Default = TRUE.
$calc.aicT/F, calculate AIC? Default = TRUE.
$printT/F, print table to console? Default = TRUE.
$save.csvT/F, save table as a CSV file? Default = FALSE.
- plot.opts
list of options for plots:
$out.typecharacter, either
'pdf'or'png'(default ='png'because I am not suresystem('pdftk')will work across platforms.)$civector of T/F, length = 1 (applied to all models) or number of models
$yearsvector, which years to plot? Default = all (model and projection years).
$whichvector, which plots to make? Default = all. See details.
$relative.tocharacter, name of "base" model to plot differences relative to.
$alphascalar, (1-alpha)% confidence intervals will be plotted. Default = 0.05 for 95% CI.
$ages.labvector, overwrite model age labels.
$kobe.yrinteger vector (length = 1 (applied to all models or number of models, which year to use in Kobe plot (relative status). Default = terminal model year(s).
$M.ageinteger, which age to use in M time-series plot. Default =
max(data$which_F_age)(max age of F to use for full total F).$return.ggplotT/F, return a list of ggplot2 objects for later modification? Default = TRUE.
$kobe.probT/F, print probabilities for each model in each quadrant of Kobe plot? Default = TRUE.
$refpt"XSPR" or "MSY", which reference point to use. Default = "XSPR".
$browseOpen html document in web browser (if $out.type =
'png') (default = TRUE).
Value
a list with the following components:
daicVector of delta-AIC by model (if
do.table=Tandtable.opts$calc.aic=T)aicVector of AIC by model (if
do.table=Tandtable.opts$calc.aic=T)rhoMatrix of Mohn's rho by model (if
do.table=Tandtable.opts$calc.rho=T)bestName of best model (lowest AIC) (if
do.table=Tandtable.opts$calc.aic=T)tabResults table of AIC and Mohn's rho (if
do.table=T)gList of ggplot2 objects for later modification (if
do.plot=Tandplot.opts$return.ggplot=T)
Details
plot.opts$which specifies which plots to make:
- 1
3-panel of SSB (spawning stock biomass), F (fully-selected fishing mortality), and Recruitment
- 2
CV (coefficient of variation) for SSB, F, and Recruitment
- 3
Fleet selectivity (by block, averaged across years)
- 4
Index selectivity (by block, averaged across years)
- 5
Selectivity tile (fleets + indices, useful for time-varying random effects)
- 6
M time series (natural mortality, can specify which age with plot.opts$M.age)
- 7
M tile (useful for time-varying random effects)
- 8
3-panel of F X% SPR, SSB at F_X%SPR, and yield at F_X%SPR
- 9
2-panel of relative status (SSB / SSB at F_X%SPR and F / F_X%SPR)
- 10
Kobe status (relative SSB vs. relative F)
If plot.opts$return.ggplot = TRUE, a list g is returned holding the above ggplot2 objects for later modification.
g[[i]] holds the plot corresponding to i above, e.g. g[[2]] is the CV plot.
Examples
if (FALSE) { # \dontrun{
base <- read_asap3_fit()
m1 <- fit_wham(input1)
m2 <- fit_wham(input2)
mods <- list(base=base, m1=m1, m2=m2)
res <- compare_wham_models(mods)
} # }