#' Summarize the output from multiple Stock Synthesis models.
#'
#' Summarize various quantities from the model output collected by
#' [SSgetoutput()] and return them in a list of tables and vectors.
#'
#'
#' @param biglist A list of lists, one for each model. The individual lists can
#' be created by [SS_output()] or the list of lists can be
#' created by [SSgetoutput()] (which iteratively calls
#' [SS_output()]).
#' @param sizeselfactor A string or vector of strings indicating which elements
#' of the selectivity at length output to summarize. Default=c("Lsel").
#' @param ageselfactor A string or vector of strings indicating which elements
#' of the selectivity at age output to summarize. Default=c("Asel").
#' @param selfleet Vector of fleets for which selectivity will be summarized.
#' NULL=all fleets. Default=NULL.
#' @param selyr String or vector of years for which selectivity will be
#' summarized. NOTE: NOT CURRENTLY WORKING. Options: NULL=all years,
#' "startyr" = first year.
#' @param selgender Deprecated. Use selsex instead.
#' @param selsex Vector of sexes (1 and/or 2) for which selectivity will
#' be summarized. NULL=all sexes. Default=NULL.
#' @param SpawnOutputUnits Optional single value or vector of "biomass" or
#' "numbers" giving units of spawning for each model.
#' @param lowerCI Quantile for lower bound on calculated intervals. Default =
#' 0.025 for 95% intervals.
#' @param upperCI Quantile for upper bound on calculated intervals. Default =
#' 0.975 for 95% intervals.
#' @template verbose
#' @author Ian Taylor
#' @export
#' @seealso [SSgetoutput()]
SSsummarize <- function(biglist,
sizeselfactor = "Lsel",
ageselfactor = "Asel",
selfleet = NULL,
selyr = "startyr",
selgender = lifecycle::deprecated(),
selsex = 1,
SpawnOutputUnits = NULL,
lowerCI = 0.025,
upperCI = 0.975,
verbose = TRUE) {
if (lifecycle::is_present(selgender)) {
lifecycle::deprecate_warn(
when = "1.45.0",
what = "SSsummarize(selgender)",
with = "SSsummarize(selsex)"
)
selsex <- selgender
}
# confirm that biglist is a list of lists, each created by SS_output()
# this could be improved with the use of S3 classes in the future
if (!is.list(biglist) | # if whole thing is not a list
!is.list(biglist[[1]]) | # or if the first element isn't also a list
!"parameters" %in% names(biglist[[1]])) { # or if 1st list seems wrong
stop(
"Input 'biglist' needs to be a list of the lists returned by ",
"SS_output(), either by grouping those lists within 'list()', or ",
"running SSgetoutput() which calls SS_output() repeatedly ",
"and returning a big list in the appropriate format."
)
}
# loop over outputs to create list of parameters, derived quantities, and years
parnames <- NULL
dernames <- NULL
likenames <- NULL
allyears <- NULL
# figure out how many models and give them names if they don't have them already
n <- length(biglist)
modelnames <- names(biglist)
if (is.null(modelnames)) {
modelnames <- paste0("model", 1:n)
}
# accumulator age for each model
accuages <- rep(NA, n)
# do the loop
for (imodel in 1:n) {
stats <- biglist[[imodel]]
parnames <- union(parnames, stats[["parameters"]][["Label"]])
dernames <- union(dernames, stats[["derived_quants"]][["Label"]])
allyears <- union(allyears, stats[["timeseries"]][["Yr"]])
likenames <- union(likenames, rownames(stats[["likelihoods_used"]]))
# accumulator age for each model
accuages[imodel] <- stats[["accuage"]]
}
allyears <- sort(allyears) # not actually getting any timeseries stuff yet
# objects to store quantities
pars <- parsSD <- parphases <- par_prior_likes <-
as.data.frame(matrix(NA, nrow = length(parnames), ncol = n))
quants <- quantsSD <- as.data.frame(matrix(NA, nrow = length(dernames), ncol = n))
maxgrad <- NULL
nsexes <- NULL
likelihoods <- likelambdas <- as.data.frame(matrix(NA, nrow = length(likenames), ncol = n))
likelihoods_by_fleet <- NULL
likelihoods_by_tag_group <- NULL
indices <- NULL
sizesel <- NULL
agesel <- NULL
accuage <- max(accuages)
if (all(accuages == accuage)) {
growth <- as.data.frame(matrix(NA, nrow = accuage + 1, ncol = n))
names(growth) <- modelnames
} else {
warning(
"problem summarizing growth due to different ",
"accumulator ages among models"
)
growth <- NULL
}
# notes about what runs were used
sim <- NULL
listnames <- NULL
npars <- NULL
startyrs <- NULL
endyrs <- NULL
SPRratioLabels <- NULL
FvalueLabels <- NULL
sprtargs <- NULL
btargs <- NULL
minbthreshs <- NULL
FleetNames <- list()
mcmc <- list()
warn <- FALSE # flag for whether filter warning has been printed or not
if (verbose) {
message("Summarizing ", n, " models:")
}
# loop over models within biglist
for (imodel in 1:n) {
stats <- biglist[[imodel]]
listname <- names(biglist)[imodel]
if (verbose) {
message("imodel=", imodel, "/", n)
}
# gradient
maxgrad <- c(maxgrad, stats[["maximum_gradient_component"]])
# nsexes
nsexes <- c(nsexes, stats[["nsexes"]])
# start and end years
startyrs <- c(startyrs, stats[["startyr"]])
endyrs <- c(endyrs, stats[["endyr"]])
# size selectivity
sizeseltemp <- stats[["sizeselex"]]
# check for non-NULL selectivity table
if (is.null(sizeseltemp)) {
if (verbose) {
message(" no selectivity-at-length output")
}
} else {
# if factor(s) not input, get all unique values from table
if (is.null(sizeselfactor) & !is.null(sizeseltemp)) {
sizeselfactor <- unique(sizeseltemp[["Factor"]])
}
# loop over factor(s) input by user or taken from table
for (iselfactor in seq_along(sizeselfactor)) {
seltemp_i <- sizeseltemp[sizeseltemp[["Factor"]] == sizeselfactor[iselfactor], ]
seltemp_i[["imodel"]] <- imodel
seltemp_i[["name"]] <- modelnames[imodel]
# if sizesel is not NULL, then check for whether columns of new addition
# match existing file
if (is.null(sizesel) || (ncol(seltemp_i) == ncol(sizesel) &&
all(names(seltemp_i) == names(sizesel)))) {
sizesel <- rbind(sizesel, seltemp_i)
} else {
warning(
"problem summarizing size selectivity due to mismatched columns ",
"(perhaps different bins)\n"
)
}
}
rownames(sizesel) <- 1:nrow(sizesel)
}
# age selectivity
ageseltemp <- stats[["ageselex"]]
# check for NULL selectivity table
if (is.null(ageseltemp)) {
if (verbose) {
message(" no selectivity-at-age output")
}
} else {
# if factor(s) not input, get all unique values from table
if (is.null(ageselfactor)) {
ageselfactor <- unique(ageseltemp[["Factor"]])
}
# loop over factor(s) input by user or taken from table
for (iselfactor in seq_along(ageselfactor)) {
seltemp_i <- ageseltemp[ageseltemp[["Factor"]] == ageselfactor[iselfactor], ]
seltemp_i[["imodel"]] <- imodel
seltemp_i[["name"]] <- modelnames[imodel]
# if agesel is not NULL, then check for whether columns of new addition
# match existing file
if (is.null(agesel) || (ncol(seltemp_i) == ncol(agesel) &&
all(names(seltemp_i) == names(agesel)))) {
agesel <- rbind(agesel, seltemp_i)
} else {
warning(
"problem summarizing age selectivity due to mismatched columns ",
"(perhaps different bins)\n"
)
}
}
rownames(agesel) <- 1:nrow(agesel)
}
## growth (females only)
if (!is.null(growth)) {
growthtemp <- stats[["growthseries"]]
# check for non-NULL growth output
if (!is.null(growthtemp)) {
# subset for the female main morph
imorphf <- stats[["morph_indexing"]][["Index"]][
stats[["morph_indexing"]][["Sex"]] == 1 &
stats[["morph_indexing"]][["Platoon"]] %in% stats[["mainmorphs"]]
]
growthtemp <- growthtemp[growthtemp[["Morph"]] == imorphf, -(1:4)]
# remove any rows with all zeros (not sure why these occur)
growthtemp <- growthtemp[apply(growthtemp, 1, sum) > 0, ]
# get last row and bind to values from previous models
if (nrow(growthtemp) > 0) {
growth[, imodel] <- as.numeric(growthtemp[nrow(growthtemp), ])
} else {
growth[, imodel] <- NA
}
}
}
## likelihoods (total by component)
liketemp <- stats[["likelihoods_used"]]
for (irow in 1:nrow(liketemp)) {
likelihoods[likenames == rownames(liketemp)[irow], imodel] <- liketemp[["values"]][irow]
likelambdas[likenames == rownames(liketemp)[irow], imodel] <- liketemp[["lambdas"]][irow]
}
## likelihoods by fleet
# add initial column with model number to table from each model
liketemp2 <- data.frame(model = imodel, stats[["likelihoods_by_fleet"]])
# test for presence of existing table to append to with matching number of columns
if (is.null(likelihoods_by_fleet) ||
(ncol(likelihoods_by_fleet) == ncol(liketemp2) &&
all(names(likelihoods_by_fleet) == names(liketemp2)))) {
likelihoods_by_fleet <- rbind(likelihoods_by_fleet, liketemp2)
} else {
likelihoods_by_fleet <- merge(likelihoods_by_fleet, liketemp2, all = TRUE)
}
## likelihoods by tag group
# add initial column with model number to table from each model
if (!is.null(stats[["likelihoods_by_tag_group"]])) {
liketemp3 <- data.frame(model = imodel, stats[["likelihoods_by_tag_group"]])
# test for presence of existing table to append to with matching number of columns
if (is.null(likelihoods_by_tag_group) ||
(ncol(likelihoods_by_tag_group) == ncol(liketemp3) &&
all(names(likelihoods_by_tag_group) == names(liketemp3)))) {
likelihoods_by_tag_group <- rbind(likelihoods_by_tag_group, liketemp3)
} else {
warning("problem summarizing likelihoods by fleet due to mismatched columns")
}
}
## compile parameters
parstemp <- stats[["parameters"]]
for (ipar in 1:nrow(parstemp)) {
pars[parnames == parstemp[["Label"]][ipar], imodel] <- parstemp[["Value"]][ipar]
parsSD[parnames == parstemp[["Label"]][ipar], imodel] <- parstemp[["Parm_StDev"]][ipar]
parphases[parnames == parstemp[["Label"]][ipar], imodel] <- parstemp[["Phase"]][ipar]
par_prior_likes[parnames == parstemp[["Label"]][ipar], imodel] <- parstemp[["Pr_Like"]][ipar]
}
if (verbose) {
message(" N active pars = ", sum(!is.na(parstemp[["Active_Cnt"]])))
}
## compile derived quantities
quantstemp <- stats[["derived_quants"]]
for (iquant in 1:nrow(quantstemp)) {
quants[dernames == quantstemp[["Label"]][iquant], imodel] <- quantstemp[["Value"]][iquant]
quantsSD[dernames == quantstemp[["Label"]][iquant], imodel] <- quantstemp[["StdDev"]][iquant]
}
SPRratioLabels <- c(SPRratioLabels, stats[["SPRratioLabel"]])
FvalueLabels <- c(FvalueLabels, stats[["F_report_basis"]])
sprtargs <- c(sprtargs, stats[["sprtarg"]])
btargs <- c(btargs, stats[["btarg"]])
minbthreshs <- c(minbthreshs, stats[["minbthresh"]])
FleetNames[[imodel]] <- stats[["FleetNames"]]
## indices
indextemp <- stats[["cpue"]]
if (is.null(indextemp) || is.na(indextemp[[1]][1])) {
if (verbose) {
message(" no index data")
}
} else {
# temporarily remove columns added in SS version 3.30.13 (March 2019)
indextemp <- indextemp[!names(indextemp) %in% c("Area", "Subseas", "Month")]
indextemp[["name"]] <- modelnames[imodel]
indextemp[["imodel"]] <- imodel
if (is.null(indices)) {
# first pass through with nothing in combined data frame
indices <- rbind(indices, indextemp)
} else {
# after indices contains output from at least one model
# check that there are equal number of columns with matching names
# Working here
if (ncol(indextemp) == ncol(indices) &&
all(names(indextemp) == names(indices))) {
indices <- rbind(indices, indextemp)
} else {
indices <- merge(indices, indextemp, all = TRUE)
}
}
}
# number of parameters
npars <- c(npars, stats[["N_estimated_parameters"]])
# 2nd fecundity parameter indicates whether spawning output is proportional to biomass
if (!is.null(SpawnOutputUnits)) {
# if 1 value is input, repeate n times
if (length(SpawnOutputUnits) == 1) SpawnOutputUnits <- rep(SpawnOutputUnits, n)
# if total doesn't currently equal n, stop everything
if (length(SpawnOutputUnits) != n) {
stop("'SpawnOutputUnits' should have length = 1 or", n)
}
} else {
# if NULL, then make vector of NA values
SpawnOutputUnits <- rep(NA, n)
}
# if NA value in vector for current model, replace with value from model
if (is.na(SpawnOutputUnits[imodel]) & !is.null(stats[["SpawnOutputUnits"]])) {
SpawnOutputUnits[imodel] <- stats[["SpawnOutputUnits"]]
}
# get mcmc values if present
if (!is.null(stats[["mcmc"]])) {
mcmc[[imodel]] <- stats[["mcmc"]]
}
} # end loop over models
### format and process info from the models
names(pars) <- names(parsSD) <- names(parphases) <- names(par_prior_likes) <-
modelnames
names(quants) <- names(quantsSD) <- modelnames
names(likelihoods) <- names(likelambdas) <- modelnames
pars[["Label"]] <- parsSD[["Label"]] <- parphases[["Label"]] <-
par_prior_likes[["Label"]] <- parnames
quants[["Label"]] <- quantsSD[["Label"]] <- dernames
likelihoods[["Label"]] <- likelambdas[["Label"]] <- likenames
# extract year values from labels for some parameters associated with years
pars[["Yr"]] <- NA
for (ipar in 1:nrow(pars)) {
substrings <- strsplit(as.character(pars[["Label"]][ipar]), "_")[[1]]
yr <- substrings[substrings %in% allyears][1]
pars[["Yr"]][ipar] <- ifelse(is.null(yr), NA, as.numeric(yr))
}
quants[["Yr"]] <- quantsSD[["Yr"]] <- NA
for (iquant in 1:nrow(quants)) {
substrings <- strsplit(as.character(quants[["Label"]][iquant]), "_")[[1]]
yr <- substrings[substrings %in% allyears][1]
quants[["Yr"]][iquant] <- ifelse(is.null(yr), NA, as.numeric(yr))
quantsSD[["Yr"]][iquant] <- ifelse(is.null(yr), NA, as.numeric(yr))
}
# rows numbers of derived quantities that start with "SSB_"
SSBrows <- grep("SSB_", quants[["Label"]])
# row numbers that start with "SSB_" but are not part of time series
SSBexclude <- c(
grep("SSB_unfished", quants[["Label"]], ignore.case = TRUE),
grep("SSB_Btgt", quants[["Label"]], ignore.case = TRUE),
grep("SSB_SPR", quants[["Label"]], ignore.case = TRUE),
grep("SSB_MSY", quants[["Label"]], ignore.case = TRUE),
grep("SSB_F01", quants[["Label"]], ignore.case = TRUE)
)
# filter rows to only include time series
SSBrows <- setdiff(SSBrows, SSBexclude)
# identify spawning biomass parameters
SpawnBio <- quants[SSBrows, ]
SpawnBioSD <- quantsSD[SSBrows, ]
# add year values for Virgin and Initial years
minyr <- min(SpawnBio[["Yr"]], na.rm = TRUE)
SpawnBio[["Yr"]][grep("SSB_Virgin", SpawnBio[["Label"]])] <- minyr - 2
SpawnBio[["Yr"]][grep("SSB_Initial", SpawnBio[["Label"]])] <- minyr - 1
SpawnBioSD[["Yr"]] <- SpawnBio[["Yr"]]
SpawnBio <- SpawnBio[order(SpawnBio[["Yr"]]), ]
SpawnBioSD <- SpawnBioSD[order(SpawnBioSD[["Yr"]]), ]
SpawnBioLower <- SpawnBioUpper <- SpawnBioSD
SpawnBioLower[, 1:n] <- qnorm(
p = lowerCI, mean = as.matrix(SpawnBio[, 1:n]),
sd = as.matrix(SpawnBioSD[, 1:n])
)
SpawnBioUpper[, 1:n] <- qnorm(
p = upperCI, mean = as.matrix(SpawnBio[, 1:n]),
sd = as.matrix(SpawnBioSD[, 1:n])
)
# identify biomass ratio parameters
Bratio <- quants[grep("^Bratio_", quants[["Label"]]), ]
BratioSD <- quantsSD[grep("^Bratio_", quantsSD[["Label"]]), ]
BratioLower <- BratioUpper <- BratioSD
BratioLower[, 1:n] <- qnorm(
p = lowerCI, mean = as.matrix(Bratio[, 1:n]),
sd = as.matrix(BratioSD[, 1:n])
)
BratioUpper[, 1:n] <- qnorm(
p = upperCI, mean = as.matrix(Bratio[, 1:n]),
sd = as.matrix(BratioSD[, 1:n])
)
# identify SPR ratio derived quantities
SPRratio <- quants[grep("^SPRratio_", quants[["Label"]]), ]
SPRratioSD <- quantsSD[grep("^SPRratio_", quantsSD[["Label"]]), ]
SPRratioLower <- SPRratioUpper <- SPRratioSD
SPRratioLower[, 1:n] <- qnorm(
p = lowerCI, mean = as.matrix(SPRratio[, 1:n]),
sd = as.matrix(SPRratioSD[, 1:n])
)
SPRratioUpper[, 1:n] <- qnorm(
p = upperCI, mean = as.matrix(SPRratio[, 1:n]),
sd = as.matrix(SPRratioSD[, 1:n])
)
# identify F derived quantities
Fvalue <- quants[grep("^F_", quants[["Label"]]), ]
FvalueSD <- quantsSD[grep("^F_", quantsSD[["Label"]]), ]
FvalueLower <- FvalueUpper <- FvalueSD
FvalueLower[, 1:n] <- qnorm(
p = lowerCI, mean = as.matrix(Fvalue[, 1:n]),
sd = as.matrix(FvalueSD[, 1:n])
)
FvalueUpper[, 1:n] <- qnorm(
p = upperCI, mean = as.matrix(Fvalue[, 1:n]),
sd = as.matrix(FvalueSD[, 1:n])
)
# identify recruitment parameters and their uncertainty
recruits <- quants[grep("^Recr_", quants[["Label"]]), ]
recruitsSD <- quantsSD[grep("^Recr_", quantsSD[["Label"]]), ]
if (length(grep("Recr_Unfished", recruits[["Label"]], ignore.case = TRUE)) > 0) {
recruits <- recruits[-grep("Recr_Unfished", recruits[["Label"]], ignore.case = TRUE), ]
recruitsSD <- recruitsSD[-grep("Recr_Unfished", recruitsSD[["Label"]], ignore.case = TRUE), ]
}
minyr <- min(recruits[["Yr"]], na.rm = TRUE)
recruits[["Yr"]][grep("Recr_Virgin", recruits[["Label"]])] <- minyr - 2
recruits[["Yr"]][grep("Recr_Initial", recruits[["Label"]])] <- minyr - 1
recruitsSD[["Yr"]][grep("Recr_Virgin", recruitsSD[["Label"]])] <- minyr - 2
recruitsSD[["Yr"]][grep("Recr_Initial", recruitsSD[["Label"]])] <- minyr - 1
recruits <- recruits[order(recruits[["Yr"]]), ]
recruitsSD <- recruitsSD[order(recruitsSD[["Yr"]]), ]
recruitsLower <- recruitsUpper <- recruitsSD
# assume recruitments have log-normal distribution
# from first principals (multiplicative survival probabilities)
# and from their basis as exponential of normal recdevs
sdlog <- sqrt(log(1 + (as.matrix(recruitsSD[, 1:n]) /
as.matrix(recruits[, 1:n]))^2))
recruitsLower[, 1:n] <- qlnorm(
p = lowerCI,
meanlog = log(as.matrix(recruits[, 1:n])),
sdlog = sdlog
)
recruitsUpper[, 1:n] <- qlnorm(
p = upperCI,
meanlog = log(as.matrix(recruits[, 1:n])),
sdlog = sdlog
)
# identify parameters that are recruitment deviations
pars[["recdev"]] <- FALSE
pars[["recdev"]][grep("RecrDev", pars[["Label"]])] <- TRUE
pars[["recdev"]][grep("InitAge", pars[["Label"]])] <- TRUE
pars[["recdev"]][grep("ForeRecr", pars[["Label"]])] <- TRUE
# if there are any initial age parameters, figure out what year they're from
InitAgeRows <- grep("InitAge", pars[["Label"]])
if (length(InitAgeRows) > 0) {
# separate out values from string
temp <- unlist(strsplit(pars[["Label"]][InitAgeRows], "InitAge_"))
# get odd entries in above separation
InitAgeVals <- as.numeric(temp[seq(2, length(temp), 2)])
# make empty matrix to store values
InitAgeYrs <- matrix(NA, nrow = length(InitAgeRows), ncol = n)
# loop over models
for (imodel in 1:n) {
# get parameters
modelpars <- pars[, imodel]
# get vector of years associated with recdevs
devyears <- pars[["Yr"]][!is.na(modelpars) & pars[["recdev"]]]
# figure out first year of recdevs that already have years figured out
if (any(!is.na(devyears))) {
minyr <- min(devyears, na.rm = TRUE)
} else {
minyr <- NA
}
# determine which parameter values are associated with InitAge and not NA
good <- !is.na(modelpars[InitAgeRows])
# if minyr was not NA, and is above 0 and there are good InitAge values,
# then compute the associated year
if (!is.na(minyr) & minyr > 0 & any(good)) {
InitAgeYrs[good, imodel] <- minyr - InitAgeVals[good]
}
}
# check for differences in assignment of initial ages
if (any(apply(InitAgeYrs, 1, max, na.rm = TRUE) -
apply(InitAgeYrs, 1, min, na.rm = TRUE) != 0)) {
warning(
"years for InitAge parameters differ between models,",
"use InitAgeYrs matrix"
)
} else {
pars[["Yr"]][InitAgeRows] <- apply(InitAgeYrs, 1, max, na.rm = TRUE)
}
} else {
# no parameters seem to be associated with initial age structure
InitAgeYrs <- NA
}
if (any(pars[["recdev"]])) {
recdevs <- pars[pars[["recdev"]], ]
recdevsSD <- parsSD[pars[["recdev"]], ]
myorder <- order(recdevs[["Yr"]]) # save order for use in both values and SDs
recdevs <- recdevs[myorder, 1:(n + 2)]
recdevsSD <- recdevsSD[myorder, 1:(n + 1)]
recdevsSD[["Yr"]] <- recdevs[["Yr"]]
recdevsLower <- recdevsUpper <- recdevsSD
recdevsLower[, 1:n] <- qnorm(
p = lowerCI, mean = as.matrix(recdevs[, 1:n]),
sd = as.matrix(recdevsSD[, 1:n])
)
recdevsUpper[, 1:n] <- qnorm(
p = upperCI, mean = as.matrix(recdevs[, 1:n]),
sd = as.matrix(recdevsSD[, 1:n])
)
} else {
recdevs <- recdevsSD <- recdevsLower <- recdevsUpper <- NULL
}
# function to merge duplicate rows caused by different parameter labels
# that are associated with the same year, such as the recdev for 2016
# being called "ForeRecr_2016", "Late_RecrDev_2016", or "Main_RecrDev_2016",
# in 3 different models depending on the ending year of each model and the
# choice of recdev vector breaks
merge.duplicates <- function(x) {
if (!is.null(x)) {
if (length(unique(x[["Yr"]])) < length(x[["Yr"]])) {
# n should be number of models
n <- sum(!names(x) %in% c("Label", "Yr"))
x2 <- NULL # alternative data.frame
for (Yr in unique(x[["Yr"]])) {
x.Yr <- x[which(x[["Yr"]] == Yr), ]
if (nrow(x.Yr) == 1) {
# if only 1 row associated with this year add to new data.frame
x2 <- rbind(x2, x.Yr)
} else {
# more than 1 row associated with this year
# create empty row with matching names
newrow <- data.frame(t(rep(NA, n)),
Label = paste0("Multiple_labels_", Yr), Yr = Yr
)
names(newrow) <- names(x)
# loop over models to pick the (hopefully) unique value among rows
for (icol in 1:n) {
good <- !is.na(x.Yr[, icol])
if (sum(good) > 1) {
# warn if more than 1 value
warning("multiple recdevs values associated with year =", Yr)
}
if (sum(good) == 1) {
# put good value into new row
newrow[, icol] <- x.Yr[good, icol]
}
# if there are no good values, this model likely ends prior to Yr
}
# add new row to new data.frame
x2 <- rbind(x2, newrow)
} # end test for duplicates for particular year
} # end loop over years
} else { # end test for duplicates in general
# if no duplicates, just return data.frame
x2 <- x
}
} else { # test for is.null(x)
return(x)
}
return(x2)
}
# helper fxn b/c name of DM parameter changed
# todo: delete when these models no longer need to be maintained
copy.dm <- function(data, oldgrep = "EffN", newgrep = "theta") {
oldrows <- grep(oldgrep, data[, "Label"])
if (length(oldrows) == 0) {
return(data)
}
newrows <- grep(newgrep, data[, "Label"])
fix <- which(apply(
X = data[newrows, grep("model", colnames(data))],
MARGIN = 2, function(vector) all(is.na(vector))
))
if (length(oldrows) != length(newrows) | length(fix) == 0) {
return(data)
}
if (get("verbose", envir = parent.frame()) &
deparse(substitute(data)) == "pars") {
message(
"For model(s) ", paste(fix, collapse = ", "),
", values in 'pars', 'parsSD', 'parphases', and 'par_prior_likes' for\n",
paste(data[oldrows, "Label"], data[newrows, "Label"],
sep = " -> ", collapse = ", "
),
"\nwere copied from x -> y."
)
}
data[newrows, fix] <- data[oldrows, fix]
return(data)
}
# function to sort by year
sort.fn <- function(x) {
if (!is.null(x)) {
return(x[order(x[["Yr"]]), ])
} else {
return()
}
}
mylist <- list()
mylist[["n"]] <- n
mylist[["npars"]] <- npars
mylist[["modelnames"]] <- modelnames
mylist[["maxgrad"]] <- maxgrad
mylist[["nsexes"]] <- nsexes
mylist[["startyrs"]] <- startyrs
mylist[["endyrs"]] <- endyrs
mylist[["pars"]] <- copy.dm(pars)
mylist[["parsSD"]] <- copy.dm(parsSD)
mylist[["parphases"]] <- copy.dm(parphases)
mylist[["par_prior_likes"]] <- copy.dm(par_prior_likes)
mylist[["quants"]] <- quants
mylist[["quantsSD"]] <- quantsSD
mylist[["likelihoods"]] <- likelihoods
mylist[["likelambdas"]] <- likelambdas
mylist[["likelihoods_by_fleet"]] <- likelihoods_by_fleet
mylist[["likelihoods_by_tag_group"]] <- likelihoods_by_tag_group
mylist[["SpawnBio"]] <- sort.fn(SpawnBio)
mylist[["SpawnBioSD"]] <- sort.fn(SpawnBioSD)
mylist[["SpawnBioLower"]] <- sort.fn(SpawnBioLower)
mylist[["SpawnBioUpper"]] <- sort.fn(SpawnBioUpper)
mylist[["Bratio"]] <- sort.fn(Bratio)
mylist[["BratioSD"]] <- sort.fn(BratioSD)
mylist[["BratioLower"]] <- sort.fn(BratioLower)
mylist[["BratioUpper"]] <- sort.fn(BratioUpper)
mylist[["SPRratio"]] <- sort.fn(SPRratio)
mylist[["SPRratioSD"]] <- sort.fn(SPRratioSD)
mylist[["SPRratioLower"]] <- sort.fn(SPRratioLower)
mylist[["SPRratioUpper"]] <- sort.fn(SPRratioUpper)
mylist[["SPRratioLabels"]] <- SPRratioLabels
mylist[["Fvalue"]] <- sort.fn(Fvalue)
mylist[["FvalueSD"]] <- sort.fn(FvalueSD)
mylist[["FvalueLower"]] <- sort.fn(FvalueLower)
mylist[["FvalueUpper"]] <- sort.fn(FvalueUpper)
mylist[["FvalueLabels"]] <- FvalueLabels
mylist[["sprtargs"]] <- sprtargs
mylist[["btargs"]] <- btargs
mylist[["minbthreshs"]] <- minbthreshs
mylist[["recruits"]] <- sort.fn(recruits)
mylist[["recruitsSD"]] <- sort.fn(recruitsSD)
mylist[["recruitsLower"]] <- sort.fn(recruitsLower)
mylist[["recruitsUpper"]] <- sort.fn(recruitsUpper)
mylist[["recdevs"]] <- merge.duplicates(sort.fn(recdevs))
mylist[["recdevsSD"]] <- merge.duplicates(sort.fn(recdevsSD))
mylist[["recdevsLower"]] <- merge.duplicates(sort.fn(recdevsLower))
mylist[["recdevsUpper"]] <- merge.duplicates(sort.fn(recdevsUpper))
mylist[["growth"]] <- growth
mylist[["sizesel"]] <- sizesel
mylist[["agesel"]] <- agesel
mylist[["indices"]] <- indices
mylist[["InitAgeYrs"]] <- InitAgeYrs
mylist[["lowerCI"]] <- lowerCI
mylist[["upperCI"]] <- upperCI
mylist[["SpawnOutputUnits"]] <- SpawnOutputUnits
mylist[["FleetNames"]] <- FleetNames
mylist[["mcmc"]] <- mcmc
# mylist[["lbinspop"]] <- as.numeric(names(stats[["sizeselex"]])[-(1:5)])
if (verbose) {
message(
"Summary finished. ",
"To avoid printing details above, use 'verbose = FALSE'."
)
}
return(invisible(mylist))
} # end function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.