R/SSplotBiology.R
In r4ss/r4ss: R Code for Stock Synthesis
Defines functions
SSplotBiology
Documented in
SSplotBiology
#' Plot biology related quantities.
#'
#' Plot biology related quantities from Stock Synthesis model output, including
#' mean weight, maturity, fecundity, and spawning output.
#'
#' @template replist
#' @template plot
#' @template print
#' @param add add to existing plot
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows:
#' \itemize{
#' \item 1 growth curve only
#' \item 2 growth curve with CV and SD
#' \item 3 growth curve with maturity and weight
#' \item 4 distribution of length at age (still in development)
#' \item 5 length or wtatage matrix
#' \item 6 maturity
#' \item 7 fecundity from model parameters
#' \item 8 fecundity at weight from BIOLOGY section
#' \item 9 fecundity at length from BIOLOGY section
#' \item 10 spawning output at length
#' \item 11 spawning output at age
#' \item 21 Natural mortality (if age-dependent)
#' \item 22 Time-varying growth persp
#' \item 23 Time-varying growth contour
#' \item 24 plot time-series of any time-varying quantities (created if the
#' MGparm_By_Year_after_adjustments table (report:7) is available in the
#' Report.sso file)
#' \item 31 hermaphroditism transition probability
#' \item 32 sex ratio in ending year
#' (only plotted when model has hermaphroditism)
#' }
#' Additional plots not created by default
#' \itemize{
#' \item 101 diagram with labels showing female growth curve
#' \item 102 diagram with labels showing female growth curve & male offsets
#' \item 103 diagram with labels showing female CV = f(A) (offset type 2)
#' \item 104 diagram with labels showing female CV = f(A) & male offset (type 2)
#' \item 105 diagram with labels showing female CV = f(A) (offset type 3)
#' \item 106 diagram with labels showing female CV = f(A) & male offset (type 3)
#' }
#' @param seas which season to plot (values other than 1 only work in
#' seasonal models but but maybe not fully implemented)
#' @param morphs Which morphs to plot (if more than 1 per sex)? By default this
#' will be `replist[["mainmorphs"]]`
#' @param forecast Include forecast years in plots of time-varying biology?
#' @param minyr optional input for minimum year to show in plots
#' @param maxyr optional input for maximum year to show in plots
#' @param colvec vector of length 3 with colors for various points/lines
#' @template areacols
#' @param ltyvec vector of length 2 with lty for females/males in growth plots
#' values can be applied to other plots in the future
#' @param shadealpha Transparency parameter used to make default shadecol
#' values (see ?rgb for more info)
#' @template legendloc
#' @template plotdir
#' @template labels
#' @template pwidth
#' @template pheight
#' @template pheight_tall
#' @template punits
#' @template res
#' @template ptsize
#' @template cex.main
#' @param imageplot_text Whether to add numerical text to the image plots
#' when using weight at age. Defaults to FALSE.
#' @param imageplot_text_round The number of significant digits to which
#' the image plot text is rounded. Defaults to 0, meaning whole numbers. If
#' all your values are small and there's no contrast in the text, you might
#' want to make this 1 or 2.
#' @template mainTitle
#' @template verbose
#' @author Ian Stewart, Ian Taylor
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotBiology <-
function(
replist,
plot = TRUE,
print = FALSE,
add = FALSE,
subplots = 1:32,
seas = 1,
morphs = NULL,
forecast = FALSE,
minyr = -Inf,
maxyr = Inf,
colvec = c("red", "blue", "grey20"),
areacols = NULL,
ltyvec = c(1, 2),
shadealpha = 0.1,
imageplot_text = FALSE,
imageplot_text_round = 0,
legendloc = "topleft",
plotdir = "default",
labels = c(
"Length (cm)", # 1
"Age (yr)", # 2
"Maturity", # 3
"Mean weight (kg) in last year", # 4
replist[["SpawnOutputLabel"]], # 5
"Length (cm, beginning of the year)", # 6
"Natural mortality", # 7
"Female weight (kg)", # 8
"Female length (cm)", # 9
"Fecundity", # 10
"Default fecundity label", # 11
"Year", # 12
"Hermaphroditism transition rate", # 13
"Fraction females by age in ending year"
), # 14
pwidth = 6.5,
pheight = 5.0,
pheight_tall = 6.5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = TRUE,
verbose = TRUE
) {
#### current (Aug 18, 2017) order of plots:
# subplot 1: growth_curve_fn - growth curve only
# subplot 2: growth_curve_plus_fn - growth curve with CV and SD
# subplot 3: growth_curve_plus_fn - growth curve with maturity and weight
# subplot 4: SSplotAgeMatrix (separate function) - distribution of length at age
# subplot 5: weight-length or wtatage matrix
# subplot 6: maturity
# subplot 7: fec_pars_fn - fecundity from model parameters
# subplot 8: fec_weight_fn - fecundity at weight from BIOLOGY section
# subplot 9: fec_len_fn - fecundity at length from BIOLOGY section
# subplot 10: spawn_output_len_fn - spawning output at length
# subplot 11: spawn_output_age_fn - spawning output at age
# subplot 21 mfunc - Natural mortality (if age-dependent)
# subplot 22 [no function] - Time-varying growth persp
# subplot 23 [no function] - Time-varying growth contour
# subplot 24 timeVaryingParmFunc - plot time-series of any time-varying quantities
#### unfinished addition on 3-Mar-16
# subplot 25 [no function] - matrix of M by age and time
#### vectors related to hermaphroditism
# subplot 31 Herma_Trans
# subplot 32 Herma_Cum
# EXTRA PLOTS NOT PRODUCED BY DEFAULT
# subplot 101: diagram with labels showing female growth curve
# subplot 102: diagram with labels showing female growth curve & male offsets
# subplot 103: diagram with labels showing female CV = f(A) (offset type 2)
# subplot 104: diagram with labels showing female CV = f(A) & male offset (type 2)
# subplot 105: diagram with labels showing female CV = f(A) (offset type 3)
# subplot 106: diagram with labels showing female CV = f(A) & male offset (type 3)
# test for presence of wtatage_switch
wtatage_switch <- replist[["wtatage_switch"]]
if (wtatage_switch) {
message(
"Note: this model uses the empirical weight-at-age input.\n",
" Plots of many quantities related to growth are skipped."
)
} else {
if (is.null(replist[["endgrowth"]])) {
message(
"Skipping biology plots because model output doesn't include\n",
" biology-at-age information (endgrowth), likely because\n",
" brief output was specified in starter.ss."
)
return()
}
}
plotinfo <- NULL
ians_blues <- c(
"white",
"grey",
"lightblue",
"skyblue",
"steelblue1",
"slateblue",
topo.colors(6),
"blue",
"blue2",
"blue3",
"blue4",
"black"
)
ians_contour <- c("white", rep("blue", 100))
# convert colvec to semi-transparent colors for shading polygons
shadecolvec <- rep(NA, length(colvec))
for (icol in seq_along(colvec)) {
tmp <- col2rgb(colvec[icol]) / 255
shadecolvec[icol] <- rgb(
red = tmp[1],
green = tmp[2],
blue = tmp[3],
alpha = shadealpha
)
}
#### plot function 1
# mean weight, maturity, fecundity, spawning output
# get objects from replist
nseasons <- replist[["nseasons"]]
nareas <- replist[["nareas"]]
growdat <- replist[["endgrowth"]][
replist[["endgrowth"]][["Seas"]] == seas,
]
growthCVtype <- replist[["growthCVtype"]]
biology <- replist[["biology"]]
startyr <- replist[["startyr"]]
FecType <- replist[["FecType"]]
FecPar1name <- replist[["FecPar1name"]]
FecPar2name <- replist[["FecPar2name"]]
FecPar1 <- replist[["FecPar1"]]
FecPar2 <- replist[["FecPar2"]]
parameters <- replist[["parameters"]]
nsexes <- replist[["nsexes"]]
accuage <- replist[["accuage"]]
startyr <- replist[["startyr"]]
endyr <- replist[["endyr"]]
growthvaries <- replist[["growthvaries"]]
growthseries <- replist[["growthseries"]]
ageselex <- replist[["ageselex"]]
MGparmAdj <- replist[["MGparmAdj"]]
wtatage <- replist[["wtatage"]]
# remove comments on each line of wtatage that were added with SSv3.30
if ("comment" %in% names(wtatage)) {
wtatage <- wtatage[, -grep("comment", names(wtatage))]
}
M_at_age <- replist[["M_at_age"]]
Growth_Parameters <- replist[["Growth_Parameters"]]
if (is.null(morphs)) {
morphs <- replist[["mainmorphs"]]
if (is.null(morphs)) {
# if morphs not supplied but growth is available, use middle
# morph within each sex
morphs <- median(growdat[["Morph"]][growdat[["Sex"]] == 1])
if (nsexes == 2) {
morphs <- c(morphs, median(growdat[["Morph"]][growdat[["Sex"]] == 2]))
}
}
}
# get any derived quantities related to growth curve uncertainty
Grow_std <- replist[["derived_quants"]][
grep("Grow_std_", replist[["derived_quants"]][["Label"]]),
]
if (nrow(Grow_std) == 0) {
Grow_std <- NULL
} else {
# convert things like "Grow_std_1_Fem_A_25" into
# in more recent 3.30.14 versions, the label appears as
# "Grow_std_GP:_1_Fem_A_25"
# so the "shift" below = 0 or 1 to adjust the position accordingly
# "pattern 1, female, age 25"
Grow_std[["pattern"]] <- NA
Grow_std[["sex_char"]] <- NA
Grow_std[["sex"]] <- NA
Grow_std[["age"]] <- NA
shift <- length(grep("GP:", Grow_std[["Label"]][1]))
for (irow in 1:nrow(Grow_std)) {
tmp <- strsplit(Grow_std[["Label"]][irow], split = "_")[[1]]
Grow_std[["pattern"]][irow] <- as.numeric(tmp[3 + shift])
Grow_std[["sex_char"]][irow] <- tmp[4 + shift]
Grow_std[["age"]][irow] <- as.numeric(tmp[6 + shift])
}
Grow_std[["sex"]][Grow_std[["sex_char"]] == "Fem"] <- 1
Grow_std[["sex"]][Grow_std[["sex_char"]] == "Mal"] <- 2
#### now it should look something like this:
## Label Value StdDev (Val-1.0)/Stddev
## Grow_std_GP:_1_Fem_A_0 Grow_std_GP:_1_Fem_A_0 27.3188 0.393286 NA
## Grow_std_GP:_1_Fem_A_5 Grow_std_GP:_1_Fem_A_5 50.6099 0.178201 NA
## CumNorm pattern sex_char sex age
## Grow_std_GP:_1_Fem_A_0 NA 1 Fem 1 0
## Grow_std_GP:_1_Fem_A_5 NA 1 Fem 1 5
}
# get any derived quantities related to M uncertainty
NatM_std <- replist[["derived_quants"]][
grep("NatM_std_", replist[["derived_quants"]][["Label"]]),
]
if (nrow(NatM_std) == 0) {
NatM_std <- NULL
} else {
# convert things like "NatM_std_GP:_1_Fem_A_1" into
# "pattern 1, female, age 25"
NatM_std[["pattern"]] <- NA
NatM_std[["sex_char"]] <- NA
NatM_std[["sex"]] <- NA
NatM_std[["age"]] <- NA
for (irow in 1:nrow(NatM_std)) {
tmp <- strsplit(NatM_std[["Label"]][irow], split = "_")[[1]]
NatM_std[["pattern"]][irow] <- as.numeric(tmp[4])
NatM_std[["sex_char"]][irow] <- tmp[5]
NatM_std[["age"]][irow] <- as.numeric(tmp[7])
}
NatM_std[["sex"]][NatM_std[["sex_char"]] == "Fem"] <- 1
NatM_std[["sex"]][NatM_std[["sex_char"]] == "Mal"] <- 2
#### now it should look something like this:
## Label Value StdDev (Val-1.0)/Stddev
## NatM_std_GP:_1_Fem_A_1 NatM_std_GP:_1_Fem_A_1 0.564563 0 NA
## NatM_std_GP:_1_Fem_A_2 NatM_std_GP:_1_Fem_A_2 0.510574 0 NA
## CumNorm pattern sex_char sex age
## NatM_std_GP:_1_Fem_A_1 NA 1 Fem 1 1
## NatM_std_GP:_1_Fem_A_2 NA 1 Fem 1 2
}
if (!seas %in% 1:nseasons) {
stop("'seas' input should be within 1:nseasons")
}
if (nseasons > 1) {
labels[6] <- gsub(
"beginning of the year",
paste("beginning of season", seas),
labels[6]
)
}
if (plotdir == "default") {
plotdir <- replist[["inputs"]][["dir"]]
}
# check dimensions
if (length(morphs) > nsexes) {
warning(
"!Error with morph indexing in SSplotBiology function.\n",
" Code is not set up to handle multiple growth patterns or birth seasons.\n"
)
}
## # stuff from selectivity that is not used
## FecundAtAge <- ageselex[ageselex[["factor"]]=="Fecund", names(ageselex)%in%0:accuage]
## WtAtAge <- ageselex[ageselex[["factor"]]=="bodywt", names(ageselex)%in%0:accuage]
# determine fecundity type
# define labels and x-variable
if (FecType == 1) {
fec_ylab <- "Eggs per kg"
fec_xlab <- labels[8]
FecX <- biology[["Wt_F"]]
FecY <- FecPar1 + FecPar2 * FecX
}
if (labels[11] != "Default fecundity label") {
fec_ylab <- labels[11]
}
# Beginning of season 1 (or specified season) mean length at age
# with 95% range of lengths (by sex if applicable)
## Ian T.: consider somehow generalizing to allow looping over growth pattern
if (!is.null(growdat)) {
# calculate CVs from SD and mean length
growdat[["CV_Beg"]] <- growdat[["SD_Beg"]] / growdat[["Len_Beg"]]
# female growth
growdatF <- growdat[
growdat[["Sex"]] == 1 &
growdat[["Morph"]] == morphs[1],
]
growdatF[["Sd_Size"]] <- growdatF[["SD_Beg"]]
if (growthCVtype == "logSD=f(A)") {
# lognormal distribution of length at age
growdatF[["high"]] <- qlnorm(
0.975,
meanlog = log(growdatF[["Len_Beg"]]),
sdlog = growdatF[["Sd_Size"]]
)
growdatF[["low"]] <- qlnorm(
0.025,
meanlog = log(growdatF[["Len_Beg"]]),
sdlog = growdatF[["Sd_Size"]]
)
} else {
# normal distribution of length at age
growdatF[["high"]] <- qnorm(
0.975,
mean = growdatF[["Len_Beg"]],
sd = growdatF[["Sd_Size"]]
)
growdatF[["low"]] <- qnorm(
0.025,
mean = growdatF[["Len_Beg"]],
sd = growdatF[["Sd_Size"]]
)
}
if (nsexes > 1) {
# do males if 2-sex model
growdatM <- growdat[
growdat[["Sex"]] == 2 & growdat[["Morph"]] == morphs[2],
]
# IAN T. this should probably be generalized
xm <- growdatM[["Age_Beg"]]
growdatM[["Sd_Size"]] <- growdatM[["SD_Beg"]]
if (growthCVtype == "logSD=f(A)") {
# lognormal distribution of length at age
growdatM[["high"]] <- qlnorm(
0.975,
meanlog = log(growdatM[["Len_Beg"]]),
sdlog = growdatM[["Sd_Size"]]
)
growdatM[["low"]] <- qlnorm(
0.025,
meanlog = log(growdatM[["Len_Beg"]]),
sdlog = growdatM[["Sd_Size"]]
)
} else {
# normal distribution of length at age
growdatM[["high"]] <- qnorm(
0.975,
mean = growdatM[["Len_Beg"]],
sd = growdatM[["Sd_Size"]]
)
growdatM[["low"]] <- qnorm(
0.025,
mean = growdatM[["Len_Beg"]],
sd = growdatM[["Sd_Size"]]
)
}
}
}
clean_wtatage <- function(x) {
## quick function to check for valid wtatage matrix and remove first
## redundant row if it's there. Used in weight_plot and maturity_plot.
if (nrow(x) < 2) {
message("Not enough rows in weight-at-age matrix to plot")
return(NULL)
}
if (all(x[1, ] == x[2, ])) {
x <- x[-1, ]
}
return(x)
}
weight_plot <- function(sex) {
# weight
## This needs to be a function of sex since it can be called
## either once for a single sex model or twice to produce plots for
## each one.
if (!wtatage_switch) {
# if empirical weight-at-age is not used
x <- biology[["Len_mean"]]
if (!add) {
ymax <- max(biology[["Wt_F"]])
if (nsexes > 1) {
ymax <- max(ymax, biology[["Wt_M"]])
}
plot(
x,
x,
ylim = c(0, 1.1 * ymax),
xlab = labels[1],
ylab = labels[4],
type = "n",
las = 1,
yaxs = "i"
)
}
lines(x, biology[["Wt_F"]], type = "o", col = colvec[1])
if (nsexes > 1) {
lines(x, biology[["Wt_M"]], type = "o", col = colvec[2])
if (!add) {
legend(
legendloc,
bty = "n",
c("Females", "Males"),
lty = 1,
col = c(colvec[1], colvec[2])
)
}
}
} else {
## if empirical weight-at-age IS used
# calculate maximum weight across all sexes
wtmax <- wtatage[
wtatage[["fleet"]] == -1 &
wtatage[["seas"]] == seas,
-(1:6)
] |>
max()
# get matrix (including initial year column)
wtmat <- wtatage[
wtatage[["fleet"]] == -1 &
wtatage[["sex"]] == sex &
wtatage[["seas"]] == seas,
-(2:6)
]
# check matrix
wtmat <- clean_wtatage(wtmat)
# plot it
if (!is.null(wtmat)) {
makeimage(wtmat, main = "", lastbin = wtmax)
}
}
}
maturity_plot <- function() {
# maturity
if (!wtatage_switch) {
# if empirical weight-at-age is not used
x <- biology[["Len_mean"]]
if (min(biology[["Mat"]]) < 1) {
# if length based
if (!add) {
plot(
x,
biology[["Mat"]],
xlab = labels[1],
ylab = labels[3],
las = 1,
yaxs = "i",
ylim = c(0, max(biology[["Mat"]])),
type = "o",
col = colvec[1]
)
}
if (add) lines(x, biology[["Mat"]], type = "o", col = colvec[1])
} else {
# else is age based
if (!add) {
plot(
growdatF[["Age_Beg"]],
growdatF[["Age_Mat"]],
xlab = labels[2],
las = 1,
yaxs = "i",
ylim = c(0, 1.1 * max(growdatF[["Age_Mat"]])),
ylab = labels[3],
type = "o",
col = colvec[1]
)
} else {
lines(
growdatF[["Age_Beg"]],
growdatF[["Age_Mat"]],
type = "o",
col = colvec[1]
)
}
}
} else {
# if empirical weight-at-age IS used
fecmat <- wtatage[wtatage[["fleet"]] == -2 & wtatage[["sex"]] == 1, ]
if (nrow(fecmat) > 1) {
# figure out which seasons have fecundity values (maybe always only one?)
fecseasons <- sort(unique(fecmat[["seas"]]))
seas_label <- NULL
for (iseas in fecseasons) {
fecmat_seas <- fecmat[fecmat[["seas"]] == iseas, -(2:6)]
# label the season only if a multi-season model
# also testing for length of fecseasons, but that's probably redundant
if (nseasons > 1 | length(fecseasons) > 1) {
seas_label <- paste("in season", iseas)
}
main <- paste("", seas_label)
fecmat_seas <- clean_wtatage(fecmat_seas)
## persp(x=abs(fecmat_seas[,1]),
## y=0:accuage,
## z=as.matrix(fecmat_seas[,-1]),
## theta=70,phi=30,xlab="Year",ylab="Age",zlab="",
## main=main)
}
makeimage(fecmat_seas, main = "")
}
}
}
makeimage <- function(mat, main = "", lastbin = NULL) {
yrvec <- abs(mat[["year"]])
##### this stuff was used to add a row of mean values
## if(is.null(meanvec)){
## meanvec <- mat[,1]
## mat <- mat[,-1]
## }
## yrvec2 <- c(-2:-1+min(yrvec), yrvec)
## mat2 <- cbind(meanvec,NA,mat)
# subset based on minyr and maxyr
subset <- yrvec >= minyr & yrvec <= maxyr
yrvec2 <- yrvec[subset]
mat2 <- mat[subset, -1]
if (is.null(lastbin)) {
lastbin <- max(mat2)
}
z <- t(mat2)
breaks <- seq(0, lastbin, length = 51)
col <- rainbow(60)[1:50]
## Make a two panel plot, for plot and then legend
layout(
matrix(c(1, 2), nrow = 1, ncol = 2),
widths = c(4, .5),
heights = c(4, 4)
)
# save current parameter settings
par_old <- par()
# change parameters
par(mar = c(4.2, 4.2, 1, .5) + .1)
image(
x = 0:accuage,
y = yrvec2,
z = z,
axes = F,
xlab = "Age",
ylab = "Year",
col = col,
breaks = breaks,
main = ifelse(mainTitle, main, "")
)
axis(1, cex.axis = .7)
axis(2, las = 1, cex.axis = .7)
box()
## add text if desired
if (imageplot_text) {
zdataframe <- expand.grid(age = 0:accuage, yr = yrvec2)
zdataframe[["z"]] <- c(t(mat2))
zdataframe[["font"]] <- 1
## Rounding should depend on how big the numbers get. This is untested
## and may need to be adjusted.
ztext <- format(round(zdataframe[["z"]], imageplot_text_round))
ztext[ztext == " NA" | ztext == " NA"] <- ""
text(
x = zdataframe[["age"]],
y = zdataframe[["Yr"]],
label = ztext,
font = zdataframe[["font"]],
cex = .7
)
}
## add a legend to the image plot
## Code adapated from online function. Taken on 11/10/2015 from:
## http://menugget.blogspot.com/2011/08/adding-scale-to-image-plot.html#more
## ------------------------------------------------------------
zlim <- range(z, na.rm = TRUE)
zlim[2] <- zlim[2] + c(zlim[2] - zlim[1]) * (1E-3) # adds a bit to the range in both directions
zlim[1] <- zlim[1] - c(zlim[2] - zlim[1]) * (1E-3)
poly <- vector(mode = "list", length(col))
for (i in seq(poly)) {
poly[[i]] <- c(breaks[i], breaks[i + 1], breaks[i + 1], breaks[i])
}
ylim <- range(breaks)
xlim <- c(0, 1)
par(mar = c(4.3, .5, 1, 3))
plot(
1,
1,
type = "n",
ylim = ylim,
xlim = xlim,
axes = FALSE,
ann = FALSE,
xaxs = "i",
yaxs = "i"
)
for (i in seq(poly)) {
polygon(c(0, 0, 1, 1), poly[[i]], col = col[i], border = NA)
}
axis(4, cex.axis = .7, las = 2)
box()
## turn off the matrix plot
layout(mat = c(1, 1))
# restore default single panel settings
par(mfcol = c(1, 1), mar = par_old[["mar"]], oma = par_old[["oma"]])
}
fec_pars_fn <- function() {
# fecundity from model parameters
ymax <- 1.1 * max(FecY)
if (!add) {
plot(
FecX,
FecY,
xlab = fec_xlab,
ylab = fec_ylab,
ylim = c(0, ymax),
las = 1,
yaxs = "i",
col = colvec[2],
pch = 19
)
lines(FecX, rep(FecPar1, length(FecX)), col = colvec[1])
text(
mean(range(FecX)),
FecPar1 - 0.05 * ymax,
"Egg output proportional to spawning biomass"
)
} else {
points(FecX, FecY, col = colvec[2], pch = 19)
}
}
fecundityOK <- all(!is.na(biology[["Fec"]]))
fec_weight_fn <- function() {
# fecundity at weight from BIOLOGY section
ymax <- 1.1 * max(biology[["Fec"]])
if (!add) {
plot(
biology[["Wt_F"]],
biology[["Fec"]],
xlab = labels[8],
ylab = labels[10],
las = 1,
yaxs = "i",
ylim = c(0, ymax),
col = colvec[1],
type = "o"
)
} else {
points(
biology[["Len_mean"]],
biology[["Fec"]],
col = colvec[1],
type = "o"
)
}
}
fec_len_fn <- function() {
# fecundity at length from BIOLOGY section
ymax <- 1.1 * max(biology[["Fec"]])
if (!add) {
plot(
biology[["Len_mean"]],
biology[["Fec"]],
xlab = labels[9],
ylab = labels[10],
las = 1,
yaxs = "i",
ylim = c(0, 1.1 * ymax),
col = colvec[1],
type = "o",
yaxs = "i"
)
} else {
points(
biology[["Len_mean"]],
biology[["Fec"]],
col = colvec[1],
type = "o"
)
}
}
spawn_output_len_fn <- function() {
# spawning output at length
x <- biology[["Len_mean"]]
y <- biology[["Mat*Fec"]]
ymax <- 1.1 * max(y)
if (!add) {
plot(
x,
y,
xlab = labels[1],
ylab = labels[5],
type = "o",
col = colvec[1],
las = 1,
yaxs = "i",
ylim = c(0, ymax)
)
} else {
lines(x, y, type = "o", col = colvec[1])
}
}
spawn_output_age_fn <- function() {
# spawning output at age
x <- growdat[["Age_Beg"]][growdat[["Sex"]] == 1]
y <- growdat$"Mat*Fecund"[growdat[["Sex"]] == 1]
ymax <- 1.1 * max(y)
if (!add) {
plot(
x,
y,
xlab = labels[2],
ylab = labels[5],
type = "o",
col = colvec[1],
las = 1,
yaxs = "i",
ylim = c(0, ymax)
)
} else {
lines(x, y, type = "o", col = colvec[1])
}
}
if (!wtatage_switch) {
ymax <- max(biology[["Len_mean"]])
x <- growdatF[["Age_Beg"]]
}
main <- "Ending year expected growth (with 95% intervals)"
# if(nseasons > 1){main <- paste(main," season 1",sep="")}
col_index1 <- 3 # default is grey for single-sex model
if (nsexes > 1) {
col_index1 <- 1 # change line for females to red
}
add_uncertainty_fn <- function(std_table, sex, age_offset = 0.5) {
# function to add uncertainty estimates to mean length at age or M at age
# previous settings for warnings
old_warn <- options()[["warn"]]
# turn off "zero-length arrow" warning
options(warn = -1)
# confirm presence of table of values and request to add to plot
if (!is.null(std_table)) {
std_table.sex <- std_table[std_table[["sex"]] == sex, ]
# confirm presence for requested sex
if (!is.null(std_table.sex)) {
for (irow in 1:nrow(std_table.sex)) {
std.age <- std_table.sex[["age"]][irow]
# this might not work for log-normal length at age
mean <- std_table.sex[["Value"]][irow]
ages <- std_table.sex[["age"]][irow] + age_offset
high <- qnorm(
0.975,
mean = mean,
sd = std_table.sex[["StdDev"]][irow]
)
low <- qnorm(
0.025,
mean = mean,
sd = std_table.sex[["StdDev"]][irow]
)
# set sex-dependent color
if (sex == 1) {
col <- colvec[col_index1]
}
if (sex == 2) {
col <- colvec[2]
}
arrows(
x0 = ages,
x1 = ages,
y0 = low,
y1 = high,
length = 0.04,
angle = 90,
code = 3,
col = col
)
}
} # end check for missing sex
} # end check for empty table
# returning to old warnings (after turning off in case of zero-length arrow)
options(warn = old_warn)
}
growth_curve_fn <- function(add_labels = TRUE, add_uncertainty = TRUE) {
# growth
x <- growdatF[["Age_Beg"]]
# make empty plot unless this is being added to existing figure
if (!add) {
plot(
x,
growdatF[["Len_Beg"]],
ylim = c(0, 1.1 * ymax),
type = "n",
xlab = "",
ylab = "",
axes = FALSE,
xaxs = "i",
yaxs = "i"
)
abline(h = 0, col = "grey")
if (add_labels) {
title(
main = ifelse(mainTitle, main, ""),
xlab = labels[2],
ylab = labels[6],
cex.main = cex.main
)
axis(1)
axis(2, las = 1)
}
}
lty <- 1
# make shaded polygon
polygon(
c(x, rev(x)),
c(growdatF[["low"]], rev(growdatF[["high"]])),
border = NA,
col = shadecolvec[col_index1]
)
# add lines for mean, low and high
lines(
x,
growdatF[["Len_Beg"]],
col = colvec[col_index1],
lwd = 2,
lty = ltyvec[1]
)
lines(
x,
growdatF[["high"]],
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
lines(x, growdatF[["low"]], col = colvec[col_index1], lwd = 1, lty = "12")
# add 95% intervals for uncertainty in growth from $derived_quants
if (add_uncertainty) {
add_uncertainty_fn(std_table = Grow_std, sex = 1, age_offset = 0.5)
}
# add males if they are present in the model
if (nsexes > 1) {
# make shaded polygon
polygon(
c(xm, rev(xm)),
c(growdatM[["low"]], rev(growdatM[["high"]])),
border = NA,
col = shadecolvec[2]
)
# add lines for mean, low and high
lines(
xm,
growdatM[["Len_Beg"]],
col = colvec[2],
lwd = 2,
lty = ltyvec[2]
)
lines(xm, growdatM[["high"]], col = colvec[2], lwd = 1, lty = "13")
lines(xm, growdatM[["low"]], col = colvec[2], lwd = 1, lty = "13")
# add 95% intervals for uncertainty in growth from $derived_quants
if (add_uncertainty) {
add_uncertainty_fn(std_table = Grow_std, sex = 2, age_offset = 0.5)
}
}
if (!add) {
grid()
box()
}
if (nsexes > 1) {
legend(
legendloc,
bty = "n",
c("Females", "Males"),
lty = ltyvec,
lwd = 2,
col = c(colvec[1], colvec[2])
)
}
}
if (plot & 1 %in% subplots & !wtatage_switch) {
growth_curve_fn()
}
if (print & 1 %in% subplots & !wtatage_switch) {
file <- "bio1_sizeatage.png"
caption <- paste(
"Length at age in the beginning of the year (or season) in the ending",
"year of the model. Shaded area indicates 95% distribution of",
"length at age around estimated growth curve."
)
if (!is.null(Grow_std)) {
caption <- paste(
caption,
"Vertical intervals around growth curve indicate",
"estimated 95% uncertainty intervals in estimated mean growth."
)
}
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_fn()
dev.off()
}
growth_curve_plus_fn <- function(add_labels = TRUE, option = 1) {
# function to add panels to growth curve with info on variability in
# length at age or info on maturity, weight, and fecundity
# save current parameter settings
par_old <- par()
# create 4-panel layout
layout(
mat = matrix(1:4, byrow = TRUE, ncol = 2),
widths = c(2, 1),
heights = c(2, 1)
)
par(mar = c(1, 1, 1, 1), oma = c(5, 5, 5, 4))
# plot growth curve
growth_curve_fn(add_labels = FALSE)
axis(2, las = 1)
# add age axis on top
axis(3)
# add some labels
mtext(side = 2, line = 2.5, labels[6], las = 0)
mtext(side = 3, line = 2.5, labels[2], las = 0)
box()
if (option == 1) {
lab1 <- "SD_Beg"
lab1long <- "SD of lengths"
lab2 <- "CV_Beg"
lab2long <- "CV of lengths"
lab1max <- 1.1 * max(growdat[[lab1]], na.rm = TRUE)
lab2max <- 1.05 * max(growdat[[lab2]], na.rm = TRUE)
lab1_axis_vec <- NULL
}
if (option == 2) {
lab1 <- "Len_Mat"
lab1long <- "Frac. mature"
lab2 <- "Wt_Beg"
lab2long <- "Mean weight"
lab1max <- 1
lab2max <- max(c(biology[["Wt_F"]], biology[["Wt_M"]]), na.rm = TRUE)
lab1_axis_vec <- c(0, 0.5, 1)
}
# calculate scaling factor between CVs and SDs to share each panel
lab2_to_lab1_scale <- lab1max / lab2max
lab2_axis_vec <- pretty(c(0, lab2max))
# make empty panel in top-right location
plot(
0,
type = "n",
xaxs = "i",
xlim = c(0, 1.0 * lab1max),
yaxs = "i",
ylim = par()[["usr"]][3:4],
axes = FALSE
)
if (option == 1) {
# if plotting SD and CV, then get it from age-based data
# add line for lab1 vs. Len
lines(
growdatF[[lab1]],
growdatF[["Len_Beg"]],
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
# add line for lab2 vs. Len
lines(
growdatF[[lab2]] * lab2_to_lab1_scale,
growdatF[["Len_Beg"]],
col = colvec[col_index1],
lwd = 3
)
if (nsexes > 1) {
# add lines for males
if (option == 1) {
# add line for lab1 vs. Age
# (unless it's maturity, which is not defined for males)
lines(
growdatM[[lab1]],
growdatM[["Len_Beg"]],
col = colvec[2],
lwd = 1,
lty = "13"
)
}
# add line for lab2 vs. Len
lines(
growdatM[[lab2]] * lab2_to_lab1_scale,
growdatM[["Len_Beg"]],
col = colvec[2],
lwd = 3,
lty = 2
)
}
}
if (option == 2) {
# if plotting maturity and fecundity, then get this panel from length-based data
lines(
biology[["Wt_F"]] * lab2_to_lab1_scale,
biology[["Len_mean"]],
col = colvec[col_index1],
lwd = 3
)
lines(
biology[["Mat"]],
biology[["Len_mean"]],
col = colvec[col_index1],
lty = "12"
)
if (nsexes > 1) {
lines(
biology[["Wt_M"]] * lab2_to_lab1_scale,
biology[["Len_mean"]],
col = colvec[2],
lwd = 3,
lty = 2
)
}
}
# add axes and labels
mtext(side = 4, line = 2.5, labels[6], las = 0)
mtext(side = 1, line = 2.5, lab1long, las = 0)
mtext(side = 3, line = 2.5, lab2long, las = 0)
axis(1, at = lab1_axis_vec, las = 1)
axis(
3,
at = lab2_axis_vec * lab2_to_lab1_scale,
labels = lab2_axis_vec,
las = 1
)
axis(4, las = 1)
box()
# make empty plot in lower-left position
plot(
0,
type = "n",
xaxs = "i",
xlim = range(growdat[["Age_Beg"]]),
yaxs = "i",
ylim = c(0, 1.0 * lab1max),
axes = FALSE
)
# add line for lab1 vs. Age
if (option == 1) {
lines(
growdatF[["Age_Beg"]],
growdatF[[lab1]],
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
}
if (option == 2) {
# maturity curve is product of age-based and length-based factors
lines(
growdatF[["Age_Beg"]],
growdatF[["Len_Mat"]] * abs(growdatF[["Age_Mat"]]),
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
}
# add line for lab2 vs. Age
lines(
growdatF[["Age_Beg"]],
growdatF[[lab2]] * lab2_to_lab1_scale,
col = colvec[col_index1],
lwd = 3
)
if (nsexes > 1) {
# add lines for males
if (option == 1) {
# add line for lab1 vs. Age (unless it's maturity, which is not defined for males)
lines(
growdatM[["Age_Beg"]],
growdatM[[lab1]],
col = colvec[2],
lwd = 1,
lty = "13"
)
}
# add line for lab2 vs. Age
lines(
growdatM[["Age_Beg"]],
growdatM[[lab2]] * lab2_to_lab1_scale,
col = colvec[2],
lwd = 3,
lty = 2
)
}
# add axes and labels
mtext(side = 1, line = 2.5, labels[2], las = 0)
mtext(side = 4, line = 2.5, lab1long, las = 0)
mtext(side = 2, line = 2.5, lab2long, las = 0)
axis(1, las = 1)
axis(
2,
at = lab2_axis_vec * lab2_to_lab1_scale,
labels = lab2_axis_vec,
las = 1
)
axis(4, at = lab1_axis_vec, las = 1)
box()
# restore default single panel settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
} # end growth_curve_plus_fn()
# make plots of growth curve with CV and SD of length
if (plot & 2 %in% subplots & !wtatage_switch) {
growth_curve_plus_fn(option = 1)
}
if (print & 2 %in% subplots & !wtatage_switch) {
file <- "bio2_sizeatage_plus_CV_and_SD.png"
caption <- paste(
"Length at age (top-left panel) with",
"CV (thick line) and SD (thin line) of",
"length at age shown in top-right and lower-left panels"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_plus_fn(option = 1)
dev.off()
}
# make plots of growth curve with weight-length curve and maturity
if (plot & 3 %in% subplots & !wtatage_switch) {
growth_curve_plus_fn(option = 2)
}
if (print & 3 %in% subplots & !wtatage_switch) {
file <- "bio3_sizeatage_plus_WT_and_MAT.png"
caption <- paste(
"Length at age (top-left panel) with",
"weight (thick line) and maturity (thin line)",
"shown in top-right and lower-left panels"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_plus_fn(option = 2)
dev.off()
}
# plot distribution of length at age (by season, sub-season, and morph)
if (4 %in% subplots & !wtatage_switch) {
if (!is.null(replist[["ALK"]])) {
plotinfo.tmp <- SSplotAgeMatrix(
replist = replist,
option = 1,
plot = plot,
print = print,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
cex.main = cex.main,
mainTitle = mainTitle
)
# SSplotAgeMatrix adds a "category" column which isn't present
# in plotinfo until the end of this SSplotBiology function.
plotinfo.tmp <- plotinfo.tmp[, c("file", "caption", "alt_text")]
plotinfo <- rbind(plotinfo, plotinfo.tmp)
} else {
message(
"Skipped some plots because AGE_LENGTH_KEY unavailable in report file\n",
"because starter file set to produce limited report detail."
)
}
}
# function for illustrating parameterization of growth curves
growth_curve_labeled_fn <- function(option = 1) {
# growth
if (is.null(Growth_Parameters)) {
message(
"Need updated SS_output function to get Growth_Parameters output\n"
)
return()
}
# save current parameter settings
par_old <- par()
par(mar = c(4, 4, 1, 4))
AminF <- Growth_Parameters[["A1"]][1]
AmaxF <- Growth_Parameters[["A2"]][1]
AminM <- Growth_Parameters[["A1"]][2]
AmaxM <- Growth_Parameters[["A2"]][2]
L_at_AminF <- Growth_Parameters[["L_a_A1"]][1]
L_at_AmaxF <- Growth_Parameters[["L_a_A2"]][1]
L_at_AminM <- Growth_Parameters[["L_a_A1"]][2]
L_at_AmaxM <- Growth_Parameters[["L_a_A2"]][2]
LinfF <- Growth_Parameters[["Linf"]][1]
LinfM <- Growth_Parameters[["Linf"]][2]
ymax <- max(biology[["Len_mean"]])
plot(
0,
type = "n",
xlim = c(0, 1 + max(growdatF[["Age_Beg"]])),
ylim = c(0, 1.1 * ymax),
xlab = "",
ylab = "",
axes = FALSE,
xaxs = "i",
yaxs = "i"
)
abline(h = 0, col = "grey")
# title(main=main, xlab=labels[2], ylab=labels[6], cex.main=cex.main)
axis(
1,
at = c(0, AminF, AmaxF, accuage),
labels = expression(0, italic(A[1]), italic(A[2]), italic(N[ages]))
)
axis(
2,
at = c(0, L_at_AminF, L_at_AmaxF, LinfF),
labels = expression(
0,
italic(L[A[1]]),
italic(L[A[2]]),
italic(L[infinity])
),
las = 1
)
# add growth curve itself
lines(
growdatF[["Age_Beg"]],
growdatF[["Len_Beg"]],
col = 1,
lwd = 3,
lty = 1
)
# add growth curve for males
if (nsexes > 1 & option == 2) {
lines(
growdatM[["Age_Beg"]],
growdatM[["Len_Beg"]],
col = 1,
lwd = 2,
lty = 2
)
}
if (option == 1) {
# lines connecting axes to curve
lines(c(AminF, AminF, 0), c(0, L_at_AminF, L_at_AminF), lty = 3)
lines(c(AmaxF, AmaxF, 0), c(0, L_at_AmaxF, L_at_AmaxF), lty = 3)
abline(h = LinfF, lty = 3)
}
if (option == 2) {
# lines connecting two curves
lines(c(0, AminF), c(L_at_AminF, L_at_AminF), lty = 3)
lines(c(0, AmaxF), c(L_at_AmaxF, L_at_AmaxF), lty = 3)
lines(c(accuage + 10, AmaxM), c(L_at_AmaxM, L_at_AmaxM), lty = 3)
lines(c(accuage + 10, AminM), c(L_at_AminM, L_at_AminM), lty = 3)
arrows(
x0 = AminF,
x1 = AminF,
y0 = L_at_AminF,
y1 = L_at_AminM,
lwd = 3,
col = 3,
length = 0.1
)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = L_at_AmaxF,
y1 = L_at_AmaxM,
lwd = 3,
col = 3,
length = 0.1
)
text(
AminF,
mean(c(L_at_AminF, L_at_AminM)),
"Male parameter 1\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
text(
AmaxF,
mean(c(L_at_AmaxF, L_at_AmaxM)),
"Male parameter 2\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
axis(
4,
at = c(0, L_at_AminM, L_at_AmaxM),
labels = expression(0, italic(L[A[1]]^male), italic(L[A[2]]^male)),
las = 1
)
abline(h = LinfF, lty = 3)
}
box()
# restore old parameter settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
}
# plots of diagrams for illustrating parameterization
if (plot & 101 %in% subplots & !wtatage_switch) {
growth_curve_labeled_fn(option = 1)
}
if (print & 101 %in% subplots & !wtatage_switch) {
file <- "bio101_growth_illustration.png"
caption <- "Illustration of growth parameters"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_labeled_fn(option = 1)
dev.off()
}
if (plot & 102 %in% subplots & !wtatage_switch) {
growth_curve_labeled_fn(option = 2)
}
if (print & 102 %in% subplots & !wtatage_switch) {
file <- "bio102_growth_illustration2.png"
caption <- "Illustration of growth parameters with male offsets"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_labeled_fn(option = 2)
dev.off()
}
# function for illustrating parameterization of CVs around growth curves
CV_values_labeled_fn <- function(option = 1) {
# growth
if (is.null(Growth_Parameters)) {
message(
"Need updated SS_output function to get Growth_Parameters output"
)
return()
}
# save current parameter settings
par_old <- par()
par(mar = c(4, 4, 1, 4))
AminF <- Growth_Parameters[["A1"]][1]
AmaxF <- Growth_Parameters[["A2"]][1]
AminM <- Growth_Parameters[["A1"]][2]
AmaxM <- Growth_Parameters[["A2"]][2]
CV_at_AminF <- Growth_Parameters[["CVmin"]][1]
CV_at_AmaxF <- Growth_Parameters[["CVmax"]][1]
CV_at_AminM <- Growth_Parameters[["CVmin"]][2]
CV_at_AmaxM <- Growth_Parameters[["CVmax"]][2]
ymax <- max(CV_at_AminF, CV_at_AmaxF, CV_at_AminM, CV_at_AmaxM)
plot(
0,
type = "n",
xlim = c(0, 1 + max(growdatF[["Age_Beg"]])),
ylim = c(0, 1.1 * ymax),
xlab = "",
ylab = "",
axes = FALSE,
xaxs = "i",
yaxs = "i"
)
abline(h = 0, col = "grey")
# title(main=main, xlab=labels[2], ylab=labels[6], cex.main=cex.main)
axis(
1,
at = c(0, AminF, AmaxF, accuage),
labels = expression(0, italic(A[1]), italic(A[2]), italic(N[ages]))
)
axis(
2,
at = c(0, CV_at_AminF, CV_at_AmaxF),
labels = expression(0, italic(CV[A[1]]), italic(CV[A[2]])),
las = 1
)
# add growth curve itself
lines(
growdatF[["Age_Beg"]],
growdatF[["CV_Beg"]],
col = 1,
lwd = 3,
lty = 1
)
# add growth curve for males
if (nsexes > 1 & option %in% c(2, 4)) {
lines(
growdatM[["Age_Beg"]],
growdatM[["CV_Beg"]],
col = 1,
lwd = 2,
lty = 2
)
}
if (option == 1) {
# lines connecting axes to curve
lines(c(AminF, AminF, 0), c(0, CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(AmaxF, AmaxF, 0), c(0, CV_at_AmaxF, CV_at_AmaxF), lty = 3)
}
if (option == 2) {
# lines connecting two curves
lines(c(0, AminF), c(CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(0, AmaxF), c(CV_at_AmaxF, CV_at_AmaxF), lty = 3)
lines(c(accuage + 10, AmaxM), c(CV_at_AmaxM, CV_at_AmaxM), lty = 3)
lines(c(accuage + 10, AminM), c(CV_at_AminM, CV_at_AminM), lty = 3)
arrows(
x0 = AminF,
x1 = AminF,
y0 = CV_at_AminF,
y1 = CV_at_AminM,
lwd = 3,
col = 3,
length = 0.1
)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AmaxF,
y1 = CV_at_AmaxM,
lwd = 3,
col = 3,
length = 0.1
)
text(
AminF,
mean(c(CV_at_AminM)),
"Male parameter 1\n(exponential offset)",
col = 3,
adj = c(0, 1.5)
)
text(
AmaxF,
mean(c(CV_at_AmaxM)),
"Male parameter 2\n(exponential offset)",
col = 3,
adj = c(0, 1.5)
)
axis(
4,
at = c(0, CV_at_AminM, CV_at_AmaxM),
labels = expression(0, italic(CV[A[1]]^male), italic(CV[A[2]]^male)),
las = 1
)
}
if (option == 3) {
# lines connecting axes to curve
lines(c(AminF, AminF, 0), c(0, CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(AmaxF, AmaxF, 0), c(0, CV_at_AmaxF, CV_at_AmaxF), lty = 3)
lines(c(0, AmaxF), c(CV_at_AminF, CV_at_AminF), lty = 3)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AminF,
y1 = CV_at_AmaxF,
lwd = 3,
col = 3,
length = 0.1
)
text(
AmaxF,
mean(c(CV_at_AminF, CV_at_AmaxF)),
"Female parameter 2\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
}
if (option == 4) {
# lines connecting two curves
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AminF,
y1 = CV_at_AmaxF,
lwd = 3,
col = 3,
length = 0.1
)
text(
AmaxF,
mean(c(CV_at_AminF, CV_at_AmaxF)),
"Female parameter 2\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
lines(c(0, AmaxF), c(CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(0, AmaxF), c(CV_at_AmaxF, CV_at_AmaxF), lty = 3)
lines(c(accuage + 10, AmaxM), c(CV_at_AmaxM, CV_at_AmaxM), lty = 3)
lines(c(accuage + 10, AminM), c(CV_at_AminM, CV_at_AminM), lty = 3)
arrows(
x0 = AminF,
x1 = AminF,
y0 = CV_at_AminF,
y1 = CV_at_AminM,
lwd = 3,
col = 3,
length = 0.1
)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AminM,
y1 = CV_at_AmaxM,
lwd = 3,
col = 3,
length = 0.1
)
text(
AminF,
mean(c(CV_at_AminF, CV_at_AminM)),
"Male parameter 1\n(exponential offset)",
col = 3,
adj = c(0, 2)
)
text(
AmaxF,
mean(c(CV_at_AminM, CV_at_AmaxM)),
"Male parameter 2\n(exponential offset)",
col = 3,
adj = c(0, 2)
)
axis(
4,
at = c(0, CV_at_AminM, CV_at_AmaxM),
labels = expression(0, italic(CV[A[1]]^male), italic(CV[A[2]]^male)),
las = 1
)
}
box()
# restore old parameter settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
}
# extra plots illustrating parameterization of CV in growth
if (plot & 103 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 1)
}
if (print & 103 %in% subplots & !wtatage_switch) {
file <- "bio103_CV_illustration.png"
caption <- "Illustration of growth variability parameters"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 1)
dev.off()
}
if (plot & 104 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 2)
}
if (print & 104 %in% subplots & !wtatage_switch) {
file <- "bio104_CV_illustration2.png"
caption <- "Illustration of growth variability parameters with male offsets"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 2)
dev.off()
}
if (plot & 105 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 3)
}
if (print & 105 %in% subplots & !wtatage_switch) {
file <- "bio105_CV_illustration.png"
caption <- "Illustration of growth variability parameters for offset type 3"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 3)
dev.off()
}
if (plot & 106 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 4)
}
if (print & 106 %in% subplots & !wtatage_switch) {
file <- "bio106_CV_illustration2.png"
caption <- "Illustration of growth variability parameters with male offsets"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 4)
dev.off()
}
x <- biology[["Len_mean"]]
## NOTE: weight plots are now a special case since they are broken down
## by whether the model is 1-sex or 2-sex. In the latter two separate
## plots need to be made.
if (plot) {
# plot to screen or to PDF file
if (5 %in% subplots & !is.null(biology)) {
weight_plot(sex = 1)
if (nsexes == 2) {
weight_plot(sex = 2)
}
}
if (6 %in% subplots) {
maturity_plot()
}
if (!wtatage_switch) {
if (7 %in% subplots & FecType == 1) {
fec_pars_fn()
}
if (8 %in% subplots & fecundityOK) {
fec_weight_fn()
}
if (9 %in% subplots & fecundityOK) {
fec_len_fn()
}
if (10 %in% subplots) {
spawn_output_len_fn()
}
if (11 %in% subplots) {
spawn_output_age_fn()
}
}
}
if (print) {
# print to PNG files
if (5 %in% subplots) {
file <- "bio5_weightatsize.png"
caption <- "Weight-length relationship"
if (wtatage_switch) {
caption <- "Weight-at-age"
if (nsexes == 1) {
caption <- paste(caption, "for females")
}
}
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = ifelse(wtatage_switch, pheight_tall, pheight), # weight-at-age matrix works better when taller
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
weight_plot(sex = 1)
dev.off()
# add separate plot for males (only for empirical weight-at-age models)
if (wtatage_switch & nsexes == 2) {
file <- "bio5_weightatsize2.png"
caption <- "Weight-at-age for males"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight_tall,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
weight_plot(sex = 2)
dev.off()
}
}
if (6 %in% subplots) {
file <- "bio6_maturity.png"
if (!wtatage_switch) {
caption <- paste(
"Maturity at",
ifelse(min(biology[["Mat"]]) < 1, "length", "age")
)
} else {
caption <- "Spawning output at age (maturity x fecundity)"
}
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = ifelse(wtatage_switch, pheight_tall, pheight), # spawning-output-at-age matrix works better when taller
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
maturity_plot()
dev.off()
}
if (!wtatage_switch) {
if (7 %in% subplots & FecType == 1) {
file <- "bio7_fecundity.png"
caption <- "Fecundity"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
fec_pars_fn()
dev.off()
}
if (8 %in% subplots & fecundityOK) {
file <- "bio8_fecundity_wt.png"
caption <- "Fecundity as a function of weight"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
fec_weight_fn()
dev.off()
}
if (9 %in% subplots & fecundityOK) {
file <- "bio9_fecundity_len.png"
caption <- "Fecundity as a function of length"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
fec_len_fn()
dev.off()
}
if (10 %in% subplots) {
file <- "bio10_spawningoutput_len.png"
caption <- paste(
"Spawning output at length.",
"This is the product of maturity and fecundity",
"unless maturity is age-based, in which case only",
"fecundity is represented."
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
spawn_output_len_fn()
dev.off()
}
if (11 %in% subplots) {
file <- "bio11_spawningoutput_age.png"
caption <- paste(
"Spawning output at age.",
"This is the product of maturity and fecundity.",
"When these processes are length-based they are",
"converted into the age dimension using the matrix",
"of length at age."
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
spawn_output_age_fn()
dev.off()
}
} # end check for wtatage model
} # end check for print=TRUE (printing to PNG files)
# Natural mortality (if age-dependent -- need to add time-varying M plot)
MatAge <- growdatF[["M"]] # female mortality in the ending year
if (nsexes > 1) {
growdatM <- growdat[growdat[["Morph"]] == morphs[2], ]
MatAge_M <- growdatM[["M"]]
}
if (min(MatAge) != max(MatAge) & 21 %in% subplots) {
ymax <- 1.1 * max(MatAge)
if (nsexes > 1) {
ymax <- 1.1 * max(MatAge, MatAge_M)
}
# function to plut natural mortality
mfunc <- function(add_uncertainty = TRUE) {
if (!add) {
plot(
growdatF[["Age_Beg"]],
MatAge,
col = colvec[col_index1],
lwd = 2,
ylim = c(0, ymax),
yaxs = "i",
type = "n",
ylab = labels[7],
xlab = labels[2],
las = 1
)
}
lines(
growdatF[["Age_Beg"]],
MatAge,
col = colvec[col_index1],
lwd = 2,
type = "o"
)
# add uncertainty in M-at-age (if available from derived_quants)
if (add_uncertainty) {
add_uncertainty_fn(std_table = NatM_std, sex = 1, age_offset = 0.0)
}
if (nsexes > 1) {
growdatM <- growdat[growdat[["Morph"]] == morphs[2], ]
lines(
growdatM[["Age_Beg"]],
growdatM[["M"]],
lty = ltyvec[2],
col = colvec[2],
lwd = 2,
type = "o"
)
legend(
"bottomleft",
bty = "n",
c("Females", "Males"),
lty = ltyvec,
lwd = 2,
col = c(colvec[col_index1], colvec[2])
)
# add uncertainty in M-at-age (if available from derived_quants)
if (add_uncertainty) {
add_uncertainty_fn(std_table = NatM_std, sex = 2, age_offset = 0.0)
}
}
}
# run function if requested to make plot
if (plot & 21 %in% subplots) {
mfunc()
}
# run function if requested to write figure to PNG file
if (print & 21 %in% subplots) {
file <- "bio21_natmort.png"
caption <- "Natural mortality at age"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
mfunc()
dev.off()
}
}
# Time-varying growth
if (is.null(growthvaries)) {
if (verbose) {
message(
"No check for time-varying growth because starter file set to produce\n",
"limited report detail."
)
}
} else {
# temporarily disable multi-season plotting of time-varying growth
if (is.null(growthseries)) {
warning(
"No time-varying growth info because starter file set to produce\n",
"limited report detail."
)
} else {
# if growth is time varying and weight-at-age not used
if (growthvaries & !wtatage_switch) {
for (i in 1:nsexes) {
growdatuse <- growthseries[
growthseries[["Yr"]] >= startyr - 2 &
growthseries[["Morph"]] == morphs[i],
]
# trim based on forecast=TRUE/FALSE
if (!forecast) {
growdatuse <- growdatuse[growdatuse[["Yr"]] <= endyr, ]
}
# trim based on minyr and maxyr
growdatuse <- growdatuse[
growdatuse[["Yr"]] >= minyr &
growdatuse[["Yr"]] <= maxyr,
]
# assemble vectors and matrix
x <- 0:accuage
y <- growdatuse[["Yr"]]
z <- as.matrix(growdatuse[, -(1:4)])
# check for time-varying growth
if (replist[["growthvaries"]]) {
z <- t(z)
if (i == 1) {
main <- "Female time-varying growth"
}
if (nsexes == 1) {
main <- "Time-varying growth"
}
if (i == 2) {
main <- "Male time-varying growth"
}
if (nseasons > 1) {
main <- paste0(main, " season 1")
}
if (plot) {
if (22 %in% subplots) {
persp(
x,
y,
z,
col = "white",
xlab = labels[2],
ylab = "",
zlab = labels[1],
expand = 0.5,
box = TRUE,
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
ticktype = "detailed",
phi = 35,
theta = -10
)
}
if (23 %in% subplots) {
contour(
x,
y,
z,
nlevels = 12,
xlab = labels[2],
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
col = ians_contour,
lwd = 2
)
}
}
if (print) {
if (22 %in% subplots) {
file <- paste0("bio22_timevarygrowthsurf_sex", i, ".png")
caption <- "Perspective plot of time-varying growth"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
persp(
x,
y,
z,
col = "white",
xlab = labels[2],
ylab = "",
zlab = labels[1],
expand = 0.5,
box = TRUE,
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
ticktype = "detailed",
phi = 35,
theta = -10
)
dev.off()
}
if (23 %in% subplots) {
file <- paste0("bio23_timevarygrowthcontour_sex", i, ".png")
caption <- "Contour plot of time-varying growth"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
contour(
x,
y,
z,
nlevels = 12,
xlab = labels[2],
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
col = ians_contour,
lwd = 2
)
dev.off()
}
} # end print
} # end if time-varying
} # end loop over sexes
} # end if growth is time varying
} # end of if data available for time varying growth
} # end disable of time-varying growth for multi-season models
# plot time-series of any time-varying quantities
if (24 %in% subplots) {
if (!is.null(MGparmAdj)) {
# general function to work for any parameter
timeVaryingParmFunc <- function(parmlabel, forecast = FALSE) {
if (forecast) {
MGparmAdj.tmp <- MGparmAdj
} else {
MGparmAdj.tmp <- MGparmAdj[MGparmAdj[["Yr"]] <= endyr, ]
}
# trim based on minyr and maxyr
MGparmAdj.tmp <- MGparmAdj.tmp[
MGparmAdj.tmp[["Yr"]] >= minyr &
MGparmAdj.tmp[["Yr"]] <= maxyr,
]
# make plot
plot(
MGparmAdj.tmp[["Yr"]],
MGparmAdj.tmp[[parmlabel]],
xlab = labels[12],
ylab = parmlabel,
type = "l",
lwd = 3,
col = colvec[2]
)
}
# check to make sure MGparmAdj looks as expected
# (maybe had different or conditional format in old SS versions)
if (!is.null(ncol(MGparmAdj)) && ncol(MGparmAdj) > 1) {
# loop over columns looking for time-varying parameters
for (icol in 2:ncol(MGparmAdj)) {
parmlabel <- names(MGparmAdj)[icol]
# exclude column indicating change added with version 3.30.06.02
if (parmlabel != "Change?") {
parmvals <- MGparmAdj[, icol]
# check for changes
if (length(unique(parmvals[MGparmAdj[["Yr"]] <= endyr])) > 1) {
# make plot
if (plot) {
timeVaryingParmFunc(parmlabel, forecast = forecast)
}
if (print) {
file <- paste0("bio24_time-varying_", parmlabel, ".png")
# replace % sign which cause problems for filename
file <- gsub(
pattern = "%",
replacement = "percent",
x = file,
fixed = TRUE
)
# replace : which which cause problems for filename
file <- gsub(
pattern = ":",
replacement = "_",
x = file,
fixed = TRUE
)
caption <- "Time-varying mortality and growth parameters"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
timeVaryingParmFunc(parmlabel, forecast = forecast)
dev.off()
}
}
}
}
}
} else {
message(
"Skipping timevarying quantity plots (subplot 24), most likely\n",
"because the MGparm_By_Year_after_adjustments table (report:7)\n",
"is not reported in the Report.sso file."
)
}
}
# plot of vectors found in models with hermaphroditism
if ("Herma_Trans" %in% names(growdatF)) {
# Hermaphro_Option from SS3 User Manual
# 1 = invoke female to male age-specific function; and
# -1 = invoke male to female age-specific function.
if (any(!is.na(growdatF[["Herma_Trans"]]))) {
Hermaphro_Option <- 1
labels[13] <- paste(labels[13], "(F to M)")
} else {
Hermaphro_Option <- -1
labels[13] <- paste(labels[13], "(M to F)")
}
herma_func1 <- function() {
# if females change to males, there will be numerical values here:
if (Hermaphro_Option == 1) {
df <- growdatF
} else {
df <- growdatM
}
# make plot (age vector should be the same for females and males)
plot(
df[["Age_Beg"]],
df[["Herma_Trans"]],
xaxs = "i",
ylim = c(0, 1),
las = 1,
xlab = labels[12],
ylab = labels[13],
type = "l",
lwd = 3,
col = colvec[2]
)
abline(h = c(0, 1), col = "grey")
}
herma_func2 <- function(area = 1) {
# new area-specific columns in 3.30.21 have format "sex_ratio_area:1"
# assuming that there will always be one column per area
colnames <- grep("sex_ratio", names(growdatF), value = TRUE)
matplot(
growdatF[["Age_Beg"]],
growdatF[, colnames],
xaxs = "i",
ylim = c(0, 1),
las = 1,
lty = 1:nareas,
xlab = labels[2],
ylab = labels[14],
type = "l",
lwd = 3,
col = areacols
)
# add grey line at 1.0
abline(h = c(0, 1), col = "grey")
# add legend if multiple lines per area
if (length(colnames) > 1) {
legend(
legendloc,
bty = "n",
col = areacols,
lty = 1:nareas,
lwd = 3,
legend = paste("area", 1:nareas)
)
}
}
# set default colors if not specified
areacols <- get_areacols(areacols, nareas)
if (31 %in% subplots) {
if (plot) {
herma_func1()
}
if (print) {
plotinfo <- save_png(
plotinfo = plotinfo,
file = "bio31_hermaphrodite_transition.png",
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = labels[13]
)
herma_func1()
dev.off()
}
}
if (32 %in% subplots) {
if (plot) {
herma_func2()
}
if (print) {
plotinfo <- save_png(
plotinfo = plotinfo,
file = "bio32_hermaphrodite_cumulative.png",
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = labels[14]
)
herma_func2()
dev.off()
}
}
}
## ### Ian T.: finish this stuff
## # Matrix of M-at-age if mortality varies with age and time
## # strip off forecast years in M-at-age matrix
## M_at_age <- M_at_age[M_at_age[["Year"]] <= endyr,]
## if(any(M_at_age[["Bio_Pattern"]]!=1)){}
## for(i in unique(M_at_age[["Bio_Pattern"]])){
## for(isex in unique(
## M_at_age_tmpYear <- M_at_age[["Year"]]
## M_at_age_ages <- 0:accuage
## M_at_age_matrix <- as.matrix(M_at_age[,-
#### this stuff can be canibalized to get the M-at-age matrix working
## if(growthvaries){ # if growth is time varying
## for(i in 1:nsexes)
## {
## growdatuse <- growthseries[growthseries[["Yr"]] >= startyr-2 &
## growthseries[["Morph"]]==morphs[i],]
## x <- 0:accuage
## y <- growdatuse[["Yr"]]
## z <- as.matrix(growdatuse[,-(1:4)])
## time <- FALSE
## for(t in 1:ncol(z)) if(max(z[,t])!=min(z[,t])) time <- TRUE
## if(time)
## {
## z <- t(z)
## if(i==1){main <- "Female time-varying growth"}
## if(nsexes==1){main <- "Time-varying growth"}
## if(i==2){main <- "Male time-varying growth"}
## if(nseasons > 1){main <- paste(main," season 1",sep="")}
## if(plot){
## if(12 %in% subplots)
## persp(x,y,z,col="white",xlab=labels[2],ylab="",zlab=labels[1],expand=0.5,
## box=TRUE,main=main,cex.main=cex.main,ticktype="detailed",
## phi=35,theta=-10)
## if(13 %in% subplots)
## contour(x,y,z,nlevels=12,xlab=labels[2],
## main=main,cex.main=cex.main,col=ians_contour,lwd=2)}
## if(print){
## if(12 %in% subplots){
## file <- paste("bio12_timevarygrowthsurf_sex",i,".png",sep="")
## caption <- "Perspective plot of time-varying growth"
## plotinfo <- save_png(file=file, caption=caption)
## persp(x,y,z,col="white",xlab=labels[2],ylab="",zlab=labels[1],expand=0.5,
## box=TRUE,main=main,cex.main=cex.main,ticktype="detailed",
## phi=35,theta=-10)
## dev.off()
## }
## if(13 %in% subplots){
## file <- paste("bio13_timevarygrowthcontour_sex",i,".png",sep="")
## caption <- "Contour plot of time-varying growth"
## plotinfo <- save_png(file=file, caption=caption)
## contour(x,y,z,nlevels=12,xlab=labels[2],
## main=main,cex.main=cex.main,col=ians_contour,lwd=2)
## dev.off()
## }
## } # end print
## } # end if time-varying
## } # end loop over sexes
## } # end of if data available for time varying growth
## }# end disable of time-varying growth for multi-season models
# add category and return plotinfo
if (!is.null(plotinfo)) {
plotinfo[["category"]] <- "Bio"
}
return(invisible(plotinfo))
}
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 SSplotBiology
Documented in SSplotBiology
#' Plot biology related quantities.
#'
#' Plot biology related quantities from Stock Synthesis model output, including
#' mean weight, maturity, fecundity, and spawning output.
#'
#' @template replist
#' @template plot
#' @template print
#' @param add add to existing plot
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows:
#' \itemize{
#' \item 1 growth curve only
#' \item 2 growth curve with CV and SD
#' \item 3 growth curve with maturity and weight
#' \item 4 distribution of length at age (still in development)
#' \item 5 length or wtatage matrix
#' \item 6 maturity
#' \item 7 fecundity from model parameters
#' \item 8 fecundity at weight from BIOLOGY section
#' \item 9 fecundity at length from BIOLOGY section
#' \item 10 spawning output at length
#' \item 11 spawning output at age
#' \item 21 Natural mortality (if age-dependent)
#' \item 22 Time-varying growth persp
#' \item 23 Time-varying growth contour
#' \item 24 plot time-series of any time-varying quantities (created if the
#' MGparm_By_Year_after_adjustments table (report:7) is available in the
#' Report.sso file)
#' \item 31 hermaphroditism transition probability
#' \item 32 sex ratio in ending year
#' (only plotted when model has hermaphroditism)
#' }
#' Additional plots not created by default
#' \itemize{
#' \item 101 diagram with labels showing female growth curve
#' \item 102 diagram with labels showing female growth curve & male offsets
#' \item 103 diagram with labels showing female CV = f(A) (offset type 2)
#' \item 104 diagram with labels showing female CV = f(A) & male offset (type 2)
#' \item 105 diagram with labels showing female CV = f(A) (offset type 3)
#' \item 106 diagram with labels showing female CV = f(A) & male offset (type 3)
#' }
#' @param seas which season to plot (values other than 1 only work in
#' seasonal models but but maybe not fully implemented)
#' @param morphs Which morphs to plot (if more than 1 per sex)? By default this
#' will be `replist[["mainmorphs"]]`
#' @param forecast Include forecast years in plots of time-varying biology?
#' @param minyr optional input for minimum year to show in plots
#' @param maxyr optional input for maximum year to show in plots
#' @param colvec vector of length 3 with colors for various points/lines
#' @template areacols
#' @param ltyvec vector of length 2 with lty for females/males in growth plots
#' values can be applied to other plots in the future
#' @param shadealpha Transparency parameter used to make default shadecol
#' values (see ?rgb for more info)
#' @template legendloc
#' @template plotdir
#' @template labels
#' @template pwidth
#' @template pheight
#' @template pheight_tall
#' @template punits
#' @template res
#' @template ptsize
#' @template cex.main
#' @param imageplot_text Whether to add numerical text to the image plots
#' when using weight at age. Defaults to FALSE.
#' @param imageplot_text_round The number of significant digits to which
#' the image plot text is rounded. Defaults to 0, meaning whole numbers. If
#' all your values are small and there's no contrast in the text, you might
#' want to make this 1 or 2.
#' @template mainTitle
#' @template verbose
#' @author Ian Stewart, Ian Taylor
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotBiology <-
function(
replist,
plot = TRUE,
print = FALSE,
add = FALSE,
subplots = 1:32,
seas = 1,
morphs = NULL,
forecast = FALSE,
minyr = -Inf,
maxyr = Inf,
colvec = c("red", "blue", "grey20"),
areacols = NULL,
ltyvec = c(1, 2),
shadealpha = 0.1,
imageplot_text = FALSE,
imageplot_text_round = 0,
legendloc = "topleft",
plotdir = "default",
labels = c(
"Length (cm)", # 1
"Age (yr)", # 2
"Maturity", # 3
"Mean weight (kg) in last year", # 4
replist[["SpawnOutputLabel"]], # 5
"Length (cm, beginning of the year)", # 6
"Natural mortality", # 7
"Female weight (kg)", # 8
"Female length (cm)", # 9
"Fecundity", # 10
"Default fecundity label", # 11
"Year", # 12
"Hermaphroditism transition rate", # 13
"Fraction females by age in ending year"
), # 14
pwidth = 6.5,
pheight = 5.0,
pheight_tall = 6.5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = TRUE,
verbose = TRUE
) {
#### current (Aug 18, 2017) order of plots:
# subplot 1: growth_curve_fn - growth curve only
# subplot 2: growth_curve_plus_fn - growth curve with CV and SD
# subplot 3: growth_curve_plus_fn - growth curve with maturity and weight
# subplot 4: SSplotAgeMatrix (separate function) - distribution of length at age
# subplot 5: weight-length or wtatage matrix
# subplot 6: maturity
# subplot 7: fec_pars_fn - fecundity from model parameters
# subplot 8: fec_weight_fn - fecundity at weight from BIOLOGY section
# subplot 9: fec_len_fn - fecundity at length from BIOLOGY section
# subplot 10: spawn_output_len_fn - spawning output at length
# subplot 11: spawn_output_age_fn - spawning output at age
# subplot 21 mfunc - Natural mortality (if age-dependent)
# subplot 22 [no function] - Time-varying growth persp
# subplot 23 [no function] - Time-varying growth contour
# subplot 24 timeVaryingParmFunc - plot time-series of any time-varying quantities
#### unfinished addition on 3-Mar-16
# subplot 25 [no function] - matrix of M by age and time
#### vectors related to hermaphroditism
# subplot 31 Herma_Trans
# subplot 32 Herma_Cum
# EXTRA PLOTS NOT PRODUCED BY DEFAULT
# subplot 101: diagram with labels showing female growth curve
# subplot 102: diagram with labels showing female growth curve & male offsets
# subplot 103: diagram with labels showing female CV = f(A) (offset type 2)
# subplot 104: diagram with labels showing female CV = f(A) & male offset (type 2)
# subplot 105: diagram with labels showing female CV = f(A) (offset type 3)
# subplot 106: diagram with labels showing female CV = f(A) & male offset (type 3)
# test for presence of wtatage_switch
wtatage_switch <- replist[["wtatage_switch"]]
if (wtatage_switch) {
message(
"Note: this model uses the empirical weight-at-age input.\n",
" Plots of many quantities related to growth are skipped."
)
} else {
if (is.null(replist[["endgrowth"]])) {
message(
"Skipping biology plots because model output doesn't include\n",
" biology-at-age information (endgrowth), likely because\n",
" brief output was specified in starter.ss."
)
return()
}
}
plotinfo <- NULL
ians_blues <- c(
"white",
"grey",
"lightblue",
"skyblue",
"steelblue1",
"slateblue",
topo.colors(6),
"blue",
"blue2",
"blue3",
"blue4",
"black"
)
ians_contour <- c("white", rep("blue", 100))
# convert colvec to semi-transparent colors for shading polygons
shadecolvec <- rep(NA, length(colvec))
for (icol in seq_along(colvec)) {
tmp <- col2rgb(colvec[icol]) / 255
shadecolvec[icol] <- rgb(
red = tmp[1],
green = tmp[2],
blue = tmp[3],
alpha = shadealpha
)
}
#### plot function 1
# mean weight, maturity, fecundity, spawning output
# get objects from replist
nseasons <- replist[["nseasons"]]
nareas <- replist[["nareas"]]
growdat <- replist[["endgrowth"]][
replist[["endgrowth"]][["Seas"]] == seas,
]
growthCVtype <- replist[["growthCVtype"]]
biology <- replist[["biology"]]
startyr <- replist[["startyr"]]
FecType <- replist[["FecType"]]
FecPar1name <- replist[["FecPar1name"]]
FecPar2name <- replist[["FecPar2name"]]
FecPar1 <- replist[["FecPar1"]]
FecPar2 <- replist[["FecPar2"]]
parameters <- replist[["parameters"]]
nsexes <- replist[["nsexes"]]
accuage <- replist[["accuage"]]
startyr <- replist[["startyr"]]
endyr <- replist[["endyr"]]
growthvaries <- replist[["growthvaries"]]
growthseries <- replist[["growthseries"]]
ageselex <- replist[["ageselex"]]
MGparmAdj <- replist[["MGparmAdj"]]
wtatage <- replist[["wtatage"]]
# remove comments on each line of wtatage that were added with SSv3.30
if ("comment" %in% names(wtatage)) {
wtatage <- wtatage[, -grep("comment", names(wtatage))]
}
M_at_age <- replist[["M_at_age"]]
Growth_Parameters <- replist[["Growth_Parameters"]]
if (is.null(morphs)) {
morphs <- replist[["mainmorphs"]]
if (is.null(morphs)) {
# if morphs not supplied but growth is available, use middle
# morph within each sex
morphs <- median(growdat[["Morph"]][growdat[["Sex"]] == 1])
if (nsexes == 2) {
morphs <- c(morphs, median(growdat[["Morph"]][growdat[["Sex"]] == 2]))
}
}
}
# get any derived quantities related to growth curve uncertainty
Grow_std <- replist[["derived_quants"]][
grep("Grow_std_", replist[["derived_quants"]][["Label"]]),
]
if (nrow(Grow_std) == 0) {
Grow_std <- NULL
} else {
# convert things like "Grow_std_1_Fem_A_25" into
# in more recent 3.30.14 versions, the label appears as
# "Grow_std_GP:_1_Fem_A_25"
# so the "shift" below = 0 or 1 to adjust the position accordingly
# "pattern 1, female, age 25"
Grow_std[["pattern"]] <- NA
Grow_std[["sex_char"]] <- NA
Grow_std[["sex"]] <- NA
Grow_std[["age"]] <- NA
shift <- length(grep("GP:", Grow_std[["Label"]][1]))
for (irow in 1:nrow(Grow_std)) {
tmp <- strsplit(Grow_std[["Label"]][irow], split = "_")[[1]]
Grow_std[["pattern"]][irow] <- as.numeric(tmp[3 + shift])
Grow_std[["sex_char"]][irow] <- tmp[4 + shift]
Grow_std[["age"]][irow] <- as.numeric(tmp[6 + shift])
}
Grow_std[["sex"]][Grow_std[["sex_char"]] == "Fem"] <- 1
Grow_std[["sex"]][Grow_std[["sex_char"]] == "Mal"] <- 2
#### now it should look something like this:
## Label Value StdDev (Val-1.0)/Stddev
## Grow_std_GP:_1_Fem_A_0 Grow_std_GP:_1_Fem_A_0 27.3188 0.393286 NA
## Grow_std_GP:_1_Fem_A_5 Grow_std_GP:_1_Fem_A_5 50.6099 0.178201 NA
## CumNorm pattern sex_char sex age
## Grow_std_GP:_1_Fem_A_0 NA 1 Fem 1 0
## Grow_std_GP:_1_Fem_A_5 NA 1 Fem 1 5
}
# get any derived quantities related to M uncertainty
NatM_std <- replist[["derived_quants"]][
grep("NatM_std_", replist[["derived_quants"]][["Label"]]),
]
if (nrow(NatM_std) == 0) {
NatM_std <- NULL
} else {
# convert things like "NatM_std_GP:_1_Fem_A_1" into
# "pattern 1, female, age 25"
NatM_std[["pattern"]] <- NA
NatM_std[["sex_char"]] <- NA
NatM_std[["sex"]] <- NA
NatM_std[["age"]] <- NA
for (irow in 1:nrow(NatM_std)) {
tmp <- strsplit(NatM_std[["Label"]][irow], split = "_")[[1]]
NatM_std[["pattern"]][irow] <- as.numeric(tmp[4])
NatM_std[["sex_char"]][irow] <- tmp[5]
NatM_std[["age"]][irow] <- as.numeric(tmp[7])
}
NatM_std[["sex"]][NatM_std[["sex_char"]] == "Fem"] <- 1
NatM_std[["sex"]][NatM_std[["sex_char"]] == "Mal"] <- 2
#### now it should look something like this:
## Label Value StdDev (Val-1.0)/Stddev
## NatM_std_GP:_1_Fem_A_1 NatM_std_GP:_1_Fem_A_1 0.564563 0 NA
## NatM_std_GP:_1_Fem_A_2 NatM_std_GP:_1_Fem_A_2 0.510574 0 NA
## CumNorm pattern sex_char sex age
## NatM_std_GP:_1_Fem_A_1 NA 1 Fem 1 1
## NatM_std_GP:_1_Fem_A_2 NA 1 Fem 1 2
}
if (!seas %in% 1:nseasons) {
stop("'seas' input should be within 1:nseasons")
}
if (nseasons > 1) {
labels[6] <- gsub(
"beginning of the year",
paste("beginning of season", seas),
labels[6]
)
}
if (plotdir == "default") {
plotdir <- replist[["inputs"]][["dir"]]
}
# check dimensions
if (length(morphs) > nsexes) {
warning(
"!Error with morph indexing in SSplotBiology function.\n",
" Code is not set up to handle multiple growth patterns or birth seasons.\n"
)
}
## # stuff from selectivity that is not used
## FecundAtAge <- ageselex[ageselex[["factor"]]=="Fecund", names(ageselex)%in%0:accuage]
## WtAtAge <- ageselex[ageselex[["factor"]]=="bodywt", names(ageselex)%in%0:accuage]
# determine fecundity type
# define labels and x-variable
if (FecType == 1) {
fec_ylab <- "Eggs per kg"
fec_xlab <- labels[8]
FecX <- biology[["Wt_F"]]
FecY <- FecPar1 + FecPar2 * FecX
}
if (labels[11] != "Default fecundity label") {
fec_ylab <- labels[11]
}
# Beginning of season 1 (or specified season) mean length at age
# with 95% range of lengths (by sex if applicable)
## Ian T.: consider somehow generalizing to allow looping over growth pattern
if (!is.null(growdat)) {
# calculate CVs from SD and mean length
growdat[["CV_Beg"]] <- growdat[["SD_Beg"]] / growdat[["Len_Beg"]]
# female growth
growdatF <- growdat[
growdat[["Sex"]] == 1 &
growdat[["Morph"]] == morphs[1],
]
growdatF[["Sd_Size"]] <- growdatF[["SD_Beg"]]
if (growthCVtype == "logSD=f(A)") {
# lognormal distribution of length at age
growdatF[["high"]] <- qlnorm(
0.975,
meanlog = log(growdatF[["Len_Beg"]]),
sdlog = growdatF[["Sd_Size"]]
)
growdatF[["low"]] <- qlnorm(
0.025,
meanlog = log(growdatF[["Len_Beg"]]),
sdlog = growdatF[["Sd_Size"]]
)
} else {
# normal distribution of length at age
growdatF[["high"]] <- qnorm(
0.975,
mean = growdatF[["Len_Beg"]],
sd = growdatF[["Sd_Size"]]
)
growdatF[["low"]] <- qnorm(
0.025,
mean = growdatF[["Len_Beg"]],
sd = growdatF[["Sd_Size"]]
)
}
if (nsexes > 1) {
# do males if 2-sex model
growdatM <- growdat[
growdat[["Sex"]] == 2 & growdat[["Morph"]] == morphs[2],
]
# IAN T. this should probably be generalized
xm <- growdatM[["Age_Beg"]]
growdatM[["Sd_Size"]] <- growdatM[["SD_Beg"]]
if (growthCVtype == "logSD=f(A)") {
# lognormal distribution of length at age
growdatM[["high"]] <- qlnorm(
0.975,
meanlog = log(growdatM[["Len_Beg"]]),
sdlog = growdatM[["Sd_Size"]]
)
growdatM[["low"]] <- qlnorm(
0.025,
meanlog = log(growdatM[["Len_Beg"]]),
sdlog = growdatM[["Sd_Size"]]
)
} else {
# normal distribution of length at age
growdatM[["high"]] <- qnorm(
0.975,
mean = growdatM[["Len_Beg"]],
sd = growdatM[["Sd_Size"]]
)
growdatM[["low"]] <- qnorm(
0.025,
mean = growdatM[["Len_Beg"]],
sd = growdatM[["Sd_Size"]]
)
}
}
}
clean_wtatage <- function(x) {
## quick function to check for valid wtatage matrix and remove first
## redundant row if it's there. Used in weight_plot and maturity_plot.
if (nrow(x) < 2) {
message("Not enough rows in weight-at-age matrix to plot")
return(NULL)
}
if (all(x[1, ] == x[2, ])) {
x <- x[-1, ]
}
return(x)
}
weight_plot <- function(sex) {
# weight
## This needs to be a function of sex since it can be called
## either once for a single sex model or twice to produce plots for
## each one.
if (!wtatage_switch) {
# if empirical weight-at-age is not used
x <- biology[["Len_mean"]]
if (!add) {
ymax <- max(biology[["Wt_F"]])
if (nsexes > 1) {
ymax <- max(ymax, biology[["Wt_M"]])
}
plot(
x,
x,
ylim = c(0, 1.1 * ymax),
xlab = labels[1],
ylab = labels[4],
type = "n",
las = 1,
yaxs = "i"
)
}
lines(x, biology[["Wt_F"]], type = "o", col = colvec[1])
if (nsexes > 1) {
lines(x, biology[["Wt_M"]], type = "o", col = colvec[2])
if (!add) {
legend(
legendloc,
bty = "n",
c("Females", "Males"),
lty = 1,
col = c(colvec[1], colvec[2])
)
}
}
} else {
## if empirical weight-at-age IS used
# calculate maximum weight across all sexes
wtmax <- wtatage[
wtatage[["fleet"]] == -1 &
wtatage[["seas"]] == seas,
-(1:6)
] |>
max()
# get matrix (including initial year column)
wtmat <- wtatage[
wtatage[["fleet"]] == -1 &
wtatage[["sex"]] == sex &
wtatage[["seas"]] == seas,
-(2:6)
]
# check matrix
wtmat <- clean_wtatage(wtmat)
# plot it
if (!is.null(wtmat)) {
makeimage(wtmat, main = "", lastbin = wtmax)
}
}
}
maturity_plot <- function() {
# maturity
if (!wtatage_switch) {
# if empirical weight-at-age is not used
x <- biology[["Len_mean"]]
if (min(biology[["Mat"]]) < 1) {
# if length based
if (!add) {
plot(
x,
biology[["Mat"]],
xlab = labels[1],
ylab = labels[3],
las = 1,
yaxs = "i",
ylim = c(0, max(biology[["Mat"]])),
type = "o",
col = colvec[1]
)
}
if (add) lines(x, biology[["Mat"]], type = "o", col = colvec[1])
} else {
# else is age based
if (!add) {
plot(
growdatF[["Age_Beg"]],
growdatF[["Age_Mat"]],
xlab = labels[2],
las = 1,
yaxs = "i",
ylim = c(0, 1.1 * max(growdatF[["Age_Mat"]])),
ylab = labels[3],
type = "o",
col = colvec[1]
)
} else {
lines(
growdatF[["Age_Beg"]],
growdatF[["Age_Mat"]],
type = "o",
col = colvec[1]
)
}
}
} else {
# if empirical weight-at-age IS used
fecmat <- wtatage[wtatage[["fleet"]] == -2 & wtatage[["sex"]] == 1, ]
if (nrow(fecmat) > 1) {
# figure out which seasons have fecundity values (maybe always only one?)
fecseasons <- sort(unique(fecmat[["seas"]]))
seas_label <- NULL
for (iseas in fecseasons) {
fecmat_seas <- fecmat[fecmat[["seas"]] == iseas, -(2:6)]
# label the season only if a multi-season model
# also testing for length of fecseasons, but that's probably redundant
if (nseasons > 1 | length(fecseasons) > 1) {
seas_label <- paste("in season", iseas)
}
main <- paste("", seas_label)
fecmat_seas <- clean_wtatage(fecmat_seas)
## persp(x=abs(fecmat_seas[,1]),
## y=0:accuage,
## z=as.matrix(fecmat_seas[,-1]),
## theta=70,phi=30,xlab="Year",ylab="Age",zlab="",
## main=main)
}
makeimage(fecmat_seas, main = "")
}
}
}
makeimage <- function(mat, main = "", lastbin = NULL) {
yrvec <- abs(mat[["year"]])
##### this stuff was used to add a row of mean values
## if(is.null(meanvec)){
## meanvec <- mat[,1]
## mat <- mat[,-1]
## }
## yrvec2 <- c(-2:-1+min(yrvec), yrvec)
## mat2 <- cbind(meanvec,NA,mat)
# subset based on minyr and maxyr
subset <- yrvec >= minyr & yrvec <= maxyr
yrvec2 <- yrvec[subset]
mat2 <- mat[subset, -1]
if (is.null(lastbin)) {
lastbin <- max(mat2)
}
z <- t(mat2)
breaks <- seq(0, lastbin, length = 51)
col <- rainbow(60)[1:50]
## Make a two panel plot, for plot and then legend
layout(
matrix(c(1, 2), nrow = 1, ncol = 2),
widths = c(4, .5),
heights = c(4, 4)
)
# save current parameter settings
par_old <- par()
# change parameters
par(mar = c(4.2, 4.2, 1, .5) + .1)
image(
x = 0:accuage,
y = yrvec2,
z = z,
axes = F,
xlab = "Age",
ylab = "Year",
col = col,
breaks = breaks,
main = ifelse(mainTitle, main, "")
)
axis(1, cex.axis = .7)
axis(2, las = 1, cex.axis = .7)
box()
## add text if desired
if (imageplot_text) {
zdataframe <- expand.grid(age = 0:accuage, yr = yrvec2)
zdataframe[["z"]] <- c(t(mat2))
zdataframe[["font"]] <- 1
## Rounding should depend on how big the numbers get. This is untested
## and may need to be adjusted.
ztext <- format(round(zdataframe[["z"]], imageplot_text_round))
ztext[ztext == " NA" | ztext == " NA"] <- ""
text(
x = zdataframe[["age"]],
y = zdataframe[["Yr"]],
label = ztext,
font = zdataframe[["font"]],
cex = .7
)
}
## add a legend to the image plot
## Code adapated from online function. Taken on 11/10/2015 from:
## http://menugget.blogspot.com/2011/08/adding-scale-to-image-plot.html#more
## ------------------------------------------------------------
zlim <- range(z, na.rm = TRUE)
zlim[2] <- zlim[2] + c(zlim[2] - zlim[1]) * (1E-3) # adds a bit to the range in both directions
zlim[1] <- zlim[1] - c(zlim[2] - zlim[1]) * (1E-3)
poly <- vector(mode = "list", length(col))
for (i in seq(poly)) {
poly[[i]] <- c(breaks[i], breaks[i + 1], breaks[i + 1], breaks[i])
}
ylim <- range(breaks)
xlim <- c(0, 1)
par(mar = c(4.3, .5, 1, 3))
plot(
1,
1,
type = "n",
ylim = ylim,
xlim = xlim,
axes = FALSE,
ann = FALSE,
xaxs = "i",
yaxs = "i"
)
for (i in seq(poly)) {
polygon(c(0, 0, 1, 1), poly[[i]], col = col[i], border = NA)
}
axis(4, cex.axis = .7, las = 2)
box()
## turn off the matrix plot
layout(mat = c(1, 1))
# restore default single panel settings
par(mfcol = c(1, 1), mar = par_old[["mar"]], oma = par_old[["oma"]])
}
fec_pars_fn <- function() {
# fecundity from model parameters
ymax <- 1.1 * max(FecY)
if (!add) {
plot(
FecX,
FecY,
xlab = fec_xlab,
ylab = fec_ylab,
ylim = c(0, ymax),
las = 1,
yaxs = "i",
col = colvec[2],
pch = 19
)
lines(FecX, rep(FecPar1, length(FecX)), col = colvec[1])
text(
mean(range(FecX)),
FecPar1 - 0.05 * ymax,
"Egg output proportional to spawning biomass"
)
} else {
points(FecX, FecY, col = colvec[2], pch = 19)
}
}
fecundityOK <- all(!is.na(biology[["Fec"]]))
fec_weight_fn <- function() {
# fecundity at weight from BIOLOGY section
ymax <- 1.1 * max(biology[["Fec"]])
if (!add) {
plot(
biology[["Wt_F"]],
biology[["Fec"]],
xlab = labels[8],
ylab = labels[10],
las = 1,
yaxs = "i",
ylim = c(0, ymax),
col = colvec[1],
type = "o"
)
} else {
points(
biology[["Len_mean"]],
biology[["Fec"]],
col = colvec[1],
type = "o"
)
}
}
fec_len_fn <- function() {
# fecundity at length from BIOLOGY section
ymax <- 1.1 * max(biology[["Fec"]])
if (!add) {
plot(
biology[["Len_mean"]],
biology[["Fec"]],
xlab = labels[9],
ylab = labels[10],
las = 1,
yaxs = "i",
ylim = c(0, 1.1 * ymax),
col = colvec[1],
type = "o",
yaxs = "i"
)
} else {
points(
biology[["Len_mean"]],
biology[["Fec"]],
col = colvec[1],
type = "o"
)
}
}
spawn_output_len_fn <- function() {
# spawning output at length
x <- biology[["Len_mean"]]
y <- biology[["Mat*Fec"]]
ymax <- 1.1 * max(y)
if (!add) {
plot(
x,
y,
xlab = labels[1],
ylab = labels[5],
type = "o",
col = colvec[1],
las = 1,
yaxs = "i",
ylim = c(0, ymax)
)
} else {
lines(x, y, type = "o", col = colvec[1])
}
}
spawn_output_age_fn <- function() {
# spawning output at age
x <- growdat[["Age_Beg"]][growdat[["Sex"]] == 1]
y <- growdat$"Mat*Fecund"[growdat[["Sex"]] == 1]
ymax <- 1.1 * max(y)
if (!add) {
plot(
x,
y,
xlab = labels[2],
ylab = labels[5],
type = "o",
col = colvec[1],
las = 1,
yaxs = "i",
ylim = c(0, ymax)
)
} else {
lines(x, y, type = "o", col = colvec[1])
}
}
if (!wtatage_switch) {
ymax <- max(biology[["Len_mean"]])
x <- growdatF[["Age_Beg"]]
}
main <- "Ending year expected growth (with 95% intervals)"
# if(nseasons > 1){main <- paste(main," season 1",sep="")}
col_index1 <- 3 # default is grey for single-sex model
if (nsexes > 1) {
col_index1 <- 1 # change line for females to red
}
add_uncertainty_fn <- function(std_table, sex, age_offset = 0.5) {
# function to add uncertainty estimates to mean length at age or M at age
# previous settings for warnings
old_warn <- options()[["warn"]]
# turn off "zero-length arrow" warning
options(warn = -1)
# confirm presence of table of values and request to add to plot
if (!is.null(std_table)) {
std_table.sex <- std_table[std_table[["sex"]] == sex, ]
# confirm presence for requested sex
if (!is.null(std_table.sex)) {
for (irow in 1:nrow(std_table.sex)) {
std.age <- std_table.sex[["age"]][irow]
# this might not work for log-normal length at age
mean <- std_table.sex[["Value"]][irow]
ages <- std_table.sex[["age"]][irow] + age_offset
high <- qnorm(
0.975,
mean = mean,
sd = std_table.sex[["StdDev"]][irow]
)
low <- qnorm(
0.025,
mean = mean,
sd = std_table.sex[["StdDev"]][irow]
)
# set sex-dependent color
if (sex == 1) {
col <- colvec[col_index1]
}
if (sex == 2) {
col <- colvec[2]
}
arrows(
x0 = ages,
x1 = ages,
y0 = low,
y1 = high,
length = 0.04,
angle = 90,
code = 3,
col = col
)
}
} # end check for missing sex
} # end check for empty table
# returning to old warnings (after turning off in case of zero-length arrow)
options(warn = old_warn)
}
growth_curve_fn <- function(add_labels = TRUE, add_uncertainty = TRUE) {
# growth
x <- growdatF[["Age_Beg"]]
# make empty plot unless this is being added to existing figure
if (!add) {
plot(
x,
growdatF[["Len_Beg"]],
ylim = c(0, 1.1 * ymax),
type = "n",
xlab = "",
ylab = "",
axes = FALSE,
xaxs = "i",
yaxs = "i"
)
abline(h = 0, col = "grey")
if (add_labels) {
title(
main = ifelse(mainTitle, main, ""),
xlab = labels[2],
ylab = labels[6],
cex.main = cex.main
)
axis(1)
axis(2, las = 1)
}
}
lty <- 1
# make shaded polygon
polygon(
c(x, rev(x)),
c(growdatF[["low"]], rev(growdatF[["high"]])),
border = NA,
col = shadecolvec[col_index1]
)
# add lines for mean, low and high
lines(
x,
growdatF[["Len_Beg"]],
col = colvec[col_index1],
lwd = 2,
lty = ltyvec[1]
)
lines(
x,
growdatF[["high"]],
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
lines(x, growdatF[["low"]], col = colvec[col_index1], lwd = 1, lty = "12")
# add 95% intervals for uncertainty in growth from $derived_quants
if (add_uncertainty) {
add_uncertainty_fn(std_table = Grow_std, sex = 1, age_offset = 0.5)
}
# add males if they are present in the model
if (nsexes > 1) {
# make shaded polygon
polygon(
c(xm, rev(xm)),
c(growdatM[["low"]], rev(growdatM[["high"]])),
border = NA,
col = shadecolvec[2]
)
# add lines for mean, low and high
lines(
xm,
growdatM[["Len_Beg"]],
col = colvec[2],
lwd = 2,
lty = ltyvec[2]
)
lines(xm, growdatM[["high"]], col = colvec[2], lwd = 1, lty = "13")
lines(xm, growdatM[["low"]], col = colvec[2], lwd = 1, lty = "13")
# add 95% intervals for uncertainty in growth from $derived_quants
if (add_uncertainty) {
add_uncertainty_fn(std_table = Grow_std, sex = 2, age_offset = 0.5)
}
}
if (!add) {
grid()
box()
}
if (nsexes > 1) {
legend(
legendloc,
bty = "n",
c("Females", "Males"),
lty = ltyvec,
lwd = 2,
col = c(colvec[1], colvec[2])
)
}
}
if (plot & 1 %in% subplots & !wtatage_switch) {
growth_curve_fn()
}
if (print & 1 %in% subplots & !wtatage_switch) {
file <- "bio1_sizeatage.png"
caption <- paste(
"Length at age in the beginning of the year (or season) in the ending",
"year of the model. Shaded area indicates 95% distribution of",
"length at age around estimated growth curve."
)
if (!is.null(Grow_std)) {
caption <- paste(
caption,
"Vertical intervals around growth curve indicate",
"estimated 95% uncertainty intervals in estimated mean growth."
)
}
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_fn()
dev.off()
}
growth_curve_plus_fn <- function(add_labels = TRUE, option = 1) {
# function to add panels to growth curve with info on variability in
# length at age or info on maturity, weight, and fecundity
# save current parameter settings
par_old <- par()
# create 4-panel layout
layout(
mat = matrix(1:4, byrow = TRUE, ncol = 2),
widths = c(2, 1),
heights = c(2, 1)
)
par(mar = c(1, 1, 1, 1), oma = c(5, 5, 5, 4))
# plot growth curve
growth_curve_fn(add_labels = FALSE)
axis(2, las = 1)
# add age axis on top
axis(3)
# add some labels
mtext(side = 2, line = 2.5, labels[6], las = 0)
mtext(side = 3, line = 2.5, labels[2], las = 0)
box()
if (option == 1) {
lab1 <- "SD_Beg"
lab1long <- "SD of lengths"
lab2 <- "CV_Beg"
lab2long <- "CV of lengths"
lab1max <- 1.1 * max(growdat[[lab1]], na.rm = TRUE)
lab2max <- 1.05 * max(growdat[[lab2]], na.rm = TRUE)
lab1_axis_vec <- NULL
}
if (option == 2) {
lab1 <- "Len_Mat"
lab1long <- "Frac. mature"
lab2 <- "Wt_Beg"
lab2long <- "Mean weight"
lab1max <- 1
lab2max <- max(c(biology[["Wt_F"]], biology[["Wt_M"]]), na.rm = TRUE)
lab1_axis_vec <- c(0, 0.5, 1)
}
# calculate scaling factor between CVs and SDs to share each panel
lab2_to_lab1_scale <- lab1max / lab2max
lab2_axis_vec <- pretty(c(0, lab2max))
# make empty panel in top-right location
plot(
0,
type = "n",
xaxs = "i",
xlim = c(0, 1.0 * lab1max),
yaxs = "i",
ylim = par()[["usr"]][3:4],
axes = FALSE
)
if (option == 1) {
# if plotting SD and CV, then get it from age-based data
# add line for lab1 vs. Len
lines(
growdatF[[lab1]],
growdatF[["Len_Beg"]],
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
# add line for lab2 vs. Len
lines(
growdatF[[lab2]] * lab2_to_lab1_scale,
growdatF[["Len_Beg"]],
col = colvec[col_index1],
lwd = 3
)
if (nsexes > 1) {
# add lines for males
if (option == 1) {
# add line for lab1 vs. Age
# (unless it's maturity, which is not defined for males)
lines(
growdatM[[lab1]],
growdatM[["Len_Beg"]],
col = colvec[2],
lwd = 1,
lty = "13"
)
}
# add line for lab2 vs. Len
lines(
growdatM[[lab2]] * lab2_to_lab1_scale,
growdatM[["Len_Beg"]],
col = colvec[2],
lwd = 3,
lty = 2
)
}
}
if (option == 2) {
# if plotting maturity and fecundity, then get this panel from length-based data
lines(
biology[["Wt_F"]] * lab2_to_lab1_scale,
biology[["Len_mean"]],
col = colvec[col_index1],
lwd = 3
)
lines(
biology[["Mat"]],
biology[["Len_mean"]],
col = colvec[col_index1],
lty = "12"
)
if (nsexes > 1) {
lines(
biology[["Wt_M"]] * lab2_to_lab1_scale,
biology[["Len_mean"]],
col = colvec[2],
lwd = 3,
lty = 2
)
}
}
# add axes and labels
mtext(side = 4, line = 2.5, labels[6], las = 0)
mtext(side = 1, line = 2.5, lab1long, las = 0)
mtext(side = 3, line = 2.5, lab2long, las = 0)
axis(1, at = lab1_axis_vec, las = 1)
axis(
3,
at = lab2_axis_vec * lab2_to_lab1_scale,
labels = lab2_axis_vec,
las = 1
)
axis(4, las = 1)
box()
# make empty plot in lower-left position
plot(
0,
type = "n",
xaxs = "i",
xlim = range(growdat[["Age_Beg"]]),
yaxs = "i",
ylim = c(0, 1.0 * lab1max),
axes = FALSE
)
# add line for lab1 vs. Age
if (option == 1) {
lines(
growdatF[["Age_Beg"]],
growdatF[[lab1]],
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
}
if (option == 2) {
# maturity curve is product of age-based and length-based factors
lines(
growdatF[["Age_Beg"]],
growdatF[["Len_Mat"]] * abs(growdatF[["Age_Mat"]]),
col = colvec[col_index1],
lwd = 1,
lty = "12"
)
}
# add line for lab2 vs. Age
lines(
growdatF[["Age_Beg"]],
growdatF[[lab2]] * lab2_to_lab1_scale,
col = colvec[col_index1],
lwd = 3
)
if (nsexes > 1) {
# add lines for males
if (option == 1) {
# add line for lab1 vs. Age (unless it's maturity, which is not defined for males)
lines(
growdatM[["Age_Beg"]],
growdatM[[lab1]],
col = colvec[2],
lwd = 1,
lty = "13"
)
}
# add line for lab2 vs. Age
lines(
growdatM[["Age_Beg"]],
growdatM[[lab2]] * lab2_to_lab1_scale,
col = colvec[2],
lwd = 3,
lty = 2
)
}
# add axes and labels
mtext(side = 1, line = 2.5, labels[2], las = 0)
mtext(side = 4, line = 2.5, lab1long, las = 0)
mtext(side = 2, line = 2.5, lab2long, las = 0)
axis(1, las = 1)
axis(
2,
at = lab2_axis_vec * lab2_to_lab1_scale,
labels = lab2_axis_vec,
las = 1
)
axis(4, at = lab1_axis_vec, las = 1)
box()
# restore default single panel settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
} # end growth_curve_plus_fn()
# make plots of growth curve with CV and SD of length
if (plot & 2 %in% subplots & !wtatage_switch) {
growth_curve_plus_fn(option = 1)
}
if (print & 2 %in% subplots & !wtatage_switch) {
file <- "bio2_sizeatage_plus_CV_and_SD.png"
caption <- paste(
"Length at age (top-left panel) with",
"CV (thick line) and SD (thin line) of",
"length at age shown in top-right and lower-left panels"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_plus_fn(option = 1)
dev.off()
}
# make plots of growth curve with weight-length curve and maturity
if (plot & 3 %in% subplots & !wtatage_switch) {
growth_curve_plus_fn(option = 2)
}
if (print & 3 %in% subplots & !wtatage_switch) {
file <- "bio3_sizeatage_plus_WT_and_MAT.png"
caption <- paste(
"Length at age (top-left panel) with",
"weight (thick line) and maturity (thin line)",
"shown in top-right and lower-left panels"
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_plus_fn(option = 2)
dev.off()
}
# plot distribution of length at age (by season, sub-season, and morph)
if (4 %in% subplots & !wtatage_switch) {
if (!is.null(replist[["ALK"]])) {
plotinfo.tmp <- SSplotAgeMatrix(
replist = replist,
option = 1,
plot = plot,
print = print,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
cex.main = cex.main,
mainTitle = mainTitle
)
# SSplotAgeMatrix adds a "category" column which isn't present
# in plotinfo until the end of this SSplotBiology function.
plotinfo.tmp <- plotinfo.tmp[, c("file", "caption", "alt_text")]
plotinfo <- rbind(plotinfo, plotinfo.tmp)
} else {
message(
"Skipped some plots because AGE_LENGTH_KEY unavailable in report file\n",
"because starter file set to produce limited report detail."
)
}
}
# function for illustrating parameterization of growth curves
growth_curve_labeled_fn <- function(option = 1) {
# growth
if (is.null(Growth_Parameters)) {
message(
"Need updated SS_output function to get Growth_Parameters output\n"
)
return()
}
# save current parameter settings
par_old <- par()
par(mar = c(4, 4, 1, 4))
AminF <- Growth_Parameters[["A1"]][1]
AmaxF <- Growth_Parameters[["A2"]][1]
AminM <- Growth_Parameters[["A1"]][2]
AmaxM <- Growth_Parameters[["A2"]][2]
L_at_AminF <- Growth_Parameters[["L_a_A1"]][1]
L_at_AmaxF <- Growth_Parameters[["L_a_A2"]][1]
L_at_AminM <- Growth_Parameters[["L_a_A1"]][2]
L_at_AmaxM <- Growth_Parameters[["L_a_A2"]][2]
LinfF <- Growth_Parameters[["Linf"]][1]
LinfM <- Growth_Parameters[["Linf"]][2]
ymax <- max(biology[["Len_mean"]])
plot(
0,
type = "n",
xlim = c(0, 1 + max(growdatF[["Age_Beg"]])),
ylim = c(0, 1.1 * ymax),
xlab = "",
ylab = "",
axes = FALSE,
xaxs = "i",
yaxs = "i"
)
abline(h = 0, col = "grey")
# title(main=main, xlab=labels[2], ylab=labels[6], cex.main=cex.main)
axis(
1,
at = c(0, AminF, AmaxF, accuage),
labels = expression(0, italic(A[1]), italic(A[2]), italic(N[ages]))
)
axis(
2,
at = c(0, L_at_AminF, L_at_AmaxF, LinfF),
labels = expression(
0,
italic(L[A[1]]),
italic(L[A[2]]),
italic(L[infinity])
),
las = 1
)
# add growth curve itself
lines(
growdatF[["Age_Beg"]],
growdatF[["Len_Beg"]],
col = 1,
lwd = 3,
lty = 1
)
# add growth curve for males
if (nsexes > 1 & option == 2) {
lines(
growdatM[["Age_Beg"]],
growdatM[["Len_Beg"]],
col = 1,
lwd = 2,
lty = 2
)
}
if (option == 1) {
# lines connecting axes to curve
lines(c(AminF, AminF, 0), c(0, L_at_AminF, L_at_AminF), lty = 3)
lines(c(AmaxF, AmaxF, 0), c(0, L_at_AmaxF, L_at_AmaxF), lty = 3)
abline(h = LinfF, lty = 3)
}
if (option == 2) {
# lines connecting two curves
lines(c(0, AminF), c(L_at_AminF, L_at_AminF), lty = 3)
lines(c(0, AmaxF), c(L_at_AmaxF, L_at_AmaxF), lty = 3)
lines(c(accuage + 10, AmaxM), c(L_at_AmaxM, L_at_AmaxM), lty = 3)
lines(c(accuage + 10, AminM), c(L_at_AminM, L_at_AminM), lty = 3)
arrows(
x0 = AminF,
x1 = AminF,
y0 = L_at_AminF,
y1 = L_at_AminM,
lwd = 3,
col = 3,
length = 0.1
)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = L_at_AmaxF,
y1 = L_at_AmaxM,
lwd = 3,
col = 3,
length = 0.1
)
text(
AminF,
mean(c(L_at_AminF, L_at_AminM)),
"Male parameter 1\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
text(
AmaxF,
mean(c(L_at_AmaxF, L_at_AmaxM)),
"Male parameter 2\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
axis(
4,
at = c(0, L_at_AminM, L_at_AmaxM),
labels = expression(0, italic(L[A[1]]^male), italic(L[A[2]]^male)),
las = 1
)
abline(h = LinfF, lty = 3)
}
box()
# restore old parameter settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
}
# plots of diagrams for illustrating parameterization
if (plot & 101 %in% subplots & !wtatage_switch) {
growth_curve_labeled_fn(option = 1)
}
if (print & 101 %in% subplots & !wtatage_switch) {
file <- "bio101_growth_illustration.png"
caption <- "Illustration of growth parameters"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_labeled_fn(option = 1)
dev.off()
}
if (plot & 102 %in% subplots & !wtatage_switch) {
growth_curve_labeled_fn(option = 2)
}
if (print & 102 %in% subplots & !wtatage_switch) {
file <- "bio102_growth_illustration2.png"
caption <- "Illustration of growth parameters with male offsets"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
growth_curve_labeled_fn(option = 2)
dev.off()
}
# function for illustrating parameterization of CVs around growth curves
CV_values_labeled_fn <- function(option = 1) {
# growth
if (is.null(Growth_Parameters)) {
message(
"Need updated SS_output function to get Growth_Parameters output"
)
return()
}
# save current parameter settings
par_old <- par()
par(mar = c(4, 4, 1, 4))
AminF <- Growth_Parameters[["A1"]][1]
AmaxF <- Growth_Parameters[["A2"]][1]
AminM <- Growth_Parameters[["A1"]][2]
AmaxM <- Growth_Parameters[["A2"]][2]
CV_at_AminF <- Growth_Parameters[["CVmin"]][1]
CV_at_AmaxF <- Growth_Parameters[["CVmax"]][1]
CV_at_AminM <- Growth_Parameters[["CVmin"]][2]
CV_at_AmaxM <- Growth_Parameters[["CVmax"]][2]
ymax <- max(CV_at_AminF, CV_at_AmaxF, CV_at_AminM, CV_at_AmaxM)
plot(
0,
type = "n",
xlim = c(0, 1 + max(growdatF[["Age_Beg"]])),
ylim = c(0, 1.1 * ymax),
xlab = "",
ylab = "",
axes = FALSE,
xaxs = "i",
yaxs = "i"
)
abline(h = 0, col = "grey")
# title(main=main, xlab=labels[2], ylab=labels[6], cex.main=cex.main)
axis(
1,
at = c(0, AminF, AmaxF, accuage),
labels = expression(0, italic(A[1]), italic(A[2]), italic(N[ages]))
)
axis(
2,
at = c(0, CV_at_AminF, CV_at_AmaxF),
labels = expression(0, italic(CV[A[1]]), italic(CV[A[2]])),
las = 1
)
# add growth curve itself
lines(
growdatF[["Age_Beg"]],
growdatF[["CV_Beg"]],
col = 1,
lwd = 3,
lty = 1
)
# add growth curve for males
if (nsexes > 1 & option %in% c(2, 4)) {
lines(
growdatM[["Age_Beg"]],
growdatM[["CV_Beg"]],
col = 1,
lwd = 2,
lty = 2
)
}
if (option == 1) {
# lines connecting axes to curve
lines(c(AminF, AminF, 0), c(0, CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(AmaxF, AmaxF, 0), c(0, CV_at_AmaxF, CV_at_AmaxF), lty = 3)
}
if (option == 2) {
# lines connecting two curves
lines(c(0, AminF), c(CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(0, AmaxF), c(CV_at_AmaxF, CV_at_AmaxF), lty = 3)
lines(c(accuage + 10, AmaxM), c(CV_at_AmaxM, CV_at_AmaxM), lty = 3)
lines(c(accuage + 10, AminM), c(CV_at_AminM, CV_at_AminM), lty = 3)
arrows(
x0 = AminF,
x1 = AminF,
y0 = CV_at_AminF,
y1 = CV_at_AminM,
lwd = 3,
col = 3,
length = 0.1
)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AmaxF,
y1 = CV_at_AmaxM,
lwd = 3,
col = 3,
length = 0.1
)
text(
AminF,
mean(c(CV_at_AminM)),
"Male parameter 1\n(exponential offset)",
col = 3,
adj = c(0, 1.5)
)
text(
AmaxF,
mean(c(CV_at_AmaxM)),
"Male parameter 2\n(exponential offset)",
col = 3,
adj = c(0, 1.5)
)
axis(
4,
at = c(0, CV_at_AminM, CV_at_AmaxM),
labels = expression(0, italic(CV[A[1]]^male), italic(CV[A[2]]^male)),
las = 1
)
}
if (option == 3) {
# lines connecting axes to curve
lines(c(AminF, AminF, 0), c(0, CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(AmaxF, AmaxF, 0), c(0, CV_at_AmaxF, CV_at_AmaxF), lty = 3)
lines(c(0, AmaxF), c(CV_at_AminF, CV_at_AminF), lty = 3)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AminF,
y1 = CV_at_AmaxF,
lwd = 3,
col = 3,
length = 0.1
)
text(
AmaxF,
mean(c(CV_at_AminF, CV_at_AmaxF)),
"Female parameter 2\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
}
if (option == 4) {
# lines connecting two curves
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AminF,
y1 = CV_at_AmaxF,
lwd = 3,
col = 3,
length = 0.1
)
text(
AmaxF,
mean(c(CV_at_AminF, CV_at_AmaxF)),
"Female parameter 2\n(exponential offset)",
col = 3,
adj = c(-.1, 0)
)
lines(c(0, AmaxF), c(CV_at_AminF, CV_at_AminF), lty = 3)
lines(c(0, AmaxF), c(CV_at_AmaxF, CV_at_AmaxF), lty = 3)
lines(c(accuage + 10, AmaxM), c(CV_at_AmaxM, CV_at_AmaxM), lty = 3)
lines(c(accuage + 10, AminM), c(CV_at_AminM, CV_at_AminM), lty = 3)
arrows(
x0 = AminF,
x1 = AminF,
y0 = CV_at_AminF,
y1 = CV_at_AminM,
lwd = 3,
col = 3,
length = 0.1
)
arrows(
x0 = AmaxF,
x1 = AmaxF,
y0 = CV_at_AminM,
y1 = CV_at_AmaxM,
lwd = 3,
col = 3,
length = 0.1
)
text(
AminF,
mean(c(CV_at_AminF, CV_at_AminM)),
"Male parameter 1\n(exponential offset)",
col = 3,
adj = c(0, 2)
)
text(
AmaxF,
mean(c(CV_at_AminM, CV_at_AmaxM)),
"Male parameter 2\n(exponential offset)",
col = 3,
adj = c(0, 2)
)
axis(
4,
at = c(0, CV_at_AminM, CV_at_AmaxM),
labels = expression(0, italic(CV[A[1]]^male), italic(CV[A[2]]^male)),
las = 1
)
}
box()
# restore old parameter settings
par(
mfcol = par_old[["mfcol"]],
mar = par_old[["mar"]],
oma = par_old[["oma"]]
)
}
# extra plots illustrating parameterization of CV in growth
if (plot & 103 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 1)
}
if (print & 103 %in% subplots & !wtatage_switch) {
file <- "bio103_CV_illustration.png"
caption <- "Illustration of growth variability parameters"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 1)
dev.off()
}
if (plot & 104 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 2)
}
if (print & 104 %in% subplots & !wtatage_switch) {
file <- "bio104_CV_illustration2.png"
caption <- "Illustration of growth variability parameters with male offsets"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 2)
dev.off()
}
if (plot & 105 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 3)
}
if (print & 105 %in% subplots & !wtatage_switch) {
file <- "bio105_CV_illustration.png"
caption <- "Illustration of growth variability parameters for offset type 3"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 3)
dev.off()
}
if (plot & 106 %in% subplots & !wtatage_switch) {
CV_values_labeled_fn(option = 4)
}
if (print & 106 %in% subplots & !wtatage_switch) {
file <- "bio106_CV_illustration2.png"
caption <- "Illustration of growth variability parameters with male offsets"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
CV_values_labeled_fn(option = 4)
dev.off()
}
x <- biology[["Len_mean"]]
## NOTE: weight plots are now a special case since they are broken down
## by whether the model is 1-sex or 2-sex. In the latter two separate
## plots need to be made.
if (plot) {
# plot to screen or to PDF file
if (5 %in% subplots & !is.null(biology)) {
weight_plot(sex = 1)
if (nsexes == 2) {
weight_plot(sex = 2)
}
}
if (6 %in% subplots) {
maturity_plot()
}
if (!wtatage_switch) {
if (7 %in% subplots & FecType == 1) {
fec_pars_fn()
}
if (8 %in% subplots & fecundityOK) {
fec_weight_fn()
}
if (9 %in% subplots & fecundityOK) {
fec_len_fn()
}
if (10 %in% subplots) {
spawn_output_len_fn()
}
if (11 %in% subplots) {
spawn_output_age_fn()
}
}
}
if (print) {
# print to PNG files
if (5 %in% subplots) {
file <- "bio5_weightatsize.png"
caption <- "Weight-length relationship"
if (wtatage_switch) {
caption <- "Weight-at-age"
if (nsexes == 1) {
caption <- paste(caption, "for females")
}
}
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = ifelse(wtatage_switch, pheight_tall, pheight), # weight-at-age matrix works better when taller
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
weight_plot(sex = 1)
dev.off()
# add separate plot for males (only for empirical weight-at-age models)
if (wtatage_switch & nsexes == 2) {
file <- "bio5_weightatsize2.png"
caption <- "Weight-at-age for males"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight_tall,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
weight_plot(sex = 2)
dev.off()
}
}
if (6 %in% subplots) {
file <- "bio6_maturity.png"
if (!wtatage_switch) {
caption <- paste(
"Maturity at",
ifelse(min(biology[["Mat"]]) < 1, "length", "age")
)
} else {
caption <- "Spawning output at age (maturity x fecundity)"
}
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = ifelse(wtatage_switch, pheight_tall, pheight), # spawning-output-at-age matrix works better when taller
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
maturity_plot()
dev.off()
}
if (!wtatage_switch) {
if (7 %in% subplots & FecType == 1) {
file <- "bio7_fecundity.png"
caption <- "Fecundity"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
fec_pars_fn()
dev.off()
}
if (8 %in% subplots & fecundityOK) {
file <- "bio8_fecundity_wt.png"
caption <- "Fecundity as a function of weight"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
fec_weight_fn()
dev.off()
}
if (9 %in% subplots & fecundityOK) {
file <- "bio9_fecundity_len.png"
caption <- "Fecundity as a function of length"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
fec_len_fn()
dev.off()
}
if (10 %in% subplots) {
file <- "bio10_spawningoutput_len.png"
caption <- paste(
"Spawning output at length.",
"This is the product of maturity and fecundity",
"unless maturity is age-based, in which case only",
"fecundity is represented."
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
spawn_output_len_fn()
dev.off()
}
if (11 %in% subplots) {
file <- "bio11_spawningoutput_age.png"
caption <- paste(
"Spawning output at age.",
"This is the product of maturity and fecundity.",
"When these processes are length-based they are",
"converted into the age dimension using the matrix",
"of length at age."
)
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
spawn_output_age_fn()
dev.off()
}
} # end check for wtatage model
} # end check for print=TRUE (printing to PNG files)
# Natural mortality (if age-dependent -- need to add time-varying M plot)
MatAge <- growdatF[["M"]] # female mortality in the ending year
if (nsexes > 1) {
growdatM <- growdat[growdat[["Morph"]] == morphs[2], ]
MatAge_M <- growdatM[["M"]]
}
if (min(MatAge) != max(MatAge) & 21 %in% subplots) {
ymax <- 1.1 * max(MatAge)
if (nsexes > 1) {
ymax <- 1.1 * max(MatAge, MatAge_M)
}
# function to plut natural mortality
mfunc <- function(add_uncertainty = TRUE) {
if (!add) {
plot(
growdatF[["Age_Beg"]],
MatAge,
col = colvec[col_index1],
lwd = 2,
ylim = c(0, ymax),
yaxs = "i",
type = "n",
ylab = labels[7],
xlab = labels[2],
las = 1
)
}
lines(
growdatF[["Age_Beg"]],
MatAge,
col = colvec[col_index1],
lwd = 2,
type = "o"
)
# add uncertainty in M-at-age (if available from derived_quants)
if (add_uncertainty) {
add_uncertainty_fn(std_table = NatM_std, sex = 1, age_offset = 0.0)
}
if (nsexes > 1) {
growdatM <- growdat[growdat[["Morph"]] == morphs[2], ]
lines(
growdatM[["Age_Beg"]],
growdatM[["M"]],
lty = ltyvec[2],
col = colvec[2],
lwd = 2,
type = "o"
)
legend(
"bottomleft",
bty = "n",
c("Females", "Males"),
lty = ltyvec,
lwd = 2,
col = c(colvec[col_index1], colvec[2])
)
# add uncertainty in M-at-age (if available from derived_quants)
if (add_uncertainty) {
add_uncertainty_fn(std_table = NatM_std, sex = 2, age_offset = 0.0)
}
}
}
# run function if requested to make plot
if (plot & 21 %in% subplots) {
mfunc()
}
# run function if requested to write figure to PNG file
if (print & 21 %in% subplots) {
file <- "bio21_natmort.png"
caption <- "Natural mortality at age"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
mfunc()
dev.off()
}
}
# Time-varying growth
if (is.null(growthvaries)) {
if (verbose) {
message(
"No check for time-varying growth because starter file set to produce\n",
"limited report detail."
)
}
} else {
# temporarily disable multi-season plotting of time-varying growth
if (is.null(growthseries)) {
warning(
"No time-varying growth info because starter file set to produce\n",
"limited report detail."
)
} else {
# if growth is time varying and weight-at-age not used
if (growthvaries & !wtatage_switch) {
for (i in 1:nsexes) {
growdatuse <- growthseries[
growthseries[["Yr"]] >= startyr - 2 &
growthseries[["Morph"]] == morphs[i],
]
# trim based on forecast=TRUE/FALSE
if (!forecast) {
growdatuse <- growdatuse[growdatuse[["Yr"]] <= endyr, ]
}
# trim based on minyr and maxyr
growdatuse <- growdatuse[
growdatuse[["Yr"]] >= minyr &
growdatuse[["Yr"]] <= maxyr,
]
# assemble vectors and matrix
x <- 0:accuage
y <- growdatuse[["Yr"]]
z <- as.matrix(growdatuse[, -(1:4)])
# check for time-varying growth
if (replist[["growthvaries"]]) {
z <- t(z)
if (i == 1) {
main <- "Female time-varying growth"
}
if (nsexes == 1) {
main <- "Time-varying growth"
}
if (i == 2) {
main <- "Male time-varying growth"
}
if (nseasons > 1) {
main <- paste0(main, " season 1")
}
if (plot) {
if (22 %in% subplots) {
persp(
x,
y,
z,
col = "white",
xlab = labels[2],
ylab = "",
zlab = labels[1],
expand = 0.5,
box = TRUE,
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
ticktype = "detailed",
phi = 35,
theta = -10
)
}
if (23 %in% subplots) {
contour(
x,
y,
z,
nlevels = 12,
xlab = labels[2],
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
col = ians_contour,
lwd = 2
)
}
}
if (print) {
if (22 %in% subplots) {
file <- paste0("bio22_timevarygrowthsurf_sex", i, ".png")
caption <- "Perspective plot of time-varying growth"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
persp(
x,
y,
z,
col = "white",
xlab = labels[2],
ylab = "",
zlab = labels[1],
expand = 0.5,
box = TRUE,
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
ticktype = "detailed",
phi = 35,
theta = -10
)
dev.off()
}
if (23 %in% subplots) {
file <- paste0("bio23_timevarygrowthcontour_sex", i, ".png")
caption <- "Contour plot of time-varying growth"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
contour(
x,
y,
z,
nlevels = 12,
xlab = labels[2],
main = ifelse(mainTitle, main, ""),
cex.main = cex.main,
col = ians_contour,
lwd = 2
)
dev.off()
}
} # end print
} # end if time-varying
} # end loop over sexes
} # end if growth is time varying
} # end of if data available for time varying growth
} # end disable of time-varying growth for multi-season models
# plot time-series of any time-varying quantities
if (24 %in% subplots) {
if (!is.null(MGparmAdj)) {
# general function to work for any parameter
timeVaryingParmFunc <- function(parmlabel, forecast = FALSE) {
if (forecast) {
MGparmAdj.tmp <- MGparmAdj
} else {
MGparmAdj.tmp <- MGparmAdj[MGparmAdj[["Yr"]] <= endyr, ]
}
# trim based on minyr and maxyr
MGparmAdj.tmp <- MGparmAdj.tmp[
MGparmAdj.tmp[["Yr"]] >= minyr &
MGparmAdj.tmp[["Yr"]] <= maxyr,
]
# make plot
plot(
MGparmAdj.tmp[["Yr"]],
MGparmAdj.tmp[[parmlabel]],
xlab = labels[12],
ylab = parmlabel,
type = "l",
lwd = 3,
col = colvec[2]
)
}
# check to make sure MGparmAdj looks as expected
# (maybe had different or conditional format in old SS versions)
if (!is.null(ncol(MGparmAdj)) && ncol(MGparmAdj) > 1) {
# loop over columns looking for time-varying parameters
for (icol in 2:ncol(MGparmAdj)) {
parmlabel <- names(MGparmAdj)[icol]
# exclude column indicating change added with version 3.30.06.02
if (parmlabel != "Change?") {
parmvals <- MGparmAdj[, icol]
# check for changes
if (length(unique(parmvals[MGparmAdj[["Yr"]] <= endyr])) > 1) {
# make plot
if (plot) {
timeVaryingParmFunc(parmlabel, forecast = forecast)
}
if (print) {
file <- paste0("bio24_time-varying_", parmlabel, ".png")
# replace % sign which cause problems for filename
file <- gsub(
pattern = "%",
replacement = "percent",
x = file,
fixed = TRUE
)
# replace : which which cause problems for filename
file <- gsub(
pattern = ":",
replacement = "_",
x = file,
fixed = TRUE
)
caption <- "Time-varying mortality and growth parameters"
plotinfo <- save_png(
plotinfo = plotinfo,
file = file,
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = caption
)
timeVaryingParmFunc(parmlabel, forecast = forecast)
dev.off()
}
}
}
}
}
} else {
message(
"Skipping timevarying quantity plots (subplot 24), most likely\n",
"because the MGparm_By_Year_after_adjustments table (report:7)\n",
"is not reported in the Report.sso file."
)
}
}
# plot of vectors found in models with hermaphroditism
if ("Herma_Trans" %in% names(growdatF)) {
# Hermaphro_Option from SS3 User Manual
# 1 = invoke female to male age-specific function; and
# -1 = invoke male to female age-specific function.
if (any(!is.na(growdatF[["Herma_Trans"]]))) {
Hermaphro_Option <- 1
labels[13] <- paste(labels[13], "(F to M)")
} else {
Hermaphro_Option <- -1
labels[13] <- paste(labels[13], "(M to F)")
}
herma_func1 <- function() {
# if females change to males, there will be numerical values here:
if (Hermaphro_Option == 1) {
df <- growdatF
} else {
df <- growdatM
}
# make plot (age vector should be the same for females and males)
plot(
df[["Age_Beg"]],
df[["Herma_Trans"]],
xaxs = "i",
ylim = c(0, 1),
las = 1,
xlab = labels[12],
ylab = labels[13],
type = "l",
lwd = 3,
col = colvec[2]
)
abline(h = c(0, 1), col = "grey")
}
herma_func2 <- function(area = 1) {
# new area-specific columns in 3.30.21 have format "sex_ratio_area:1"
# assuming that there will always be one column per area
colnames <- grep("sex_ratio", names(growdatF), value = TRUE)
matplot(
growdatF[["Age_Beg"]],
growdatF[, colnames],
xaxs = "i",
ylim = c(0, 1),
las = 1,
lty = 1:nareas,
xlab = labels[2],
ylab = labels[14],
type = "l",
lwd = 3,
col = areacols
)
# add grey line at 1.0
abline(h = c(0, 1), col = "grey")
# add legend if multiple lines per area
if (length(colnames) > 1) {
legend(
legendloc,
bty = "n",
col = areacols,
lty = 1:nareas,
lwd = 3,
legend = paste("area", 1:nareas)
)
}
}
# set default colors if not specified
areacols <- get_areacols(areacols, nareas)
if (31 %in% subplots) {
if (plot) {
herma_func1()
}
if (print) {
plotinfo <- save_png(
plotinfo = plotinfo,
file = "bio31_hermaphrodite_transition.png",
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = labels[13]
)
herma_func1()
dev.off()
}
}
if (32 %in% subplots) {
if (plot) {
herma_func2()
}
if (print) {
plotinfo <- save_png(
plotinfo = plotinfo,
file = "bio32_hermaphrodite_cumulative.png",
plotdir = plotdir,
pwidth = pwidth,
pheight = pheight,
punits = punits,
res = res,
ptsize = ptsize,
caption = labels[14]
)
herma_func2()
dev.off()
}
}
}
## ### Ian T.: finish this stuff
## # Matrix of M-at-age if mortality varies with age and time
## # strip off forecast years in M-at-age matrix
## M_at_age <- M_at_age[M_at_age[["Year"]] <= endyr,]
## if(any(M_at_age[["Bio_Pattern"]]!=1)){}
## for(i in unique(M_at_age[["Bio_Pattern"]])){
## for(isex in unique(
## M_at_age_tmpYear <- M_at_age[["Year"]]
## M_at_age_ages <- 0:accuage
## M_at_age_matrix <- as.matrix(M_at_age[,-
#### this stuff can be canibalized to get the M-at-age matrix working
## if(growthvaries){ # if growth is time varying
## for(i in 1:nsexes)
## {
## growdatuse <- growthseries[growthseries[["Yr"]] >= startyr-2 &
## growthseries[["Morph"]]==morphs[i],]
## x <- 0:accuage
## y <- growdatuse[["Yr"]]
## z <- as.matrix(growdatuse[,-(1:4)])
## time <- FALSE
## for(t in 1:ncol(z)) if(max(z[,t])!=min(z[,t])) time <- TRUE
## if(time)
## {
## z <- t(z)
## if(i==1){main <- "Female time-varying growth"}
## if(nsexes==1){main <- "Time-varying growth"}
## if(i==2){main <- "Male time-varying growth"}
## if(nseasons > 1){main <- paste(main," season 1",sep="")}
## if(plot){
## if(12 %in% subplots)
## persp(x,y,z,col="white",xlab=labels[2],ylab="",zlab=labels[1],expand=0.5,
## box=TRUE,main=main,cex.main=cex.main,ticktype="detailed",
## phi=35,theta=-10)
## if(13 %in% subplots)
## contour(x,y,z,nlevels=12,xlab=labels[2],
## main=main,cex.main=cex.main,col=ians_contour,lwd=2)}
## if(print){
## if(12 %in% subplots){
## file <- paste("bio12_timevarygrowthsurf_sex",i,".png",sep="")
## caption <- "Perspective plot of time-varying growth"
## plotinfo <- save_png(file=file, caption=caption)
## persp(x,y,z,col="white",xlab=labels[2],ylab="",zlab=labels[1],expand=0.5,
## box=TRUE,main=main,cex.main=cex.main,ticktype="detailed",
## phi=35,theta=-10)
## dev.off()
## }
## if(13 %in% subplots){
## file <- paste("bio13_timevarygrowthcontour_sex",i,".png",sep="")
## caption <- "Contour plot of time-varying growth"
## plotinfo <- save_png(file=file, caption=caption)
## contour(x,y,z,nlevels=12,xlab=labels[2],
## main=main,cex.main=cex.main,col=ians_contour,lwd=2)
## dev.off()
## }
## } # end print
## } # end if time-varying
## } # end loop over sexes
## } # end of if data available for time varying growth
## }# end disable of time-varying growth for multi-season models
# add category and return plotinfo
if (!is.null(plotinfo)) {
plotinfo[["category"]] <- "Bio"
}
return(invisible(plotinfo))
}
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.