R/sam.R
In einarhjorleifsson/fishvice: A collection of various personal functions used in my work on fish stock assessment
Defines functions
sam_ibya
sam_get_fit
Documented in
sam_get_fit
sam_ibya
# A: sam rbx -------------------------------------------------------------------
#' sam_get_fit
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return A "sam"-object
#'
#' @export
sam_get_fit <- function(assessment) {
stockassessment::fitfromweb(assessment, character.only = TRUE)
}
#' @title sam_ibya
#'
#' @description sam input data as a tibble. Note that "tuning" variables are
#' not returned.
#'
#' @param fit A "sam"-object
#' @param scale A scaler (default 1)
#' @param long A boolean indicating if returned table wide (default TRUE) variables
#' as names within column 'var'. Alternative (FALSE) not yet active.
#' @param run Name of the run
#'
#' @return A tibble containing the following variables:
#' \itemize{
#' \item year:
#' \item age:
#' \item cW: Catch weights
#' \item dW: Discards weights
#' \item lF: Fraction of fish landed
#' \item lW: Landed weights
#' \item m: Assumed natural mortality (M)
#' \item mat: Maturity ogive
#' \item oC: Observed catch
#' \item pF: Proportional F prior to spawning
#' \item pM: Proportional M prior to spawning
#' \item sW: Spawning stock weight
#' }
#'
#' @export
#'
#'
sam_ibya <- function(fit, scale = 1, long = TRUE, run) {
if(class(fit)[[1]] != "sam")
stop('Object has to be of class "sam"')
lh <- function(x, cn) {
x <-
x %>%
as.data.frame()
names(x) <- stringr::str_replace(names(x), ".Residual catch", "")
x %>%
dplyr::mutate(year = row.names(.) %>% as.integer()) %>%
tidyr::gather(age, {{cn}}, -year, convert = TRUE) %>%
tibble::as_tibble()
}
data <- fit$data
obs <-
data$aux %>%
as.data.frame() %>%
tibble::as_tibble() %>%
dplyr::mutate(obs = exp(data$logobs))
d <-
obs %>%
dplyr::filter(fleet == 1) %>%
dplyr::select(year, age, oC = obs) %>%
dplyr::mutate(oC = oC / scale) %>%
dplyr::full_join(lh(data$catchMeanWeight, cW), by = c("year", "age")) %>%
dplyr::full_join(lh(data$disMeanWeight, dW), by = c("year", "age")) %>%
dplyr::full_join(lh(data$landFrac, lF), by = c("year", "age")) %>%
dplyr::full_join(lh(data$landMeanWeight, lW), by = c("year", "age")) %>%
dplyr::full_join(lh(data$propMat, mat), by = c("year", "age")) %>%
dplyr::full_join(lh(data$natMor, m), by = c("year", "age")) %>%
dplyr::full_join(lh(data$propF, pF), by = c("year", "age")) %>%
dplyr::full_join(lh(data$propM, pM), by = c("year", "age")) %>%
dplyr::full_join(lh(data$stockMeanWeight, sW), by = c("year", "age"))
if(!missing(run)) d$run <- run
return(d)
}
#' @title sam_rbya
#'
#' @description Extracts stock in numbers (Nay) and fishing mortality (Fay) at
#' age from object "sam". If a tibble ibya (e.g. generated by {sam_ibya} is
#' included in the argument, these estimates get added to that tibble.
#'
#' @param fit A "sam" object
#' @param data A boolean (default TRUE) specifying if input data should also
#' be returned.
#' @param scale A scaler (default 1)
#' @param long A boolean indicating if returned table wide (default TRUE) variables
#' as names within column 'var'. Alternative (FALSE) not yet active.
#' @param run Name of the run
#'
#' @return A tibble, containing at minimum:
#' \itemize{
#' \item year:
#' \item age:
#' \item n: stock in numbers
#' \item f: fishing mortality
#' }
#'
#' @export
#'
sam_rbya <- function(fit, data = TRUE, scale = 1, long = TRUE, run) {
nay <-
stockassessment::ntable(fit) %>%
as.data.frame() %>%
dplyr::mutate(year = rownames(.) %>% as.integer()) %>%
tidyr::gather(age, n, -year, convert = TRUE) %>%
dplyr::mutate(n = n / scale)
fay <-
stockassessment::faytable(fit) %>%
as.data.frame() %>%
dplyr::mutate(year = rownames(.) %>% as.integer()) %>%
tidyr::gather(age, f, -year, convert = TRUE)
res <-
nay %>%
dplyr::full_join(fay, by = c("year", "age")) %>%
tibble::as_tibble()
if(data) {
res <-
fit %>%
sam_ibya(scale = scale, long = FALSE) %>%
dplyr::full_join(res, by = c("year", "age"))
}
if(long) {
res <-
res %>%
tidyr::gather(variable, val, -c(year, age))
}
if(!missing(run)) res$run <- run
return(res)
}
#' sam_rby
#'
#' @description Get assessment summary data from "sam"-object
#'
#' @param fit XXX
#' @param scale A scaler (default 1)
#' @param run Name of the run
#'
#' @return A tibble containing the following variables:
#' \itemize{
#' \item year:
#' \item est: Medium value
#' \item low: Lower 2.5% quantile??
#' \item high: Upper 2.5% quantile??
#' \item variable: Name of variable (catch, recruitment, ssb, tsb and fbar)
#' }
#' @export
#'
sam_rby <- function(fit, scale = 1, run) {
if(class(fit) == "samset") {
base <- attributes(fit)$fit |> fv_sam_rby(scale = scale, run = run)
max.year <- max(base$year)
base$assyear <- max.year
res <- purrr::map(fit, sam_rby, scale = scale, run = run)
names(res) <- max.year - 1:length(res)
ret <-
dplyr::bind_rows(base,
res |>
dplyr::bind_rows(.id = "assyear") |>
dplyr::mutate(assyear = as.numeric(assyear)))
}
if(class(fit) == "sam") {
ret <-
fit |>
fv_sam_rby(scale = scale, run = run)
}
return(ret)
}
fv_sam_rby <- function(fit, scale = 1, run) {
lh <- function(x, variable, scale = 1) {
x %>%
as.data.frame() %>%
tibble::rownames_to_column(var = "year") %>%
dplyr::mutate(year = as.integer(year),
variable = variable) %>%
tibble::as_tibble() %>%
dplyr::mutate(Estimate = Estimate / scale,
Low = Low / scale,
High = High / scale)
}
d <-
dplyr::bind_rows(stockassessment::catchtable(fit) %>% lh("catch", scale = scale),
stockassessment::rectable(fit) %>% lh("rec", scale = scale),
stockassessment::ssbtable(fit) %>% lh("ssb", scale = scale),
stockassessment::tsbtable(fit) %>% lh("tsb", scale = scale),
stockassessment::fbartable(fit) %>% lh("fbar")) %>%
dplyr::rename(est = Estimate,
low = Low,
high = High)
if(!missing(run)) d$run <- run
return(d)
}
#' Sam fit
#'
#' Observed and predicted values
#'
#' @param fit A "sam" object
#' @param scale A scale for the scaleable variables
#' @param run Name of the run
#'
#' @return A list of length 2 containing the following:
#' \itemize{
#' \item data: A tibble with the following variables:
#' \itemize{
#' \item year:
#' \item age:
#' \item fleet: Name of the fleets
#' \item o: Observed value, default log scale
#' \item p: Predicted value, default log scale
#' \item r: Residuals, not yet returned
#' }
#' }
#'
#' @export
#'
sam_opr <- function(fit, scale = 1, run) {
# code from stockassessment
fleets <- unique(fit$data$aux[,"fleet"])
log <- TRUE
idx <- fit$data$aux[,"fleet"] %in% fleets
trans <- function(x) if(log){x}else{exp(x)}
#p <- trans(fit$obj$report(c(fit$sdrep$par.fixed,fit$sdrep$par.random))$predObs[idx])
p <- fit$rep$predObs
o <- trans(fit$data$logobs[idx])
aa <- fit$data$aux[idx,"age"]
neg.age <- (aa < -1.0e-6)
aa[neg.age] <- NA
a <- paste0("a=",aa," ")
f <- paste0(" f=",strtrim(attr(fit$data,"fleetNames")[fit$data$aux[idx,"fleet"]],50))
Year <- fit$data$aux[idx,"year"]
if(length(fleets)==1){
myby <- paste(a, ":", f)
}else{
myby <- cbind(a,f)
}
fleet <- strtrim(attr(fit$data,"fleetNames")[fit$data$aux[idx,"fleet"]],50)
d <-
tibble::tibble(year = Year,
age = aa,
fleet = fleet,
o = o,
p = p) %>%
dplyr::mutate(fleet = ifelse(fleet == "Residual catch", "catch", fleet))
if(!missing(run)) d$run <- run
return(d)
}
sam_fleets <- function(fit) {
tibble::tibble(fleet_nr = sort(unique(fit$data$aux[,"fleet"])),
fleet_name = attr(fit$data,"fleetNames"),
keyBiomassTreat = fit$conf$keyBiomassTreat,
obsLikelihoodFlag = as.character(fit$conf$obsLikelihoodFlag),
obsCorStruct = as.character(fit$conf$obsCorStruct)) %>%
dplyr::mutate(fleet_name = ifelse(fleet_nr == 1, "catch", fleet_name))
}
sam_partable <- function(fit, run) {
lu <-
tibble::tribble(#~name, ~par_conf,
~out_name, ~in_name,
"logFpar", "keyLogFpar",
"logQpow", "keyQpow",
"logSdLogFsta","keyVarF", # sd random walk F
"logSdLogN", "keyVarLogN",
"logSdLogObs", "keyVarObs",
"transfIRARdist", "",
"logitReleaseSurvival", "",
"logitReleaseSurvival", "",
"logitRecapturePhi", "")
d <-
fit %>%
stockassessment:::partable() %>%
as.data.frame() %>%
tibble::as_tibble(rownames = "name") %>%
janitor::clean_names() %>%
dplyr::mutate(key = as.integer(stringr::str_replace(name, "^.+_", "")),
out_name = stringr::str_sub(name, 1, nchar(name) - nchar(key) - 1)) %>%
dplyr::left_join(lu, by = "out_name") %>%
dplyr::rename(sd = sd_par, est = exp_par) %>%
dplyr::left_join(sam_conf_tbl(fit) %>% dplyr::rename(in_name = par_conf),
by = c("key", "in_name")) %>%
dplyr::left_join(sam_fleets(fit),
by = "fleet_nr") %>%
dplyr::mutate(what = dplyr::case_when(out_name == "logFpar" ~ "catchabilities",
out_name == "logQpow" ~ "powers",
# NOTE: need to check output name
out_name == "logSdLogObs" ~ "obsvar",
out_name == "logSdLogN" ~ "process",
out_name == "logSdLogFsta" ~ "rwalkF",
TRUE ~ "rest")) %>%
dplyr::select(fleet = fleet_name, age, m = par, cv = sd, est, low, high, what, in_name, out_name)
d <-
d %>%
dplyr::mutate(age2 = stringr::str_pad(age, 2, pad = "0"),
fleet2 = stringr::str_sub(fleet, 1, 10),
fleetage = dplyr::case_when(!is.na(fleet2) & !is.na(age2) ~ paste0(fleet2, "_", age2),
is.na(fleet2) & !is.na(age2) ~ age2,
is.na(fleet2) & is.na(age2) ~ out_name,
TRUE ~ NA_character_)) %>%
dplyr::select(-c(age2, fleet2))
if(!missing(run)) d$run <- run
return(d)
}
sam_conf_tbl <- function(fit) {
lh <- function(conf, what) {
x <- conf[names(conf) == what]
x <- x[[1]]
colnames(x) <- conf$minAge:conf$maxAge
x %>%
as.data.frame() %>%
tibble::as_tibble(rownames = "fleet_nr") %>%
tidyr::gather(age, key, -fleet_nr, convert = TRUE) %>%
dplyr::mutate(fleet_nr = as.integer(fleet_nr),
par_conf = what) %>%
dplyr::arrange(fleet_nr, age) %>%
dplyr::filter(key != -1)
}
dplyr::bind_rows(lh(fit$conf, "keyLogFsta"),
lh(fit$conf, "keyLogFpar"),
lh(fit$conf, "keyVarObs"),
lh(fit$conf, "keyQpow"),
lh(fit$conf, "keyVarF")) %>%
dplyr::select(par_conf, dplyr::everything()) %>%
dplyr::bind_rows(tibble::tibble(par_conf = "keyVarLogN",
fleet_nr = -1,
age = fit$conf$minAge:fit$conf$maxAge,
key = fit$conf$keyVarLogN))
}
#' rbx tibbles from object class "sam"
#'
#' @param fit A "sam"-object
#' @param scale A scaler (default 1)
#'
#' @return A list containing tibbles "rby", "rbya", "opr" and "par"
#'
#' @export
#'
sam_rbx <- function(fit, scale = 1) {
list(rby = sam_rby(fit, scale = scale),
rbya = sam_rbya(fit, data = TRUE, scale = scale),
opr = sam_opr(fit),
par = sam_partable(fit))
}
# B: FLRSAM --------------------------------------------------------------------
# B: residuals, work in progress -----------------------------------------------
#' One-observation-ahead residuals
#'
#' @description Note, this is normally just plotted using plot(res)
#'
#' @param res An object returned from {stockassessment::residuals}
#'
#' @return A tibble, the "one-observation-ahead-residuals" is in
#' variable residual.
sam_one_observation_ahead_residuals <- function(res) {
if(class(res) != "samres") stop("Object has to be of class 'samres'")
# must be an easier way for the following:
d <- tibble::tibble(year = res$year,
fleet = res$fleet,
age = res$age,
observation = res$observation,
nll = res$nll,
grad = res$grad,
mean = res$mean,
residual = res$residual)
key <- tibble::tibble(fleet = d$fleet %>% unique(),
what = attributes(res)$fleetNames)
d %>%
dplyr::left_join(key, by = "fleet") %>%
return()
}
#empirobscorrplot(RES)
#setcap("One-observation-ahead residual correlations", "Residual correlations.")
# stockassessment:::empirobscorrplot.samres(RES)
#stockassessment:::empirobscorrplot.samres
#function (res, ...)
#{
# res <- RES
# dat <- data.frame(resid = res$residual, age = res$age, year = res$year,
# fleet = res$fleet)
# fleets <- unique(dat$fleet)
# fn <- attr(res, "fleetNames")
# str(dat)
#
# x <- list()
# for (i in 1:length(fleets)) {
# tmp <- xtabs(resid ~ age + year, data = dat[dat$fleet ==
# fleets[i], ])
# xx <- cor(t(tmp))
# x[[length(x) + 1]] <- xx
# }
#
# corplotcommon(x, fn)
# Process residuals ----------------------------------
# if(exists("RESP")){
# plot(RESP)
# setcap("Process residuals", "Standardized single-joint-sample residuals of process increments")
# stampit(fit)
# par(mfrow=c(1,1))
# }
#' Proces residuals
#'
#' @param resp An object returned from {stockassessment::procres}
#'
#' @return A tibble
sam_process_residuals <- function(resp) {
d <-
tibble::tibble(year = resp$year,
age = resp$age,
fleet = resp$fleet,
residual = resp$residual)
key <-
tibble::tibble(fleet = d$fleet %>% unique(),
fleetn = attributes(resp)$fleetNames)
d %>%
dplyr::left_join(key, by = "fleet") %>%
return()
}
# C: sam various ---------------------------------------------------------------
#' @title sam process error
#'
#' @description XXX
#'
#' @export
#'
#' @param rbya XXX
#' @param plot_it XXX
#' @param plot_catch XXX
#' @param plus_group XXX
#'
sam_process_error <- function(rbya, plot_it=FALSE, plot_catch = FALSE, plus_group=TRUE) {
last.true.age <- max(rbya$age) - plus_group
last.year <- max(rbya$year)
# # align the year-classes
# x <- rbya[,c("year","age","n")]
# x$year <- x$year - 1
# x$age <- x$age - 1
# names(x)[3] <- "n.end"
# d <- plyr::join(rbya[,c("year","age","n","m","f","oC", "cW")],x, by=c("year","age"))
# d <- d[!is.na(d$n.end),]
# # Note: could have a problem for the terminal year
# d$f[is.na(d$f)] <- 0
# # exclude plus-group
# if(plus_group) d <- d[d$age < max(d$age),]
# # process error expressed as mortality
# d$z.n <- log(d$n/d$n.end)
# d$z.f <- d$f + d$m
# d$z.d <- d$z.n - d$z.f
#
# # process error expressed as numbers
# d$n.end2 <- d$n * exp(-(d$f + d$m))
# # Calculate the difference
# d$n.d <- d$n.end - d$n.end2
d <-
rbya %>%
dplyr::mutate(yc = year - age,
f = tidyr::replace_na(f, 0)) %>%
dplyr::arrange(yc, age) %>%
dplyr::group_by(yc) %>%
# process error expressed as mortality
dplyr::mutate(n.end = dplyr::lead(n),
z.d = log(n / dplyr::lead(n)) - (f + m),
# process error expressed as numbers
n.d = dplyr::lead(n) - (n * exp(-(f + m))),
# process errror expressed as biomass
b.d = n.d * cW) %>%
dplyr::ungroup() %>%
dplyr::mutate(z.d = ifelse(age >= last.true.age, NA, z.d),
n.d = ifelse(age >= last.true.age, NA, n.d),
b.d = ifelse(age >= last.true.age, NA, b.d)) %>%
dplyr::arrange(year, age)
#
# process errror expressed as biomass
x <-
d %>%
dplyr::group_by(year) %>%
dplyr::summarise(b = sum(b.d, na.rm = TRUE),
y = sum(oC * cW, na.rm = TRUE))
if(plot_it) {
mort <-
ggplot2::ggplot(d, ggplot2::aes(year, z.d)) +
ggplot2::theme_bw() +
ggplot2::geom_text(ggplot2::aes(label = age)) +
ggplot2::stat_smooth(span = 0.1) +
ggplot2::labs(x = NULL, y = NULL,
title = "Process error expressed as deviations in mortality")
abun <-
ggplot2::ggplot(d, ggplot2::aes(year, n.d)) +
ggplot2::theme_bw() +
ggplot2::geom_col() +
ggplot2::facet_wrap(~ age, scales = "free_y") +
ggplot2::labs(x = NULL, y = NULL,
title = "Process error expressed as deviations in number of fish")
mass <- ggplot2::ggplot(x, ggplot2::aes(year, b)) +
ggplot2::theme_bw() +
ggplot2::geom_bar(stat="identity") +
ggplot2::labs(x = NULL, y = "Mass",
title = "Process error expressed as deviation in mass")
if(plot_catch) {
abun <-
abun +
ggplot2::geom_point(ggplot2::aes(y = oC),
colour = "red",
size = 0.5)
mass <-
mass +
ggplot2::geom_point(ggplot2::aes(y = y),
colour = "red",
size = 0.5)
}
return(list(rbya = d, rby = x, mort = mort, abun = abun, mass = mass))
}
return(list(rbya=d, rby=x))
}
sam_ypr <- function(fit) {
# NOTE: May want to pass arguements to ypr, like number of years, etc.
tmp <- stockassessment::ypr(fit)
res <- tibble::tibble(fbar = tmp$fbar,
yield = tmp$yield,
ssb = tmp$ssb)
ref <- tibble::tibble(ref = c("Fmax", "F01", "F35"),
fbar = c(tmp$fmax, tmp$f01, tmp$f35),
ssb = c(res$ssb[tmp$fmaxIdx],
res$ssb[tmp$f01Idx],
res$ssb[tmp$f35Idx]),
yield = c(res$yield[tmp$fmaxIdx],
res$yield[tmp$f01Idx],
res$yield[tmp$f35Idx]))
return(list(ypr = res, ref = ref))
}
#if(!all(fit$conf$obsCorStruct=="ID")){
# corplot(fit)
# setcap("Estimated correlations", "Estimates correlations between age groups for each fleet")
# stampit(fit)
#}
# CHECK OUT: stockassessment:::obscorrplot.sam
# obscov(fit, TRUE)
# Some experimental stuff
#' Sam catcability
#'
#' @param fit A "sam" object
#'
#' @return A tibble
#'
sam_catchability <- function(fit) {
Q <- fit$pl$logFpar
Qsd <- fit$plsd$logFpar
key <- fit$conf$keyLogFpar
fun <- function(x)if(x<0){NA}else{Q[x+1]}
FF <- Vectorize(fun)
ages <- fit$conf$minAge:fit$conf$maxAge
Qage <- exp(t(matrix(FF(key), nrow = nrow(key))))
d <-
Qage %>%
tidyr::as_tibble()
names(d) <- attr(fit$data, "fleetNames")[1:nrow(key)]
d$age <- ages
d %>%
tidyr::gather(fleet, value, -age) %>%
tidyr::drop_na()
}
# ------------------------------------------------------------------------------
# Get stuff from www.stockassessment.org
#' @title Get a sam directory from stockassessment.org
#'
#' @description The function copies the whole directory of an assessment run from
#' stockassessment.org to a local directory
#'
#' @param assessment Name of the assessment
#' @param user Name of the user (default "user3")
#'
sam_get_directory <- function(assessment, user = "user3") {
path <- paste("https://www.stockassessment.org/datadisk/stockassessment/userdirs",user,assessment,sep="/")
cmd <- paste0("wget --recursive --reject=png,html --level=0 --no-parent ",path,"/")
system(cmd)
# cleanup
Path <- "www.stockassessment.org/datadisk/stockassessment/userdirs"
Path <- paste(Path,user,assessment,sep="/")
cmd <- paste("mv", Path, ".")
system(cmd)
system("rm -r www.stockassessment.org")
}
#' sam_get_data
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return A list
#'
sam_get_data <- function(assessment) {
dat <- NA
load(url(sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run/data.RData")))
return(dat)
}
#' sam_get_residuals
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return XXXX
#'
sam_get_residuals <- function(assessment) {
fil <- sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run/residuals.RData")
if(!RCurl::url.exists(fil)) {
stop(paste0("File: 'residuals.RData for ",
assessment,
" does not exist at: ",
sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run"))
)
}
RES <- NA
load(url(fil))
return(RES)
}
#' sam_get_retro
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return XXXX
sam_get_retro <- function(assessment) {
RETRO <- NA
load(url(sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run/retro.RData")))
return(RETRO)
}
einarhjorleifsson/fishvice documentation built on Jan. 4, 2024, 8:43 p.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 sam_ibya sam_get_fit
Documented in sam_get_fit sam_ibya
# A: sam rbx -------------------------------------------------------------------
#' sam_get_fit
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return A "sam"-object
#'
#' @export
sam_get_fit <- function(assessment) {
stockassessment::fitfromweb(assessment, character.only = TRUE)
}
#' @title sam_ibya
#'
#' @description sam input data as a tibble. Note that "tuning" variables are
#' not returned.
#'
#' @param fit A "sam"-object
#' @param scale A scaler (default 1)
#' @param long A boolean indicating if returned table wide (default TRUE) variables
#' as names within column 'var'. Alternative (FALSE) not yet active.
#' @param run Name of the run
#'
#' @return A tibble containing the following variables:
#' \itemize{
#' \item year:
#' \item age:
#' \item cW: Catch weights
#' \item dW: Discards weights
#' \item lF: Fraction of fish landed
#' \item lW: Landed weights
#' \item m: Assumed natural mortality (M)
#' \item mat: Maturity ogive
#' \item oC: Observed catch
#' \item pF: Proportional F prior to spawning
#' \item pM: Proportional M prior to spawning
#' \item sW: Spawning stock weight
#' }
#'
#' @export
#'
#'
sam_ibya <- function(fit, scale = 1, long = TRUE, run) {
if(class(fit)[[1]] != "sam")
stop('Object has to be of class "sam"')
lh <- function(x, cn) {
x <-
x %>%
as.data.frame()
names(x) <- stringr::str_replace(names(x), ".Residual catch", "")
x %>%
dplyr::mutate(year = row.names(.) %>% as.integer()) %>%
tidyr::gather(age, {{cn}}, -year, convert = TRUE) %>%
tibble::as_tibble()
}
data <- fit$data
obs <-
data$aux %>%
as.data.frame() %>%
tibble::as_tibble() %>%
dplyr::mutate(obs = exp(data$logobs))
d <-
obs %>%
dplyr::filter(fleet == 1) %>%
dplyr::select(year, age, oC = obs) %>%
dplyr::mutate(oC = oC / scale) %>%
dplyr::full_join(lh(data$catchMeanWeight, cW), by = c("year", "age")) %>%
dplyr::full_join(lh(data$disMeanWeight, dW), by = c("year", "age")) %>%
dplyr::full_join(lh(data$landFrac, lF), by = c("year", "age")) %>%
dplyr::full_join(lh(data$landMeanWeight, lW), by = c("year", "age")) %>%
dplyr::full_join(lh(data$propMat, mat), by = c("year", "age")) %>%
dplyr::full_join(lh(data$natMor, m), by = c("year", "age")) %>%
dplyr::full_join(lh(data$propF, pF), by = c("year", "age")) %>%
dplyr::full_join(lh(data$propM, pM), by = c("year", "age")) %>%
dplyr::full_join(lh(data$stockMeanWeight, sW), by = c("year", "age"))
if(!missing(run)) d$run <- run
return(d)
}
#' @title sam_rbya
#'
#' @description Extracts stock in numbers (Nay) and fishing mortality (Fay) at
#' age from object "sam". If a tibble ibya (e.g. generated by {sam_ibya} is
#' included in the argument, these estimates get added to that tibble.
#'
#' @param fit A "sam" object
#' @param data A boolean (default TRUE) specifying if input data should also
#' be returned.
#' @param scale A scaler (default 1)
#' @param long A boolean indicating if returned table wide (default TRUE) variables
#' as names within column 'var'. Alternative (FALSE) not yet active.
#' @param run Name of the run
#'
#' @return A tibble, containing at minimum:
#' \itemize{
#' \item year:
#' \item age:
#' \item n: stock in numbers
#' \item f: fishing mortality
#' }
#'
#' @export
#'
sam_rbya <- function(fit, data = TRUE, scale = 1, long = TRUE, run) {
nay <-
stockassessment::ntable(fit) %>%
as.data.frame() %>%
dplyr::mutate(year = rownames(.) %>% as.integer()) %>%
tidyr::gather(age, n, -year, convert = TRUE) %>%
dplyr::mutate(n = n / scale)
fay <-
stockassessment::faytable(fit) %>%
as.data.frame() %>%
dplyr::mutate(year = rownames(.) %>% as.integer()) %>%
tidyr::gather(age, f, -year, convert = TRUE)
res <-
nay %>%
dplyr::full_join(fay, by = c("year", "age")) %>%
tibble::as_tibble()
if(data) {
res <-
fit %>%
sam_ibya(scale = scale, long = FALSE) %>%
dplyr::full_join(res, by = c("year", "age"))
}
if(long) {
res <-
res %>%
tidyr::gather(variable, val, -c(year, age))
}
if(!missing(run)) res$run <- run
return(res)
}
#' sam_rby
#'
#' @description Get assessment summary data from "sam"-object
#'
#' @param fit XXX
#' @param scale A scaler (default 1)
#' @param run Name of the run
#'
#' @return A tibble containing the following variables:
#' \itemize{
#' \item year:
#' \item est: Medium value
#' \item low: Lower 2.5% quantile??
#' \item high: Upper 2.5% quantile??
#' \item variable: Name of variable (catch, recruitment, ssb, tsb and fbar)
#' }
#' @export
#'
sam_rby <- function(fit, scale = 1, run) {
if(class(fit) == "samset") {
base <- attributes(fit)$fit |> fv_sam_rby(scale = scale, run = run)
max.year <- max(base$year)
base$assyear <- max.year
res <- purrr::map(fit, sam_rby, scale = scale, run = run)
names(res) <- max.year - 1:length(res)
ret <-
dplyr::bind_rows(base,
res |>
dplyr::bind_rows(.id = "assyear") |>
dplyr::mutate(assyear = as.numeric(assyear)))
}
if(class(fit) == "sam") {
ret <-
fit |>
fv_sam_rby(scale = scale, run = run)
}
return(ret)
}
fv_sam_rby <- function(fit, scale = 1, run) {
lh <- function(x, variable, scale = 1) {
x %>%
as.data.frame() %>%
tibble::rownames_to_column(var = "year") %>%
dplyr::mutate(year = as.integer(year),
variable = variable) %>%
tibble::as_tibble() %>%
dplyr::mutate(Estimate = Estimate / scale,
Low = Low / scale,
High = High / scale)
}
d <-
dplyr::bind_rows(stockassessment::catchtable(fit) %>% lh("catch", scale = scale),
stockassessment::rectable(fit) %>% lh("rec", scale = scale),
stockassessment::ssbtable(fit) %>% lh("ssb", scale = scale),
stockassessment::tsbtable(fit) %>% lh("tsb", scale = scale),
stockassessment::fbartable(fit) %>% lh("fbar")) %>%
dplyr::rename(est = Estimate,
low = Low,
high = High)
if(!missing(run)) d$run <- run
return(d)
}
#' Sam fit
#'
#' Observed and predicted values
#'
#' @param fit A "sam" object
#' @param scale A scale for the scaleable variables
#' @param run Name of the run
#'
#' @return A list of length 2 containing the following:
#' \itemize{
#' \item data: A tibble with the following variables:
#' \itemize{
#' \item year:
#' \item age:
#' \item fleet: Name of the fleets
#' \item o: Observed value, default log scale
#' \item p: Predicted value, default log scale
#' \item r: Residuals, not yet returned
#' }
#' }
#'
#' @export
#'
sam_opr <- function(fit, scale = 1, run) {
# code from stockassessment
fleets <- unique(fit$data$aux[,"fleet"])
log <- TRUE
idx <- fit$data$aux[,"fleet"] %in% fleets
trans <- function(x) if(log){x}else{exp(x)}
#p <- trans(fit$obj$report(c(fit$sdrep$par.fixed,fit$sdrep$par.random))$predObs[idx])
p <- fit$rep$predObs
o <- trans(fit$data$logobs[idx])
aa <- fit$data$aux[idx,"age"]
neg.age <- (aa < -1.0e-6)
aa[neg.age] <- NA
a <- paste0("a=",aa," ")
f <- paste0(" f=",strtrim(attr(fit$data,"fleetNames")[fit$data$aux[idx,"fleet"]],50))
Year <- fit$data$aux[idx,"year"]
if(length(fleets)==1){
myby <- paste(a, ":", f)
}else{
myby <- cbind(a,f)
}
fleet <- strtrim(attr(fit$data,"fleetNames")[fit$data$aux[idx,"fleet"]],50)
d <-
tibble::tibble(year = Year,
age = aa,
fleet = fleet,
o = o,
p = p) %>%
dplyr::mutate(fleet = ifelse(fleet == "Residual catch", "catch", fleet))
if(!missing(run)) d$run <- run
return(d)
}
sam_fleets <- function(fit) {
tibble::tibble(fleet_nr = sort(unique(fit$data$aux[,"fleet"])),
fleet_name = attr(fit$data,"fleetNames"),
keyBiomassTreat = fit$conf$keyBiomassTreat,
obsLikelihoodFlag = as.character(fit$conf$obsLikelihoodFlag),
obsCorStruct = as.character(fit$conf$obsCorStruct)) %>%
dplyr::mutate(fleet_name = ifelse(fleet_nr == 1, "catch", fleet_name))
}
sam_partable <- function(fit, run) {
lu <-
tibble::tribble(#~name, ~par_conf,
~out_name, ~in_name,
"logFpar", "keyLogFpar",
"logQpow", "keyQpow",
"logSdLogFsta","keyVarF", # sd random walk F
"logSdLogN", "keyVarLogN",
"logSdLogObs", "keyVarObs",
"transfIRARdist", "",
"logitReleaseSurvival", "",
"logitReleaseSurvival", "",
"logitRecapturePhi", "")
d <-
fit %>%
stockassessment:::partable() %>%
as.data.frame() %>%
tibble::as_tibble(rownames = "name") %>%
janitor::clean_names() %>%
dplyr::mutate(key = as.integer(stringr::str_replace(name, "^.+_", "")),
out_name = stringr::str_sub(name, 1, nchar(name) - nchar(key) - 1)) %>%
dplyr::left_join(lu, by = "out_name") %>%
dplyr::rename(sd = sd_par, est = exp_par) %>%
dplyr::left_join(sam_conf_tbl(fit) %>% dplyr::rename(in_name = par_conf),
by = c("key", "in_name")) %>%
dplyr::left_join(sam_fleets(fit),
by = "fleet_nr") %>%
dplyr::mutate(what = dplyr::case_when(out_name == "logFpar" ~ "catchabilities",
out_name == "logQpow" ~ "powers",
# NOTE: need to check output name
out_name == "logSdLogObs" ~ "obsvar",
out_name == "logSdLogN" ~ "process",
out_name == "logSdLogFsta" ~ "rwalkF",
TRUE ~ "rest")) %>%
dplyr::select(fleet = fleet_name, age, m = par, cv = sd, est, low, high, what, in_name, out_name)
d <-
d %>%
dplyr::mutate(age2 = stringr::str_pad(age, 2, pad = "0"),
fleet2 = stringr::str_sub(fleet, 1, 10),
fleetage = dplyr::case_when(!is.na(fleet2) & !is.na(age2) ~ paste0(fleet2, "_", age2),
is.na(fleet2) & !is.na(age2) ~ age2,
is.na(fleet2) & is.na(age2) ~ out_name,
TRUE ~ NA_character_)) %>%
dplyr::select(-c(age2, fleet2))
if(!missing(run)) d$run <- run
return(d)
}
sam_conf_tbl <- function(fit) {
lh <- function(conf, what) {
x <- conf[names(conf) == what]
x <- x[[1]]
colnames(x) <- conf$minAge:conf$maxAge
x %>%
as.data.frame() %>%
tibble::as_tibble(rownames = "fleet_nr") %>%
tidyr::gather(age, key, -fleet_nr, convert = TRUE) %>%
dplyr::mutate(fleet_nr = as.integer(fleet_nr),
par_conf = what) %>%
dplyr::arrange(fleet_nr, age) %>%
dplyr::filter(key != -1)
}
dplyr::bind_rows(lh(fit$conf, "keyLogFsta"),
lh(fit$conf, "keyLogFpar"),
lh(fit$conf, "keyVarObs"),
lh(fit$conf, "keyQpow"),
lh(fit$conf, "keyVarF")) %>%
dplyr::select(par_conf, dplyr::everything()) %>%
dplyr::bind_rows(tibble::tibble(par_conf = "keyVarLogN",
fleet_nr = -1,
age = fit$conf$minAge:fit$conf$maxAge,
key = fit$conf$keyVarLogN))
}
#' rbx tibbles from object class "sam"
#'
#' @param fit A "sam"-object
#' @param scale A scaler (default 1)
#'
#' @return A list containing tibbles "rby", "rbya", "opr" and "par"
#'
#' @export
#'
sam_rbx <- function(fit, scale = 1) {
list(rby = sam_rby(fit, scale = scale),
rbya = sam_rbya(fit, data = TRUE, scale = scale),
opr = sam_opr(fit),
par = sam_partable(fit))
}
# B: FLRSAM --------------------------------------------------------------------
# B: residuals, work in progress -----------------------------------------------
#' One-observation-ahead residuals
#'
#' @description Note, this is normally just plotted using plot(res)
#'
#' @param res An object returned from {stockassessment::residuals}
#'
#' @return A tibble, the "one-observation-ahead-residuals" is in
#' variable residual.
sam_one_observation_ahead_residuals <- function(res) {
if(class(res) != "samres") stop("Object has to be of class 'samres'")
# must be an easier way for the following:
d <- tibble::tibble(year = res$year,
fleet = res$fleet,
age = res$age,
observation = res$observation,
nll = res$nll,
grad = res$grad,
mean = res$mean,
residual = res$residual)
key <- tibble::tibble(fleet = d$fleet %>% unique(),
what = attributes(res)$fleetNames)
d %>%
dplyr::left_join(key, by = "fleet") %>%
return()
}
#empirobscorrplot(RES)
#setcap("One-observation-ahead residual correlations", "Residual correlations.")
# stockassessment:::empirobscorrplot.samres(RES)
#stockassessment:::empirobscorrplot.samres
#function (res, ...)
#{
# res <- RES
# dat <- data.frame(resid = res$residual, age = res$age, year = res$year,
# fleet = res$fleet)
# fleets <- unique(dat$fleet)
# fn <- attr(res, "fleetNames")
# str(dat)
#
# x <- list()
# for (i in 1:length(fleets)) {
# tmp <- xtabs(resid ~ age + year, data = dat[dat$fleet ==
# fleets[i], ])
# xx <- cor(t(tmp))
# x[[length(x) + 1]] <- xx
# }
#
# corplotcommon(x, fn)
# Process residuals ----------------------------------
# if(exists("RESP")){
# plot(RESP)
# setcap("Process residuals", "Standardized single-joint-sample residuals of process increments")
# stampit(fit)
# par(mfrow=c(1,1))
# }
#' Proces residuals
#'
#' @param resp An object returned from {stockassessment::procres}
#'
#' @return A tibble
sam_process_residuals <- function(resp) {
d <-
tibble::tibble(year = resp$year,
age = resp$age,
fleet = resp$fleet,
residual = resp$residual)
key <-
tibble::tibble(fleet = d$fleet %>% unique(),
fleetn = attributes(resp)$fleetNames)
d %>%
dplyr::left_join(key, by = "fleet") %>%
return()
}
# C: sam various ---------------------------------------------------------------
#' @title sam process error
#'
#' @description XXX
#'
#' @export
#'
#' @param rbya XXX
#' @param plot_it XXX
#' @param plot_catch XXX
#' @param plus_group XXX
#'
sam_process_error <- function(rbya, plot_it=FALSE, plot_catch = FALSE, plus_group=TRUE) {
last.true.age <- max(rbya$age) - plus_group
last.year <- max(rbya$year)
# # align the year-classes
# x <- rbya[,c("year","age","n")]
# x$year <- x$year - 1
# x$age <- x$age - 1
# names(x)[3] <- "n.end"
# d <- plyr::join(rbya[,c("year","age","n","m","f","oC", "cW")],x, by=c("year","age"))
# d <- d[!is.na(d$n.end),]
# # Note: could have a problem for the terminal year
# d$f[is.na(d$f)] <- 0
# # exclude plus-group
# if(plus_group) d <- d[d$age < max(d$age),]
# # process error expressed as mortality
# d$z.n <- log(d$n/d$n.end)
# d$z.f <- d$f + d$m
# d$z.d <- d$z.n - d$z.f
#
# # process error expressed as numbers
# d$n.end2 <- d$n * exp(-(d$f + d$m))
# # Calculate the difference
# d$n.d <- d$n.end - d$n.end2
d <-
rbya %>%
dplyr::mutate(yc = year - age,
f = tidyr::replace_na(f, 0)) %>%
dplyr::arrange(yc, age) %>%
dplyr::group_by(yc) %>%
# process error expressed as mortality
dplyr::mutate(n.end = dplyr::lead(n),
z.d = log(n / dplyr::lead(n)) - (f + m),
# process error expressed as numbers
n.d = dplyr::lead(n) - (n * exp(-(f + m))),
# process errror expressed as biomass
b.d = n.d * cW) %>%
dplyr::ungroup() %>%
dplyr::mutate(z.d = ifelse(age >= last.true.age, NA, z.d),
n.d = ifelse(age >= last.true.age, NA, n.d),
b.d = ifelse(age >= last.true.age, NA, b.d)) %>%
dplyr::arrange(year, age)
#
# process errror expressed as biomass
x <-
d %>%
dplyr::group_by(year) %>%
dplyr::summarise(b = sum(b.d, na.rm = TRUE),
y = sum(oC * cW, na.rm = TRUE))
if(plot_it) {
mort <-
ggplot2::ggplot(d, ggplot2::aes(year, z.d)) +
ggplot2::theme_bw() +
ggplot2::geom_text(ggplot2::aes(label = age)) +
ggplot2::stat_smooth(span = 0.1) +
ggplot2::labs(x = NULL, y = NULL,
title = "Process error expressed as deviations in mortality")
abun <-
ggplot2::ggplot(d, ggplot2::aes(year, n.d)) +
ggplot2::theme_bw() +
ggplot2::geom_col() +
ggplot2::facet_wrap(~ age, scales = "free_y") +
ggplot2::labs(x = NULL, y = NULL,
title = "Process error expressed as deviations in number of fish")
mass <- ggplot2::ggplot(x, ggplot2::aes(year, b)) +
ggplot2::theme_bw() +
ggplot2::geom_bar(stat="identity") +
ggplot2::labs(x = NULL, y = "Mass",
title = "Process error expressed as deviation in mass")
if(plot_catch) {
abun <-
abun +
ggplot2::geom_point(ggplot2::aes(y = oC),
colour = "red",
size = 0.5)
mass <-
mass +
ggplot2::geom_point(ggplot2::aes(y = y),
colour = "red",
size = 0.5)
}
return(list(rbya = d, rby = x, mort = mort, abun = abun, mass = mass))
}
return(list(rbya=d, rby=x))
}
sam_ypr <- function(fit) {
# NOTE: May want to pass arguements to ypr, like number of years, etc.
tmp <- stockassessment::ypr(fit)
res <- tibble::tibble(fbar = tmp$fbar,
yield = tmp$yield,
ssb = tmp$ssb)
ref <- tibble::tibble(ref = c("Fmax", "F01", "F35"),
fbar = c(tmp$fmax, tmp$f01, tmp$f35),
ssb = c(res$ssb[tmp$fmaxIdx],
res$ssb[tmp$f01Idx],
res$ssb[tmp$f35Idx]),
yield = c(res$yield[tmp$fmaxIdx],
res$yield[tmp$f01Idx],
res$yield[tmp$f35Idx]))
return(list(ypr = res, ref = ref))
}
#if(!all(fit$conf$obsCorStruct=="ID")){
# corplot(fit)
# setcap("Estimated correlations", "Estimates correlations between age groups for each fleet")
# stampit(fit)
#}
# CHECK OUT: stockassessment:::obscorrplot.sam
# obscov(fit, TRUE)
# Some experimental stuff
#' Sam catcability
#'
#' @param fit A "sam" object
#'
#' @return A tibble
#'
sam_catchability <- function(fit) {
Q <- fit$pl$logFpar
Qsd <- fit$plsd$logFpar
key <- fit$conf$keyLogFpar
fun <- function(x)if(x<0){NA}else{Q[x+1]}
FF <- Vectorize(fun)
ages <- fit$conf$minAge:fit$conf$maxAge
Qage <- exp(t(matrix(FF(key), nrow = nrow(key))))
d <-
Qage %>%
tidyr::as_tibble()
names(d) <- attr(fit$data, "fleetNames")[1:nrow(key)]
d$age <- ages
d %>%
tidyr::gather(fleet, value, -age) %>%
tidyr::drop_na()
}
# ------------------------------------------------------------------------------
# Get stuff from www.stockassessment.org
#' @title Get a sam directory from stockassessment.org
#'
#' @description The function copies the whole directory of an assessment run from
#' stockassessment.org to a local directory
#'
#' @param assessment Name of the assessment
#' @param user Name of the user (default "user3")
#'
sam_get_directory <- function(assessment, user = "user3") {
path <- paste("https://www.stockassessment.org/datadisk/stockassessment/userdirs",user,assessment,sep="/")
cmd <- paste0("wget --recursive --reject=png,html --level=0 --no-parent ",path,"/")
system(cmd)
# cleanup
Path <- "www.stockassessment.org/datadisk/stockassessment/userdirs"
Path <- paste(Path,user,assessment,sep="/")
cmd <- paste("mv", Path, ".")
system(cmd)
system("rm -r www.stockassessment.org")
}
#' sam_get_data
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return A list
#'
sam_get_data <- function(assessment) {
dat <- NA
load(url(sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run/data.RData")))
return(dat)
}
#' sam_get_residuals
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return XXXX
#'
sam_get_residuals <- function(assessment) {
fil <- sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run/residuals.RData")
if(!RCurl::url.exists(fil)) {
stop(paste0("File: 'residuals.RData for ",
assessment,
" does not exist at: ",
sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run"))
)
}
RES <- NA
load(url(fil))
return(RES)
}
#' sam_get_retro
#'
#' @param assessment Name of run on stockassessment.org
#'
#' @return XXXX
sam_get_retro <- function(assessment) {
RETRO <- NA
load(url(sub("SN", assessment, "https://stockassessment.org/datadisk/stockassessment/userdirs/user3/SN/run/retro.RData")))
return(RETRO)
}
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.