R/SSplotComps.R
In r4ss/r4ss: R Code for Stock Synthesis
Defines functions
SSplotComps
Documented in
SSplotComps
#' Plot composition data and fits.
#'
#' Plot composition data and fits from Stock Synthesis output. Multi-figure
#' plots depend on `make_multifig`.
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows, where subplots 21 to 24
#' (aggregated across years) are provided first, and subplots
#' 1 to 10 are all repeated for each fleet
#' \itemize{
#' \item 1 index data by fleet
#' \item 1 multi-panel composition plot
#' \item 2 single panel bubble plot for numbers at length or age
#' \item 3 multi-panel bubble plots for conditional age-at-length
#' \item 4 multi-panel plot of fit to conditional age-at-length for specific years
#' \item 5 Pearson residuals for A-L key
#' \item 6 multi-panel plot of point and line fit to conditional
#' age-at-length for specific length bins
#' \item 7 sample size plot
#' \item 8 TA1.8 Francis plot for marginal data with
#' Dirichlet-Multinomial and no Francis adjustment
#' \item 9 TA1.8 Francis weighting plot for marginal data
#' \item 10 TA1.8 Francis plot for conditional data with
#' Dirichlet-Multinomial and no Francis adjustment
#' \item 11 TA1.8 Francis weighting plot for conditional data
#' \item 12 Andre's mean age and std. dev. in conditional AAL
#' \item 21 composition by fleet aggregating across years
#' \item 22 composition by fleet aggregating across years within each season
#' \item 23 composition by fleet aggregating across seasons within a year
#' \item 24 bubble plot comparison of length or age residuals
# across fleets within partition
#' }
#' @param kind indicator of type of plot can be "LEN", "SIZE", "AGE", "cond",
#' "GSTAGE", "GSTLEN", "L@A", or "W@A".
#' @param sizemethod if kind = "SIZE" then this switch chooses which of the
#' generalized size bin methods will be plotted.
#' @param aalyear Years to plot multi-panel conditional age-at-length fits for
#' all length bins; must be in a "c(YYYY,YYYY)" format. Useful for checking the
#' fit of a dominant year class, critical time period, etc. Default=-1.
#' @param aalbin The length bin for which multi-panel plots of the fit to
#' conditional age-at-length data will be produced for all years. Useful to
#' see if growth curves are ok, or to see the information on year classes move
#' through the conditional data. Default=-1.
#' @template plot
#' @template print
#' @template fleets
#' @template fleetnames
#' @param sexes which sexes to show plots for. Default="all" which will include
#' males, females, and unsexed. This option is not fully implemented for all
#' plots.
#' @param yupper upper limit on ymax for polygon/histogram composition plots
#' @param datonly make plots of data without fits?
#' @param samplesizeplots make sample size plots?
#' @param compresidplots make plots of residuals for fit to composition data?
#' @param bub make bubble plot for numbers at age or size?
#' @param showyears Add labels for years to sample size plots?
#' @param showsampsize add sample sizes to plot
#' @param showeffN add effective sample sizes to plot
#' @param aggregates_by_mkt separate plots of aggregates across years
#' into different plots for each market category (retained, discarded)?
#' @param sampsizeline show line for input sample sizes on top of conditional
#' age-at-length plots (TRUE/FALSE, still in development)
#' @param effNline show line for effective sample sizes on top of conditional
#' age-at-length plots (TRUE/FALSE, still in development)
#' @param minnbubble number of unique x values before adding buffer. see
#' ?bubble3 for more info.
#' @param pntscalar This scalar defines the maximum bubble size for bubble
#' plots. This option is still available but a better choice is to use cexZ1
#' which allow the same scaling throughout all plots.
#' @param scalebubbles scale data-only bubbles by sample size, not just
#' proportion within sample? Default=FALSE.
#' @param cexZ1 Character expansion (cex) for point associated with value of 1.
#' @param bublegend Add legend with example bubble sizes to bubble plots.
#' @param colvec Vector of length 3 with colors for females, males, unsexed fish
#' @param linescol Color for lines on top of polygons
#' @param xlas label style (las) input for x-axis. Default 0 has horizontal
#' labels, input 2 would provide vertical labels.
#' @param ylas label style (las) input for y-axis. Default NULL has horizontal
#' labels when all labels have fewer than 6 characters and vertical otherwise.
#' Input 0 would force vertical labels, and 1 would force horizontal.
#' @param axis1 optional position of bottom axis values
#' @param axis2 optional position of left size axis values
#' @param axis1labs optional vector of labels for axis1 (either NULL or needs to
#' match length of axis1)
#' @param sizebinlabs Vector of size bin labels corresponding to the generalized
#' size frequency method
#' @param red What color to use for females in bubble plots (default is slightly
#' transparent red)
#' @param blue What color to use for males in bubble plots (default is slightly
#' transparent blue)
#' @template pwidth
#' @template pheight
#' @template punits
#' @template ptsize
#' @template res
#' @template plotdir
#' @template cex.main
#' @param linepos should lines be added before points (linepos=1) or after
#' (linepos=2)?
#' @param fitbar show fit to bars instead of points
#' @param do.sqrt scale bubbles based on sqrt of size vector. see ?bubble3 for
#' more info.
#' @param smooth add loess smoother to observed vs. expected index plots and
#' input vs. effective sample size?
#' @param cohortlines optional vector of birth years for cohorts for which to
#' add growth curves to numbers at length bubble plots
#' @template labels
#' @param printmkt show market categories in plot titles?
#' @param printsex show sex in plot titles?
#' @param maxrows maximum (or fixed) number or rows of panels in the plot
#' @param maxcols maximum (or fixed) number or columns of panels in the plot
#' @param maxrows2 maximum number of rows for conditional age at length plots
#' @param maxcols2 maximum number of columns for conditional age at length
#' plots
#' @param rows number or rows to return to as default for next plots to come or
#' for single plots
#' @param cols number or cols to return to as default for next plots to come or
#' for single plots
#' @param andre_oma Outer margins passed to Andre's multi-panel conditional
#' age-at-length plots.
#' @param andrerows Number of rows of Andre's conditional age-at-length plots
#' within each page. Default=3.
#' @param fixdims fix the dimensions at maxrows by maxcols or resize based on
#' number of years of data
#' @param fixdims2 fix the dimensions at maxrows by maxcols in aggregate plots
#' or resize based on number of fleets
#' @param maxneff the maximum value to include on plots of input and effective
#' sample size. Occasionally a calculation of effective N blows up to very
#' large numbers, rendering it impossible to observe the relationship for other
#' data. Default=5000.
#' @template verbose
#' @param scalebins Rescale expected and observed proportions by dividing by
#' bin width for models where bins have different widths? Caution!: May not
#' work correctly in all cases.
#' @param addMeans Add parameter means in addition to medians for MCMC
#' posterior distributions in which the median and mean differ.
#' @template mainTitle
#' @param \dots additional arguments that will be passed to
#' the `par` command in the [make_multifig()] function.
#' @author Ian Taylor
#' @export
#' @seealso [SS_plots()], [make_multifig()]
SSplotComps <-
function(
replist,
subplots = c(1:10, 21, 24),
kind = "LEN",
sizemethod = 1,
aalyear = -1,
aalbin = -1,
plot = TRUE,
print = FALSE,
fleets = "all",
fleetnames = "default",
sexes = "all",
yupper = 0.4,
datonly = FALSE,
samplesizeplots = TRUE,
compresidplots = TRUE,
bub = FALSE,
showyears = TRUE,
showsampsize = TRUE,
showeffN = TRUE,
aggregates_by_mkt = FALSE,
sampsizeline = FALSE,
effNline = FALSE,
minnbubble = 3,
pntscalar = NULL,
scalebubbles = FALSE,
cexZ1 = 1.5,
bublegend = TRUE,
colvec = c(rgb(1, 0, 0, .7), rgb(0, 0, 1, .7), rgb(.1, .1, .1, .7)),
linescol = c(rgb(0, .5, 0, .7), rgb(.8, 0, 0, .7), rgb(0, 0, .8, .7)),
xlas = 0,
ylas = NULL,
axis1 = NULL,
axis2 = NULL,
axis1labs = NULL,
sizebinlabs = NULL,
blue = rgb(0, 0, 1, 0.7),
red = rgb(1, 0, 0, 0.7),
pwidth = 6.5,
pheight = 6.5,
punits = "in",
ptsize = 10,
res = 300,
plotdir = "default",
cex.main = 1,
linepos = 1,
fitbar = FALSE,
do.sqrt = TRUE,
smooth = TRUE,
cohortlines = c(),
labels = c(
"Length (cm)", # 1
"Age (yr)", # 2
"Year", # 3
"Observed sample size", # 4
"Effective sample size", # 5
"Proportion", # 6
"cm", # 7
"Frequency", # 8
"Weight", # 9
"Length", # 10
"(t)", # 11
"(numbers x1000)", # 12
"Stdev (Age)", # 13
"Conditional AAL plot, ", # 14
"Size bin"
), # 15
printmkt = TRUE,
printsex = TRUE,
maxrows = 6,
maxcols = 4,
maxrows2 = 4,
maxcols2 = 4,
rows = 1,
cols = 1,
andre_oma = c(3, 0, 3, 0),
andrerows = 4,
fixdims = TRUE,
fixdims2 = FALSE,
maxneff = 5000,
verbose = TRUE,
scalebins = FALSE,
addMeans = TRUE,
mainTitle = FALSE,
...
) {
################################################################################
# SSplotComps
################################################################################
# table to store information on each plot
plotinfo <- NULL
SS_versionNumeric <- replist[["SS_versionNumeric"]]
lendbase <- replist[["lendbase"]]
sizedbase <- replist[["sizedbase"]]
agedbase <- replist[["agedbase"]]
condbase <- replist[["condbase"]]
ghostagedbase <- replist[["ghostagedbase"]]
ghostlendbase <- replist[["ghostlendbase"]]
ladbase <- replist[["ladbase"]]
wadbase <- replist[["wadbase"]]
tagdbase1 <- replist[["tagdbase1"]]
tagdbase2 <- replist[["tagdbase2"]]
nfleets <- replist[["nfleets"]]
nseasons <- replist[["nseasons"]]
seasfracs <- replist[["seasfracs"]]
FleetNames <- replist[["FleetNames"]]
nsexes <- replist[["nsexes"]]
accuage <- replist[["accuage"]]
Age_tuning <- replist[["Age_comp_Eff_N_tuning_check"]]
# define a variety of titles and labels
titles <- NULL
titlemkt <- ""
if (plotdir == "default") {
plotdir <- replist[["inputs"]][["dir"]]
}
# sort out which fleets will be included
if (fleets[1] == "all") {
fleets <- 1:nfleets
} else {
if (length(intersect(fleets, 1:nfleets)) != length(fleets)) {
stop(
"Input 'fleets' should be 'all' or a vector of values between 1 and nfleets."
)
}
}
if (fleetnames[1] == "default") {
fleetnames <- FleetNames
}
# sort out which sexes will be included, and associated labels
if (sexes[1] == "all") {
sexes <- 0:nsexes # this can be used to subset stuff below
}
if (nsexes == 1) {
sexes <- 0:nsexes
}
if (nsexes == 1 | length(sexes) > 1) {
# single-sex models, or models with all sexes shown
# on the same plot, don't need sex-specific title
titlesex <- ""
filesex <- ""
}
if (nsexes > 1 & length(sexes) == 1) {
# multi-sex models with only 1 sex shown need
# title and filename info
if (sexes == 0) {
titlesex <- "sexes combined, "
filesex <- "sex0"
}
if (sexes == 1) {
titlesex <- "female, "
filesex <- "sex1"
}
if (sexes == 2) {
titlesex <- "male, "
filesex <- "sex2"
}
}
titlesex <- ifelse(printsex, titlesex, "")
# a few quantities related to data type and plot number
tail_warning <- FALSE
if (kind == "LEN") {
dbase_kind <- lendbase
kindlab <- labels[1]
if (datonly) {
filenamestart <- "comp_lendat_"
titledata <- "Length comp data, "
} else {
filenamestart <- "comp_lenfit_"
titledata <- "Length comps, "
}
if (
!is.null(replist[["Length_comp_error_controls"]]) &&
any(replist[["Length_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "GSTLEN") {
dbase_kind <- ghostlendbase
kindlab <- labels[1]
if (datonly) {
filenamestart <- "comp_gstlendat_"
titledata <- "Excluded length comp data, "
} else {
filenamestart <- "comp_gstlenfit_"
titledata <- "Excluded length comps, "
}
if (
!is.null(replist[["Length_comp_error_controls"]]) &&
any(replist[["Length_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "SIZE") {
dbase_kind <- sizedbase[sizedbase[["method"]] == sizemethod, ]
# if user provided sizebinlabs, then add tick and label
# associated with each bin
if (!is.null(sizebinlabs)) {
kindlab <- labels[15]
axis1 <- sort(unique(dbase_kind[["Bin"]]))
# confirm that vector of labels is the correct length
if (length(sizebinlabs) == length(axis1)) {
axis1labs <- sizebinlabs
} else {
axis1labs <- axis1
warning(
"Input 'sizebinlabs' differs in length from the unique Bin\n",
" values associated with sizemethod = ",
sizemethod,
". Using bin values instead."
)
}
} else {
# evenly spaced units in length or weight rather than sizebins
sizeunits <- unique(dbase_kind[["units"]])
if (length(sizeunits) > 1) {
stop(
"!error with size units in generalized size comp plots:\n",
" more than one unit value per method.\n"
)
}
if (sizeunits %in% c("in", "cm")) {
kindlab <- paste(labels[10], " (", sizeunits, ")", sep = "")
}
if (sizeunits %in% c("lb", "kg")) {
kindlab <- paste(labels[9], " (", sizeunits, ")", sep = "")
}
}
if (datonly) {
filenamestart <- "comp_sizedat_"
titledata <- "Size comp data, "
} else {
filenamestart <- "comp_sizefit_"
titledata <- "Size comps, "
}
# add text noting which size method is represented
if (length(unique(sizedbase[["method"]])) > 1) {
filenamestart <- paste0(filenamestart, "method", sizemethod, "_")
titledata <- paste0(titledata, " size method ", sizemethod, ", ")
}
}
if (kind == "AGE") {
dbase_kind <- agedbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_agedat_"
titledata <- "Age comp data, "
} else {
filenamestart <- "comp_agefit_"
titledata <- "Age comps, "
}
if (
!is.null(replist[["Age_comp_error_controls"]]) &&
any(replist[["Age_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "cond") {
dbase_kind <- condbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_condAALdat_"
titledata <- "Conditional age-at-length data, "
} else {
filenamestart <- "comp_condAALfit_"
titledata <- "Conditional age-at-length, "
}
}
if (kind == "GSTAGE") {
dbase_kind <- ghostagedbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_gstagedat_"
titledata <- "Excluded age comp data, "
} else {
filenamestart <- "comp_gstagefit_"
titledata <- "Excluded age comps, "
}
if (
!is.null(replist[["Age_comp_error_controls"]]) &&
any(replist[["Age_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "GSTcond") {
dbase_kind <- ghostagedbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_gstCAALdat_"
titledata <- "Excluded conditional age-at-length data, "
} else {
filenamestart <- "comp_gstCAALfit_"
titledata <- "Excluded conditional age-at-length comps, "
}
}
if (kind == "L@A") {
dbase_kind <- ladbase[ladbase[["Nsamp_adj"]] != 0, ] # remove values with 0 sample size
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_LAAdat_"
titledata <- "Mean length at age data, "
} else {
filenamestart <- "comp_LAAfit_"
titledata <- "Mean length at age fit, "
}
dbase_kind[["SD"]] <- dbase_kind[["Lbin_lo"]] / dbase_kind[["Nsamp_adj"]]
}
if (kind == "W@A") {
dbase_kind <- wadbase[wadbase[["Nsamp_adj"]] != 0, ] # remove values with 0 sample size
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_WAAdat_"
titledata <- "Mean weight at age data, "
} else {
filenamestart <- "comp_WAAfit_"
titledata <- "Mean weight at age fit, "
}
dbase_kind[["SD"]] <- dbase_kind[["Lbin_lo"]] / dbase_kind[["Nsamp_adj"]]
}
if (
!(kind %in%
c("LEN", "SIZE", "AGE", "cond", "GSTAGE", "GSTLEN", "L@A", "W@A"))
) {
stop(
"Input 'kind' to SSplotComps needs to be one of the following:\n ",
"'LEN','SIZE','AGE','cond','GSTAGE','GSTLEN','L@A','W@A'."
)
}
if (datonly) {
titletype <- titledata
} else {
titletype <- "Pearson residuals, "
}
# partition group is used by some aggregate plots (subplot 21+)
# either equal to partition or constant across all samples
if (nrow(dbase_kind) > 0) {
if (aggregates_by_mkt) {
dbase_kind[["Part_group"]] <- dbase_kind[["Part"]]
} else {
dbase_kind[["Part_group"]] <- -1 # code for all partitions combined
}
}
# Add asterix to indicate super periods and then remove rows labeled "skip".
# It would be better to somehow show the range of years, but that seems difficult.
if (
any(dbase_kind[["SuprPer"]] == "Sup" & dbase_kind[["Used"]] == "skip")
) {
message(
"Removing super-period composition values labeled 'skip'\n",
"and designating super-period values with a '*'"
)
dbase_kind <- dbase_kind[
dbase_kind[["SuprPer"]] == "No" | dbase_kind[["Used"]] != "skip",
]
dbase_kind[["YrSeasName"]] <- paste(
dbase_kind[["YrSeasName"]],
ifelse(dbase_kind[["SuprPer"]] == "Sup", "*", ""),
sep = ""
)
}
ageerr_warning <- TRUE
# subset data based on requested range of fleets and sexes
dbase_kind <- dbase_kind[
dbase_kind[["Fleet"]] %in% fleets & dbase_kind[["sex"]] %in% sexes,
]
### subplot 21: by fleet aggregating across years
if (21 %in% subplots & kind != "cond") {
# for age or length comps, but not conditional AAL
# check for the presence of data
if (nrow(dbase_kind) > 0) {
# loop over partitions (discard, retain, total)
for (j in unique(dbase_kind[["Part_group"]])) {
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
dbase <- dbase_kind[dbase_kind[["Part_group"]] == j, ]
if (nrow(dbase) > 0) {
# market category
if (j == -1) {
titlemkt <- ""
}
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
## aggregating identifiers for plot titles and filenames
## note: titlesex is set at the top of this function
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste(filesex)
if (j > -1) {
# add market category to filename if it's not a mix
filename_fltsexmkt <- paste0(filename_fltsexmkt, "mkt", j)
}
caption1 <- paste(
titledata,
title_sexmkt,
"aggregated across time by fleet",
sep = ""
) # total title
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
Bins <- sort(unique(dbase[["Bin"]]))
nbins <- length(Bins)
df <- data.frame(
Nsamp_adj = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
obs = dbase[["Obs"]] * dbase[["Nsamp_adj"]],
exp = dbase[["Exp"]] * dbase[["Nsamp_adj"]]
)
if (
"Nsamp_DM" %in% names(dbase) && any(!is.na(dbase[["Nsamp_DM"]]))
) {
df[["Nsamp_DM"]] <- dbase[["Nsamp_DM"]]
}
agg <- aggregate(
x = df,
by = list(
bin = dbase[["Bin"]],
f = dbase[["Fleet"]],
sex = dbase[["sex"]],
mkt = dbase[["Part"]]
),
FUN = sum
)
agg <- agg[agg[["f"]] %in% fleets, ]
agg[["obs"]] <- agg[["obs"]] / agg[["Nsamp_adj"]]
agg[["exp"]] <- agg[["exp"]] / agg[["Nsamp_adj"]]
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
# loop over fleets
for (f in unique(agg[["f"]])) {
# loop over fleets within market
for (mkt in unique(agg[["mkt"]][agg[["f"]] == f])) {
sub <- agg[["f"]] == f & agg[["mkt"]] == mkt
agg[["Nsamp_adj"]][sub] <- max(agg[["Nsamp_adj"]][sub])
if (
"Nsamp_DM" %in% names(agg) && any(!is.na(agg[["Nsamp_DM"]]))
) {
agg[["Nsamp_DM"]][sub] <- max(
agg[["Nsamp_DM"]][sub],
na.rm = TRUE
)
} else {
if (any(!is.na(agg[["effN"]][sub]))) {
agg[["effN"]][sub] <- max(agg[["effN"]][sub], na.rm = TRUE)
} else {
agg[["effN"]][sub] <- NA
}
}
}
}
namesvec <- fleetnames[agg[["f"]]]
# check for multiple market categories in a fleet to plot separately
max_n_mkt <- max(apply(table(agg[["f"]], agg[["mkt"]]) > 0, 1, sum))
if (max_n_mkt > 0) {
mktnames <- c("", "(discards)", "(retained)")
namesvec <- paste(
fleetnames[agg[["f"]]],
mktnames[agg[["mkt"]] + 1]
)
}
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
# group remaining calculations as a function
tempfun7 <- function(ipage, ...) {
# test for Dirichlet-Multinomial likelihood
if (
"Nsamp_DM" %in% names(agg) && any(!is.na(agg[["Nsamp_DM"]]))
) {
# Dirichlet-Multinomial likelihood
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["f"]] + agg[["mkt"]] / 10,
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["Nsamp_DM"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "Sum of N input=",
effN_label = "Sum of N adj.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(namesvec, "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
} else {
# standard multinomial likelihood
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["f"]] + agg[["mkt"]] / 10,
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "Sum of N adj.=",
effN_label = "Sum of N eff.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(namesvec, "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
}
}
if (plot) {
tempfun7(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(agg[["f"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
# add warning about tail compression
if (ipage == 1 & tail_warning) {
caption <- paste0(
caption,
".\n <br> ",
"Note: tail compression in use makes aggregated plots more difficult to interpret."
)
}
# add warning about partitions
if (max_n_mkt > 0 & ipage == 1) {
caption <-
paste0(
caption,
".\n <br> ",
"Labels 'retained' and 'discard' indicate",
" discarded or retained sampled for each fleet.",
" Panels without this designation represent",
" the whole catch.\n"
)
}
# remove double period if present
caption <- gsub("..", ".", caption, fixed = TRUE)
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
"<br> (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
"_aggregated_across_time.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun7(ipage = ipage, ...)
dev.off()
}
} # end print function
} else {
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=dbase[["effN"]],showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase[["YrSeasName"]],"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase[["YrSeasName"]]),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
} # end test for presence of observations in this partition group
} # end loop over partitions group
# } # end loop over combined/not-combined sex
} # end if data
} # end subplot 21
### subplot 22: by fleet aggregating across years within each season
if (22 %in% subplots & kind != "cond" & nseasons > 1) {
# for age or length comps, but not conditional AAL
dbasef <- dbase_kind[dbase_kind[["Fleet"]] %in% fleets, ]
if ("Nsamp_DM" %in% names(dbasef) && any(!is.na(dbasef[["Nsamp_DM"]]))) {
warning(
"Sample sizes in plots by fleet aggregating across years within each season have not yet been updated to reflect Dirichlet-Multinomial likelihood"
)
}
# check for the presence of data
if (nrow(dbasef) > 0) {
testor <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0
]) >
0
testor[2] <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3)
]) >
0
testor[3] <- length(dbasef[["sex"]][dbasef[["sex"]] == 2]) > 0
# loop over sex combinations
for (k in (1:3)[testor]) {
if (k == 1) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0,
]
}
if (k == 2) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3),
]
}
if (k == 3) {
dbase_k <- dbasef[dbasef[["sex"]] == 2, ]
}
sex <- ifelse(k == 3, 2, 1)
# loop over partitions (discard, retain, total)
for (j in unique(dbase_k[["Part"]])) {
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
dbase <- dbase_k[dbase_k[["Part"]] == j, ]
if (nrow(dbase) > 0) {
## assemble pieces of plot title
# sex
if (k == 1) {
titlesex <- "sexes combined, "
}
if (k == 2) {
titlesex <- "female, "
}
if (k == 3) {
titlesex <- "male, "
}
titlesex <- ifelse(printsex, titlesex, "")
# market category
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste("sex", k, "mkt", j, sep = "")
caption1 <- paste0(
titledata,
title_sexmkt,
"\naggregated within season by fleet"
) # total title
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
Bins <- sort(unique(dbase[["Bin"]]))
nbins <- length(Bins)
df <- data.frame(
Nsamp_adj = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
obs = dbase[["Obs"]] * dbase[["Nsamp_adj"]],
exp = dbase[["Exp"]] * dbase[["Nsamp_adj"]]
)
agg <- aggregate(
x = df,
by = list(
bin = dbase[["Bin"]],
f = dbase[["Fleet"]],
s = dbase[["Seas"]]
),
FUN = sum
)
agg <- agg[agg[["f"]] %in% fleets, ]
if (any(agg[["s"]] <= 0)) {
message(
"super-periods may not work correctly in plots of aggregated comps"
)
agg <- agg[agg[["s"]] > 0, ]
}
agg[["obs"]] <- agg[["obs"]] / agg[["Nsamp_adj"]]
agg[["exp"]] <- agg[["exp"]] / agg[["Nsamp_adj"]]
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for (f in unique(agg[["f"]])) {
# loop over fleets
for (s in unique(agg[["s"]][agg[["f"]] == f])) {
# loop over seasons within fleet
infleetseas <- agg[["f"]] == f & agg[["s"]] == s
agg[["Nsamp_adj"]][infleetseas] <- max(agg[["Nsamp_adj"]][
infleetseas
])
agg[["effN"]][infleetseas] <- max(agg[["effN"]][infleetseas])
}
}
agg[["fseas"]] <- agg[["f"]] + seasfracs[agg[["s"]]]
namesvec <- paste(
fleetnames[agg[["f"]]],
" s",
agg[["s"]],
sep = ""
)
# group remaining calculations as a function
tempfun8 <- function(ipage, ...) {
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["fseas"]],
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(namesvec, "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=dbase[["effN"]],showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase[["YrSeasName"]],"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase[["YrSeasName"]]),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
if (plot) {
tempfun8(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(agg[["fseas"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
caption <- caption1
pagetext <- ""
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
"_aggregated_within_season.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun8(ipage = ipage, ...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined sex
} # end if data
} # end subplot 22
### subplot 23: by fleet aggregating across seasons within a year
if (23 %in% subplots & kind != "cond" & nseasons > 1) {
# for age or length comps, but not conditional AAL
# loop over fleets
for (f in fleets) {
dbasef <- dbase_kind[dbase_kind[["Fleet"]] == f, ]
if (
"Nsamp_DM" %in% names(dbasef) && any(!is.na(dbasef[["Nsamp_DM"]]))
) {
warning(
"Sample sizes in plots by fleet aggregating across seasons within a year have not yet been updated to reflect Dirichlet-Multinomial likelihood"
)
}
# check for the presence of data
if (nrow(dbasef) > 0) {
testor <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0
]) >
0
testor[2] <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3)
]) >
0
testor[3] <- length(dbasef[["sex"]][dbasef[["sex"]] == 2]) > 0
# loop over sex combinations
for (k in (1:3)[testor]) {
if (k == 1) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0,
]
}
if (k == 2) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3),
]
}
if (k == 3) {
dbase_k <- dbasef[dbasef[["sex"]] == 2, ]
}
sex <- ifelse(k == 3, 2, 1)
# loop over partitions (discard, retain, total)
for (j in unique(dbase_k[["Part"]])) {
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
dbase <- dbase_k[dbase_k[["Part"]] == j, ]
if (nrow(dbase) > 0) {
## assemble pieces of plot title
# sex
if (k == 1) {
titlesex <- "sexes combined, "
}
if (k == 2) {
titlesex <- "female, "
}
if (k == 3) {
titlesex <- "male, "
}
titlesex <- ifelse(printsex, titlesex, "")
# market category
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste(
"flt",
f,
"sex",
k,
"mkt",
j,
sep = ""
)
Bins <- sort(unique(dbase[["Bin"]]))
nbins <- length(Bins)
df <- data.frame(
Nsamp_adj = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
obs = dbase[["Obs"]] * dbase[["Nsamp_adj"]],
exp = dbase[["Exp"]] * dbase[["Nsamp_adj"]]
)
agg <- aggregate(
x = df,
by = list(
bin = dbase[["Bin"]],
f = dbase[["Fleet"]],
y = floor(dbase[["Yr.S"]])
),
FUN = sum
)
agg <- agg[agg[["f"]] %in% fleets, ]
agg[["obs"]] <- agg[["obs"]] / agg[["Nsamp_adj"]]
agg[["exp"]] <- agg[["exp"]] / agg[["Nsamp_adj"]]
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for (f in unique(agg[["f"]])) {
for (y in unique(agg[["y"]][agg[["f"]] == f])) {
infleetyr <- agg[["f"]] == f & agg[["y"]] == y
agg[["Nsamp_adj"]][infleetyr] <- max(agg[["Nsamp_adj"]][
infleetyr
])
agg[["effN"]][infleetyr] <- max(agg[["effN"]][infleetyr])
}
}
agg[["fy"]] <- agg[["f"]] + agg[["y"]] / 10000
# total title
caption <- paste(
titledata,
title_sexmkt,
fleetnames[f],
"\naggregated across seasons within year",
sep = ""
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
# group remaining calculations as a function
tempfun9 <- function(ipage, ...) {
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["fy"]],
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(agg[["y"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=dbase[["effN"]],showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase[["YrSeasName"]],"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase[["YrSeasName"]]),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
} # end tempfun
if (plot) {
tempfun9(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(agg[["fy"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
"_aggregated_across_seasons_within_year.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun9(ipage = ipage, ...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined sex
} # end if data
} # end loop over fleets
} # end subplot 23
### subplot 24: bubble plot comparison of length or age residuals
### across fleets within sex/partition
if (24 %in% subplots & kind %in% c("LEN", "AGE")) {
# loop over partitions groups (everything, or separate discard, retain, total)
for (j in unique(dbase_kind[["Part_group"]])) {
# subset data.frame for this partition group and subset of fleets of interest
dbase_parts <- dbase_kind[dbase_kind[["Part_group"]] == j, ]
# new column combining fleet and partition
# where 3.1 is fleet=3, partition=1
dbase_parts[["FleetPart"]] <- dbase_parts[["Fleet"]] +
0.1 * dbase_parts[["Part"]]
# table of info on each panel
panel_table <- data.frame(
FleetPart = sort(unique(dbase_parts[["FleetPart"]]))
)
# separate the pieces out again
panel_table[["Fleet"]] <- floor(panel_table[["FleetPart"]])
# round below is necessary because values were coming out as
# 1.9999999999999996 instead of 2
panel_table[["Part"]] <- round(
10 * (panel_table[["FleetPart"]] - panel_table[["Fleet"]])
)
# fleet name to use for each panel
panel_table[["Name"]] <- fleetnames[panel_table[["Fleet"]]]
# check for multiple market categories in a fleet to plot separately
max_n_mkt <- max(apply(
table(
panel_table[["Fleet"]],
panel_table[["Part"]]
) >
0,
1,
sum
))
if (max_n_mkt > 1) {
# if multiple categories within a fleet, append to name
mktnames <- c("", "(discards)", "(retained)")
panel_table[["Name"]] <- paste(
panel_table[["Name"]],
mktnames[panel_table[["Part"]] + 1]
)
}
npanels <- nrow(panel_table)
panelvec <- 1:npanels
xlim <- range(dbase_parts[["Yr.S"]]) # set xlim based on range across all fleets
xaxislab <- sort(unique(floor(dbase_parts[["Yr.S"]]))) # label with all years
# get growth curves if requested
if (length(cohortlines) > 0) {
growdat <- replist[["endgrowth"]]
growdatF <- growdat[
growdat[["Sex"]] == 1 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 1]),
]
if (nsexes > 1) {
growdatM <- growdat[
growdat[["Sex"]] == 2 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 2]),
]
}
}
# market category
if (j == -1) {
titlemkt <- ""
}
if (j == 0) {
titlemkt <- "whole catch"
}
if (j == 1) {
titlemkt <- "discard"
}
if (j == 2) {
titlemkt <- "retained"
}
titlemkt <- ifelse(printmkt, titlemkt, "")
## title_sexmkt <- paste(titlesex,titlemkt,sep="")
caption_base <- paste0(titletype, titlemkt, ", comparing across fleets")
# hack to remove ", ," from caption
caption_base <- gsub(", ,", ", ", caption_base)
if (mainTitle) {
ptitle <- caption_base
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
# if not all market categories, append mkt value to filename
filenamemkt <- ifelse(j > -1, paste("mkt", j, sep = ""), "")
multifleet.bubble.fun <- function(ipage = 0) {
# a function to wrap up multi-fleet bubble plots
# old graphics parameter settings
par_old <- par()
# multi-figure plot with as many rows as fleets, or the maxrows value
par(
mfrow = c(min(npanels, maxrows), 1),
mar = c(0.5, 0, 0, 0),
oma = c(4, 6, ifelse(mainTitle, 3, 1), 1)
)
# set up some stuff for cases where there are more fleets than panels in one plot
panelrange <- 1:npanels
npages <- ceiling(npanels / maxrows) # how many pages of plots
if (npages > 1 & ipage != 0) {
# range of which panels to print for each page
panelrange <- intersect(
panelrange,
1:maxrows + maxrows * (ipage - 1)
)
}
# loop over panels (fleet x partition)
for (ipanel in panelvec[panelrange]) {
flt <- panel_table[["Fleet"]][ipanel] # fleet number
mkt <- panel_table[["Part"]][ipanel] # market category
# subset database of values
dbase <- dbase_parts[
dbase_parts[["Fleet"]] == flt &
dbase_parts[["Part"]] == mkt,
]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), sex, and partition,
### not sure if multiple ageing error methods is supported at the moment,
### haven't tested -Ian 6/7/17
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(
table(dbase[["Yr.S"]], dbase[["Ageerr"]]) > 0,
1,
sum
))
if (max_n_ageerr > 1) {
if (ageerr_warning) {
message(
"Multiple samples with different ageing error types within fleet/year.\n",
"Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
"Bubble plots may be misleading with overlapping bubbles.\n"
)
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish
# between types within a year (may not work well for this plot)
dbase[["Yr.S"]] <- dbase[["Yr.S"]] +
dbase[["Ageerr"]] / (1000 * max_n_ageerr)
dbase[["YrSeasName"]] <- paste(
dbase[["YrSeasName"]],
"a",
dbase[["Ageerr"]],
sep = ""
)
}
# calculate smallest difference among years
# which is used to adjust offsets males and females
xdiff <- 0.1 *
sort(unique(diff(sort(unique(dbase[["Yr.S"]])), na.rm = TRUE)))[1]
# not sure what cases would have missing xdiff
# from above calculation, but it definitely happens
# with only one year of data, so setting default
# to work for that case
if (is.na(xdiff)) {
xdiff <- 0.1
}
# define colors
xvals <- dbase[["Yr.S"]]
cols <- rep(colvec[3], nrow(dbase))
if (nsexes > 1) {
xvals[dbase[["sex"]] > 0] <- dbase[["Yr.S"]][dbase[["sex"]] > 0] -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
if (length(unique(dbase[["Yr.S"]])) == 1) {
# if only one year, don't bother showing points
# as mid-year values
# this may not be ideal for seasonal models
xvals[dbase[["sex"]] > 0] <- floor(dbase[["Yr.S"]][
dbase[["sex"]] > 0
]) -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
}
cols[dbase[["sex"]] > 0] <- colvec[dbase[["sex"]][
dbase[["sex"]] > 0
]]
}
# determine bubble size and colors
if (datonly) {
z <- dbase[["Obs"]]
if (scalebubbles) {
z <- dbase[["Nsamp_adj"]] * dbase[["Obs"]]
} # if requested, scale by sample sizes
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
} else {
z <- dbase[["Pearson"]]
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
# make bubbles for a single fleet
# this section is a modified version of tempfun2 above
ylim <- range(dbase[["Bin"]])
ylim[2] <- ylim[2] + 0.2 * diff(ylim) # add buffer of 10% at the top for fleet name
bubble3(
x = xvals,
y = dbase[["Bin"]],
z = z,
col = cols,
cexZ1 = cexZ1,
legend = bublegend,
las = 1,
main = "",
cex.main = cex.main,
maxsize = pntscalar,
allopen = allopen,
xlim = xlim,
ylim = ylim,
axis1 = FALSE
)
### add label at top left of each panel
legend(
"topleft",
title = panel_table[["Name"]][ipanel],
legend = NA,
bty = "n",
cex = 1.5
)
### alternative way with legends on the side
### (probably not as good once the partition has been added)
# mtext(namesvec[ipanel],side=2,line=4.5,cex=par()[["cex"]])
# add lines for growth of individual cohorts if requested
if (length(cohortlines) > 0) {
for (icohort in seq_along(cohortlines)) {
message(" Adding line for", cohortlines[icohort], "cohort\n")
if (kind == "LEN") {
lines(
growdatF[["Age"]] + cohortlines[icohort],
growdatF[["Len_Mid"]],
col = colvec[1]
) # females
if (nsexes > 1) {
lines(
growdatM[["Age"]] + cohortlines[icohort],
growdatM[["Len_Mid"]],
col = colvec[2]
) # males
}
}
if (kind == "AGE") {
lines(
0.5 +
c(cohortlines[icohort], cohortlines[icohort] + accuage),
c(0, accuage),
col = "red"
)
}
}
}
if (
par()[["mfg"]][1] == par()[["mfg"]][3] |
ipanel == tail(panelvec, 1)
) {
# label all years on x-axis of last panel
axis(1, at = xaxislab)
} else {
# or just tick marks for other panels
axis(1, at = xaxislab, labels = rep("", length(xaxislab)))
}
if (par()[["mfg"]][1] == 1) {
# add title after making first panel
title(
main = ptitle,
outer = TRUE,
xlab = labels[3],
ylab = kindlab
)
}
} # end loop over fleets
# restore previous graphics parameter settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
} # end function wrapping up a single page of the residual comparison plot
# make plots or write to PNG file
if (length(panelvec) > 0) {
if (plot) {
multifleet.bubble.fun(ipage = 0)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(nrow(panel_table) / maxrows)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption_base
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste0(
caption,
" (plot ",
ipage,
" of ",
npages,
")"
)
}
if (ipage == 1 & length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
#### current scaling allows comparison across panels, so warning below
#### has been turned off
## caption <- paste(caption,
## "<br>Note: bubble sizes are scaled to maximum within each panel.",
## "<br>Thus, comparisons across panels should focus on patterns, not bubble sizes.")
file <- paste(
filenamestart,
filenamemkt,
pagetext,
"_multi-fleet_comparison.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
multifleet.bubble.fun(ipage = ipage)
dev.off()
} # end loop over pages within printing PNG
} # end printing to PNG files
} # end test for non-zero number of fleets
} # end loop over partitions
## } # end loop over sexes
## } # end loop over sex combinations
# restore default single panel settings
par(mfcol = c(rows, cols), mar = c(5, 4, 4, 2) + .1, oma = rep(0, 4))
} # end subplot 24
# loop over fleets
for (f in fleets) {
# check for the presence of data
if (length(dbase_kind[["Obs"]][dbase_kind[["Fleet"]] == f]) > 0) {
dbasef <- dbase_kind[dbase_kind[["Fleet"]] == f, ]
# get mean sample quantities to show in conditional age-at-length figures
if (kind %in% c("cond", "GSTcond") && f %in% Age_tuning[["Fleet"]]) {
#### these values not to be trusted in the presence of excluded data:
## HarmEffNage <- Age_tuning$"HarMean(effN)"[Age_tuning[["Fleet"]]==f]
## MeanNage <- Age_tuning$"mean(inputN*Adj)"[Age_tuning[["Fleet"]]==f]
HarmEffNage <- NULL
MeanNage <- NULL
} else {
HarmEffNage <- NULL
MeanNage <- NULL
}
# get table of info about the comp data of this kind/fleet
if (kind %in% c("AGE", "GSTAGE", "cond", "GSTcond")) {
data_info <- replist[["Age_comp_error_controls"]]
if (is.null(data_info)) {
# versions prior to 3.30.21
data_info <- replist[["age_data_info"]]
}
}
if (kind %in% c("LEN", "GSTLEN")) {
data_info <- replist[["Length_comp_error_controls"]]
if (is.null(data_info)) {
# versions prior to 3.30.21
data_info <- replist[["len_data_info"]]
}
}
if (kind %in% c("SIZE")) {
data_info <- replist[["Size_comp_error_controls"]]
} # output wasn't available in earlier versions of SS3
# loop over partitions (discard, retain, total)
for (j in unique(dbasef[["Part"]])) {
dbase <- dbasef[dbasef[["Part"]] == j, ]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
## # starting with SSv3.24a, the Yr.S column is already in the output, otherwise fill it in
## if(!"Yr.S" %in% names(dbase)){
## # add fraction of season to distinguish between samples
## dbase[["Yr.S"]] <- dbase[["Yr"]] + (0.5/nseasons)*dbase[["Seas"]]
## }
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(
table(dbase[["Yr.S"]], dbase[["Ageerr"]]) > 0,
1,
sum
))
if (max_n_ageerr > 1) {
if (ageerr_warning) {
message(
"Note: multiple samples with different ageing error types within fleet/year.\n",
" Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
" Bubble plots may be misleading with overlapping bubbles.\n"
)
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish between types within a year
dbase[["Yr.S"]] <- dbase[["Yr.S"]] + dbase[["Ageerr"]] / 1000
dbase[["YrSeasName"]] <- paste(
dbase[["YrSeasName"]],
"a",
dbase[["Ageerr"]],
sep = ""
)
}
## assemble pieces of plot title
# market category
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste("flt", f, filesex, "mkt", j, sep = "")
### subplot 1: multi-panel composition plot
if (1 %in% subplots & kind != "cond") {
# for age or length comps, but not conditional AAL
caption <- paste(titledata, title_sexmkt, fleetnames[f], sep = "") # total title
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun <- function(ipage, ...) {
sexvec <- dbase[["sex"]]
# a function to combine a bunch of repeated commands
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
# test for Dirichlet-Multinomial likelihood
if (
"Nsamp_DM" %in%
names(dbase) &&
any(!is.na(dbase[["Nsamp_DM"]]))
) {
# Dirichlet-Multinomial likelihood
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["Nsamp_DM"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "N input=",
effN_label = "N adj.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
} else {
# standard multinomial likelihood
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "N adj.=",
effN_label = "N eff.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
}
}
if (kind == "GSTAGE") {
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
showsampsize = FALSE,
showeffN = FALSE,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
}
if (kind == "GSTLEN") {
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
showsampsize = FALSE,
showeffN = FALSE,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
...
)
}
if (kind %in% c("L@A", "W@A")) {
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
ptsSD = dbase[["SD"]],
sampsize = dbase[["Nsamp_adj"]],
effN = 0,
showsampsize = FALSE,
showeffN = FALSE,
nlegends = 1,
legtext = list(dbase[["YrSeasName"]]),
bars = bars,
linepos = (1 - datonly) * linepos,
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = ifelse(kind == "W@A", labels[9], labels[1]),
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
...
)
}
} # end tempfun
if (plot) {
tempfun(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(dbase[["Yr.S"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption_count <- ""
if (npages > 1) {
pagetext <- paste0("_page", ipage)
caption_count <- paste0(" (plot ", ipage, " of ", npages, ")")
}
caption_extra <- ""
if (ipage == 1) {
if (
"Nsamp_DM" %in%
names(dbase) &&
any(!is.na(dbase[["Nsamp_DM"]]))
) {
# if using DM likelihood, get Theta value for this fleet & partition
# partition weighting added in 3.30.21
if ("partition" %in% names(data_info)) {
# added in 3.30.21
sub <- data_info[["Fleet"]] == f &
data_info[["partition"]] == j
} else {
if ("Fleet" %in% names(data_info)) {
# partition weighting not available for all comp types,
# so can't rely on "partition" to distinguish
sub <- data_info[["Fleet"]] == f
} else {
# earlier version without partition weighting
# has one row per fleet so no "fleet" column
sub <- f
}
}
if (kind != "SIZE") {
ParmSelect <- data_info[sub, "ParmSelect"]
CompError <- data_info[sub, "CompError"]
} else {
# Size_comp_error_controls have different format
sub <- sizemethod
ParmSelect <- data_info[sub, "error_parm_ID"]
CompError <- data_info[sub, "error_type"]
}
# D-M option 1 (linear)
if (CompError == 1) {
Theta <- as.numeric(replist[[
"Dirichlet_Multinomial_pars"
]][["Theta"]][ParmSelect])
# note: in caption below Θ = Theta
caption_extra <-
paste0(
".<br>'N input' is the input sample size. ",
"'N adj.' is the sample size after adjustment by the ",
"Dirichlet-Multinomial <i>Θ</i> parameter based on the ",
"formula N adj. = 1 / (1+<i>Θ</i>) + N * <i>Θ</i> / (1+<i>Θ</i>). ",
"<br>For this fleet, <i>Θ</i> = ",
round(Theta, 3),
" and the sample size multiplier is approximately ",
"<i>Θ</i> / (1+<i>Θ</i>) = ",
round(Theta / (1 + Theta), 3),
"."
)
}
# D-M option 2 (saturating)
if (CompError == 2) {
beta <- as.numeric(replist[[
"Dirichlet_Multinomial_pars"
]][["Theta"]][ParmSelect])
# note: in captions below β = beta
caption_extra <-
paste0(
".<br>'N input' is the input sample size. ",
"'N adj.' is the sample size after adjustment by the ",
"Dirichlet-Multinomial <i>β</i> parameter based on the ",
"formula N adj. = (N + N<i>β</i>) / (N + <i>β</i>). ",
"<br><br>For this fleet, <i>β</i> = ",
round(beta, 1),
". But due to the saturating functional form, there is no single ",
"approximate sample size multiplier."
)
}
caption_extra <- paste0(
caption_extra,
" For more info, see ",
"<a href=https://doi.org/10.1016/j.fishres.2016.06.005>",
"Thorson et al. (2017)</a>."
)
} else {
# if not using Dirichlet-Multinomial likelihood
caption_extra <-
paste0(
".<br><br>'N adj.' is the input sample size ",
"after data-weighting adjustment. ",
"N eff. is the calculated effective sample size used ",
"in the McAllister-Ianelli tuning method."
)
}
} # end check for ipage == 1
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
caption <- paste0(caption, caption_count, caption_extra)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun(ipage = ipage, ...)
dev.off()
}
}
} # end subplot 1
# some things related to the next two bubble plots (single or multi-panel)
if (datonly) {
z <- dbase[["Obs"]]
if (scalebubbles) {
z <- dbase[["Nsamp_adj"]] * dbase[["Obs"]] # if requested, scale by sample sizes
}
col <- rep("black", 2)
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
} else {
z <- dbase[["Pearson"]]
col <- rep(colvec[3], 2)
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
### subplot 2: single panel bubble plot for numbers at length or age
if (2 %in% subplots & bub & kind != "cond") {
# get growth curves if requested
if (length(cohortlines) > 0) {
growdat <- replist[["endgrowth"]]
growdatF <- growdat[
growdat[["Sex"]] == 1 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 1]),
]
if (nsexes > 1) {
growdatM <- growdat[
growdat[["Sex"]] == 2 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 2]),
]
}
}
# assemble caption that may also be used for plot title
caption <- paste(titletype, title_sexmkt, fleetnames[f], sep = "")
caption <- paste(
caption,
" (max=",
round(max(z), digits = 2),
")",
sep = ""
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun2 <- function() {
xvals <- dbase[["Yr.S"]]
# calculate smallest difference among years
# which is used to adjust offsets males and females
xdiff <- 0.1 *
sort(unique(diff(sort(unique(dbase[["Yr.S"]])), na.rm = TRUE)))[
1
]
# not sure what cases would have missing xdiff
# from above calculation, but it definitely happens
# with only one year of data, so setting default
# to work for that case
if (is.na(xdiff)) {
xdiff <- 0.1
}
cols <- rep(colvec[3], nrow(dbase))
if (nsexes > 1) {
xvals[dbase[["sex"]] > 0] <- dbase[["Yr.S"]][
dbase[["sex"]] > 0
] -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
if (length(unique(dbase[["Yr.S"]])) == 1) {
# if only one year, don't bother showing points
# as mid-year values
# this may not be ideal for seasonal models
xvals[dbase[["sex"]] > 0] <- floor(dbase[["Yr.S"]][
dbase[["sex"]] > 0
]) -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
}
cols[dbase[["sex"]] > 0] <- colvec[dbase[["sex"]][
dbase[["sex"]] > 0
]]
}
bubble3(
x = xvals,
y = dbase[["Bin"]],
z = z,
xlab = labels[3],
ylab = kindlab,
col = cols,
cexZ1 = cexZ1,
legend = bublegend,
las = 1,
main = ptitle,
cex.main = cex.main,
maxsize = pntscalar,
allopen = allopen,
minnbubble = minnbubble
)
# add lines for growth of individual cohorts if requested
if (length(cohortlines) > 0) {
for (icohort in seq_along(cohortlines)) {
message(" Adding line for", cohortlines[icohort], "cohort\n")
if (kind == "LEN") {
if (nsexes > 1) {
lines(
growdatF[["Age_Mid"]] + cohortlines[icohort],
growdatF[["Len_Mid"]],
col = colvec[1]
) # females
lines(
growdatM[["Age_Mid"]] + cohortlines[icohort],
growdatM[["Len_Mid"]],
col = colvec[2]
) # males
} else {
lines(
growdatF[["Age_Mid"]] + cohortlines[icohort],
growdatF[["Len_Mid"]],
col = colvec[3]
) # single-sex growth
}
}
if (kind %in% c("AGE", "GSTAGE")) {
lines(
0.5 +
c(cohortlines[icohort], cohortlines[icohort] + accuage),
c(0, accuage),
col = colvec[3],
lty = 3
) # one-one line for age
}
}
}
}
if (plot) {
tempfun2()
}
if (print) {
# set up plotting to png file if required
pagetext <- ""
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
if (length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
file <- paste(
filenamestart,
filetype,
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun2()
dev.off() # close device if png
}
} # end bubble plot
### subplot 3: multi-panel bubble plots for conditional age-at-length
if (3 %in% subplots & kind == "cond") {
# assemble caption that may also be used for plot title
caption1 <- paste(titletype, title_sexmkt, fleetnames[f], sep = "")
caption1 <- paste(
caption1,
" (max=",
round(max(z), digits = 2),
")",
sep = ""
)
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
# calculate scaling of lines showing effect and input sample size
sampsizeline.old <- sampsizeline
effNline.old <- effNline
if (is.logical(sampsizeline) && sampsizeline) {
# scaling when displaying only adjusted input sample size
sampsizeline <- max(dbase[["Bin"]]) /
max(dbase[["Nsamp_adj"]], na.rm = TRUE)
if (!datonly && is.logical(effNline) && effNline) {
# scaling when displaying both input and effective
sampsizeline <- effNline <- max(dbase[["Bin"]]) /
max(dbase[["Nsamp_adj"]], dbase[["effN"]], na.rm = TRUE)
message(
" Fleet ",
f,
" ",
titlesex,
"adj. input & effective N in red & green scaled by ",
effNline,
"\n",
sep = ""
)
} else {
message(
" Fleet ",
f,
" ",
titlesex,
"adj. input N in red scaled by ",
sampsizeline,
"\n",
sep = ""
)
}
}
# function to make plots
tempfun3 <- function(ipage, ...) {
sexvec <- dbase[["sex"]]
col.index <- sexvec
col.index[col.index == 0] <- 3
cols <- colvec[col.index]
yrvec <- dbase[["Yr.S"]] + dbase[["sex"]] * 1e-6
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Lbin_mid"]],
yr = yrvec,
size = z,
sampsize = dbase[["Nsamp_adj"]],
showsampsize = showsampsize,
effN = dbase[["effN"]],
showeffN = FALSE,
cexZ1 = cexZ1,
bublegend = bublegend,
nlegends = 1,
legtext = list(dbase[["YrSeasName"]]),
bars = FALSE,
linepos = 0,
main = ptitle,
cex.main = cex.main,
xlab = labels[2],
ylab = labels[1],
ymin0 = FALSE,
maxrows = maxrows2,
maxcols = maxcols2,
fixdims = fixdims,
allopen = allopen,
minnbubble = minnbubble,
# ptscol=col[1],ptscol2=col[2],
ptscol = cols,
ipage = ipage,
scalebins = scalebins,
sampsizeline = sampsizeline,
effNline = effNline,
sampsizemean = MeanNage,
effNmean = HarmEffNage,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
...
)
}
if (plot) {
tempfun3(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(dbase[["Yr.S"]])) *
length(unique(dbase[["sex"]])) /
maxrows2 /
maxcols2
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filetype,
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun3(ipage = ipage, ...)
dev.off() # close device if png
}
}
sampsizeline <- sampsizeline.old
effNline <- effNline.old
} # end conditional bubble plot
#
### subplots 4 and 5: multi-panel plot of point and line fit to
### conditional age-at-length
### and Pearson residuals of A-L key for specific years
if (
(4 %in% subplots | 5 %in% subplots) &
aalyear[1] > 0 &
kind == "cond"
) {
for (y in seq_along(aalyear)) {
aalyr <- aalyear[y]
if (length(dbase[["Obs"]][dbase[["Yr"]] == aalyr]) > 0) {
ydbase <- dbase[dbase[["Yr"]] == aalyr, ]
sexvec <- ydbase[["sex"]]
if (4 %in% subplots) {
### subplot 4: multi-panel plot of fit to conditional age-at-length for specific years
caption1 <- paste(
aalyr,
" age-at-length bin, ",
title_sexmkt,
fleetnames[f],
sep = ""
)
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
lenbinlegend <- paste(
ydbase[["Lbin_lo"]],
labels[7],
sep = ""
)
lenbinlegend[ydbase[["Lbin_range"]] > 0] <- paste(
ydbase[["Lbin_lo"]],
"-",
ydbase[["Lbin_hi"]],
labels[7],
sep = ""
)
tempfun4 <- function(ipage, ...) {
# temporary function to aid repeating the big function call
make_multifig(
ptsx = ydbase[["Bin"]],
ptsy = ydbase[["Obs"]],
yr = ydbase[["Lbin_lo"]],
linesx = ydbase[["Bin"]],
linesy = ydbase[["Exp"]],
sampsize = ydbase[["Nsamp_adj"]],
effN = ydbase[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
nlegends = 3,
legtext = list(lenbinlegend, "sampsize", "effN"),
bars = FALSE,
linepos = linepos,
main = ptitle,
cex.main = cex.main,
xlab = labels[2],
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
}
if (plot) {
tempfun4(ipage = 0, ...)
}
if (print) {
npages <- ceiling(
length(unique(ydbase[["Yr.S"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
if (length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
file <- paste0(
filenamestart,
filename_fltsexmkt,
"_",
aalyr,
pagetext,
".png"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun4(ipage = ipage, ...)
dev.off() # close device if print
}
}
} # end if 4 in subplots
if (5 %in% subplots) {
### subplot 5: Pearson residuals for A-L key
z <- ydbase[["Pearson"]]
col.index <- sexvec
col.index[col.index == 0] <- 3
cols <- colvec[col.index]
x.vec <- ydbase[["Bin"]] + ydbase[["sex"]] * 1e-6
caption <- paste(
aalyr,
" Pearson residuals for A-L key, ",
title_sexmkt,
fleetnames[f],
sep = ""
)
caption <- paste(
caption,
" (max=",
round(abs(max(z)), digits = 2),
")",
sep = ""
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun5 <- function() {
bubble3(
x = x.vec,
y = ydbase[["Lbin_lo"]],
z = z,
xlab = labels[2],
ylab = labels[1],
col = cols,
las = 1,
main = ptitle,
cex.main = cex.main,
maxsize = pntscalar,
cexZ1 = cexZ1,
legend = bublegend,
allopen = FALSE,
minnbubble = minnbubble
)
}
if (plot) {
tempfun5()
}
if (print) {
if (length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
file <- paste0(
filenamestart,
"yearresids_",
filename_fltsexmkt,
"_",
aalyr,
pagetext,
".png"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun5()
dev.off() # close device if print
}
} # end if 5 in subplots
}
}
}
### subplot 6: multi-panel plot of point and line fit to conditional
### age-at-length for specific length bins
if (6 %in% subplots & aalbin[1] > 0) {
badbins <- setdiff(aalbin, dbase[["Lbin_hi"]])
goodbins <- intersect(aalbin, dbase[["Lbin_hi"]])
if (length(goodbins) > 0) {
if (length(badbins) > 0) {
message(
"Error! the following inputs for 'aalbin' do not match the Lbin_hi values for the conditional age-at-length data:",
badbins,
"\n",
" the following inputs for 'aalbin' are fine:",
goodbins,
"\n"
)
}
for (ibin in seq_along(goodbins)) {
# loop over good bins
ilenbin <- goodbins[ibin]
abindbase <- dbase[dbase[["Lbin_hi"]] == ilenbin, ]
if (nrow(abindbase) > 0) {
# check for data associated with this bin
sexvec <- abindbase[["sex"]]
caption1 <- paste0(
"Age-at-length ",
ilenbin,
labels[7],
", ",
title_sexmkt,
fleetnames[f]
)
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun6 <- function(ipage, ...) {
# temporary function to aid repeating the big function call
make_multifig(
ptsx = abindbase[["Bin"]],
ptsy = abindbase[["Obs"]],
yr = abindbase[["Yr.S"]],
linesx = abindbase[["Bin"]],
linesy = abindbase[["Exp"]],
sampsize = abindbase[["Nsamp_adj"]],
effN = abindbase[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
nlegends = 3,
legtext = list(
abindbase[["YrSeasName"]],
"sampsize",
"effN"
),
bars = bars,
linepos = (1 - datonly) * linepos,
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
sexvec = sexvec,
...
)
}
if (plot) {
tempfun6(ipage = 0, ...)
}
if (print) {
npages <- ceiling(
length(unique(abindbase[["Yr.S"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste0("_page", ipage)
caption <- paste0(
caption,
" (plot ",
ipage,
" of ",
npages,
")"
)
}
file <- paste0(
filenamestart,
filename_fltsexmkt,
"_length",
ilenbin,
labels[7],
pagetext,
".png"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun6(ipage = ipage, ...)
dev.off() # close device if print
}
} # end print
} # end if data
} # end loop over length bins
} # end if length(goodbins)>0
} # end if plot requested
### subplot 7: sample size plot
if (
7 %in%
subplots &
samplesizeplots &
!datonly &
!("Nsamp_DM" %in%
names(dbase) &&
any(!is.na(dbase[["Nsamp_DM"]]))) &
!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))
) {
caption <- paste0(
"N-EffN comparison, ",
titledata,
title_sexmkt,
fleetnames[f]
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
lfitfunc <- function() {
if (kind == "cond") {
# trap nonrobust effective n's
# should this only be for conditional age-at-length or all plots?
dbasegood <- dbase[
dbase[["Obs"]] >= 0.0001 &
dbase[["Exp"]] < 0.99 &
!is.na(dbase[["effN"]]) &
dbase[["effN"]] < maxneff,
]
} else {
dbasegood <- dbase
}
if (nrow(dbasegood) > 0) {
# thinning out columns and removing rows with redundant information
# (for the purposes of this function)
dbasegood2 <- dbasegood[, c("YrSeasName", "Nsamp_adj", "effN")]
dbasegood2 <- unique(dbasegood2)
plot(
dbasegood2[["Nsamp_adj"]],
dbasegood2[["effN"]],
xlab = labels[4],
main = ptitle,
cex.main = cex.main,
ylim = c(0, 1.15 * max(dbasegood2[["effN"]])),
xlim = c(0, 1.15 * max(dbasegood2[["Nsamp_adj"]])),
col = colvec[3],
pch = 19,
ylab = labels[5],
xaxs = "i",
yaxs = "i"
)
# add labels for the years if requested
if (showyears) {
par(xpd = TRUE) # allows the label to go over plot boundary
text(
x = dbasegood2[["Nsamp_adj"]],
y = dbasegood2[["effN"]],
dbasegood2[["YrSeasName"]],
adj = c(-0.2, 0.5)
)
par(xpd = FALSE) # restores default clipping
}
abline(0, 1, col = "black", lty = 1)
# add loess smoother if there's at least 6 points with a range greater than 2
if (
smooth &
length(unique(dbasegood2[["Nsamp_adj"]])) > 6 &
diff(range(dbasegood2[["Nsamp_adj"]])) > 2
) {
old_warn <- options()[["warn"]] # previous warnings setting
options(warn = -1) # turn off loess warnings
psmooth <- loess(
dbasegood2[["effN"]] ~ dbasegood2[["Nsamp_adj"]],
degree = 1
)
options(warn = old_warn) # returning to old value
lines(
psmooth[["x"]][order(psmooth[["x"]])],
psmooth[["fit"]][order(psmooth[["x"]])],
lwd = 1.2,
col = "red",
lty = "dashed"
)
}
if (addMeans) {
# vertical line with label for mean input sample size
abline(
v = mean(dbasegood2[["Nsamp_adj"]]),
lty = "22",
col = "green3"
)
text(
x = mean(dbasegood2[["Nsamp_adj"]]),
y = 0,
col = "green3",
"arithmetic mean",
srt = 90,
adj = c(-0.1, -0.3)
)
# horizontal line with label for harmonic effective sample size
abline(
h = 1 / mean(1 / dbasegood2[["effN"]]),
lty = "22",
col = "green3"
)
text(
x = 0,
y = 1 / mean(1 / dbasegood2[["effN"]]),
col = "green3",
"harmonic mean",
adj = c(-0.1, -0.3)
)
}
}
}
if (plot) {
lfitfunc()
}
if (print) {
# set up plotting to png file if required
file <- paste(
filenamestart,
"sampsize_",
filename_fltsexmkt,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
lfitfunc()
dev.off()
}
} # end subplot 7
### subplots 8 and 9: Francis plot for non-conditional data
if (nrow(replist[["Dirichlet_Multinomial_pars"]]) == 0) {
# only include subplots 8 and 10 if using DM likelihood
subplots <- subplots[!subplots %in% c(8, 10)]
}
for (whichplot in 8:9) {
# don't include Francis adjustment in subplot 8, just 9
# and plot only values by partition in subplot 8
if (whichplot == 8) {
plotadj <- FALSE
part <- j
} else {
plotadj <- TRUE
part <- 0:2
}
if (whichplot %in% subplots & kind %in% c("LEN", "SIZE", "AGE")) {
# convert "AGE" to "age" so that SSMethod.TA1.8 can find "agedbase", etc.
kind2 <- tolower(kind)
if (plot) {
tmp <- SSMethod.TA1.8(
fit = replist,
type = kind2,
part = part,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj,
printit = FALSE
)
}
if (print) {
# set up plotting to png file if required
file <- paste0(
filenamestart,
ifelse(
whichplot == 8,
"data_weighting_timeseries_",
"data_weighting_TA1-8_"
),
fleetnames[f],
".png"
)
if (length(part) == 1) {
file <- gsub(
pattern = "timeseries",
replacement = paste0("timeseries_part", part),
x = file
)
}
# not using save_png because caption isn't available until after
# plot is created
# old command: plotinfo <- save_png(file=file, caption=caption)
png(
filename = file.path(plotdir, file),
width = pwidth,
height = pheight,
units = punits,
res = res,
pointsize = ptsize
)
# run function
tmp <- SSMethod.TA1.8(
fit = replist,
type = kind2,
part = part,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj,
printit = FALSE
)
# create caption
caption <- paste0(
"Mean ",
gsub("len", "length", tolower(kind)),
" for ",
fleetnames[f],
" with 95% confidence intervals",
" based on current sample sizes."
)
if (!datonly & whichplot == 8) {
caption <- gsub(
pattern = "sizes.",
replacement = "sizes (including any D-M weighting).",
x = caption
)
}
if (!datonly & whichplot == 9) {
caption <- paste0(
caption,
"<br>Francis data weighting method TA1.8:"
)
if (!is.null(tmp[1])) {
vals <- paste0(
"thinner intervals (with capped ends) show ",
"result of further adjusting sample sizes ",
"based on suggested multiplier ",
"(with 95% interval) for ",
kind2,
" data from ",
fleetnames[f],
":<br>",
round(tmp[1], 4),
" (",
round(tmp[2], 4),
"-",
round(tmp[3], 4),
")"
)
} else {
vals <- "too few points to calculate adjustments."
}
caption <- paste(
caption,
vals,
"<br>For more info, see ",
"<a href=https://doi.org/10.1139/f2011-025>",
"Francis (2011)</a>."
)
} # end additional caption text about Francis weighting
# add caption to the plotinfo table (normally done by save_png)
plotinfo <- rbind(
plotinfo,
data.frame(
file = file,
caption = caption,
alt_text = NA
)
)
dev.off() # close device if png
} # end test for print to PNG option
} # end subplot 8 or 9
} # end loop over 8 and 9
### subplots 10 and 11: Francis weighting for conditional data
for (whichplot in 10:11) {
# don't include Francis adjustment in subplot 10, just 11
plotadj <- (whichplot == 11)
if (
whichplot %in%
subplots &
kind == "cond" &
(f %in% condbase[["Fleet"]])
) {
if (plot) {
SSMethod.Cond.TA1.8(
fit = replist,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj
)
}
if (print) {
# set up plotting to png file if required
file <- paste0(
filenamestart,
ifelse(
whichplot == 10,
"data_weighting_timeseries_condAge",
"data_weighting_TA1-8_condAge"
),
gsub(" ", "", fleetnames[f]),
".png"
)
# not using save_png because caption isn't available until after
# plot is created
png(
filename = file.path(plotdir, file),
width = pwidth,
height = pheight,
units = punits,
res = res,
pointsize = ptsize
)
# run function
tmp <- SSMethod.Cond.TA1.8(
fit = replist,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj
)
# create caption
caption <- paste0(
"Mean age from conditional data",
" (aggregated across length bins) for ",
fleetnames[f],
" with 95% confidence intervals ",
" based on current samples sizes."
)
if (!datonly & whichplot == 10) {
caption <- gsub(
pattern = "sizes.",
replacement = "sizes (including any D-M weighting).",
x = caption
)
}
if (!datonly & whichplot == 11) {
caption <- paste0(
caption,
"<br>Francis data weighting method TA1.8:"
)
if (!is.null(tmp[1])) {
vals <- paste0(
"thinner intervals (with capped ends) show ",
"result of further adjusting sample sizes ",
"based on suggested multiplier ",
"(with 95% interval) for ",
"conditional age-at-length data from ",
fleetnames[f],
":<br>",
round(tmp[1], 4),
" (",
round(tmp[2], 4),
"-",
round(tmp[3], 4),
")",
sep = ""
)
} else {
vals <- "too few points to calculate adjustments."
}
caption <- paste(
caption,
vals,
"<br>For more info, see ",
"<a href=https://doi.org/10.1139/f2011-025>",
"Francis (2011)</a> and ",
"<a href=https://doi.org/10.1016/j.fishres.2015.12.006>",
"Punt (2017)</a>."
)
} # end test for datonly
# add caption to the plotinfo table (normally done by save_png)
plotinfo <- rbind(
plotinfo,
data.frame(
file = file,
caption = caption,
alt_text = NA
)
)
dev.off() # close device if png
} # end test for print to PNG option
} # end plot 10 or 11
} # end loop over 10 and 11
### subplot 12: Andre's mean age and std. dev. in conditional AAL
if (
12 %in%
subplots &
kind == "cond" &
length(unique(dbase[["Bin"]])) > 1
) {
caption1 <- paste(labels[14], title_sexmkt, fleetnames[f], sep = "")
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
andrefun <- function(ipage = 0) {
Lens <- sort(unique(dbase[["Lbin_lo"]]))
Yrs <- sort(unique(dbase[["Yr.S"]]))
ymax <- 1.1 * max(dbase[["Bin"]], na.rm = TRUE)
xmax <- max(condbase[["Lbin_hi"]], na.rm = TRUE)
xmin <- min(condbase[["Lbin_lo"]], na.rm = TRUE)
# do some stuff so that figures that span multiple pages can be output as separate PNG files
npanels <- length(Yrs)
npages <- npanels / andrerows
panelrange <- 1:npanels
if (npages > 1 & ipage != 0) {
panelrange <- intersect(
panelrange,
1:andrerows + andrerows * (ipage - 1)
)
}
Yrs2 <- Yrs[panelrange]
par(mfrow = c(andrerows, 2), mar = c(2, 4, 1, 1), oma = andre_oma)
for (Yr in Yrs2) {
y <- dbase[dbase[["Yr.S"]] == Yr, ]
Size <- NULL
Size2 <- NULL
Obs <- NULL
Obs2 <- NULL
Pred <- NULL
Pred2 <- NULL
Upp <- NULL
Low <- NULL
Upp2 <- NULL
Low2 <- NULL
for (Ilen in Lens) {
z <- y[y[["Lbin_lo"]] == Ilen, ]
if (length(z[, 1]) > 0) {
weightsPred <- z[["Exp"]] / sum(z[["Exp"]])
weightsObs <- z[["Obs"]] / sum(z[["Obs"]])
ObsV <- sum(z[["Bin"]] * weightsObs)
ObsV2 <- sum(z[["Bin"]] * z[["Bin"]] * weightsObs)
PredV <- sum(z[["Bin"]] * weightsPred)
PredV2 <- sum(z[["Bin"]] * z[["Bin"]] * weightsPred)
# Overdispersion on N
# NN <- z[["Nsamp_adj"]][1]*0.01 # Andre did this for reasons unknown
NN <- z[["Nsamp_adj"]][1]
if (
max(z[["Obs"]], na.rm = TRUE) > 1.0e-4 &
!is.na(NN) &&
NN > 0
) {
Size <- c(Size, Ilen)
Obs <- c(Obs, ObsV)
Pred <- c(Pred, PredV)
varn <- sqrt(PredV2 - PredV * PredV) / sqrt(NN)
Pred2 <- c(Pred2, varn)
varn <- sqrt(max(0, ObsV2 - ObsV * ObsV, na.rm = TRUE)) /
sqrt(NN)
Obs2 <- c(Obs2, varn)
Low <- c(Low, ObsV - 1.64 * varn)
Upp <- c(Upp, ObsV + 1.64 * varn)
if (NN > 1) {
Size2 <- c(Size2, Ilen)
Low2 <- c(
Low2,
varn * sqrt((NN - 1) / qchisq(0.95, NN))
)
Upp2 <- c(
Upp2,
varn * sqrt((NN - 1) / qchisq(0.05, NN))
)
}
}
}
}
if (length(Obs) > 0) {
## next line was replaced with setting at the top,
## for consistency across years
# ymax <- max(Pred,Obs,Upp)*1.1
plot(
Size,
Obs,
type = "n",
xlab = "",
ylab = "Age",
xlim = c(xmin, xmax),
ylim = c(0, ymax),
yaxs = "i"
)
label <- ifelse(nseasons == 1, floor(Yr), Yr)
text(
x = par("usr")[1],
y = .9 * ymax,
labels = label,
adj = c(-.5, 0),
font = 2,
cex = 1.2
)
if (length(Low) > 1) {
polygon(
c(Size, rev(Size)),
c(Low, rev(Upp)),
col = "grey95",
border = NA
)
}
if (!datonly) {
lines(Size, Pred, col = 4, lwd = 3)
}
points(Size, Obs, pch = 16)
lines(Size, Low, lty = 3)
lines(Size, Upp, lty = 3)
if (par("mfg")[1] == 1) {
title(
main = ptitle,
xlab = labels[1],
outer = TRUE,
line = 1
)
}
box()
ymax2 <- max(Obs2, Pred2, na.rm = TRUE) * 1.1
plot(
Size,
Obs2,
type = "n",
xlab = labels[1],
ylab = labels[13],
xlim = c(xmin, xmax),
ylim = c(0, ymax2),
yaxs = "i"
)
if (length(Low2) > 1) {
polygon(
c(Size2, rev(Size2)),
c(Low2, rev(Upp2)),
col = "grey95",
border = NA
)
}
if (!datonly) {
lines(Size, Pred2, col = 4, lwd = 3)
}
points(Size, Obs2, pch = 16)
lines(Size2, Low2, lty = 3)
lines(Size2, Upp2, lty = 3)
if (!datonly & par("mfg")[1] == 1) {
legend(
"topleft",
legend = c("Observed (with 90% interval)", "Expected"),
bty = "n",
col = c(1, 4),
pch = c(16, NA),
lty = c(NA, 1),
lwd = 3
)
}
box()
} # end if data exist
} # end loop over years
} # end andrefun
if (plot) {
andrefun()
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(length(unique(dbase[["Yr.S"]])) / andrerows)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
if (ipage == 1) {
# add more information only to first page of plots
caption <- paste(
caption,
"\nThese plots show mean age and std. dev. in conditional A@L.<br>",
"Left plots are mean A@L by size-class (obs. and exp.) ",
"with 90% CIs based on adding 1.64 SE of mean to the data.<br>",
"Right plots in each pair are SE of mean A@L (obs. and exp.) ",
"with 90% CIs based on the chi-square distribution.",
sep = ""
)
}
file <- paste(
filenamestart,
"Andre_plots",
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
andrefun(ipage = ipage)
dev.off() # close device if png
} # end loop over pages
} # end test for print to PNG option
} # end subplot 12
} # end loop over partitions (index j)
# } # end test for whether sex in vector of requested sexes
# } # end loop over combined/not-combined sex
} # end if data
} # end loop over fleets
if (!is.null(plotinfo)) {
plotinfo[["category"]] <- "Comp"
}
return(invisible(plotinfo))
} # end SSplotComps function
###########################
r4ss/r4ss documentation built on July 5, 2025, 3:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions SSplotComps
Documented in SSplotComps
#' Plot composition data and fits.
#'
#' Plot composition data and fits from Stock Synthesis output. Multi-figure
#' plots depend on `make_multifig`.
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows, where subplots 21 to 24
#' (aggregated across years) are provided first, and subplots
#' 1 to 10 are all repeated for each fleet
#' \itemize{
#' \item 1 index data by fleet
#' \item 1 multi-panel composition plot
#' \item 2 single panel bubble plot for numbers at length or age
#' \item 3 multi-panel bubble plots for conditional age-at-length
#' \item 4 multi-panel plot of fit to conditional age-at-length for specific years
#' \item 5 Pearson residuals for A-L key
#' \item 6 multi-panel plot of point and line fit to conditional
#' age-at-length for specific length bins
#' \item 7 sample size plot
#' \item 8 TA1.8 Francis plot for marginal data with
#' Dirichlet-Multinomial and no Francis adjustment
#' \item 9 TA1.8 Francis weighting plot for marginal data
#' \item 10 TA1.8 Francis plot for conditional data with
#' Dirichlet-Multinomial and no Francis adjustment
#' \item 11 TA1.8 Francis weighting plot for conditional data
#' \item 12 Andre's mean age and std. dev. in conditional AAL
#' \item 21 composition by fleet aggregating across years
#' \item 22 composition by fleet aggregating across years within each season
#' \item 23 composition by fleet aggregating across seasons within a year
#' \item 24 bubble plot comparison of length or age residuals
# across fleets within partition
#' }
#' @param kind indicator of type of plot can be "LEN", "SIZE", "AGE", "cond",
#' "GSTAGE", "GSTLEN", "L@A", or "W@A".
#' @param sizemethod if kind = "SIZE" then this switch chooses which of the
#' generalized size bin methods will be plotted.
#' @param aalyear Years to plot multi-panel conditional age-at-length fits for
#' all length bins; must be in a "c(YYYY,YYYY)" format. Useful for checking the
#' fit of a dominant year class, critical time period, etc. Default=-1.
#' @param aalbin The length bin for which multi-panel plots of the fit to
#' conditional age-at-length data will be produced for all years. Useful to
#' see if growth curves are ok, or to see the information on year classes move
#' through the conditional data. Default=-1.
#' @template plot
#' @template print
#' @template fleets
#' @template fleetnames
#' @param sexes which sexes to show plots for. Default="all" which will include
#' males, females, and unsexed. This option is not fully implemented for all
#' plots.
#' @param yupper upper limit on ymax for polygon/histogram composition plots
#' @param datonly make plots of data without fits?
#' @param samplesizeplots make sample size plots?
#' @param compresidplots make plots of residuals for fit to composition data?
#' @param bub make bubble plot for numbers at age or size?
#' @param showyears Add labels for years to sample size plots?
#' @param showsampsize add sample sizes to plot
#' @param showeffN add effective sample sizes to plot
#' @param aggregates_by_mkt separate plots of aggregates across years
#' into different plots for each market category (retained, discarded)?
#' @param sampsizeline show line for input sample sizes on top of conditional
#' age-at-length plots (TRUE/FALSE, still in development)
#' @param effNline show line for effective sample sizes on top of conditional
#' age-at-length plots (TRUE/FALSE, still in development)
#' @param minnbubble number of unique x values before adding buffer. see
#' ?bubble3 for more info.
#' @param pntscalar This scalar defines the maximum bubble size for bubble
#' plots. This option is still available but a better choice is to use cexZ1
#' which allow the same scaling throughout all plots.
#' @param scalebubbles scale data-only bubbles by sample size, not just
#' proportion within sample? Default=FALSE.
#' @param cexZ1 Character expansion (cex) for point associated with value of 1.
#' @param bublegend Add legend with example bubble sizes to bubble plots.
#' @param colvec Vector of length 3 with colors for females, males, unsexed fish
#' @param linescol Color for lines on top of polygons
#' @param xlas label style (las) input for x-axis. Default 0 has horizontal
#' labels, input 2 would provide vertical labels.
#' @param ylas label style (las) input for y-axis. Default NULL has horizontal
#' labels when all labels have fewer than 6 characters and vertical otherwise.
#' Input 0 would force vertical labels, and 1 would force horizontal.
#' @param axis1 optional position of bottom axis values
#' @param axis2 optional position of left size axis values
#' @param axis1labs optional vector of labels for axis1 (either NULL or needs to
#' match length of axis1)
#' @param sizebinlabs Vector of size bin labels corresponding to the generalized
#' size frequency method
#' @param red What color to use for females in bubble plots (default is slightly
#' transparent red)
#' @param blue What color to use for males in bubble plots (default is slightly
#' transparent blue)
#' @template pwidth
#' @template pheight
#' @template punits
#' @template ptsize
#' @template res
#' @template plotdir
#' @template cex.main
#' @param linepos should lines be added before points (linepos=1) or after
#' (linepos=2)?
#' @param fitbar show fit to bars instead of points
#' @param do.sqrt scale bubbles based on sqrt of size vector. see ?bubble3 for
#' more info.
#' @param smooth add loess smoother to observed vs. expected index plots and
#' input vs. effective sample size?
#' @param cohortlines optional vector of birth years for cohorts for which to
#' add growth curves to numbers at length bubble plots
#' @template labels
#' @param printmkt show market categories in plot titles?
#' @param printsex show sex in plot titles?
#' @param maxrows maximum (or fixed) number or rows of panels in the plot
#' @param maxcols maximum (or fixed) number or columns of panels in the plot
#' @param maxrows2 maximum number of rows for conditional age at length plots
#' @param maxcols2 maximum number of columns for conditional age at length
#' plots
#' @param rows number or rows to return to as default for next plots to come or
#' for single plots
#' @param cols number or cols to return to as default for next plots to come or
#' for single plots
#' @param andre_oma Outer margins passed to Andre's multi-panel conditional
#' age-at-length plots.
#' @param andrerows Number of rows of Andre's conditional age-at-length plots
#' within each page. Default=3.
#' @param fixdims fix the dimensions at maxrows by maxcols or resize based on
#' number of years of data
#' @param fixdims2 fix the dimensions at maxrows by maxcols in aggregate plots
#' or resize based on number of fleets
#' @param maxneff the maximum value to include on plots of input and effective
#' sample size. Occasionally a calculation of effective N blows up to very
#' large numbers, rendering it impossible to observe the relationship for other
#' data. Default=5000.
#' @template verbose
#' @param scalebins Rescale expected and observed proportions by dividing by
#' bin width for models where bins have different widths? Caution!: May not
#' work correctly in all cases.
#' @param addMeans Add parameter means in addition to medians for MCMC
#' posterior distributions in which the median and mean differ.
#' @template mainTitle
#' @param \dots additional arguments that will be passed to
#' the `par` command in the [make_multifig()] function.
#' @author Ian Taylor
#' @export
#' @seealso [SS_plots()], [make_multifig()]
SSplotComps <-
function(
replist,
subplots = c(1:10, 21, 24),
kind = "LEN",
sizemethod = 1,
aalyear = -1,
aalbin = -1,
plot = TRUE,
print = FALSE,
fleets = "all",
fleetnames = "default",
sexes = "all",
yupper = 0.4,
datonly = FALSE,
samplesizeplots = TRUE,
compresidplots = TRUE,
bub = FALSE,
showyears = TRUE,
showsampsize = TRUE,
showeffN = TRUE,
aggregates_by_mkt = FALSE,
sampsizeline = FALSE,
effNline = FALSE,
minnbubble = 3,
pntscalar = NULL,
scalebubbles = FALSE,
cexZ1 = 1.5,
bublegend = TRUE,
colvec = c(rgb(1, 0, 0, .7), rgb(0, 0, 1, .7), rgb(.1, .1, .1, .7)),
linescol = c(rgb(0, .5, 0, .7), rgb(.8, 0, 0, .7), rgb(0, 0, .8, .7)),
xlas = 0,
ylas = NULL,
axis1 = NULL,
axis2 = NULL,
axis1labs = NULL,
sizebinlabs = NULL,
blue = rgb(0, 0, 1, 0.7),
red = rgb(1, 0, 0, 0.7),
pwidth = 6.5,
pheight = 6.5,
punits = "in",
ptsize = 10,
res = 300,
plotdir = "default",
cex.main = 1,
linepos = 1,
fitbar = FALSE,
do.sqrt = TRUE,
smooth = TRUE,
cohortlines = c(),
labels = c(
"Length (cm)", # 1
"Age (yr)", # 2
"Year", # 3
"Observed sample size", # 4
"Effective sample size", # 5
"Proportion", # 6
"cm", # 7
"Frequency", # 8
"Weight", # 9
"Length", # 10
"(t)", # 11
"(numbers x1000)", # 12
"Stdev (Age)", # 13
"Conditional AAL plot, ", # 14
"Size bin"
), # 15
printmkt = TRUE,
printsex = TRUE,
maxrows = 6,
maxcols = 4,
maxrows2 = 4,
maxcols2 = 4,
rows = 1,
cols = 1,
andre_oma = c(3, 0, 3, 0),
andrerows = 4,
fixdims = TRUE,
fixdims2 = FALSE,
maxneff = 5000,
verbose = TRUE,
scalebins = FALSE,
addMeans = TRUE,
mainTitle = FALSE,
...
) {
################################################################################
# SSplotComps
################################################################################
# table to store information on each plot
plotinfo <- NULL
SS_versionNumeric <- replist[["SS_versionNumeric"]]
lendbase <- replist[["lendbase"]]
sizedbase <- replist[["sizedbase"]]
agedbase <- replist[["agedbase"]]
condbase <- replist[["condbase"]]
ghostagedbase <- replist[["ghostagedbase"]]
ghostlendbase <- replist[["ghostlendbase"]]
ladbase <- replist[["ladbase"]]
wadbase <- replist[["wadbase"]]
tagdbase1 <- replist[["tagdbase1"]]
tagdbase2 <- replist[["tagdbase2"]]
nfleets <- replist[["nfleets"]]
nseasons <- replist[["nseasons"]]
seasfracs <- replist[["seasfracs"]]
FleetNames <- replist[["FleetNames"]]
nsexes <- replist[["nsexes"]]
accuage <- replist[["accuage"]]
Age_tuning <- replist[["Age_comp_Eff_N_tuning_check"]]
# define a variety of titles and labels
titles <- NULL
titlemkt <- ""
if (plotdir == "default") {
plotdir <- replist[["inputs"]][["dir"]]
}
# sort out which fleets will be included
if (fleets[1] == "all") {
fleets <- 1:nfleets
} else {
if (length(intersect(fleets, 1:nfleets)) != length(fleets)) {
stop(
"Input 'fleets' should be 'all' or a vector of values between 1 and nfleets."
)
}
}
if (fleetnames[1] == "default") {
fleetnames <- FleetNames
}
# sort out which sexes will be included, and associated labels
if (sexes[1] == "all") {
sexes <- 0:nsexes # this can be used to subset stuff below
}
if (nsexes == 1) {
sexes <- 0:nsexes
}
if (nsexes == 1 | length(sexes) > 1) {
# single-sex models, or models with all sexes shown
# on the same plot, don't need sex-specific title
titlesex <- ""
filesex <- ""
}
if (nsexes > 1 & length(sexes) == 1) {
# multi-sex models with only 1 sex shown need
# title and filename info
if (sexes == 0) {
titlesex <- "sexes combined, "
filesex <- "sex0"
}
if (sexes == 1) {
titlesex <- "female, "
filesex <- "sex1"
}
if (sexes == 2) {
titlesex <- "male, "
filesex <- "sex2"
}
}
titlesex <- ifelse(printsex, titlesex, "")
# a few quantities related to data type and plot number
tail_warning <- FALSE
if (kind == "LEN") {
dbase_kind <- lendbase
kindlab <- labels[1]
if (datonly) {
filenamestart <- "comp_lendat_"
titledata <- "Length comp data, "
} else {
filenamestart <- "comp_lenfit_"
titledata <- "Length comps, "
}
if (
!is.null(replist[["Length_comp_error_controls"]]) &&
any(replist[["Length_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "GSTLEN") {
dbase_kind <- ghostlendbase
kindlab <- labels[1]
if (datonly) {
filenamestart <- "comp_gstlendat_"
titledata <- "Excluded length comp data, "
} else {
filenamestart <- "comp_gstlenfit_"
titledata <- "Excluded length comps, "
}
if (
!is.null(replist[["Length_comp_error_controls"]]) &&
any(replist[["Length_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "SIZE") {
dbase_kind <- sizedbase[sizedbase[["method"]] == sizemethod, ]
# if user provided sizebinlabs, then add tick and label
# associated with each bin
if (!is.null(sizebinlabs)) {
kindlab <- labels[15]
axis1 <- sort(unique(dbase_kind[["Bin"]]))
# confirm that vector of labels is the correct length
if (length(sizebinlabs) == length(axis1)) {
axis1labs <- sizebinlabs
} else {
axis1labs <- axis1
warning(
"Input 'sizebinlabs' differs in length from the unique Bin\n",
" values associated with sizemethod = ",
sizemethod,
". Using bin values instead."
)
}
} else {
# evenly spaced units in length or weight rather than sizebins
sizeunits <- unique(dbase_kind[["units"]])
if (length(sizeunits) > 1) {
stop(
"!error with size units in generalized size comp plots:\n",
" more than one unit value per method.\n"
)
}
if (sizeunits %in% c("in", "cm")) {
kindlab <- paste(labels[10], " (", sizeunits, ")", sep = "")
}
if (sizeunits %in% c("lb", "kg")) {
kindlab <- paste(labels[9], " (", sizeunits, ")", sep = "")
}
}
if (datonly) {
filenamestart <- "comp_sizedat_"
titledata <- "Size comp data, "
} else {
filenamestart <- "comp_sizefit_"
titledata <- "Size comps, "
}
# add text noting which size method is represented
if (length(unique(sizedbase[["method"]])) > 1) {
filenamestart <- paste0(filenamestart, "method", sizemethod, "_")
titledata <- paste0(titledata, " size method ", sizemethod, ", ")
}
}
if (kind == "AGE") {
dbase_kind <- agedbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_agedat_"
titledata <- "Age comp data, "
} else {
filenamestart <- "comp_agefit_"
titledata <- "Age comps, "
}
if (
!is.null(replist[["Age_comp_error_controls"]]) &&
any(replist[["Age_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "cond") {
dbase_kind <- condbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_condAALdat_"
titledata <- "Conditional age-at-length data, "
} else {
filenamestart <- "comp_condAALfit_"
titledata <- "Conditional age-at-length, "
}
}
if (kind == "GSTAGE") {
dbase_kind <- ghostagedbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_gstagedat_"
titledata <- "Excluded age comp data, "
} else {
filenamestart <- "comp_gstagefit_"
titledata <- "Excluded age comps, "
}
if (
!is.null(replist[["Age_comp_error_controls"]]) &&
any(replist[["Age_comp_error_controls"]][["mintailcomp"]] >= 0)
) {
tail_warning <- TRUE
}
}
if (kind == "GSTcond") {
dbase_kind <- ghostagedbase
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_gstCAALdat_"
titledata <- "Excluded conditional age-at-length data, "
} else {
filenamestart <- "comp_gstCAALfit_"
titledata <- "Excluded conditional age-at-length comps, "
}
}
if (kind == "L@A") {
dbase_kind <- ladbase[ladbase[["Nsamp_adj"]] != 0, ] # remove values with 0 sample size
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_LAAdat_"
titledata <- "Mean length at age data, "
} else {
filenamestart <- "comp_LAAfit_"
titledata <- "Mean length at age fit, "
}
dbase_kind[["SD"]] <- dbase_kind[["Lbin_lo"]] / dbase_kind[["Nsamp_adj"]]
}
if (kind == "W@A") {
dbase_kind <- wadbase[wadbase[["Nsamp_adj"]] != 0, ] # remove values with 0 sample size
kindlab <- labels[2]
if (datonly) {
filenamestart <- "comp_WAAdat_"
titledata <- "Mean weight at age data, "
} else {
filenamestart <- "comp_WAAfit_"
titledata <- "Mean weight at age fit, "
}
dbase_kind[["SD"]] <- dbase_kind[["Lbin_lo"]] / dbase_kind[["Nsamp_adj"]]
}
if (
!(kind %in%
c("LEN", "SIZE", "AGE", "cond", "GSTAGE", "GSTLEN", "L@A", "W@A"))
) {
stop(
"Input 'kind' to SSplotComps needs to be one of the following:\n ",
"'LEN','SIZE','AGE','cond','GSTAGE','GSTLEN','L@A','W@A'."
)
}
if (datonly) {
titletype <- titledata
} else {
titletype <- "Pearson residuals, "
}
# partition group is used by some aggregate plots (subplot 21+)
# either equal to partition or constant across all samples
if (nrow(dbase_kind) > 0) {
if (aggregates_by_mkt) {
dbase_kind[["Part_group"]] <- dbase_kind[["Part"]]
} else {
dbase_kind[["Part_group"]] <- -1 # code for all partitions combined
}
}
# Add asterix to indicate super periods and then remove rows labeled "skip".
# It would be better to somehow show the range of years, but that seems difficult.
if (
any(dbase_kind[["SuprPer"]] == "Sup" & dbase_kind[["Used"]] == "skip")
) {
message(
"Removing super-period composition values labeled 'skip'\n",
"and designating super-period values with a '*'"
)
dbase_kind <- dbase_kind[
dbase_kind[["SuprPer"]] == "No" | dbase_kind[["Used"]] != "skip",
]
dbase_kind[["YrSeasName"]] <- paste(
dbase_kind[["YrSeasName"]],
ifelse(dbase_kind[["SuprPer"]] == "Sup", "*", ""),
sep = ""
)
}
ageerr_warning <- TRUE
# subset data based on requested range of fleets and sexes
dbase_kind <- dbase_kind[
dbase_kind[["Fleet"]] %in% fleets & dbase_kind[["sex"]] %in% sexes,
]
### subplot 21: by fleet aggregating across years
if (21 %in% subplots & kind != "cond") {
# for age or length comps, but not conditional AAL
# check for the presence of data
if (nrow(dbase_kind) > 0) {
# loop over partitions (discard, retain, total)
for (j in unique(dbase_kind[["Part_group"]])) {
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
dbase <- dbase_kind[dbase_kind[["Part_group"]] == j, ]
if (nrow(dbase) > 0) {
# market category
if (j == -1) {
titlemkt <- ""
}
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
## aggregating identifiers for plot titles and filenames
## note: titlesex is set at the top of this function
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste(filesex)
if (j > -1) {
# add market category to filename if it's not a mix
filename_fltsexmkt <- paste0(filename_fltsexmkt, "mkt", j)
}
caption1 <- paste(
titledata,
title_sexmkt,
"aggregated across time by fleet",
sep = ""
) # total title
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
Bins <- sort(unique(dbase[["Bin"]]))
nbins <- length(Bins)
df <- data.frame(
Nsamp_adj = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
obs = dbase[["Obs"]] * dbase[["Nsamp_adj"]],
exp = dbase[["Exp"]] * dbase[["Nsamp_adj"]]
)
if (
"Nsamp_DM" %in% names(dbase) && any(!is.na(dbase[["Nsamp_DM"]]))
) {
df[["Nsamp_DM"]] <- dbase[["Nsamp_DM"]]
}
agg <- aggregate(
x = df,
by = list(
bin = dbase[["Bin"]],
f = dbase[["Fleet"]],
sex = dbase[["sex"]],
mkt = dbase[["Part"]]
),
FUN = sum
)
agg <- agg[agg[["f"]] %in% fleets, ]
agg[["obs"]] <- agg[["obs"]] / agg[["Nsamp_adj"]]
agg[["exp"]] <- agg[["exp"]] / agg[["Nsamp_adj"]]
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
# loop over fleets
for (f in unique(agg[["f"]])) {
# loop over fleets within market
for (mkt in unique(agg[["mkt"]][agg[["f"]] == f])) {
sub <- agg[["f"]] == f & agg[["mkt"]] == mkt
agg[["Nsamp_adj"]][sub] <- max(agg[["Nsamp_adj"]][sub])
if (
"Nsamp_DM" %in% names(agg) && any(!is.na(agg[["Nsamp_DM"]]))
) {
agg[["Nsamp_DM"]][sub] <- max(
agg[["Nsamp_DM"]][sub],
na.rm = TRUE
)
} else {
if (any(!is.na(agg[["effN"]][sub]))) {
agg[["effN"]][sub] <- max(agg[["effN"]][sub], na.rm = TRUE)
} else {
agg[["effN"]][sub] <- NA
}
}
}
}
namesvec <- fleetnames[agg[["f"]]]
# check for multiple market categories in a fleet to plot separately
max_n_mkt <- max(apply(table(agg[["f"]], agg[["mkt"]]) > 0, 1, sum))
if (max_n_mkt > 0) {
mktnames <- c("", "(discards)", "(retained)")
namesvec <- paste(
fleetnames[agg[["f"]]],
mktnames[agg[["mkt"]] + 1]
)
}
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
# group remaining calculations as a function
tempfun7 <- function(ipage, ...) {
# test for Dirichlet-Multinomial likelihood
if (
"Nsamp_DM" %in% names(agg) && any(!is.na(agg[["Nsamp_DM"]]))
) {
# Dirichlet-Multinomial likelihood
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["f"]] + agg[["mkt"]] / 10,
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["Nsamp_DM"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "Sum of N input=",
effN_label = "Sum of N adj.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(namesvec, "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
} else {
# standard multinomial likelihood
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["f"]] + agg[["mkt"]] / 10,
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "Sum of N adj.=",
effN_label = "Sum of N eff.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(namesvec, "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
}
}
if (plot) {
tempfun7(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(agg[["f"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
# add warning about tail compression
if (ipage == 1 & tail_warning) {
caption <- paste0(
caption,
".\n <br> ",
"Note: tail compression in use makes aggregated plots more difficult to interpret."
)
}
# add warning about partitions
if (max_n_mkt > 0 & ipage == 1) {
caption <-
paste0(
caption,
".\n <br> ",
"Labels 'retained' and 'discard' indicate",
" discarded or retained sampled for each fleet.",
" Panels without this designation represent",
" the whole catch.\n"
)
}
# remove double period if present
caption <- gsub("..", ".", caption, fixed = TRUE)
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
"<br> (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
"_aggregated_across_time.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun7(ipage = ipage, ...)
dev.off()
}
} # end print function
} else {
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=dbase[["effN"]],showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase[["YrSeasName"]],"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase[["YrSeasName"]]),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
} # end test for presence of observations in this partition group
} # end loop over partitions group
# } # end loop over combined/not-combined sex
} # end if data
} # end subplot 21
### subplot 22: by fleet aggregating across years within each season
if (22 %in% subplots & kind != "cond" & nseasons > 1) {
# for age or length comps, but not conditional AAL
dbasef <- dbase_kind[dbase_kind[["Fleet"]] %in% fleets, ]
if ("Nsamp_DM" %in% names(dbasef) && any(!is.na(dbasef[["Nsamp_DM"]]))) {
warning(
"Sample sizes in plots by fleet aggregating across years within each season have not yet been updated to reflect Dirichlet-Multinomial likelihood"
)
}
# check for the presence of data
if (nrow(dbasef) > 0) {
testor <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0
]) >
0
testor[2] <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3)
]) >
0
testor[3] <- length(dbasef[["sex"]][dbasef[["sex"]] == 2]) > 0
# loop over sex combinations
for (k in (1:3)[testor]) {
if (k == 1) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0,
]
}
if (k == 2) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3),
]
}
if (k == 3) {
dbase_k <- dbasef[dbasef[["sex"]] == 2, ]
}
sex <- ifelse(k == 3, 2, 1)
# loop over partitions (discard, retain, total)
for (j in unique(dbase_k[["Part"]])) {
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
dbase <- dbase_k[dbase_k[["Part"]] == j, ]
if (nrow(dbase) > 0) {
## assemble pieces of plot title
# sex
if (k == 1) {
titlesex <- "sexes combined, "
}
if (k == 2) {
titlesex <- "female, "
}
if (k == 3) {
titlesex <- "male, "
}
titlesex <- ifelse(printsex, titlesex, "")
# market category
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste("sex", k, "mkt", j, sep = "")
caption1 <- paste0(
titledata,
title_sexmkt,
"\naggregated within season by fleet"
) # total title
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
Bins <- sort(unique(dbase[["Bin"]]))
nbins <- length(Bins)
df <- data.frame(
Nsamp_adj = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
obs = dbase[["Obs"]] * dbase[["Nsamp_adj"]],
exp = dbase[["Exp"]] * dbase[["Nsamp_adj"]]
)
agg <- aggregate(
x = df,
by = list(
bin = dbase[["Bin"]],
f = dbase[["Fleet"]],
s = dbase[["Seas"]]
),
FUN = sum
)
agg <- agg[agg[["f"]] %in% fleets, ]
if (any(agg[["s"]] <= 0)) {
message(
"super-periods may not work correctly in plots of aggregated comps"
)
agg <- agg[agg[["s"]] > 0, ]
}
agg[["obs"]] <- agg[["obs"]] / agg[["Nsamp_adj"]]
agg[["exp"]] <- agg[["exp"]] / agg[["Nsamp_adj"]]
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for (f in unique(agg[["f"]])) {
# loop over fleets
for (s in unique(agg[["s"]][agg[["f"]] == f])) {
# loop over seasons within fleet
infleetseas <- agg[["f"]] == f & agg[["s"]] == s
agg[["Nsamp_adj"]][infleetseas] <- max(agg[["Nsamp_adj"]][
infleetseas
])
agg[["effN"]][infleetseas] <- max(agg[["effN"]][infleetseas])
}
}
agg[["fseas"]] <- agg[["f"]] + seasfracs[agg[["s"]]]
namesvec <- paste(
fleetnames[agg[["f"]]],
" s",
agg[["s"]],
sep = ""
)
# group remaining calculations as a function
tempfun8 <- function(ipage, ...) {
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["fseas"]],
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(namesvec, "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=dbase[["effN"]],showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase[["YrSeasName"]],"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase[["YrSeasName"]]),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
if (plot) {
tempfun8(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(agg[["fseas"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
caption <- caption1
pagetext <- ""
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
"_aggregated_within_season.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun8(ipage = ipage, ...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined sex
} # end if data
} # end subplot 22
### subplot 23: by fleet aggregating across seasons within a year
if (23 %in% subplots & kind != "cond" & nseasons > 1) {
# for age or length comps, but not conditional AAL
# loop over fleets
for (f in fleets) {
dbasef <- dbase_kind[dbase_kind[["Fleet"]] == f, ]
if (
"Nsamp_DM" %in% names(dbasef) && any(!is.na(dbasef[["Nsamp_DM"]]))
) {
warning(
"Sample sizes in plots by fleet aggregating across seasons within a year have not yet been updated to reflect Dirichlet-Multinomial likelihood"
)
}
# check for the presence of data
if (nrow(dbasef) > 0) {
testor <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0
]) >
0
testor[2] <- length(dbasef[["sex"]][
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3)
]) >
0
testor[3] <- length(dbasef[["sex"]][dbasef[["sex"]] == 2]) > 0
# loop over sex combinations
for (k in (1:3)[testor]) {
if (k == 1) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] == 0,
]
}
if (k == 2) {
dbase_k <- dbasef[
dbasef[["sex"]] == 1 & dbasef[["Pick_sex"]] %in% c(1, 3),
]
}
if (k == 3) {
dbase_k <- dbasef[dbasef[["sex"]] == 2, ]
}
sex <- ifelse(k == 3, 2, 1)
# loop over partitions (discard, retain, total)
for (j in unique(dbase_k[["Part"]])) {
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
dbase <- dbase_k[dbase_k[["Part"]] == j, ]
if (nrow(dbase) > 0) {
## assemble pieces of plot title
# sex
if (k == 1) {
titlesex <- "sexes combined, "
}
if (k == 2) {
titlesex <- "female, "
}
if (k == 3) {
titlesex <- "male, "
}
titlesex <- ifelse(printsex, titlesex, "")
# market category
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste(
"flt",
f,
"sex",
k,
"mkt",
j,
sep = ""
)
Bins <- sort(unique(dbase[["Bin"]]))
nbins <- length(Bins)
df <- data.frame(
Nsamp_adj = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
obs = dbase[["Obs"]] * dbase[["Nsamp_adj"]],
exp = dbase[["Exp"]] * dbase[["Nsamp_adj"]]
)
agg <- aggregate(
x = df,
by = list(
bin = dbase[["Bin"]],
f = dbase[["Fleet"]],
y = floor(dbase[["Yr.S"]])
),
FUN = sum
)
agg <- agg[agg[["f"]] %in% fleets, ]
agg[["obs"]] <- agg[["obs"]] / agg[["Nsamp_adj"]]
agg[["exp"]] <- agg[["exp"]] / agg[["Nsamp_adj"]]
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for (f in unique(agg[["f"]])) {
for (y in unique(agg[["y"]][agg[["f"]] == f])) {
infleetyr <- agg[["f"]] == f & agg[["y"]] == y
agg[["Nsamp_adj"]][infleetyr] <- max(agg[["Nsamp_adj"]][
infleetyr
])
agg[["effN"]][infleetyr] <- max(agg[["effN"]][infleetyr])
}
}
agg[["fy"]] <- agg[["f"]] + agg[["y"]] / 10000
# total title
caption <- paste(
titledata,
title_sexmkt,
fleetnames[f],
"\naggregated across seasons within year",
sep = ""
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
# group remaining calculations as a function
tempfun9 <- function(ipage, ...) {
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
make_multifig(
ptsx = agg[["bin"]],
ptsy = agg[["obs"]],
yr = agg[["fy"]],
linesx = agg[["bin"]],
linesy = agg[["exp"]],
sampsize = agg[["Nsamp_adj"]],
effN = agg[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(agg[["y"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims2,
ipage = ipage,
lwd = 2,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = agg[["sex"]],
yupper = yupper,
...
)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=dbase[["effN"]],showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase[["YrSeasName"]],"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase[["Bin"]],ptsy=dbase[["Obs"]],yr=dbase[["Yr.S"]],linesx=dbase[["Bin"]],linesy=dbase[["Exp"]],
## sampsize=dbase[["Nsamp_adj"]],effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase[["YrSeasName"]]),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
} # end tempfun
if (plot) {
tempfun9(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(agg[["fy"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
"_aggregated_across_seasons_within_year.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun9(ipage = ipage, ...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined sex
} # end if data
} # end loop over fleets
} # end subplot 23
### subplot 24: bubble plot comparison of length or age residuals
### across fleets within sex/partition
if (24 %in% subplots & kind %in% c("LEN", "AGE")) {
# loop over partitions groups (everything, or separate discard, retain, total)
for (j in unique(dbase_kind[["Part_group"]])) {
# subset data.frame for this partition group and subset of fleets of interest
dbase_parts <- dbase_kind[dbase_kind[["Part_group"]] == j, ]
# new column combining fleet and partition
# where 3.1 is fleet=3, partition=1
dbase_parts[["FleetPart"]] <- dbase_parts[["Fleet"]] +
0.1 * dbase_parts[["Part"]]
# table of info on each panel
panel_table <- data.frame(
FleetPart = sort(unique(dbase_parts[["FleetPart"]]))
)
# separate the pieces out again
panel_table[["Fleet"]] <- floor(panel_table[["FleetPart"]])
# round below is necessary because values were coming out as
# 1.9999999999999996 instead of 2
panel_table[["Part"]] <- round(
10 * (panel_table[["FleetPart"]] - panel_table[["Fleet"]])
)
# fleet name to use for each panel
panel_table[["Name"]] <- fleetnames[panel_table[["Fleet"]]]
# check for multiple market categories in a fleet to plot separately
max_n_mkt <- max(apply(
table(
panel_table[["Fleet"]],
panel_table[["Part"]]
) >
0,
1,
sum
))
if (max_n_mkt > 1) {
# if multiple categories within a fleet, append to name
mktnames <- c("", "(discards)", "(retained)")
panel_table[["Name"]] <- paste(
panel_table[["Name"]],
mktnames[panel_table[["Part"]] + 1]
)
}
npanels <- nrow(panel_table)
panelvec <- 1:npanels
xlim <- range(dbase_parts[["Yr.S"]]) # set xlim based on range across all fleets
xaxislab <- sort(unique(floor(dbase_parts[["Yr.S"]]))) # label with all years
# get growth curves if requested
if (length(cohortlines) > 0) {
growdat <- replist[["endgrowth"]]
growdatF <- growdat[
growdat[["Sex"]] == 1 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 1]),
]
if (nsexes > 1) {
growdatM <- growdat[
growdat[["Sex"]] == 2 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 2]),
]
}
}
# market category
if (j == -1) {
titlemkt <- ""
}
if (j == 0) {
titlemkt <- "whole catch"
}
if (j == 1) {
titlemkt <- "discard"
}
if (j == 2) {
titlemkt <- "retained"
}
titlemkt <- ifelse(printmkt, titlemkt, "")
## title_sexmkt <- paste(titlesex,titlemkt,sep="")
caption_base <- paste0(titletype, titlemkt, ", comparing across fleets")
# hack to remove ", ," from caption
caption_base <- gsub(", ,", ", ", caption_base)
if (mainTitle) {
ptitle <- caption_base
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
# if not all market categories, append mkt value to filename
filenamemkt <- ifelse(j > -1, paste("mkt", j, sep = ""), "")
multifleet.bubble.fun <- function(ipage = 0) {
# a function to wrap up multi-fleet bubble plots
# old graphics parameter settings
par_old <- par()
# multi-figure plot with as many rows as fleets, or the maxrows value
par(
mfrow = c(min(npanels, maxrows), 1),
mar = c(0.5, 0, 0, 0),
oma = c(4, 6, ifelse(mainTitle, 3, 1), 1)
)
# set up some stuff for cases where there are more fleets than panels in one plot
panelrange <- 1:npanels
npages <- ceiling(npanels / maxrows) # how many pages of plots
if (npages > 1 & ipage != 0) {
# range of which panels to print for each page
panelrange <- intersect(
panelrange,
1:maxrows + maxrows * (ipage - 1)
)
}
# loop over panels (fleet x partition)
for (ipanel in panelvec[panelrange]) {
flt <- panel_table[["Fleet"]][ipanel] # fleet number
mkt <- panel_table[["Part"]][ipanel] # market category
# subset database of values
dbase <- dbase_parts[
dbase_parts[["Fleet"]] == flt &
dbase_parts[["Part"]] == mkt,
]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), sex, and partition,
### not sure if multiple ageing error methods is supported at the moment,
### haven't tested -Ian 6/7/17
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(
table(dbase[["Yr.S"]], dbase[["Ageerr"]]) > 0,
1,
sum
))
if (max_n_ageerr > 1) {
if (ageerr_warning) {
message(
"Multiple samples with different ageing error types within fleet/year.\n",
"Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
"Bubble plots may be misleading with overlapping bubbles.\n"
)
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish
# between types within a year (may not work well for this plot)
dbase[["Yr.S"]] <- dbase[["Yr.S"]] +
dbase[["Ageerr"]] / (1000 * max_n_ageerr)
dbase[["YrSeasName"]] <- paste(
dbase[["YrSeasName"]],
"a",
dbase[["Ageerr"]],
sep = ""
)
}
# calculate smallest difference among years
# which is used to adjust offsets males and females
xdiff <- 0.1 *
sort(unique(diff(sort(unique(dbase[["Yr.S"]])), na.rm = TRUE)))[1]
# not sure what cases would have missing xdiff
# from above calculation, but it definitely happens
# with only one year of data, so setting default
# to work for that case
if (is.na(xdiff)) {
xdiff <- 0.1
}
# define colors
xvals <- dbase[["Yr.S"]]
cols <- rep(colvec[3], nrow(dbase))
if (nsexes > 1) {
xvals[dbase[["sex"]] > 0] <- dbase[["Yr.S"]][dbase[["sex"]] > 0] -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
if (length(unique(dbase[["Yr.S"]])) == 1) {
# if only one year, don't bother showing points
# as mid-year values
# this may not be ideal for seasonal models
xvals[dbase[["sex"]] > 0] <- floor(dbase[["Yr.S"]][
dbase[["sex"]] > 0
]) -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
}
cols[dbase[["sex"]] > 0] <- colvec[dbase[["sex"]][
dbase[["sex"]] > 0
]]
}
# determine bubble size and colors
if (datonly) {
z <- dbase[["Obs"]]
if (scalebubbles) {
z <- dbase[["Nsamp_adj"]] * dbase[["Obs"]]
} # if requested, scale by sample sizes
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
} else {
z <- dbase[["Pearson"]]
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
# make bubbles for a single fleet
# this section is a modified version of tempfun2 above
ylim <- range(dbase[["Bin"]])
ylim[2] <- ylim[2] + 0.2 * diff(ylim) # add buffer of 10% at the top for fleet name
bubble3(
x = xvals,
y = dbase[["Bin"]],
z = z,
col = cols,
cexZ1 = cexZ1,
legend = bublegend,
las = 1,
main = "",
cex.main = cex.main,
maxsize = pntscalar,
allopen = allopen,
xlim = xlim,
ylim = ylim,
axis1 = FALSE
)
### add label at top left of each panel
legend(
"topleft",
title = panel_table[["Name"]][ipanel],
legend = NA,
bty = "n",
cex = 1.5
)
### alternative way with legends on the side
### (probably not as good once the partition has been added)
# mtext(namesvec[ipanel],side=2,line=4.5,cex=par()[["cex"]])
# add lines for growth of individual cohorts if requested
if (length(cohortlines) > 0) {
for (icohort in seq_along(cohortlines)) {
message(" Adding line for", cohortlines[icohort], "cohort\n")
if (kind == "LEN") {
lines(
growdatF[["Age"]] + cohortlines[icohort],
growdatF[["Len_Mid"]],
col = colvec[1]
) # females
if (nsexes > 1) {
lines(
growdatM[["Age"]] + cohortlines[icohort],
growdatM[["Len_Mid"]],
col = colvec[2]
) # males
}
}
if (kind == "AGE") {
lines(
0.5 +
c(cohortlines[icohort], cohortlines[icohort] + accuage),
c(0, accuage),
col = "red"
)
}
}
}
if (
par()[["mfg"]][1] == par()[["mfg"]][3] |
ipanel == tail(panelvec, 1)
) {
# label all years on x-axis of last panel
axis(1, at = xaxislab)
} else {
# or just tick marks for other panels
axis(1, at = xaxislab, labels = rep("", length(xaxislab)))
}
if (par()[["mfg"]][1] == 1) {
# add title after making first panel
title(
main = ptitle,
outer = TRUE,
xlab = labels[3],
ylab = kindlab
)
}
} # end loop over fleets
# restore previous graphics parameter settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
} # end function wrapping up a single page of the residual comparison plot
# make plots or write to PNG file
if (length(panelvec) > 0) {
if (plot) {
multifleet.bubble.fun(ipage = 0)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(nrow(panel_table) / maxrows)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption_base
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste0(
caption,
" (plot ",
ipage,
" of ",
npages,
")"
)
}
if (ipage == 1 & length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
#### current scaling allows comparison across panels, so warning below
#### has been turned off
## caption <- paste(caption,
## "<br>Note: bubble sizes are scaled to maximum within each panel.",
## "<br>Thus, comparisons across panels should focus on patterns, not bubble sizes.")
file <- paste(
filenamestart,
filenamemkt,
pagetext,
"_multi-fleet_comparison.png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
multifleet.bubble.fun(ipage = ipage)
dev.off()
} # end loop over pages within printing PNG
} # end printing to PNG files
} # end test for non-zero number of fleets
} # end loop over partitions
## } # end loop over sexes
## } # end loop over sex combinations
# restore default single panel settings
par(mfcol = c(rows, cols), mar = c(5, 4, 4, 2) + .1, oma = rep(0, 4))
} # end subplot 24
# loop over fleets
for (f in fleets) {
# check for the presence of data
if (length(dbase_kind[["Obs"]][dbase_kind[["Fleet"]] == f]) > 0) {
dbasef <- dbase_kind[dbase_kind[["Fleet"]] == f, ]
# get mean sample quantities to show in conditional age-at-length figures
if (kind %in% c("cond", "GSTcond") && f %in% Age_tuning[["Fleet"]]) {
#### these values not to be trusted in the presence of excluded data:
## HarmEffNage <- Age_tuning$"HarMean(effN)"[Age_tuning[["Fleet"]]==f]
## MeanNage <- Age_tuning$"mean(inputN*Adj)"[Age_tuning[["Fleet"]]==f]
HarmEffNage <- NULL
MeanNage <- NULL
} else {
HarmEffNage <- NULL
MeanNage <- NULL
}
# get table of info about the comp data of this kind/fleet
if (kind %in% c("AGE", "GSTAGE", "cond", "GSTcond")) {
data_info <- replist[["Age_comp_error_controls"]]
if (is.null(data_info)) {
# versions prior to 3.30.21
data_info <- replist[["age_data_info"]]
}
}
if (kind %in% c("LEN", "GSTLEN")) {
data_info <- replist[["Length_comp_error_controls"]]
if (is.null(data_info)) {
# versions prior to 3.30.21
data_info <- replist[["len_data_info"]]
}
}
if (kind %in% c("SIZE")) {
data_info <- replist[["Size_comp_error_controls"]]
} # output wasn't available in earlier versions of SS3
# loop over partitions (discard, retain, total)
for (j in unique(dbasef[["Part"]])) {
dbase <- dbasef[dbasef[["Part"]] == j, ]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, sex, and partition
## # starting with SSv3.24a, the Yr.S column is already in the output, otherwise fill it in
## if(!"Yr.S" %in% names(dbase)){
## # add fraction of season to distinguish between samples
## dbase[["Yr.S"]] <- dbase[["Yr"]] + (0.5/nseasons)*dbase[["Seas"]]
## }
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(
table(dbase[["Yr.S"]], dbase[["Ageerr"]]) > 0,
1,
sum
))
if (max_n_ageerr > 1) {
if (ageerr_warning) {
message(
"Note: multiple samples with different ageing error types within fleet/year.\n",
" Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
" Bubble plots may be misleading with overlapping bubbles.\n"
)
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish between types within a year
dbase[["Yr.S"]] <- dbase[["Yr.S"]] + dbase[["Ageerr"]] / 1000
dbase[["YrSeasName"]] <- paste(
dbase[["YrSeasName"]],
"a",
dbase[["Ageerr"]],
sep = ""
)
}
## assemble pieces of plot title
# market category
if (j == 0) {
titlemkt <- "whole catch, "
}
if (j == 1) {
titlemkt <- "discard, "
}
if (j == 2) {
titlemkt <- "retained, "
}
titlemkt <- ifelse(printmkt, titlemkt, "")
# plot bars for data only or if input 'fitbar=TRUE'
if (datonly | fitbar) {
bars <- TRUE
} else {
bars <- FALSE
}
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex, titlemkt, sep = "")
filename_fltsexmkt <- paste("flt", f, filesex, "mkt", j, sep = "")
### subplot 1: multi-panel composition plot
if (1 %in% subplots & kind != "cond") {
# for age or length comps, but not conditional AAL
caption <- paste(titledata, title_sexmkt, fleetnames[f], sep = "") # total title
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun <- function(ipage, ...) {
sexvec <- dbase[["sex"]]
# a function to combine a bunch of repeated commands
if (!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))) {
# test for Dirichlet-Multinomial likelihood
if (
"Nsamp_DM" %in%
names(dbase) &&
any(!is.na(dbase[["Nsamp_DM"]]))
) {
# Dirichlet-Multinomial likelihood
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["Nsamp_DM"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "N input=",
effN_label = "N adj.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
} else {
# standard multinomial likelihood
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
sampsize_label = "N adj.=",
effN_label = "N eff.=",
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
}
}
if (kind == "GSTAGE") {
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
showsampsize = FALSE,
showeffN = FALSE,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
}
if (kind == "GSTLEN") {
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
sampsize = dbase[["Nsamp_adj"]],
effN = dbase[["effN"]],
showsampsize = FALSE,
showeffN = FALSE,
bars = bars,
linepos = (1 - datonly) * linepos,
nlegends = 3,
legtext = list(dbase[["YrSeasName"]], "sampsize", "effN"),
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
...
)
}
if (kind %in% c("L@A", "W@A")) {
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Obs"]],
yr = dbase[["Yr.S"]],
linesx = dbase[["Bin"]],
linesy = dbase[["Exp"]],
ptsSD = dbase[["SD"]],
sampsize = dbase[["Nsamp_adj"]],
effN = 0,
showsampsize = FALSE,
showeffN = FALSE,
nlegends = 1,
legtext = list(dbase[["YrSeasName"]]),
bars = bars,
linepos = (1 - datonly) * linepos,
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = ifelse(kind == "W@A", labels[9], labels[1]),
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
...
)
}
} # end tempfun
if (plot) {
tempfun(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(dbase[["Yr.S"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption_count <- ""
if (npages > 1) {
pagetext <- paste0("_page", ipage)
caption_count <- paste0(" (plot ", ipage, " of ", npages, ")")
}
caption_extra <- ""
if (ipage == 1) {
if (
"Nsamp_DM" %in%
names(dbase) &&
any(!is.na(dbase[["Nsamp_DM"]]))
) {
# if using DM likelihood, get Theta value for this fleet & partition
# partition weighting added in 3.30.21
if ("partition" %in% names(data_info)) {
# added in 3.30.21
sub <- data_info[["Fleet"]] == f &
data_info[["partition"]] == j
} else {
if ("Fleet" %in% names(data_info)) {
# partition weighting not available for all comp types,
# so can't rely on "partition" to distinguish
sub <- data_info[["Fleet"]] == f
} else {
# earlier version without partition weighting
# has one row per fleet so no "fleet" column
sub <- f
}
}
if (kind != "SIZE") {
ParmSelect <- data_info[sub, "ParmSelect"]
CompError <- data_info[sub, "CompError"]
} else {
# Size_comp_error_controls have different format
sub <- sizemethod
ParmSelect <- data_info[sub, "error_parm_ID"]
CompError <- data_info[sub, "error_type"]
}
# D-M option 1 (linear)
if (CompError == 1) {
Theta <- as.numeric(replist[[
"Dirichlet_Multinomial_pars"
]][["Theta"]][ParmSelect])
# note: in caption below Θ = Theta
caption_extra <-
paste0(
".<br>'N input' is the input sample size. ",
"'N adj.' is the sample size after adjustment by the ",
"Dirichlet-Multinomial <i>Θ</i> parameter based on the ",
"formula N adj. = 1 / (1+<i>Θ</i>) + N * <i>Θ</i> / (1+<i>Θ</i>). ",
"<br>For this fleet, <i>Θ</i> = ",
round(Theta, 3),
" and the sample size multiplier is approximately ",
"<i>Θ</i> / (1+<i>Θ</i>) = ",
round(Theta / (1 + Theta), 3),
"."
)
}
# D-M option 2 (saturating)
if (CompError == 2) {
beta <- as.numeric(replist[[
"Dirichlet_Multinomial_pars"
]][["Theta"]][ParmSelect])
# note: in captions below β = beta
caption_extra <-
paste0(
".<br>'N input' is the input sample size. ",
"'N adj.' is the sample size after adjustment by the ",
"Dirichlet-Multinomial <i>β</i> parameter based on the ",
"formula N adj. = (N + N<i>β</i>) / (N + <i>β</i>). ",
"<br><br>For this fleet, <i>β</i> = ",
round(beta, 1),
". But due to the saturating functional form, there is no single ",
"approximate sample size multiplier."
)
}
caption_extra <- paste0(
caption_extra,
" For more info, see ",
"<a href=https://doi.org/10.1016/j.fishres.2016.06.005>",
"Thorson et al. (2017)</a>."
)
} else {
# if not using Dirichlet-Multinomial likelihood
caption_extra <-
paste0(
".<br><br>'N adj.' is the input sample size ",
"after data-weighting adjustment. ",
"N eff. is the calculated effective sample size used ",
"in the McAllister-Ianelli tuning method."
)
}
} # end check for ipage == 1
file <- paste(
filenamestart,
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
caption <- paste0(caption, caption_count, caption_extra)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun(ipage = ipage, ...)
dev.off()
}
}
} # end subplot 1
# some things related to the next two bubble plots (single or multi-panel)
if (datonly) {
z <- dbase[["Obs"]]
if (scalebubbles) {
z <- dbase[["Nsamp_adj"]] * dbase[["Obs"]] # if requested, scale by sample sizes
}
col <- rep("black", 2)
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
} else {
z <- dbase[["Pearson"]]
col <- rep(colvec[3], 2)
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
### subplot 2: single panel bubble plot for numbers at length or age
if (2 %in% subplots & bub & kind != "cond") {
# get growth curves if requested
if (length(cohortlines) > 0) {
growdat <- replist[["endgrowth"]]
growdatF <- growdat[
growdat[["Sex"]] == 1 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 1]),
]
if (nsexes > 1) {
growdatM <- growdat[
growdat[["Sex"]] == 2 &
growdat[["Morph"]] ==
min(growdat[["Morph"]][growdat[["Sex"]] == 2]),
]
}
}
# assemble caption that may also be used for plot title
caption <- paste(titletype, title_sexmkt, fleetnames[f], sep = "")
caption <- paste(
caption,
" (max=",
round(max(z), digits = 2),
")",
sep = ""
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun2 <- function() {
xvals <- dbase[["Yr.S"]]
# calculate smallest difference among years
# which is used to adjust offsets males and females
xdiff <- 0.1 *
sort(unique(diff(sort(unique(dbase[["Yr.S"]])), na.rm = TRUE)))[
1
]
# not sure what cases would have missing xdiff
# from above calculation, but it definitely happens
# with only one year of data, so setting default
# to work for that case
if (is.na(xdiff)) {
xdiff <- 0.1
}
cols <- rep(colvec[3], nrow(dbase))
if (nsexes > 1) {
xvals[dbase[["sex"]] > 0] <- dbase[["Yr.S"]][
dbase[["sex"]] > 0
] -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
if (length(unique(dbase[["Yr.S"]])) == 1) {
# if only one year, don't bother showing points
# as mid-year values
# this may not be ideal for seasonal models
xvals[dbase[["sex"]] > 0] <- floor(dbase[["Yr.S"]][
dbase[["sex"]] > 0
]) -
(dbase[["sex"]][dbase[["sex"]] > 0] - 1.5) * xdiff
}
cols[dbase[["sex"]] > 0] <- colvec[dbase[["sex"]][
dbase[["sex"]] > 0
]]
}
bubble3(
x = xvals,
y = dbase[["Bin"]],
z = z,
xlab = labels[3],
ylab = kindlab,
col = cols,
cexZ1 = cexZ1,
legend = bublegend,
las = 1,
main = ptitle,
cex.main = cex.main,
maxsize = pntscalar,
allopen = allopen,
minnbubble = minnbubble
)
# add lines for growth of individual cohorts if requested
if (length(cohortlines) > 0) {
for (icohort in seq_along(cohortlines)) {
message(" Adding line for", cohortlines[icohort], "cohort\n")
if (kind == "LEN") {
if (nsexes > 1) {
lines(
growdatF[["Age_Mid"]] + cohortlines[icohort],
growdatF[["Len_Mid"]],
col = colvec[1]
) # females
lines(
growdatM[["Age_Mid"]] + cohortlines[icohort],
growdatM[["Len_Mid"]],
col = colvec[2]
) # males
} else {
lines(
growdatF[["Age_Mid"]] + cohortlines[icohort],
growdatF[["Len_Mid"]],
col = colvec[3]
) # single-sex growth
}
}
if (kind %in% c("AGE", "GSTAGE")) {
lines(
0.5 +
c(cohortlines[icohort], cohortlines[icohort] + accuage),
c(0, accuage),
col = colvec[3],
lty = 3
) # one-one line for age
}
}
}
}
if (plot) {
tempfun2()
}
if (print) {
# set up plotting to png file if required
pagetext <- ""
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
if (length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
file <- paste(
filenamestart,
filetype,
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun2()
dev.off() # close device if png
}
} # end bubble plot
### subplot 3: multi-panel bubble plots for conditional age-at-length
if (3 %in% subplots & kind == "cond") {
# assemble caption that may also be used for plot title
caption1 <- paste(titletype, title_sexmkt, fleetnames[f], sep = "")
caption1 <- paste(
caption1,
" (max=",
round(max(z), digits = 2),
")",
sep = ""
)
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
# calculate scaling of lines showing effect and input sample size
sampsizeline.old <- sampsizeline
effNline.old <- effNline
if (is.logical(sampsizeline) && sampsizeline) {
# scaling when displaying only adjusted input sample size
sampsizeline <- max(dbase[["Bin"]]) /
max(dbase[["Nsamp_adj"]], na.rm = TRUE)
if (!datonly && is.logical(effNline) && effNline) {
# scaling when displaying both input and effective
sampsizeline <- effNline <- max(dbase[["Bin"]]) /
max(dbase[["Nsamp_adj"]], dbase[["effN"]], na.rm = TRUE)
message(
" Fleet ",
f,
" ",
titlesex,
"adj. input & effective N in red & green scaled by ",
effNline,
"\n",
sep = ""
)
} else {
message(
" Fleet ",
f,
" ",
titlesex,
"adj. input N in red scaled by ",
sampsizeline,
"\n",
sep = ""
)
}
}
# function to make plots
tempfun3 <- function(ipage, ...) {
sexvec <- dbase[["sex"]]
col.index <- sexvec
col.index[col.index == 0] <- 3
cols <- colvec[col.index]
yrvec <- dbase[["Yr.S"]] + dbase[["sex"]] * 1e-6
make_multifig(
ptsx = dbase[["Bin"]],
ptsy = dbase[["Lbin_mid"]],
yr = yrvec,
size = z,
sampsize = dbase[["Nsamp_adj"]],
showsampsize = showsampsize,
effN = dbase[["effN"]],
showeffN = FALSE,
cexZ1 = cexZ1,
bublegend = bublegend,
nlegends = 1,
legtext = list(dbase[["YrSeasName"]]),
bars = FALSE,
linepos = 0,
main = ptitle,
cex.main = cex.main,
xlab = labels[2],
ylab = labels[1],
ymin0 = FALSE,
maxrows = maxrows2,
maxcols = maxcols2,
fixdims = fixdims,
allopen = allopen,
minnbubble = minnbubble,
# ptscol=col[1],ptscol2=col[2],
ptscol = cols,
ipage = ipage,
scalebins = scalebins,
sampsizeline = sampsizeline,
effNline = effNline,
sampsizemean = MeanNage,
effNmean = HarmEffNage,
colvec = colvec,
linescol = linescol,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
...
)
}
if (plot) {
tempfun3(ipage = 0, ...)
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(
length(unique(dbase[["Yr.S"]])) *
length(unique(dbase[["sex"]])) /
maxrows2 /
maxcols2
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
file <- paste(
filenamestart,
filetype,
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun3(ipage = ipage, ...)
dev.off() # close device if png
}
}
sampsizeline <- sampsizeline.old
effNline <- effNline.old
} # end conditional bubble plot
#
### subplots 4 and 5: multi-panel plot of point and line fit to
### conditional age-at-length
### and Pearson residuals of A-L key for specific years
if (
(4 %in% subplots | 5 %in% subplots) &
aalyear[1] > 0 &
kind == "cond"
) {
for (y in seq_along(aalyear)) {
aalyr <- aalyear[y]
if (length(dbase[["Obs"]][dbase[["Yr"]] == aalyr]) > 0) {
ydbase <- dbase[dbase[["Yr"]] == aalyr, ]
sexvec <- ydbase[["sex"]]
if (4 %in% subplots) {
### subplot 4: multi-panel plot of fit to conditional age-at-length for specific years
caption1 <- paste(
aalyr,
" age-at-length bin, ",
title_sexmkt,
fleetnames[f],
sep = ""
)
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
lenbinlegend <- paste(
ydbase[["Lbin_lo"]],
labels[7],
sep = ""
)
lenbinlegend[ydbase[["Lbin_range"]] > 0] <- paste(
ydbase[["Lbin_lo"]],
"-",
ydbase[["Lbin_hi"]],
labels[7],
sep = ""
)
tempfun4 <- function(ipage, ...) {
# temporary function to aid repeating the big function call
make_multifig(
ptsx = ydbase[["Bin"]],
ptsy = ydbase[["Obs"]],
yr = ydbase[["Lbin_lo"]],
linesx = ydbase[["Bin"]],
linesy = ydbase[["Exp"]],
sampsize = ydbase[["Nsamp_adj"]],
effN = ydbase[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
nlegends = 3,
legtext = list(lenbinlegend, "sampsize", "effN"),
bars = FALSE,
linepos = linepos,
main = ptitle,
cex.main = cex.main,
xlab = labels[2],
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
xlas = xlas,
ylas = ylas,
axis1 = axis1,
axis2 = axis2,
axis1labs = axis1labs,
sexvec = sexvec,
yupper = yupper,
...
)
}
if (plot) {
tempfun4(ipage = 0, ...)
}
if (print) {
npages <- ceiling(
length(unique(ydbase[["Yr.S"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
if (length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
file <- paste0(
filenamestart,
filename_fltsexmkt,
"_",
aalyr,
pagetext,
".png"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun4(ipage = ipage, ...)
dev.off() # close device if print
}
}
} # end if 4 in subplots
if (5 %in% subplots) {
### subplot 5: Pearson residuals for A-L key
z <- ydbase[["Pearson"]]
col.index <- sexvec
col.index[col.index == 0] <- 3
cols <- colvec[col.index]
x.vec <- ydbase[["Bin"]] + ydbase[["sex"]] * 1e-6
caption <- paste(
aalyr,
" Pearson residuals for A-L key, ",
title_sexmkt,
fleetnames[f],
sep = ""
)
caption <- paste(
caption,
" (max=",
round(abs(max(z)), digits = 2),
")",
sep = ""
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun5 <- function() {
bubble3(
x = x.vec,
y = ydbase[["Lbin_lo"]],
z = z,
xlab = labels[2],
ylab = labels[1],
col = cols,
las = 1,
main = ptitle,
cex.main = cex.main,
maxsize = pntscalar,
cexZ1 = cexZ1,
legend = bublegend,
allopen = FALSE,
minnbubble = minnbubble
)
}
if (plot) {
tempfun5()
}
if (print) {
if (length(grep("Pearson", caption)) > 0) {
caption <- paste(
caption,
"<br> \nClosed bubbles are positive residuals",
"(observed > expected)",
"and open bubbles are negative residuals",
"(observed < expected)."
)
}
file <- paste0(
filenamestart,
"yearresids_",
filename_fltsexmkt,
"_",
aalyr,
pagetext,
".png"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun5()
dev.off() # close device if print
}
} # end if 5 in subplots
}
}
}
### subplot 6: multi-panel plot of point and line fit to conditional
### age-at-length for specific length bins
if (6 %in% subplots & aalbin[1] > 0) {
badbins <- setdiff(aalbin, dbase[["Lbin_hi"]])
goodbins <- intersect(aalbin, dbase[["Lbin_hi"]])
if (length(goodbins) > 0) {
if (length(badbins) > 0) {
message(
"Error! the following inputs for 'aalbin' do not match the Lbin_hi values for the conditional age-at-length data:",
badbins,
"\n",
" the following inputs for 'aalbin' are fine:",
goodbins,
"\n"
)
}
for (ibin in seq_along(goodbins)) {
# loop over good bins
ilenbin <- goodbins[ibin]
abindbase <- dbase[dbase[["Lbin_hi"]] == ilenbin, ]
if (nrow(abindbase) > 0) {
# check for data associated with this bin
sexvec <- abindbase[["sex"]]
caption1 <- paste0(
"Age-at-length ",
ilenbin,
labels[7],
", ",
title_sexmkt,
fleetnames[f]
)
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
tempfun6 <- function(ipage, ...) {
# temporary function to aid repeating the big function call
make_multifig(
ptsx = abindbase[["Bin"]],
ptsy = abindbase[["Obs"]],
yr = abindbase[["Yr.S"]],
linesx = abindbase[["Bin"]],
linesy = abindbase[["Exp"]],
sampsize = abindbase[["Nsamp_adj"]],
effN = abindbase[["effN"]],
showsampsize = showsampsize,
showeffN = showeffN,
nlegends = 3,
legtext = list(
abindbase[["YrSeasName"]],
"sampsize",
"effN"
),
bars = bars,
linepos = (1 - datonly) * linepos,
main = ptitle,
cex.main = cex.main,
xlab = kindlab,
ylab = labels[6],
maxrows = maxrows,
maxcols = maxcols,
rows = rows,
cols = cols,
fixdims = fixdims,
ipage = ipage,
scalebins = scalebins,
sexvec = sexvec,
...
)
}
if (plot) {
tempfun6(ipage = 0, ...)
}
if (print) {
npages <- ceiling(
length(unique(abindbase[["Yr.S"]])) / maxrows / maxcols
)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste0("_page", ipage)
caption <- paste0(
caption,
" (plot ",
ipage,
" of ",
npages,
")"
)
}
file <- paste0(
filenamestart,
filename_fltsexmkt,
"_length",
ilenbin,
labels[7],
pagetext,
".png"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
tempfun6(ipage = ipage, ...)
dev.off() # close device if print
}
} # end print
} # end if data
} # end loop over length bins
} # end if length(goodbins)>0
} # end if plot requested
### subplot 7: sample size plot
if (
7 %in%
subplots &
samplesizeplots &
!datonly &
!("Nsamp_DM" %in%
names(dbase) &&
any(!is.na(dbase[["Nsamp_DM"]]))) &
!(kind %in% c("GSTAGE", "GSTLEN", "L@A", "W@A"))
) {
caption <- paste0(
"N-EffN comparison, ",
titledata,
title_sexmkt,
fleetnames[f]
)
if (mainTitle) {
ptitle <- caption
} else {
ptitle <- ""
}
titles <- c(ptitle, titles) # compiling list of all plot titles
lfitfunc <- function() {
if (kind == "cond") {
# trap nonrobust effective n's
# should this only be for conditional age-at-length or all plots?
dbasegood <- dbase[
dbase[["Obs"]] >= 0.0001 &
dbase[["Exp"]] < 0.99 &
!is.na(dbase[["effN"]]) &
dbase[["effN"]] < maxneff,
]
} else {
dbasegood <- dbase
}
if (nrow(dbasegood) > 0) {
# thinning out columns and removing rows with redundant information
# (for the purposes of this function)
dbasegood2 <- dbasegood[, c("YrSeasName", "Nsamp_adj", "effN")]
dbasegood2 <- unique(dbasegood2)
plot(
dbasegood2[["Nsamp_adj"]],
dbasegood2[["effN"]],
xlab = labels[4],
main = ptitle,
cex.main = cex.main,
ylim = c(0, 1.15 * max(dbasegood2[["effN"]])),
xlim = c(0, 1.15 * max(dbasegood2[["Nsamp_adj"]])),
col = colvec[3],
pch = 19,
ylab = labels[5],
xaxs = "i",
yaxs = "i"
)
# add labels for the years if requested
if (showyears) {
par(xpd = TRUE) # allows the label to go over plot boundary
text(
x = dbasegood2[["Nsamp_adj"]],
y = dbasegood2[["effN"]],
dbasegood2[["YrSeasName"]],
adj = c(-0.2, 0.5)
)
par(xpd = FALSE) # restores default clipping
}
abline(0, 1, col = "black", lty = 1)
# add loess smoother if there's at least 6 points with a range greater than 2
if (
smooth &
length(unique(dbasegood2[["Nsamp_adj"]])) > 6 &
diff(range(dbasegood2[["Nsamp_adj"]])) > 2
) {
old_warn <- options()[["warn"]] # previous warnings setting
options(warn = -1) # turn off loess warnings
psmooth <- loess(
dbasegood2[["effN"]] ~ dbasegood2[["Nsamp_adj"]],
degree = 1
)
options(warn = old_warn) # returning to old value
lines(
psmooth[["x"]][order(psmooth[["x"]])],
psmooth[["fit"]][order(psmooth[["x"]])],
lwd = 1.2,
col = "red",
lty = "dashed"
)
}
if (addMeans) {
# vertical line with label for mean input sample size
abline(
v = mean(dbasegood2[["Nsamp_adj"]]),
lty = "22",
col = "green3"
)
text(
x = mean(dbasegood2[["Nsamp_adj"]]),
y = 0,
col = "green3",
"arithmetic mean",
srt = 90,
adj = c(-0.1, -0.3)
)
# horizontal line with label for harmonic effective sample size
abline(
h = 1 / mean(1 / dbasegood2[["effN"]]),
lty = "22",
col = "green3"
)
text(
x = 0,
y = 1 / mean(1 / dbasegood2[["effN"]]),
col = "green3",
"harmonic mean",
adj = c(-0.1, -0.3)
)
}
}
}
if (plot) {
lfitfunc()
}
if (print) {
# set up plotting to png file if required
file <- paste(
filenamestart,
"sampsize_",
filename_fltsexmkt,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
lfitfunc()
dev.off()
}
} # end subplot 7
### subplots 8 and 9: Francis plot for non-conditional data
if (nrow(replist[["Dirichlet_Multinomial_pars"]]) == 0) {
# only include subplots 8 and 10 if using DM likelihood
subplots <- subplots[!subplots %in% c(8, 10)]
}
for (whichplot in 8:9) {
# don't include Francis adjustment in subplot 8, just 9
# and plot only values by partition in subplot 8
if (whichplot == 8) {
plotadj <- FALSE
part <- j
} else {
plotadj <- TRUE
part <- 0:2
}
if (whichplot %in% subplots & kind %in% c("LEN", "SIZE", "AGE")) {
# convert "AGE" to "age" so that SSMethod.TA1.8 can find "agedbase", etc.
kind2 <- tolower(kind)
if (plot) {
tmp <- SSMethod.TA1.8(
fit = replist,
type = kind2,
part = part,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj,
printit = FALSE
)
}
if (print) {
# set up plotting to png file if required
file <- paste0(
filenamestart,
ifelse(
whichplot == 8,
"data_weighting_timeseries_",
"data_weighting_TA1-8_"
),
fleetnames[f],
".png"
)
if (length(part) == 1) {
file <- gsub(
pattern = "timeseries",
replacement = paste0("timeseries_part", part),
x = file
)
}
# not using save_png because caption isn't available until after
# plot is created
# old command: plotinfo <- save_png(file=file, caption=caption)
png(
filename = file.path(plotdir, file),
width = pwidth,
height = pheight,
units = punits,
res = res,
pointsize = ptsize
)
# run function
tmp <- SSMethod.TA1.8(
fit = replist,
type = kind2,
part = part,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj,
printit = FALSE
)
# create caption
caption <- paste0(
"Mean ",
gsub("len", "length", tolower(kind)),
" for ",
fleetnames[f],
" with 95% confidence intervals",
" based on current sample sizes."
)
if (!datonly & whichplot == 8) {
caption <- gsub(
pattern = "sizes.",
replacement = "sizes (including any D-M weighting).",
x = caption
)
}
if (!datonly & whichplot == 9) {
caption <- paste0(
caption,
"<br>Francis data weighting method TA1.8:"
)
if (!is.null(tmp[1])) {
vals <- paste0(
"thinner intervals (with capped ends) show ",
"result of further adjusting sample sizes ",
"based on suggested multiplier ",
"(with 95% interval) for ",
kind2,
" data from ",
fleetnames[f],
":<br>",
round(tmp[1], 4),
" (",
round(tmp[2], 4),
"-",
round(tmp[3], 4),
")"
)
} else {
vals <- "too few points to calculate adjustments."
}
caption <- paste(
caption,
vals,
"<br>For more info, see ",
"<a href=https://doi.org/10.1139/f2011-025>",
"Francis (2011)</a>."
)
} # end additional caption text about Francis weighting
# add caption to the plotinfo table (normally done by save_png)
plotinfo <- rbind(
plotinfo,
data.frame(
file = file,
caption = caption,
alt_text = NA
)
)
dev.off() # close device if png
} # end test for print to PNG option
} # end subplot 8 or 9
} # end loop over 8 and 9
### subplots 10 and 11: Francis weighting for conditional data
for (whichplot in 10:11) {
# don't include Francis adjustment in subplot 10, just 11
plotadj <- (whichplot == 11)
if (
whichplot %in%
subplots &
kind == "cond" &
(f %in% condbase[["Fleet"]])
) {
if (plot) {
SSMethod.Cond.TA1.8(
fit = replist,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj
)
}
if (print) {
# set up plotting to png file if required
file <- paste0(
filenamestart,
ifelse(
whichplot == 10,
"data_weighting_timeseries_condAge",
"data_weighting_TA1-8_condAge"
),
gsub(" ", "", fleetnames[f]),
".png"
)
# not using save_png because caption isn't available until after
# plot is created
png(
filename = file.path(plotdir, file),
width = pwidth,
height = pheight,
units = punits,
res = res,
pointsize = ptsize
)
# run function
tmp <- SSMethod.Cond.TA1.8(
fit = replist,
fleet = f,
fleetnames = fleetnames,
datonly = datonly,
plotadj = plotadj
)
# create caption
caption <- paste0(
"Mean age from conditional data",
" (aggregated across length bins) for ",
fleetnames[f],
" with 95% confidence intervals ",
" based on current samples sizes."
)
if (!datonly & whichplot == 10) {
caption <- gsub(
pattern = "sizes.",
replacement = "sizes (including any D-M weighting).",
x = caption
)
}
if (!datonly & whichplot == 11) {
caption <- paste0(
caption,
"<br>Francis data weighting method TA1.8:"
)
if (!is.null(tmp[1])) {
vals <- paste0(
"thinner intervals (with capped ends) show ",
"result of further adjusting sample sizes ",
"based on suggested multiplier ",
"(with 95% interval) for ",
"conditional age-at-length data from ",
fleetnames[f],
":<br>",
round(tmp[1], 4),
" (",
round(tmp[2], 4),
"-",
round(tmp[3], 4),
")",
sep = ""
)
} else {
vals <- "too few points to calculate adjustments."
}
caption <- paste(
caption,
vals,
"<br>For more info, see ",
"<a href=https://doi.org/10.1139/f2011-025>",
"Francis (2011)</a> and ",
"<a href=https://doi.org/10.1016/j.fishres.2015.12.006>",
"Punt (2017)</a>."
)
} # end test for datonly
# add caption to the plotinfo table (normally done by save_png)
plotinfo <- rbind(
plotinfo,
data.frame(
file = file,
caption = caption,
alt_text = NA
)
)
dev.off() # close device if png
} # end test for print to PNG option
} # end plot 10 or 11
} # end loop over 10 and 11
### subplot 12: Andre's mean age and std. dev. in conditional AAL
if (
12 %in%
subplots &
kind == "cond" &
length(unique(dbase[["Bin"]])) > 1
) {
caption1 <- paste(labels[14], title_sexmkt, fleetnames[f], sep = "")
if (mainTitle) {
ptitle <- caption1
} else {
ptitle <- ""
}
andrefun <- function(ipage = 0) {
Lens <- sort(unique(dbase[["Lbin_lo"]]))
Yrs <- sort(unique(dbase[["Yr.S"]]))
ymax <- 1.1 * max(dbase[["Bin"]], na.rm = TRUE)
xmax <- max(condbase[["Lbin_hi"]], na.rm = TRUE)
xmin <- min(condbase[["Lbin_lo"]], na.rm = TRUE)
# do some stuff so that figures that span multiple pages can be output as separate PNG files
npanels <- length(Yrs)
npages <- npanels / andrerows
panelrange <- 1:npanels
if (npages > 1 & ipage != 0) {
panelrange <- intersect(
panelrange,
1:andrerows + andrerows * (ipage - 1)
)
}
Yrs2 <- Yrs[panelrange]
par(mfrow = c(andrerows, 2), mar = c(2, 4, 1, 1), oma = andre_oma)
for (Yr in Yrs2) {
y <- dbase[dbase[["Yr.S"]] == Yr, ]
Size <- NULL
Size2 <- NULL
Obs <- NULL
Obs2 <- NULL
Pred <- NULL
Pred2 <- NULL
Upp <- NULL
Low <- NULL
Upp2 <- NULL
Low2 <- NULL
for (Ilen in Lens) {
z <- y[y[["Lbin_lo"]] == Ilen, ]
if (length(z[, 1]) > 0) {
weightsPred <- z[["Exp"]] / sum(z[["Exp"]])
weightsObs <- z[["Obs"]] / sum(z[["Obs"]])
ObsV <- sum(z[["Bin"]] * weightsObs)
ObsV2 <- sum(z[["Bin"]] * z[["Bin"]] * weightsObs)
PredV <- sum(z[["Bin"]] * weightsPred)
PredV2 <- sum(z[["Bin"]] * z[["Bin"]] * weightsPred)
# Overdispersion on N
# NN <- z[["Nsamp_adj"]][1]*0.01 # Andre did this for reasons unknown
NN <- z[["Nsamp_adj"]][1]
if (
max(z[["Obs"]], na.rm = TRUE) > 1.0e-4 &
!is.na(NN) &&
NN > 0
) {
Size <- c(Size, Ilen)
Obs <- c(Obs, ObsV)
Pred <- c(Pred, PredV)
varn <- sqrt(PredV2 - PredV * PredV) / sqrt(NN)
Pred2 <- c(Pred2, varn)
varn <- sqrt(max(0, ObsV2 - ObsV * ObsV, na.rm = TRUE)) /
sqrt(NN)
Obs2 <- c(Obs2, varn)
Low <- c(Low, ObsV - 1.64 * varn)
Upp <- c(Upp, ObsV + 1.64 * varn)
if (NN > 1) {
Size2 <- c(Size2, Ilen)
Low2 <- c(
Low2,
varn * sqrt((NN - 1) / qchisq(0.95, NN))
)
Upp2 <- c(
Upp2,
varn * sqrt((NN - 1) / qchisq(0.05, NN))
)
}
}
}
}
if (length(Obs) > 0) {
## next line was replaced with setting at the top,
## for consistency across years
# ymax <- max(Pred,Obs,Upp)*1.1
plot(
Size,
Obs,
type = "n",
xlab = "",
ylab = "Age",
xlim = c(xmin, xmax),
ylim = c(0, ymax),
yaxs = "i"
)
label <- ifelse(nseasons == 1, floor(Yr), Yr)
text(
x = par("usr")[1],
y = .9 * ymax,
labels = label,
adj = c(-.5, 0),
font = 2,
cex = 1.2
)
if (length(Low) > 1) {
polygon(
c(Size, rev(Size)),
c(Low, rev(Upp)),
col = "grey95",
border = NA
)
}
if (!datonly) {
lines(Size, Pred, col = 4, lwd = 3)
}
points(Size, Obs, pch = 16)
lines(Size, Low, lty = 3)
lines(Size, Upp, lty = 3)
if (par("mfg")[1] == 1) {
title(
main = ptitle,
xlab = labels[1],
outer = TRUE,
line = 1
)
}
box()
ymax2 <- max(Obs2, Pred2, na.rm = TRUE) * 1.1
plot(
Size,
Obs2,
type = "n",
xlab = labels[1],
ylab = labels[13],
xlim = c(xmin, xmax),
ylim = c(0, ymax2),
yaxs = "i"
)
if (length(Low2) > 1) {
polygon(
c(Size2, rev(Size2)),
c(Low2, rev(Upp2)),
col = "grey95",
border = NA
)
}
if (!datonly) {
lines(Size, Pred2, col = 4, lwd = 3)
}
points(Size, Obs2, pch = 16)
lines(Size2, Low2, lty = 3)
lines(Size2, Upp2, lty = 3)
if (!datonly & par("mfg")[1] == 1) {
legend(
"topleft",
legend = c("Observed (with 90% interval)", "Expected"),
bty = "n",
col = c(1, 4),
pch = c(16, NA),
lty = c(NA, 1),
lwd = 3
)
}
box()
} # end if data exist
} # end loop over years
} # end andrefun
if (plot) {
andrefun()
}
if (print) {
# set up plotting to png file if required
npages <- ceiling(length(unique(dbase[["Yr.S"]])) / andrerows)
for (ipage in 1:npages) {
pagetext <- ""
caption <- caption1
if (npages > 1) {
pagetext <- paste("_page", ipage, sep = "")
caption <- paste(
caption,
" (plot ",
ipage,
" of ",
npages,
")",
sep = ""
)
}
if (ipage == 1) {
# add more information only to first page of plots
caption <- paste(
caption,
"\nThese plots show mean age and std. dev. in conditional A@L.<br>",
"Left plots are mean A@L by size-class (obs. and exp.) ",
"with 90% CIs based on adding 1.64 SE of mean to the data.<br>",
"Right plots in each pair are SE of mean A@L (obs. and exp.) ",
"with 90% CIs based on the chi-square distribution.",
sep = ""
)
}
file <- paste(
filenamestart,
"Andre_plots",
filename_fltsexmkt,
pagetext,
".png",
sep = ""
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
andrefun(ipage = ipage)
dev.off() # close device if png
} # end loop over pages
} # end test for print to PNG option
} # end subplot 12
} # end loop over partitions (index j)
# } # end test for whether sex in vector of requested sexes
# } # end loop over combined/not-combined sex
} # end if data
} # end loop over fleets
if (!is.null(plotinfo)) {
plotinfo[["category"]] <- "Comp"
}
return(invisible(plotinfo))
} # end SSplotComps function
###########################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.