R/additional_functions.R
In MichaelFolkes/VRAP2: A Fork of the NMFS's VRAP Package for RER and populaton viability computations
#' @title Calculate confidence intervals on samples from a binomial distribution
#'
#' @param k A integer of length one reprsenting the 'successes'
#' @param n A integer of length one reprsenting the number samples taken
#' @param CIlevel A numeric of length one reprsenting the confidence level.
#' Default is 0.95
#'
#' @return A data frame.
#' @export
#'
#' @examples
calcBinomialCI <- function(k, n, CIlevel=0.95){
p <- k/n
CIlevel.upper <- CIlevel+(1-CIlevel)/2
CI <- p + c(-qnorm(CIlevel.upper),qnorm(CIlevel.upper))*sqrt((1/n)*p*(1-p))
res <- data.frame(p, CI[1], CI[2])
colnames(res) <- c("proportion", paste0(c("lower", "upper"), CIlevel, "CI"))
res
}#END calcBinomialCI
#' @title Calculate RER values from VRAP output (original method)
#'
#' @param vrap.output A list. Output from the \code{\link{lapply}} application
#' on \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several matrices. The elemented named \code{rer}
#' has the rer estimates.
#'
#' @export
#'
#' @examples
calcRER <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the LEL
#lower esc proportions:
data.quantile.lower <- aggregate(escapement~er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){ sum(x<value)/length(x) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#upper esc proportions:
#geometric mean by rep (run):
data.esc.df$escapement.positive <- data.esc.df$escapement
data.esc.df$escapement.positive[data.esc.df$escapement.positive<1] <- 1
#!! these proportions values are for ABOVE uel
data.quantile.upper <- aggregate(escapement.positive~er.target+rep, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x) exp(mean(log(x)))>data.esc.element$uel)
data.quantile.upper <- aggregate(escapement.positive~er.target, data=data.quantile.upper,FUN = function(x) sum(x)/length(x))
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement.positive'] <- 'proportion.upper'
#this is an alternate, experimental approch to calculating the upper proportion. It matches the method used for the lower proportion in that it pools the values for the final X years and takes the proportion using all runs in one step:
data.quantile.upper2 <- aggregate(escapement.positive~er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){ sum(x>value)/length(x) })
colnames(data.quantile.upper2)[colnames(data.quantile.upper2)=='escapement.positive'] <- 'proportion.upper.alternate'
data.quantile.upper <- merge(data.quantile.upper, data.quantile.upper2, by="er.target")
#data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower, data.quantile.upper, by="er.target")
#estimate the proportion, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$proportion.lower<= probabilityAtzeroER.adj & data.quantile.stats$proportion.upper>=(1-probs.row["prop.upper"]),]
results.rer.alt <- data.quantile.stats[data.quantile.stats$proportion.lower<= probabilityAtzeroER.adj & data.quantile.stats$proportion.upper.alternate>=(1-probs.row["prop.upper"]),]
#trap cases lacking any rer outcome:
if(nrow(results.rer)>0){
results.rer$rer.alternate <- max(results.rer.alt$er.target, na.rm = TRUE)
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
}
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) {
#trap cases lacking any rer outcome:
if(nrow(x)>0){
data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target, na.rm = TRUE), rer.alternate=max(x$rer.alternate, na.rm = TRUE))
}
})#END lapply
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
return(list(rer=rer, rer.detailed=results.rer, data.quantile.stats=data.quantile.stats))
}, probs.df)
return(results.rer.list)
}#END calcRER
#' @title Calculate RER values from VRAP output (modified method)
#'
#' @param vrap.output A list. Output from the \code{lapply} application on
#' \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several data frames. The elemented named \code{rer}
#' has the rer estimates. Other data frames include columns of mean and SD of
#' the proportions, by target ER. These statistics should not be trusted as
#' the frequency distributions of proportions, by targe ER, is never normal.
#' (This was discovered after the fact.) This can be confirmed by creating a
#' historgram, by target ER of either
#' \code{data.quantile.raw$proportion.lower} or
#' \code{data.quantile.raw$proportion.upper}.
#'
#' @export
#'
#' @examples
calcRER2 <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
#test if it has 'rep' variable, if not then create it:
if(!"rep" %in% names(sim.ind$output$rawoutput[[1]])){
for(i in 1:length(sim.ind$output$rawoutput)){
sim.ind$output$rawoutput[[i]]$rep <- i
}
}
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the thresholds
#lower esc probabilities:
data.quantile.lower <- aggregate(escapement~rep+er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#the distribution of probabilities, by ER, may be lognormal
#data.quantile.lower$proportion.lower[data.quantile.lower$proportion.lower==0] <- 1e-10
#data.quantile.lower$proportion.lower <- log(data.quantile.lower$proportion.lower)
data.quantile.lower.stats <- aggregate(proportion.lower~er.target, data=data.quantile.lower, FUN=function(x){c(prop.lower.mean=(mean(x)), prop.lower.sd=(sd(x)))})
data.quantile.lower.stats <- data.frame(er.target=data.quantile.lower.stats$er.target, data.quantile.lower.stats$proportion.lower)
#upper esc probabilities:
data.quantile.upper <- aggregate(escapement~rep+er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement'] <- 'proportion.upper'
#the distribution of probabilities, by ER, may be lognormal
# data.quantile.upper$proportion.upper[data.quantile.upper$proportion.upper==0] <- 1e-6
# data.quantile.upper$proportion.upper <- log(data.quantile.upper$proportion.upper)
data.quantile.upper.stats <- aggregate(proportion.upper~er.target, data=data.quantile.upper, FUN=function(x){c(prop.upper.mean=(mean(x)), prop.upper.sd=(sd(x)))})
data.quantile.upper.stats <- data.frame(er.target=data.quantile.upper.stats$er.target, data.quantile.upper.stats$proportion.upper)
data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower.stats, data.quantile.upper.stats, by="er.target")
#estimate the probability, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
VRAP2:::quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$prop.lower.mean<= probabilityAtzeroER.adj & data.quantile.stats$prop.upper.mean<=probs.row["prop.upper"],]
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target, na.rm = TRUE)))
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
data.quantile.stats$probabilityAtzeroER <- probs.df$probabilityAtzeroER[1]
data.quantile.stats$prop.upper <- probs.df$prop.upper[1]
RER.empirical <- data.quantile.raw[data.quantile.raw$proportion.lower<=rer$prob.lower & data.quantile.raw$proportion.upper<=rer$prob.upper,]
return(list(rer=rer, rer.detailed=results.rer, data.quantile.stats=data.quantile.stats, data.quantile.raw=data.quantile.raw, RER.empirical=RER.empirical))
}, probs.df)
return(results.rer.list)
}#END calcRER2
#' @title Calculate RER values from VRAP output (binomial method)
#'
#' @param x not sure about this.
#'
#' @return A data frame.
#'
#' @description Not sure about the degree of completeness of this function.
#'
#' @export
#'
#' @examples
#' \dontrun{
#' filepaths <- list.files(pattern = "*\\.rds$", full.names = TRUE, ignore.case = TRUE)
#' res.allscenarios.list <- lapply(filepaths, function(x){
#' vrap.output <- readRDS(x)
#' scenario <- names(vrap.output)
#' res.list <- extractQualifiedRuns(vrap.output)
#' return(list(qualified.reps=res.list[[1]], scenario=scenario, nruns=vrap.output[[1]]$inputs$NRuns))
#' })
#'
#' res <- calcRERbinomial(res.allscenarios.list)
#' }
calcRERbinomial <- function(x){
qualified.reps.list <- lapply(x, function(singlescenario) singlescenario$qualified.reps)
qualified.reps.df <- do.call("rbind", qualified.reps.list)
proportions.success <- aggregate(rep~er.target, data = qualified.reps.df[,c("rep","er.target")], length)
colnames(proportions.success)[colnames(proportions.success)=="rep"] <- "success.n"
nruns.total <- sum(sapply(x, function(singlescenario) singlescenario$nruns))
proportions.success$nruns.total <- nruns.total
rer.binomial <- apply(proportions.success, 1, function(x){
calcBinomialCI(k = x["success.n"], n = x["nruns.total"])
})
rer.binomial <- do.call("rbind", rer.binomial)
rer.binomial <- cbind(proportions.success, rer.binomial)
row.names(rer.binomial) <- NULL
return(rer.binomial)
}#END calcRERbinomial
#' @title Calculate RER values from VRAP output (binomial method)
#'
#' @param vrap.output A list. Output from the \code{lapply} application on
#' \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several data frames. The elemented named \code{rer}
#' has the rer estimates. Other data frames include columns of mean and SD of
#' the proportions, by target ER. These statistics should not be trusted as
#' the frequency distributions of proportions, by targe ER, is never normal.
#' (This was discovered after the fact.) This can be confirmed by creating a
#' historgram, by target ER of either
#' \code{data.quantile.raw$proportion.lower} or
#' \code{data.quantile.raw$proportion.upper}.
#'
#'
#' @examples
calcRERbinomial.old <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
#I think this was the first version and is a combo of the functions:
#extractQaulifiedRuns and calcRERbinomial
# updates to extractQaulifiedRuns have happended such that this is now behind.
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
#test if it has 'rep' variable, if not then create it:
if(!"rep" %in% names(sim.ind$output$rawoutput[[1]])){
for(i in 1:length(sim.ind$output$rawoutput)){
sim.ind$output$rawoutput[[i]]$rep <- i
}
}
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the thresholds
#lower esc probabilities:
data.quantile.lower <- aggregate(escapement~rep+er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#the distribution of probabilities, by ER, may be lognormal
#data.quantile.lower$proportion.lower[data.quantile.lower$proportion.lower==0] <- 1e-10
#data.quantile.lower$proportion.lower <- log(data.quantile.lower$proportion.lower)
data.quantile.lower.stats <- aggregate(proportion.lower~er.target, data=data.quantile.lower, FUN=function(x){c(prop.lower.mean=(mean(x)), prop.lower.sd=(sd(x)))})
data.quantile.lower.stats <- data.frame(er.target=data.quantile.lower.stats$er.target, data.quantile.lower.stats$proportion.lower)
#upper esc probabilities:
data.quantile.upper <- aggregate(escapement~rep+er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement'] <- 'proportion.upper'
#the distribution of probabilities, by ER, may be lognormal
# data.quantile.upper$proportion.upper[data.quantile.upper$proportion.upper==0] <- 1e-6
# data.quantile.upper$proportion.upper <- log(data.quantile.upper$proportion.upper)
data.quantile.upper.stats <- aggregate(proportion.upper~er.target, data=data.quantile.upper, FUN=function(x){c(prop.upper.mean=(mean(x)), prop.upper.sd=(sd(x)))})
data.quantile.upper.stats <- data.frame(er.target=data.quantile.upper.stats$er.target, data.quantile.upper.stats$proportion.upper)
data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower.stats, data.quantile.upper.stats, by="er.target")
#estimate the probability, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
VRAP2:::quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
# continue from here.
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$prop.lower.mean<= probabilityAtzeroER.adj & data.quantile.stats$prop.upper.mean<=probs.row["prop.upper"],]
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target)))
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
data.quantile.stats$probabilityAtzeroER <- probs.df$probabilityAtzeroER[1]
data.quantile.stats$prop.upper <- probs.df$prop.upper[1]
RER.empirical <- data.quantile.raw[data.quantile.raw$proportion.lower<=rer$prob.lower & data.quantile.raw$proportion.upper<=rer$prob.upper,]
proportions.success <- aggregate(rep~er.target, data = RER.empirical[,c("rep","er.target")], length)
colnames(proportions.success)[colnames(proportions.success)=="rep"] <- "success.n"
proportions.success$nruns <- vrap.output[[1]]$inputs$NRuns
#proportions.success$success.prop <- proportions.success$success.n/proportions.success$nruns
rer.binomial <- apply(proportions.success, 1, function(x){
calcBinomialCI(k = x["success.n"], n = x["nruns"])
})
rer.binomial <- do.call("rbind", rer.binomial)
rer.binomial <- cbind(proportions.success, rer.binomial)
return(list(rer.binomial=rer.binomial, rer.detailed=results.rer, data.quantile.stats=data.quantile.stats, data.quantile.raw=data.quantile.raw))
}, probs.df)
return(results.rer.list)
}#END calcRERbinomial2
#' @title Calculate RER values from summary table
#'
#' @param vrap.output
#' @param prop.lower.tolerance
#' @param prop.upper
#'
#' @return
#' @export
#'
#' @examples
calcRERfromsummary <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance=100*prop.lower.tolerance, prop.upper=100*prop.upper)
results <- lapply(vrap.output, FUN = function(x){
output.tabulated <- x$data.summary$output.tabulated
EndAv <- x$inputs$EndAv
NYears <- x$inputs$NYears
NRuns <- x$inputs$NRuns
metadata <- data.frame(NRuns=NRuns, NYears=NYears, EndAv=EndAv)
#this apply will return a list with length=#probability levels assessed. (usually=1)
apply.res <- apply(probs.df,1, FUN=function(probs.single, output.tabulated, metadata){
prop.lower <- probs.single["prop.lower.tolerance"]+output.tabulated$percent.years.belowLEL[output.tabulated$TgtER==0]
qualifyingER <- output.tabulated[output.tabulated$percent.years.belowLEL< prop.lower & output.tabulated$percent.years.aboveUEL > 100- probs.single["prop.upper"],]
if(nrow(qualifyingER)==0){
#no ER levels qualify
rer <- data.frame(metadata, lower.percent=prop.lower, upper.percent=100- probs.single["prop.upper"], rer=NA)
}else{
rer <- data.frame(metadata, lower.percent=prop.lower, upper.percent=100- probs.single["prop.upper"], rer=max(qualifyingER$TgtER, na.rm = TRUE))
}
row.names(rer) <- NULL
return(list(rer=rer, qualifyingER=qualifyingER))
}, output.tabulated, metadata)
})
#results <- results[[1]]
names(results) <- names(vrap.output)
return(results)
}#END calcRERfromsummary
#' @title Extract qualified runs from VRAP output (for binomial method)
#'
#' @param vrap.output A list. Output from the \code{lapply} application on
#' \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several data frames. The elemented named \code{rer}
#' has the rer estimates. Other data frames include columns of mean and SD of
#' the proportions, by target ER. These statistics should not be trusted as
#' the frequency distributions of proportions, by targe ER, is never normal.
#' (This was discovered after the fact.) This can be confirmed by creating a
#' historgram, by target ER of either
#' \code{data.quantile.raw$proportion.lower} or
#' \code{data.quantile.raw$proportion.upper}.
#'
#' @export
#'
#' @examples
extractQualifiedRuns <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
#test if it has 'rep' variable, if not then create it:
if(!"rep" %in% names(sim.ind$output$rawoutput[[1]])){
for(i in 1:length(sim.ind$output$rawoutput)){
sim.ind$output$rawoutput[[i]]$rep <- i
}
}
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(scenario=sim.ind$inputs$filename, lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the thresholds
#lower esc probabilities:
data.quantile.lower <- aggregate(escapement~rep+er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#the distribution of probabilities, by ER, may be lognormal
#data.quantile.lower$proportion.lower[data.quantile.lower$proportion.lower==0] <- 1e-10
#data.quantile.lower$proportion.lower <- log(data.quantile.lower$proportion.lower)
data.quantile.lower.stats <- aggregate(proportion.lower~er.target, data=data.quantile.lower, FUN=function(x){c(prop.lower.mean=(mean(x)), prop.lower.sd=(sd(x)))})
data.quantile.lower.stats <- data.frame(er.target=data.quantile.lower.stats$er.target, data.quantile.lower.stats$proportion.lower)
#upper esc probabilities:
data.quantile.upper <- aggregate(escapement~rep+er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement'] <- 'proportion.upper'
#the distribution of probabilities, by ER, may be lognormal
# data.quantile.upper$proportion.upper[data.quantile.upper$proportion.upper==0] <- 1e-6
# data.quantile.upper$proportion.upper <- log(data.quantile.upper$proportion.upper)
data.quantile.upper.stats <- aggregate(proportion.upper~er.target, data=data.quantile.upper, FUN=function(x){c(prop.upper.mean=(mean(x)), prop.upper.sd=(sd(x)))})
data.quantile.upper.stats <- data.frame(er.target=data.quantile.upper.stats$er.target, data.quantile.upper.stats$proportion.upper)
data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower.stats, data.quantile.upper.stats, by="er.target")
#estimate the probability, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
VRAP2:::quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$prop.lower.mean<= probabilityAtzeroER.adj & data.quantile.stats$prop.upper.mean<=probs.row["prop.upper"],]
#only do this if there is an rer for this simulation:
if(nrow(results.rer)>0){
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
} else {
results.rer$prob.lower <- numeric(0)
results.rer$prob.upper <- numeric(0)
}
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) {
if(nrow(x)>0) {return(data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target, na.rm = TRUE)))}
})
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
data.quantile.stats$probabilityAtzeroER <- probs.df$probabilityAtzeroER[1]
data.quantile.stats$prop.upper <- probs.df$prop.upper[1]
RER.empirical <- data.quantile.raw[data.quantile.raw$proportion.lower<=rer$prob.lower & data.quantile.raw$proportion.upper<=rer$prob.upper,]
#this traps cases when no ER achieves criteria and nrow=0
if(nrow(RER.empirical)>0){
RER.empirical.sums <- aggregate(rep~er.target, data=RER.empirical, length)
} else {
RER.empirical.sums <- data.frame(rep=integer(), er.target=numeric())
}
colnames(RER.empirical.sums)[colnames(RER.empirical.sums)=="rep"] <- "success.n"
dat.samplesize <- aggregate(rep~er.target,data= data.quantile.raw, length)
RER.empirical.sums <- merge(dat.samplesize, RER.empirical.sums, all.x = TRUE)
RER.empirical.sums$success.n[is.na(RER.empirical.sums$success.n)] <- 0
colnames(RER.empirical.sums)[colnames(RER.empirical.sums)=="rep"] <- "reps.n"
return(list(reps.qualified =RER.empirical, reps.summary=RER.empirical.sums, scenario=data.esc.element$scenario))
}, probs.df)
return(results.rer.list)
}#END extractQualifiedRuns
#' @title Read in a single .rav2 file
#'
#' @param filename A character vector of length one. The name of the .rav2 file.
#'
#' @return A list of two elements named \code{inputs.vars} and \code{inputs}.
#' The element \code{inputs.vars} is a list of just the variables and their
#' values. The element named \code{inputs} is a list with the variable names,
#' values, and description. The length of each list (\code{inputs.vars} and
#' \code{inputs}) is equal to the number of rows in the .rav2 file.
#' @export
#'
#' @details The .rav2 file is a csv file. The first column is the variable name,
#' the second column is the value for that variable, the third column is a
#' description of the variable. If commas are included in either the value or
#' descripton then the full column should be inside double quotes. If the
#' value is a numeric vector of several values, the values should be separated
#' by a space (not a comma) and not wrapped in quotes. Additional columns may
#' be added after the description column. They could include additional
#' comments and information. They will be represented, by variable, in the
#' list named \code{inputs}.
#'
#' @examples
read_rav2 <- function(filename){
x <- scan(filename, what="",quote = "", sep="\n")
x <- x[!grepl(pattern = "^#", x)]
# Separate elements by comma
#y <- strsplit(x, ",")
#this allows for input of strings that have commas
#if using commas, the full string must be in double quotes:
y <- strsplit(x, '\\"[^"]+,(*SKIP)(*FAIL)|,\\s*', perl=TRUE)
# Extract the first vector element and set it as the list element name
names(y) <- trimws(sapply(y, `[[`, 1))
# Remove the first vector element from each list element
y <- lapply(y, `[`, -1)
inputs <- lapply(y, FUN = function(x){
value <- trimws(x[1])
#split on one or more spaces:
value <- strsplit(x = value, "\\s+")[[1]]
value <- type.convert(value, as.is = TRUE)
#if the split data are characters then use the original form (unsplit):
if(any(class(value)=="character")) {
value <- trimws(x[1])
value <- gsub(pattern = "\\\"", replacement = "", x = value)
}
res <- list(value=value)
#if there's a third column in the csv file, then it's the description
#(index=2 as first column removed):
x.count <- length(x)
if(x.count>1) res$description=trimws(gsub(pattern = "\\\"", replacement = "", x = x[2]))
if(x.count>2){
#additional user-specific columns
for(i in 3:x.count){
element.name <- paste0("additional",i-2)
strsplit(x, '\\"[^"]+,(*SKIP)(*FAIL)|,\\s*', perl=TRUE)
res[element.name]=trimws(gsub(pattern = "\\\"", replacement = "", x = x[i]))
}
}
return(res)
})
#integer of the number of steps of ER or Pop Capacity to do # is 1:BufMax
inputs$BufMax <- list(value = round((inputs$BufferEnd$value - inputs$BufferStart$value) / inputs$BufferStep$value + 1), description= "integer of the number of steps of ER or Pop Capacity to do" )
inputs.vars <- lapply(inputs, "[[","value")
names(inputs.vars) <- names(inputs)
return(list(inputs.vars=inputs.vars, inputs=inputs))
}#END read_rav2
printInput <- function(filename){
input <- GetInput(filename)
}#END printInput
quantInv <- function(distr, value) {ecdf(distr)(value)}
#' @title Build multiple copies of same .rav file for sensitivity tests
#'
#' @param rav.filenames A charactdata.quantile.upper vector. The input .rav files that will be
#' copied.
#' @param sensitivity.runs A integer of length one. The number of copies to make
#' of each .rav file.
#'
#' @return Each input .rav file is copied \code{sensitivity.runs} number of
#' times. The filename of each of file is appended with a numerical sequence
#' from 1 to \code{sensitivity.runs}. The original files are deleted as all
#' copies are exactly same as the original.
#' @export
#'
#' @examples
sensitivityBuild <- function(rav2.filenames, sensitivity.runs){
invisible(
lapply(rav2.filenames, FUN = function(filename, sensitivity.runs){
run.vec <- (1:sensitivity.runs) / 10^(1+floor(log10(sensitivity.runs)))
run.vec <- formatC(run.vec,format='f', digits = 1+floor(log10(sensitivity.runs)))
run.vec <- substr(run.vec, 3, nchar(run.vec))
lapply(run.vec, FUN = function(run.single, filename){
file.copy(filename, paste0(tools::file_path_sans_ext(filename), "sensitivity.run_", run.single, ".", tools::file_ext(filename)))
}, filename)
}, sensitivity.runs)
)
#remove originals:
invisible(lapply(rav2.filenames, FUN = function(x) file.remove(x)))
}#END sensitivityBuild
tabulateSummary <- function(vrap.output){
vrap.output <- lapply(vrap.output, FUN=function(x){
# browser()
x$data.summary[(1+grep(pattern = "__Escapement__", x = x$data.summary)):length(x$data.summary)]
})
}#END tabulateSummary
#' @title Write one or many .rav2 files, which hold VRAP input variables.
#'
#' @param indexfilename A Boolean, default is FALSE. The rav2 filename can
#' sometimes get too long for windows. If so, use indexfilename=TRUE, which
#' forces a date-time stamp (to 1/10,000 of a second) to be used for a
#' filename. Additionally a translation table is written that will show both the time
#' stamp filename and the longer filename, which includes arguments.
#' @param filename The name of the .rav2 file to be written. If NA (default)
#' then the function will create a file name based on the argument values
#' choosen.
#' @param defaultargs.filename The name of a .rav2 file that includes default
#' values. If this is chosen then the default values in this function are
#' ignored. But values set in the function arguments are not ignored.
#' @param Title
#' @param RanSeed
#' @param NRuns
#' @param NYears
#' @param MinAge
#' @param MaxAge
#' @param ConvergeCrit
#' @param Debugg
#' @param SRType
#' @param BSRa
#' @param BSRb
#' @param BSRc
#' @param BSRd
#' @param AveEnv
#' @param depen
#' @param DL1
#' @param DL2
#' @param DR
#' @param EscChoice
#' @param SurvScale
#' @param SRError.FUN
#' @param SRErrorA
#' @param SRErrorB
#' @param ResCorParameter
#' @param MarSurv
#' @param BetaMarA
#' @param BetaMarB
#' @param CorrMS
#' @param FWError.FUN
#' @param GammaFlowA
#' @param GammaFlowB
#' @param MarError.FUN
#' @param GammaMarA
#' @param GammaMarB
#' @param NumBreakPoints
#' @param BaseRegime
#' @param EscpmntBreakPoint
#' @param TargetU
#' @param MgmtError
#' @param MgmtError.FUN
#' @param GammaMgmtA
#' @param GammaMgmtB
#' @param ECrit
#' @param ERecovery
#' @param EndAv
#' @param StepFunc
#' @param BufferStep
#' @param BufferStart
#' @param BufferEnd
#' @param CohortStart
#' @param NatMort
#' @param MatRate
#' @param PTU
#' @param MatU
#' @param output.raw A character vector "AEQMort Escpmnt TotAEQMort
#' TotAdultEscpmnt TotEscpmnt RanFlow RanMarine Recruits TargetUError
#' ptu.scaled MatU.scaled Cohort tmpdat TempCohort CalendarHR TargetU" These
#' are the variables that will be included in the output.
#' @param ...
#'
#' @return Writes one or more .rav2 (text/csv) files depending on number of
#' unique combinations created from the input arguments.
#' @export
#'
#' @examples
writeRav2file <- function(
indexfilename=FALSE,
filename=NA,
defaultargs.filename=NA,
Title="ER example of Beverton-Holt VRAP with no marine survival or stream flow variability",
RanSeed=0,
NRuns=1000,
NYears=25,
MinAge=2,
MaxAge=5,
ConvergeCrit=0.001,
Debugg="NO",
SRType="BEV2",
BSRa=9.27599153996063,
BSRb=3702.14272977426,
BSRc=NA,
BSRd=NA,
AveEnv=NA,
depen="YES",
DL1=400,
DL2=63,
DR=0,
EscChoice="YES",
SurvScale="YES",
SRError.FUN="rgamma(1, shape=inputs$SRErrorA, scale=inputs$SRErrorB)",
SRErrorA=0.967637142385205,
SRErrorB=1.09375572895311,
ResCorParameter=NA,
MarSurv="NO",
BetaMarA=NA,
BetaMarB=NA,
CorrMS=NA,
FWError.FUN="rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)",
GammaFlowA=0,
GammaFlowB=0,
MarError.FUN="rgamma(1, inputs$GammaMarA, scale=inputs$GammaMarB)",
GammaMarA=0,
GammaMarB=0,
NumBreakPoints=0,
BaseRegime=1,
EscpmntBreakPoint=NA,
TargetU=0.37,
MgmtError="YES",
MgmtError.FUN="rgamma(1, shape=inputs$GammaMgmtA, scale=inputs$GammaMgmtB)",
GammaMgmtA=65.3946,
GammaMgmtB=0.0158,
ECrit=416,
ERecovery=1050,
EndAv=5,
StepFunc="ER",
BufferStep=6.75675675675676E-02,
BufferStart=0,
BufferEnd=1.89189189189189,
CohortStart="5108 2554 1478 667 112",
NatMort="0.5 0.4 0.3 0.2 0.1",
MatRate="0 1.43432324005842E-02 0.226067824725145 0.685044912939359 1",
PTU="0 7.97982424808335E-02 0.234143029720483 0.352056542423367 0.33478625941973",
MatU="0 0.315672687185273 0.331197585913413 0.339907561109324 0.32870153694989",
output.raw= "AEQMort Escpmnt TotAEQMort TotAdultEscpmnt TotEscpmnt RanFlow RanMarine Recruits TargetUError ptu.scaled MatU.scaled Cohort tmpdat TempCohort CalendarHR TargetU",
...)
{
rav.df <- as.data.frame(matrix(ncol=3, byrow=TRUE,
data=c(
"Title", "argList$Title", "Title",
"RanSeed", "argList$RanSeed", "Random seed; 0 gives random seed; numbers give fixed seed",
"NRuns", "argList$NRuns", "Number of runs",
"NYears", "argList$NYears", "Number of years",
"MinAge", "argList$MinAge", "Minimum and maximum age",
"MaxAge", "argList$MaxAge", "Minimum and maximum age (for now this is fixed; do not change)",
"ConvergeCrit", "argList$ConvergeCrit", "Convergence criterion (proportional error) for target ER. converged when: abs((estimated ER- target ER)/target ER) < ConvergeCrit",
"Debugg", "argList$Debugg", "Debug file flag",
"SRType", "argList$SRType", "Spawner Recruit function (RIC2;RIC3;RIC4; BEV2;BEV3;BEV4; HOC2;HOC3;HOC4)",
"BSRa", "argList$BSRa", "S/R a parameter",
"BSRb", "argList$BSRb", "S/R b parameter",
"BSRc", "argList$BSRc", "marine survival parameter",
"BSRd", "argList$BSRd", "freshwater survival parameter",
"AveEnv", "argList$AveEnv",NA,
"#Mean and CV for M marine survival index (M^c)",NA,NA,
"#Trend; Cycle; or Autoc(orrelation) ?",NA,NA,
"#Trend/Cycle parameters: rate and 0 or increment and 1 for trend- amplitude- period & starting pt & scalar for cycle; correl and starting point (if zero is random) for autocorrelation",NA,NA,
"#Mean and CV for F flow (or other fw) index (exp(dF))",NA,NA,
"#Trend; Cycle; or Autoc(orrelation) ?",NA,NA,
"#Trend/Cycle parameters: rate and 0 or increment and 1 for trend- amplitude- period & starting pt & scalar for cycle; correl and starting point (if zero is random) for autocorrelation",NA,NA,
"FWError.FUN", "argList$FWError.FUN", "The R probability density function with arguments",
"GammaFlowA", "argList$GammaFlowA","FW (flow) error pdf parameter",
"GammaFlowB", "argList$GammaFlowB","FW (flow) error pdf parameter",
"MarError.FUN", "argList$MarError.FUN", "The R probability density function with arguments",
"GammaMarA", "argList$GammaMarA","Marine error pdf parameter",
"GammaMarB", "argList$GammaMarB","Marine error pdf parameter",
"depen", "argList$depen", "Depensation? (\"YES\" or \"NO\")",
"DL1", "argList$DL1", "1) Esc. level for depensation to start",
"DL2", "argList$DL2", "2) QET",
"DR", "argList$DR", "3)% predicted return at QET (or for r/s=1 this parameter = 1)",
"EscChoice", "argList$EscChoice", "Determine recruits from adult spawners (not total)? (\"YES\" or \"NO\")",
"SurvScale", "argList$SurvScale", "Stock-recruit variation (\"YES\" or \"NO\")",
"SRError.FUN", "argList$SRError.FUN", "The R probability density function with arguments",
"SRErrorA", "argList$SRErrorA", "error of S/R a parameter",
"SRErrorB", "argList$SRErrorB", "error of S/R b parameter",
"ResCorParameter", "argList$ResCorParameter", "S/R parameter error autocorrelation",
"MarSurv", "argList$MarSurv", "Smolt to adult survival w/variation (\"YES\" or \"NO\")",
"BetaMarA", "argList$BetaMarA", "Beta distribution a parameter",
"BetaMarB", "argList$BetaMarB", " Beta distribution b parameter",
"CorrMS", "argList$CorrMS", "Beta distribution parameter error autocorrelation",
"NumBreakPoints", "argList$NumBreakPoints", "Number of breakpoints",
"BaseRegime", "argList$BaseRegime", "Level to use as base regime",
"TargetU", "argList$TargetU", "base exploitation rate",
"MgmtError", "argList$MgmtError", "Include error (\"YES\" or \"NO\") in ER management",
"MgmtError.FUN", "argList$MgmtError.FUN", "The R probability density function with arguments",
"GammaMgmtA", "argList$GammaMgmtA", "Management error gamma pdf shape parameter",
"GammaMgmtB", "argList$GammaMgmtB", "Management error gamma pdf scale parameter",
"ECrit", "argList$ECrit", "Lower escapement threshold",
"ERecovery", "argList$ERecovery", "Upper escapement threshold (MSY)",
"EndAv", "argList$EndAv", "# yrs to average for UEL",
"StepFunc", "argList$StepFunc", "Step ER (ER) or Pop Capacity (Pop)?",
"BufferStep", "argList$BufferStep", "Buffer step size as percent of base ER or Pop capacity",
"BufferStart", "argList$BufferStart", "ER Min as proportion of base ER",
"BufferEnd", "argList$BufferEnd", "ER max as proportion of base ER",
"CohortStart", "argList$CohortStart", "Initial population size, by Age, of start year (prior to nat mort)",
"NatMort", "argList$NatMort", "Age specific natural mortality",
"MatRate", "argList$MatRate", "Age specific average maturation rate",
"PTU", "argList$PTU", "Age specific average mixed-maturity fishery fishing rates",
"MatU", "argList$MatU", "Age specific average mature fishery fishing rates")
), stringsAsFactors = FALSE)
argList <- as.list(sys.call())[-1]
if(!is.na(defaultargs.filename)){
#import the default args:
argList <- argList[names(argList) !="defaultargs.filename"]
argList.default <- read_rav2(defaultargs.filename)$inputs.vars
} else {
argList.default <- formals(writeRav2file)
argList.default <- argList.default[names(argList.default) !="..."]
}
#The argList needs to have something in it for following code to work. Otherwise no rav files are created:
if(length(argList)==0) argList$NRuns <- 1000
argLength <- lapply(argList, length)
testing.var.names <- names(which(unlist(argLength)>1))
if(length(testing.var.names)>1) cat("\n!!more than one variable being tested!!\n")
sensitivity.args <- do.call("expand.grid", argList)
#sensitivity.args is a data frame representing all combinations of the args used in the function call
apply(sensitivity.args, 1, FUN=function(x){
if(nrow(sensitivity.args)>1 | is.na(filename)){
#can't use default filename as creating more than one file,
#or no filename supplied
filename <- paste(names(x),x, sep = "_", collapse = "__")
filename <- paste0(filename, ".rav2")
}
#the filename can sometimes get too long for windows. if so, use
#indexfilename=TRUE, then a date-time stamp is used for a filename, and a
#translation table is written that will show both the time stamp filename
#and the longer filename, which includes arguments.
#must set filename here so the correct filename is included in rav2 contents
if(indexfilename){
filename.long <- filename
options(digits.secs=4)
filename <- paste0(format(Sys.time(), "%Y-%m-%d-%H-%M-%OS"), ".rav2")
translation.tbl.filename <- "rav2_filename_translation.csv"
if(file.exists(translation.tbl.filename)){
write.table(x = cbind(filename, filename.long), file = translation.tbl.filename, quote = FALSE, append = TRUE, sep=", ", row.names = FALSE, col.names = FALSE)
}else{
write.table(x = cbind(filename, filename.long) , col.names = c("index.filename", "long.filename"), file = translation.tbl.filename, quote = FALSE, sep=", ", row.names = FALSE)
}
}#END indexfilename
argList <- c(list(filename=filename), x, argList.default[!names(argList.default) %in% names(x)])
rav.str <- rep(NA, length(argList))
for(i in 1:length(argList)){
variable <- names(argList[i])
value <- paste(argList[[i]], collapse = " ")
if(grepl(pattern = ",", x = value)) value <- paste0("\"",value,"\"")
description <- rav.df[rav.df[,1]==names(argList[i]),3]
if(isTRUE(grepl(pattern = ",", x = description))) description <- paste0("\"",description,"\"")
rav.str[i] <- paste(variable, value, description, sep=", ")
}
write(rav.str, file=filename)
} )
cat(nrow(sensitivity.args), "rav2 files written")
}#END writeRav2file
#' @title Build, save, and open an R script to help run VRAP.
#'
#' @description This creates and opens a script named "vrap_script.R". This is a
#' template for the user to work with when running vrap without shiny. This is
#' intended to help the user understand the work flow, but due to file path
#' differences it may not work as is. Some object values will need updating
#' (for example the path).
#'
#' @return Opens an R script that includes the template of functions.
#' @export
#'
#' @examples
#' writeScriptVRAP()
writeScriptVRAP <- function(){
script.str <- c("
#delete old rav files:
rav2.filenames <- list.files(pattern = \"*.rav2$\")
file.remove(rav2.filenames)
writeRav2file(NYears = c(10,25) )
writeRav2file(NRuns = c(3000))
#this makes copies of .rav2 files to allow for testing on what is an appropriate value to get stable estimates of rer. sensitivity.runs is the number of copies to make of each rav file:
rav2.filenames <- list.files(pattern = \"*.rav2$\")
sensitivityBuild(rav2.filenames = rav2.filenames, sensitivity.runs = 10)
#list rav files available:
rav2.filenames <- list.files(pattern = \"*.rav2$\")
#run vrap
#a .rds file will be saved for each .rav2 file tested
#the .rds file prefix will match that of the .rav2 file.
lapply(rav2.filenames, FUN=function(x) VRAP2:::Main(x))
###run vrap in parallel:
require(parallel)
no_cores <- detectCores() - 1
# Initiate cluster
cl <- makeCluster(no_cores)
obj.names <- ls()
clusterExport(cl, obj.names)
clusterEvalQ(cl, library(VRAP2))
rav2.filenames <- list.files(pattern = \"*.rav2$\")
parLapply(cl, rav2.filenames, fun=function(x) VRAP2:::Main(x))
stopCluster(cl)
###End of parallel sims
vrap.output1 <- readRDS(\"NRuns_3000sensitivity.run_1.rds\")
vrap.output2 <- readRDS(\"NRuns_3000sensitivity.run_2.rds\")
vrap.output3 <- readRDS(\"NRuns_3000sensitivity.run_3.rds\")
vrap.output <- c(vrap.output1, vrap.output2, vrap.output3)
names(vrap.output)
rer.orig <- calcRERfromsummary(vrap.output = vrap.output, prop.lower.tolerance = 0.05)
cat(vrap.output[[1]]$data.summary$output)
rer.orig[[1]]$rer
rer.modified <- calcRER2(vrap.output = vrap.output, percent.lower = c(.05, 0.1))
rer.modified[[1]]$rer
hist(rer.modified[[1]]$RER.empirical$er.target)
#obtain the rer's for all simulations:
lapply(vrap.output, \"[[\", \"rer\")
")
write(script.str, file="vrap_script.R")
file.edit("vrap_script.R" )
}#writeScriptVRAP
MichaelFolkes/VRAP2 documentation built on June 9, 2019, 6:27 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.
#' @title Calculate confidence intervals on samples from a binomial distribution
#'
#' @param k A integer of length one reprsenting the 'successes'
#' @param n A integer of length one reprsenting the number samples taken
#' @param CIlevel A numeric of length one reprsenting the confidence level.
#' Default is 0.95
#'
#' @return A data frame.
#' @export
#'
#' @examples
calcBinomialCI <- function(k, n, CIlevel=0.95){
p <- k/n
CIlevel.upper <- CIlevel+(1-CIlevel)/2
CI <- p + c(-qnorm(CIlevel.upper),qnorm(CIlevel.upper))*sqrt((1/n)*p*(1-p))
res <- data.frame(p, CI[1], CI[2])
colnames(res) <- c("proportion", paste0(c("lower", "upper"), CIlevel, "CI"))
res
}#END calcBinomialCI
#' @title Calculate RER values from VRAP output (original method)
#'
#' @param vrap.output A list. Output from the \code{\link{lapply}} application
#' on \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several matrices. The elemented named \code{rer}
#' has the rer estimates.
#'
#' @export
#'
#' @examples
calcRER <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the LEL
#lower esc proportions:
data.quantile.lower <- aggregate(escapement~er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){ sum(x<value)/length(x) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#upper esc proportions:
#geometric mean by rep (run):
data.esc.df$escapement.positive <- data.esc.df$escapement
data.esc.df$escapement.positive[data.esc.df$escapement.positive<1] <- 1
#!! these proportions values are for ABOVE uel
data.quantile.upper <- aggregate(escapement.positive~er.target+rep, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x) exp(mean(log(x)))>data.esc.element$uel)
data.quantile.upper <- aggregate(escapement.positive~er.target, data=data.quantile.upper,FUN = function(x) sum(x)/length(x))
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement.positive'] <- 'proportion.upper'
#this is an alternate, experimental approch to calculating the upper proportion. It matches the method used for the lower proportion in that it pools the values for the final X years and takes the proportion using all runs in one step:
data.quantile.upper2 <- aggregate(escapement.positive~er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){ sum(x>value)/length(x) })
colnames(data.quantile.upper2)[colnames(data.quantile.upper2)=='escapement.positive'] <- 'proportion.upper.alternate'
data.quantile.upper <- merge(data.quantile.upper, data.quantile.upper2, by="er.target")
#data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower, data.quantile.upper, by="er.target")
#estimate the proportion, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$proportion.lower<= probabilityAtzeroER.adj & data.quantile.stats$proportion.upper>=(1-probs.row["prop.upper"]),]
results.rer.alt <- data.quantile.stats[data.quantile.stats$proportion.lower<= probabilityAtzeroER.adj & data.quantile.stats$proportion.upper.alternate>=(1-probs.row["prop.upper"]),]
#trap cases lacking any rer outcome:
if(nrow(results.rer)>0){
results.rer$rer.alternate <- max(results.rer.alt$er.target, na.rm = TRUE)
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
}
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) {
#trap cases lacking any rer outcome:
if(nrow(x)>0){
data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target, na.rm = TRUE), rer.alternate=max(x$rer.alternate, na.rm = TRUE))
}
})#END lapply
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
return(list(rer=rer, rer.detailed=results.rer, data.quantile.stats=data.quantile.stats))
}, probs.df)
return(results.rer.list)
}#END calcRER
#' @title Calculate RER values from VRAP output (modified method)
#'
#' @param vrap.output A list. Output from the \code{lapply} application on
#' \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several data frames. The elemented named \code{rer}
#' has the rer estimates. Other data frames include columns of mean and SD of
#' the proportions, by target ER. These statistics should not be trusted as
#' the frequency distributions of proportions, by targe ER, is never normal.
#' (This was discovered after the fact.) This can be confirmed by creating a
#' historgram, by target ER of either
#' \code{data.quantile.raw$proportion.lower} or
#' \code{data.quantile.raw$proportion.upper}.
#'
#' @export
#'
#' @examples
calcRER2 <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
#test if it has 'rep' variable, if not then create it:
if(!"rep" %in% names(sim.ind$output$rawoutput[[1]])){
for(i in 1:length(sim.ind$output$rawoutput)){
sim.ind$output$rawoutput[[i]]$rep <- i
}
}
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the thresholds
#lower esc probabilities:
data.quantile.lower <- aggregate(escapement~rep+er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#the distribution of probabilities, by ER, may be lognormal
#data.quantile.lower$proportion.lower[data.quantile.lower$proportion.lower==0] <- 1e-10
#data.quantile.lower$proportion.lower <- log(data.quantile.lower$proportion.lower)
data.quantile.lower.stats <- aggregate(proportion.lower~er.target, data=data.quantile.lower, FUN=function(x){c(prop.lower.mean=(mean(x)), prop.lower.sd=(sd(x)))})
data.quantile.lower.stats <- data.frame(er.target=data.quantile.lower.stats$er.target, data.quantile.lower.stats$proportion.lower)
#upper esc probabilities:
data.quantile.upper <- aggregate(escapement~rep+er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement'] <- 'proportion.upper'
#the distribution of probabilities, by ER, may be lognormal
# data.quantile.upper$proportion.upper[data.quantile.upper$proportion.upper==0] <- 1e-6
# data.quantile.upper$proportion.upper <- log(data.quantile.upper$proportion.upper)
data.quantile.upper.stats <- aggregate(proportion.upper~er.target, data=data.quantile.upper, FUN=function(x){c(prop.upper.mean=(mean(x)), prop.upper.sd=(sd(x)))})
data.quantile.upper.stats <- data.frame(er.target=data.quantile.upper.stats$er.target, data.quantile.upper.stats$proportion.upper)
data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower.stats, data.quantile.upper.stats, by="er.target")
#estimate the probability, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
VRAP2:::quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$prop.lower.mean<= probabilityAtzeroER.adj & data.quantile.stats$prop.upper.mean<=probs.row["prop.upper"],]
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target, na.rm = TRUE)))
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
data.quantile.stats$probabilityAtzeroER <- probs.df$probabilityAtzeroER[1]
data.quantile.stats$prop.upper <- probs.df$prop.upper[1]
RER.empirical <- data.quantile.raw[data.quantile.raw$proportion.lower<=rer$prob.lower & data.quantile.raw$proportion.upper<=rer$prob.upper,]
return(list(rer=rer, rer.detailed=results.rer, data.quantile.stats=data.quantile.stats, data.quantile.raw=data.quantile.raw, RER.empirical=RER.empirical))
}, probs.df)
return(results.rer.list)
}#END calcRER2
#' @title Calculate RER values from VRAP output (binomial method)
#'
#' @param x not sure about this.
#'
#' @return A data frame.
#'
#' @description Not sure about the degree of completeness of this function.
#'
#' @export
#'
#' @examples
#' \dontrun{
#' filepaths <- list.files(pattern = "*\\.rds$", full.names = TRUE, ignore.case = TRUE)
#' res.allscenarios.list <- lapply(filepaths, function(x){
#' vrap.output <- readRDS(x)
#' scenario <- names(vrap.output)
#' res.list <- extractQualifiedRuns(vrap.output)
#' return(list(qualified.reps=res.list[[1]], scenario=scenario, nruns=vrap.output[[1]]$inputs$NRuns))
#' })
#'
#' res <- calcRERbinomial(res.allscenarios.list)
#' }
calcRERbinomial <- function(x){
qualified.reps.list <- lapply(x, function(singlescenario) singlescenario$qualified.reps)
qualified.reps.df <- do.call("rbind", qualified.reps.list)
proportions.success <- aggregate(rep~er.target, data = qualified.reps.df[,c("rep","er.target")], length)
colnames(proportions.success)[colnames(proportions.success)=="rep"] <- "success.n"
nruns.total <- sum(sapply(x, function(singlescenario) singlescenario$nruns))
proportions.success$nruns.total <- nruns.total
rer.binomial <- apply(proportions.success, 1, function(x){
calcBinomialCI(k = x["success.n"], n = x["nruns.total"])
})
rer.binomial <- do.call("rbind", rer.binomial)
rer.binomial <- cbind(proportions.success, rer.binomial)
row.names(rer.binomial) <- NULL
return(rer.binomial)
}#END calcRERbinomial
#' @title Calculate RER values from VRAP output (binomial method)
#'
#' @param vrap.output A list. Output from the \code{lapply} application on
#' \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several data frames. The elemented named \code{rer}
#' has the rer estimates. Other data frames include columns of mean and SD of
#' the proportions, by target ER. These statistics should not be trusted as
#' the frequency distributions of proportions, by targe ER, is never normal.
#' (This was discovered after the fact.) This can be confirmed by creating a
#' historgram, by target ER of either
#' \code{data.quantile.raw$proportion.lower} or
#' \code{data.quantile.raw$proportion.upper}.
#'
#'
#' @examples
calcRERbinomial.old <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
#I think this was the first version and is a combo of the functions:
#extractQaulifiedRuns and calcRERbinomial
# updates to extractQaulifiedRuns have happended such that this is now behind.
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
#test if it has 'rep' variable, if not then create it:
if(!"rep" %in% names(sim.ind$output$rawoutput[[1]])){
for(i in 1:length(sim.ind$output$rawoutput)){
sim.ind$output$rawoutput[[i]]$rep <- i
}
}
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the thresholds
#lower esc probabilities:
data.quantile.lower <- aggregate(escapement~rep+er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#the distribution of probabilities, by ER, may be lognormal
#data.quantile.lower$proportion.lower[data.quantile.lower$proportion.lower==0] <- 1e-10
#data.quantile.lower$proportion.lower <- log(data.quantile.lower$proportion.lower)
data.quantile.lower.stats <- aggregate(proportion.lower~er.target, data=data.quantile.lower, FUN=function(x){c(prop.lower.mean=(mean(x)), prop.lower.sd=(sd(x)))})
data.quantile.lower.stats <- data.frame(er.target=data.quantile.lower.stats$er.target, data.quantile.lower.stats$proportion.lower)
#upper esc probabilities:
data.quantile.upper <- aggregate(escapement~rep+er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement'] <- 'proportion.upper'
#the distribution of probabilities, by ER, may be lognormal
# data.quantile.upper$proportion.upper[data.quantile.upper$proportion.upper==0] <- 1e-6
# data.quantile.upper$proportion.upper <- log(data.quantile.upper$proportion.upper)
data.quantile.upper.stats <- aggregate(proportion.upper~er.target, data=data.quantile.upper, FUN=function(x){c(prop.upper.mean=(mean(x)), prop.upper.sd=(sd(x)))})
data.quantile.upper.stats <- data.frame(er.target=data.quantile.upper.stats$er.target, data.quantile.upper.stats$proportion.upper)
data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower.stats, data.quantile.upper.stats, by="er.target")
#estimate the probability, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
VRAP2:::quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
# continue from here.
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$prop.lower.mean<= probabilityAtzeroER.adj & data.quantile.stats$prop.upper.mean<=probs.row["prop.upper"],]
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target)))
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
data.quantile.stats$probabilityAtzeroER <- probs.df$probabilityAtzeroER[1]
data.quantile.stats$prop.upper <- probs.df$prop.upper[1]
RER.empirical <- data.quantile.raw[data.quantile.raw$proportion.lower<=rer$prob.lower & data.quantile.raw$proportion.upper<=rer$prob.upper,]
proportions.success <- aggregate(rep~er.target, data = RER.empirical[,c("rep","er.target")], length)
colnames(proportions.success)[colnames(proportions.success)=="rep"] <- "success.n"
proportions.success$nruns <- vrap.output[[1]]$inputs$NRuns
#proportions.success$success.prop <- proportions.success$success.n/proportions.success$nruns
rer.binomial <- apply(proportions.success, 1, function(x){
calcBinomialCI(k = x["success.n"], n = x["nruns"])
})
rer.binomial <- do.call("rbind", rer.binomial)
rer.binomial <- cbind(proportions.success, rer.binomial)
return(list(rer.binomial=rer.binomial, rer.detailed=results.rer, data.quantile.stats=data.quantile.stats, data.quantile.raw=data.quantile.raw))
}, probs.df)
return(results.rer.list)
}#END calcRERbinomial2
#' @title Calculate RER values from summary table
#'
#' @param vrap.output
#' @param prop.lower.tolerance
#' @param prop.upper
#'
#' @return
#' @export
#'
#' @examples
calcRERfromsummary <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance=100*prop.lower.tolerance, prop.upper=100*prop.upper)
results <- lapply(vrap.output, FUN = function(x){
output.tabulated <- x$data.summary$output.tabulated
EndAv <- x$inputs$EndAv
NYears <- x$inputs$NYears
NRuns <- x$inputs$NRuns
metadata <- data.frame(NRuns=NRuns, NYears=NYears, EndAv=EndAv)
#this apply will return a list with length=#probability levels assessed. (usually=1)
apply.res <- apply(probs.df,1, FUN=function(probs.single, output.tabulated, metadata){
prop.lower <- probs.single["prop.lower.tolerance"]+output.tabulated$percent.years.belowLEL[output.tabulated$TgtER==0]
qualifyingER <- output.tabulated[output.tabulated$percent.years.belowLEL< prop.lower & output.tabulated$percent.years.aboveUEL > 100- probs.single["prop.upper"],]
if(nrow(qualifyingER)==0){
#no ER levels qualify
rer <- data.frame(metadata, lower.percent=prop.lower, upper.percent=100- probs.single["prop.upper"], rer=NA)
}else{
rer <- data.frame(metadata, lower.percent=prop.lower, upper.percent=100- probs.single["prop.upper"], rer=max(qualifyingER$TgtER, na.rm = TRUE))
}
row.names(rer) <- NULL
return(list(rer=rer, qualifyingER=qualifyingER))
}, output.tabulated, metadata)
})
#results <- results[[1]]
names(results) <- names(vrap.output)
return(results)
}#END calcRERfromsummary
#' @title Extract qualified runs from VRAP output (for binomial method)
#'
#' @param vrap.output A list. Output from the \code{lapply} application on
#' \code{\link{Main}}.
#' @param prop.lower.tolerance A numeric vector. The proportion values to be
#' added to the estimated proportion of years below LEL under zero ER.
#' Default is 0.05
#' @param prop.upper A numeric vector. The percent values to appraise frequency
#' of years with escapement above the upper escapement threshold. Default is
#' 0.20, which equates to the 80\% requirement.
#'
#' @return A list comprising several data frames. The elemented named \code{rer}
#' has the rer estimates. Other data frames include columns of mean and SD of
#' the proportions, by target ER. These statistics should not be trusted as
#' the frequency distributions of proportions, by targe ER, is never normal.
#' (This was discovered after the fact.) This can be confirmed by creating a
#' historgram, by target ER of either
#' \code{data.quantile.raw$proportion.lower} or
#' \code{data.quantile.raw$proportion.upper}.
#'
#' @export
#'
#' @examples
extractQualifiedRuns <- function(vrap.output,prop.lower.tolerance=0.05, prop.upper=0.20 ){
probs.df <- data.frame(prop.lower.tolerance, prop.upper)
data.esc.list <- lapply(vrap.output, FUN = function(sim.ind){
#test if it has 'rep' variable, if not then create it:
if(!"rep" %in% names(sim.ind$output$rawoutput[[1]])){
for(i in 1:length(sim.ind$output$rawoutput)){
sim.ind$output$rawoutput[[i]]$rep <- i
}
}
data.esc <- lapply(sim.ind$output$rawoutput, FUN=function(run.ind){
matrix(c(rep(run.ind$TargetU, length(run.ind$CalendarHR)), run.ind$CalendarHR, rep(run.ind$rep, length(run.ind$CalendarHR)), 1:length(run.ind$CalendarHR), run.ind$TotEscpmnt, run.ind$Recruits), nrow=length(run.ind$CalendarHR))
})
data.esc.df <- do.call("rbind", data.esc)
data.esc.df <- as.data.frame(data.esc.df)
colnames(data.esc.df) <- c("er.target", "er.realized", "rep", "year", "escapement", "recruits")
return(list(scenario=sim.ind$inputs$filename, lel=sim.ind$inputs$ECrit, uel=sim.ind$inputs$ERecovery, NYears=sim.ind$inputs$NYears, EndAv=sim.ind$inputs$EndAv, data.esc.df=data.esc.df))
})
names(data.esc.list) <- names(vrap.output)
results.rer.list <- lapply(data.esc.list, FUN = function(data.esc.element, probs.df){
data.esc.df <- data.esc.element$data.esc.df
#!! these proportions values are for BELOW the thresholds
#lower esc probabilities:
data.quantile.lower <- aggregate(escapement~rep+er.target, data=data.esc.df, FUN = function(x, value=data.esc.element$lel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.lower)[colnames(data.quantile.lower)=='escapement'] <- 'proportion.lower'
#the distribution of probabilities, by ER, may be lognormal
#data.quantile.lower$proportion.lower[data.quantile.lower$proportion.lower==0] <- 1e-10
#data.quantile.lower$proportion.lower <- log(data.quantile.lower$proportion.lower)
data.quantile.lower.stats <- aggregate(proportion.lower~er.target, data=data.quantile.lower, FUN=function(x){c(prop.lower.mean=(mean(x)), prop.lower.sd=(sd(x)))})
data.quantile.lower.stats <- data.frame(er.target=data.quantile.lower.stats$er.target, data.quantile.lower.stats$proportion.lower)
#upper esc probabilities:
data.quantile.upper <- aggregate(escapement~rep+er.target, data=data.esc.df[data.esc.df$year> (data.esc.element$NYears - data.esc.element$EndAv) ,], FUN = function(x, value=data.esc.element$uel){VRAP2:::quantInv(distr = x, value = value) })
colnames(data.quantile.upper)[colnames(data.quantile.upper)=='escapement'] <- 'proportion.upper'
#the distribution of probabilities, by ER, may be lognormal
# data.quantile.upper$proportion.upper[data.quantile.upper$proportion.upper==0] <- 1e-6
# data.quantile.upper$proportion.upper <- log(data.quantile.upper$proportion.upper)
data.quantile.upper.stats <- aggregate(proportion.upper~er.target, data=data.quantile.upper, FUN=function(x){c(prop.upper.mean=(mean(x)), prop.upper.sd=(sd(x)))})
data.quantile.upper.stats <- data.frame(er.target=data.quantile.upper.stats$er.target, data.quantile.upper.stats$proportion.upper)
data.quantile.raw <- merge(data.quantile.lower, data.quantile.upper, by=c("rep", "er.target"))
#combine lower and upper proportion estimates:
data.quantile.stats <- merge(data.quantile.lower.stats, data.quantile.upper.stats, by="er.target")
#estimate the probability, in each run, of being below LEL under no fishing:
data.quantile.lower.zeroER <- aggregate(escapement~rep, data=data.esc.df[data.esc.df$er.target==0,], FUN = function(x, value=data.esc.element$lel){
VRAP2:::quantInv(distr = x, value = value)
}, data.esc.element$lel)
#obtain the mean across all runs:
#the column is called 'escapement' as it's the output from the aggregate fn above that uses escapement data.
probs.df$probabilityAtzeroER <- mean(data.quantile.lower.zeroER$escapement)
results.rer <- apply(probs.df, 1, FUN=function(probs.row, data.quantile.stats){
probabilityAtzeroER.adj <- probs.row["probabilityAtzeroER"]+ probs.row["prop.lower.tolerance"]
results.rer <- data.quantile.stats[data.quantile.stats$prop.lower.mean<= probabilityAtzeroER.adj & data.quantile.stats$prop.upper.mean<=probs.row["prop.upper"],]
#only do this if there is an rer for this simulation:
if(nrow(results.rer)>0){
results.rer$prob.lower <- probabilityAtzeroER.adj
results.rer$prob.upper <- probs.row["prop.upper"]
} else {
results.rer$prob.lower <- numeric(0)
results.rer$prob.upper <- numeric(0)
}
return(results.rer)
}, data.quantile.stats)
rer <- lapply(results.rer, FUN = function(x) {
if(nrow(x)>0) {return(data.frame(prob.lower=unique(x$prob.lower), prob.upper=unique(x$prob.upper), rer=max(x$er.target, na.rm = TRUE)))}
})
rer <- do.call("rbind", rer)
results.rer <- do.call("rbind", results.rer)
data.quantile.stats$probabilityAtzeroER <- probs.df$probabilityAtzeroER[1]
data.quantile.stats$prop.upper <- probs.df$prop.upper[1]
RER.empirical <- data.quantile.raw[data.quantile.raw$proportion.lower<=rer$prob.lower & data.quantile.raw$proportion.upper<=rer$prob.upper,]
#this traps cases when no ER achieves criteria and nrow=0
if(nrow(RER.empirical)>0){
RER.empirical.sums <- aggregate(rep~er.target, data=RER.empirical, length)
} else {
RER.empirical.sums <- data.frame(rep=integer(), er.target=numeric())
}
colnames(RER.empirical.sums)[colnames(RER.empirical.sums)=="rep"] <- "success.n"
dat.samplesize <- aggregate(rep~er.target,data= data.quantile.raw, length)
RER.empirical.sums <- merge(dat.samplesize, RER.empirical.sums, all.x = TRUE)
RER.empirical.sums$success.n[is.na(RER.empirical.sums$success.n)] <- 0
colnames(RER.empirical.sums)[colnames(RER.empirical.sums)=="rep"] <- "reps.n"
return(list(reps.qualified =RER.empirical, reps.summary=RER.empirical.sums, scenario=data.esc.element$scenario))
}, probs.df)
return(results.rer.list)
}#END extractQualifiedRuns
#' @title Read in a single .rav2 file
#'
#' @param filename A character vector of length one. The name of the .rav2 file.
#'
#' @return A list of two elements named \code{inputs.vars} and \code{inputs}.
#' The element \code{inputs.vars} is a list of just the variables and their
#' values. The element named \code{inputs} is a list with the variable names,
#' values, and description. The length of each list (\code{inputs.vars} and
#' \code{inputs}) is equal to the number of rows in the .rav2 file.
#' @export
#'
#' @details The .rav2 file is a csv file. The first column is the variable name,
#' the second column is the value for that variable, the third column is a
#' description of the variable. If commas are included in either the value or
#' descripton then the full column should be inside double quotes. If the
#' value is a numeric vector of several values, the values should be separated
#' by a space (not a comma) and not wrapped in quotes. Additional columns may
#' be added after the description column. They could include additional
#' comments and information. They will be represented, by variable, in the
#' list named \code{inputs}.
#'
#' @examples
read_rav2 <- function(filename){
x <- scan(filename, what="",quote = "", sep="\n")
x <- x[!grepl(pattern = "^#", x)]
# Separate elements by comma
#y <- strsplit(x, ",")
#this allows for input of strings that have commas
#if using commas, the full string must be in double quotes:
y <- strsplit(x, '\\"[^"]+,(*SKIP)(*FAIL)|,\\s*', perl=TRUE)
# Extract the first vector element and set it as the list element name
names(y) <- trimws(sapply(y, `[[`, 1))
# Remove the first vector element from each list element
y <- lapply(y, `[`, -1)
inputs <- lapply(y, FUN = function(x){
value <- trimws(x[1])
#split on one or more spaces:
value <- strsplit(x = value, "\\s+")[[1]]
value <- type.convert(value, as.is = TRUE)
#if the split data are characters then use the original form (unsplit):
if(any(class(value)=="character")) {
value <- trimws(x[1])
value <- gsub(pattern = "\\\"", replacement = "", x = value)
}
res <- list(value=value)
#if there's a third column in the csv file, then it's the description
#(index=2 as first column removed):
x.count <- length(x)
if(x.count>1) res$description=trimws(gsub(pattern = "\\\"", replacement = "", x = x[2]))
if(x.count>2){
#additional user-specific columns
for(i in 3:x.count){
element.name <- paste0("additional",i-2)
strsplit(x, '\\"[^"]+,(*SKIP)(*FAIL)|,\\s*', perl=TRUE)
res[element.name]=trimws(gsub(pattern = "\\\"", replacement = "", x = x[i]))
}
}
return(res)
})
#integer of the number of steps of ER or Pop Capacity to do # is 1:BufMax
inputs$BufMax <- list(value = round((inputs$BufferEnd$value - inputs$BufferStart$value) / inputs$BufferStep$value + 1), description= "integer of the number of steps of ER or Pop Capacity to do" )
inputs.vars <- lapply(inputs, "[[","value")
names(inputs.vars) <- names(inputs)
return(list(inputs.vars=inputs.vars, inputs=inputs))
}#END read_rav2
printInput <- function(filename){
input <- GetInput(filename)
}#END printInput
quantInv <- function(distr, value) {ecdf(distr)(value)}
#' @title Build multiple copies of same .rav file for sensitivity tests
#'
#' @param rav.filenames A charactdata.quantile.upper vector. The input .rav files that will be
#' copied.
#' @param sensitivity.runs A integer of length one. The number of copies to make
#' of each .rav file.
#'
#' @return Each input .rav file is copied \code{sensitivity.runs} number of
#' times. The filename of each of file is appended with a numerical sequence
#' from 1 to \code{sensitivity.runs}. The original files are deleted as all
#' copies are exactly same as the original.
#' @export
#'
#' @examples
sensitivityBuild <- function(rav2.filenames, sensitivity.runs){
invisible(
lapply(rav2.filenames, FUN = function(filename, sensitivity.runs){
run.vec <- (1:sensitivity.runs) / 10^(1+floor(log10(sensitivity.runs)))
run.vec <- formatC(run.vec,format='f', digits = 1+floor(log10(sensitivity.runs)))
run.vec <- substr(run.vec, 3, nchar(run.vec))
lapply(run.vec, FUN = function(run.single, filename){
file.copy(filename, paste0(tools::file_path_sans_ext(filename), "sensitivity.run_", run.single, ".", tools::file_ext(filename)))
}, filename)
}, sensitivity.runs)
)
#remove originals:
invisible(lapply(rav2.filenames, FUN = function(x) file.remove(x)))
}#END sensitivityBuild
tabulateSummary <- function(vrap.output){
vrap.output <- lapply(vrap.output, FUN=function(x){
# browser()
x$data.summary[(1+grep(pattern = "__Escapement__", x = x$data.summary)):length(x$data.summary)]
})
}#END tabulateSummary
#' @title Write one or many .rav2 files, which hold VRAP input variables.
#'
#' @param indexfilename A Boolean, default is FALSE. The rav2 filename can
#' sometimes get too long for windows. If so, use indexfilename=TRUE, which
#' forces a date-time stamp (to 1/10,000 of a second) to be used for a
#' filename. Additionally a translation table is written that will show both the time
#' stamp filename and the longer filename, which includes arguments.
#' @param filename The name of the .rav2 file to be written. If NA (default)
#' then the function will create a file name based on the argument values
#' choosen.
#' @param defaultargs.filename The name of a .rav2 file that includes default
#' values. If this is chosen then the default values in this function are
#' ignored. But values set in the function arguments are not ignored.
#' @param Title
#' @param RanSeed
#' @param NRuns
#' @param NYears
#' @param MinAge
#' @param MaxAge
#' @param ConvergeCrit
#' @param Debugg
#' @param SRType
#' @param BSRa
#' @param BSRb
#' @param BSRc
#' @param BSRd
#' @param AveEnv
#' @param depen
#' @param DL1
#' @param DL2
#' @param DR
#' @param EscChoice
#' @param SurvScale
#' @param SRError.FUN
#' @param SRErrorA
#' @param SRErrorB
#' @param ResCorParameter
#' @param MarSurv
#' @param BetaMarA
#' @param BetaMarB
#' @param CorrMS
#' @param FWError.FUN
#' @param GammaFlowA
#' @param GammaFlowB
#' @param MarError.FUN
#' @param GammaMarA
#' @param GammaMarB
#' @param NumBreakPoints
#' @param BaseRegime
#' @param EscpmntBreakPoint
#' @param TargetU
#' @param MgmtError
#' @param MgmtError.FUN
#' @param GammaMgmtA
#' @param GammaMgmtB
#' @param ECrit
#' @param ERecovery
#' @param EndAv
#' @param StepFunc
#' @param BufferStep
#' @param BufferStart
#' @param BufferEnd
#' @param CohortStart
#' @param NatMort
#' @param MatRate
#' @param PTU
#' @param MatU
#' @param output.raw A character vector "AEQMort Escpmnt TotAEQMort
#' TotAdultEscpmnt TotEscpmnt RanFlow RanMarine Recruits TargetUError
#' ptu.scaled MatU.scaled Cohort tmpdat TempCohort CalendarHR TargetU" These
#' are the variables that will be included in the output.
#' @param ...
#'
#' @return Writes one or more .rav2 (text/csv) files depending on number of
#' unique combinations created from the input arguments.
#' @export
#'
#' @examples
writeRav2file <- function(
indexfilename=FALSE,
filename=NA,
defaultargs.filename=NA,
Title="ER example of Beverton-Holt VRAP with no marine survival or stream flow variability",
RanSeed=0,
NRuns=1000,
NYears=25,
MinAge=2,
MaxAge=5,
ConvergeCrit=0.001,
Debugg="NO",
SRType="BEV2",
BSRa=9.27599153996063,
BSRb=3702.14272977426,
BSRc=NA,
BSRd=NA,
AveEnv=NA,
depen="YES",
DL1=400,
DL2=63,
DR=0,
EscChoice="YES",
SurvScale="YES",
SRError.FUN="rgamma(1, shape=inputs$SRErrorA, scale=inputs$SRErrorB)",
SRErrorA=0.967637142385205,
SRErrorB=1.09375572895311,
ResCorParameter=NA,
MarSurv="NO",
BetaMarA=NA,
BetaMarB=NA,
CorrMS=NA,
FWError.FUN="rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)",
GammaFlowA=0,
GammaFlowB=0,
MarError.FUN="rgamma(1, inputs$GammaMarA, scale=inputs$GammaMarB)",
GammaMarA=0,
GammaMarB=0,
NumBreakPoints=0,
BaseRegime=1,
EscpmntBreakPoint=NA,
TargetU=0.37,
MgmtError="YES",
MgmtError.FUN="rgamma(1, shape=inputs$GammaMgmtA, scale=inputs$GammaMgmtB)",
GammaMgmtA=65.3946,
GammaMgmtB=0.0158,
ECrit=416,
ERecovery=1050,
EndAv=5,
StepFunc="ER",
BufferStep=6.75675675675676E-02,
BufferStart=0,
BufferEnd=1.89189189189189,
CohortStart="5108 2554 1478 667 112",
NatMort="0.5 0.4 0.3 0.2 0.1",
MatRate="0 1.43432324005842E-02 0.226067824725145 0.685044912939359 1",
PTU="0 7.97982424808335E-02 0.234143029720483 0.352056542423367 0.33478625941973",
MatU="0 0.315672687185273 0.331197585913413 0.339907561109324 0.32870153694989",
output.raw= "AEQMort Escpmnt TotAEQMort TotAdultEscpmnt TotEscpmnt RanFlow RanMarine Recruits TargetUError ptu.scaled MatU.scaled Cohort tmpdat TempCohort CalendarHR TargetU",
...)
{
rav.df <- as.data.frame(matrix(ncol=3, byrow=TRUE,
data=c(
"Title", "argList$Title", "Title",
"RanSeed", "argList$RanSeed", "Random seed; 0 gives random seed; numbers give fixed seed",
"NRuns", "argList$NRuns", "Number of runs",
"NYears", "argList$NYears", "Number of years",
"MinAge", "argList$MinAge", "Minimum and maximum age",
"MaxAge", "argList$MaxAge", "Minimum and maximum age (for now this is fixed; do not change)",
"ConvergeCrit", "argList$ConvergeCrit", "Convergence criterion (proportional error) for target ER. converged when: abs((estimated ER- target ER)/target ER) < ConvergeCrit",
"Debugg", "argList$Debugg", "Debug file flag",
"SRType", "argList$SRType", "Spawner Recruit function (RIC2;RIC3;RIC4; BEV2;BEV3;BEV4; HOC2;HOC3;HOC4)",
"BSRa", "argList$BSRa", "S/R a parameter",
"BSRb", "argList$BSRb", "S/R b parameter",
"BSRc", "argList$BSRc", "marine survival parameter",
"BSRd", "argList$BSRd", "freshwater survival parameter",
"AveEnv", "argList$AveEnv",NA,
"#Mean and CV for M marine survival index (M^c)",NA,NA,
"#Trend; Cycle; or Autoc(orrelation) ?",NA,NA,
"#Trend/Cycle parameters: rate and 0 or increment and 1 for trend- amplitude- period & starting pt & scalar for cycle; correl and starting point (if zero is random) for autocorrelation",NA,NA,
"#Mean and CV for F flow (or other fw) index (exp(dF))",NA,NA,
"#Trend; Cycle; or Autoc(orrelation) ?",NA,NA,
"#Trend/Cycle parameters: rate and 0 or increment and 1 for trend- amplitude- period & starting pt & scalar for cycle; correl and starting point (if zero is random) for autocorrelation",NA,NA,
"FWError.FUN", "argList$FWError.FUN", "The R probability density function with arguments",
"GammaFlowA", "argList$GammaFlowA","FW (flow) error pdf parameter",
"GammaFlowB", "argList$GammaFlowB","FW (flow) error pdf parameter",
"MarError.FUN", "argList$MarError.FUN", "The R probability density function with arguments",
"GammaMarA", "argList$GammaMarA","Marine error pdf parameter",
"GammaMarB", "argList$GammaMarB","Marine error pdf parameter",
"depen", "argList$depen", "Depensation? (\"YES\" or \"NO\")",
"DL1", "argList$DL1", "1) Esc. level for depensation to start",
"DL2", "argList$DL2", "2) QET",
"DR", "argList$DR", "3)% predicted return at QET (or for r/s=1 this parameter = 1)",
"EscChoice", "argList$EscChoice", "Determine recruits from adult spawners (not total)? (\"YES\" or \"NO\")",
"SurvScale", "argList$SurvScale", "Stock-recruit variation (\"YES\" or \"NO\")",
"SRError.FUN", "argList$SRError.FUN", "The R probability density function with arguments",
"SRErrorA", "argList$SRErrorA", "error of S/R a parameter",
"SRErrorB", "argList$SRErrorB", "error of S/R b parameter",
"ResCorParameter", "argList$ResCorParameter", "S/R parameter error autocorrelation",
"MarSurv", "argList$MarSurv", "Smolt to adult survival w/variation (\"YES\" or \"NO\")",
"BetaMarA", "argList$BetaMarA", "Beta distribution a parameter",
"BetaMarB", "argList$BetaMarB", " Beta distribution b parameter",
"CorrMS", "argList$CorrMS", "Beta distribution parameter error autocorrelation",
"NumBreakPoints", "argList$NumBreakPoints", "Number of breakpoints",
"BaseRegime", "argList$BaseRegime", "Level to use as base regime",
"TargetU", "argList$TargetU", "base exploitation rate",
"MgmtError", "argList$MgmtError", "Include error (\"YES\" or \"NO\") in ER management",
"MgmtError.FUN", "argList$MgmtError.FUN", "The R probability density function with arguments",
"GammaMgmtA", "argList$GammaMgmtA", "Management error gamma pdf shape parameter",
"GammaMgmtB", "argList$GammaMgmtB", "Management error gamma pdf scale parameter",
"ECrit", "argList$ECrit", "Lower escapement threshold",
"ERecovery", "argList$ERecovery", "Upper escapement threshold (MSY)",
"EndAv", "argList$EndAv", "# yrs to average for UEL",
"StepFunc", "argList$StepFunc", "Step ER (ER) or Pop Capacity (Pop)?",
"BufferStep", "argList$BufferStep", "Buffer step size as percent of base ER or Pop capacity",
"BufferStart", "argList$BufferStart", "ER Min as proportion of base ER",
"BufferEnd", "argList$BufferEnd", "ER max as proportion of base ER",
"CohortStart", "argList$CohortStart", "Initial population size, by Age, of start year (prior to nat mort)",
"NatMort", "argList$NatMort", "Age specific natural mortality",
"MatRate", "argList$MatRate", "Age specific average maturation rate",
"PTU", "argList$PTU", "Age specific average mixed-maturity fishery fishing rates",
"MatU", "argList$MatU", "Age specific average mature fishery fishing rates")
), stringsAsFactors = FALSE)
argList <- as.list(sys.call())[-1]
if(!is.na(defaultargs.filename)){
#import the default args:
argList <- argList[names(argList) !="defaultargs.filename"]
argList.default <- read_rav2(defaultargs.filename)$inputs.vars
} else {
argList.default <- formals(writeRav2file)
argList.default <- argList.default[names(argList.default) !="..."]
}
#The argList needs to have something in it for following code to work. Otherwise no rav files are created:
if(length(argList)==0) argList$NRuns <- 1000
argLength <- lapply(argList, length)
testing.var.names <- names(which(unlist(argLength)>1))
if(length(testing.var.names)>1) cat("\n!!more than one variable being tested!!\n")
sensitivity.args <- do.call("expand.grid", argList)
#sensitivity.args is a data frame representing all combinations of the args used in the function call
apply(sensitivity.args, 1, FUN=function(x){
if(nrow(sensitivity.args)>1 | is.na(filename)){
#can't use default filename as creating more than one file,
#or no filename supplied
filename <- paste(names(x),x, sep = "_", collapse = "__")
filename <- paste0(filename, ".rav2")
}
#the filename can sometimes get too long for windows. if so, use
#indexfilename=TRUE, then a date-time stamp is used for a filename, and a
#translation table is written that will show both the time stamp filename
#and the longer filename, which includes arguments.
#must set filename here so the correct filename is included in rav2 contents
if(indexfilename){
filename.long <- filename
options(digits.secs=4)
filename <- paste0(format(Sys.time(), "%Y-%m-%d-%H-%M-%OS"), ".rav2")
translation.tbl.filename <- "rav2_filename_translation.csv"
if(file.exists(translation.tbl.filename)){
write.table(x = cbind(filename, filename.long), file = translation.tbl.filename, quote = FALSE, append = TRUE, sep=", ", row.names = FALSE, col.names = FALSE)
}else{
write.table(x = cbind(filename, filename.long) , col.names = c("index.filename", "long.filename"), file = translation.tbl.filename, quote = FALSE, sep=", ", row.names = FALSE)
}
}#END indexfilename
argList <- c(list(filename=filename), x, argList.default[!names(argList.default) %in% names(x)])
rav.str <- rep(NA, length(argList))
for(i in 1:length(argList)){
variable <- names(argList[i])
value <- paste(argList[[i]], collapse = " ")
if(grepl(pattern = ",", x = value)) value <- paste0("\"",value,"\"")
description <- rav.df[rav.df[,1]==names(argList[i]),3]
if(isTRUE(grepl(pattern = ",", x = description))) description <- paste0("\"",description,"\"")
rav.str[i] <- paste(variable, value, description, sep=", ")
}
write(rav.str, file=filename)
} )
cat(nrow(sensitivity.args), "rav2 files written")
}#END writeRav2file
#' @title Build, save, and open an R script to help run VRAP.
#'
#' @description This creates and opens a script named "vrap_script.R". This is a
#' template for the user to work with when running vrap without shiny. This is
#' intended to help the user understand the work flow, but due to file path
#' differences it may not work as is. Some object values will need updating
#' (for example the path).
#'
#' @return Opens an R script that includes the template of functions.
#' @export
#'
#' @examples
#' writeScriptVRAP()
writeScriptVRAP <- function(){
script.str <- c("
#delete old rav files:
rav2.filenames <- list.files(pattern = \"*.rav2$\")
file.remove(rav2.filenames)
writeRav2file(NYears = c(10,25) )
writeRav2file(NRuns = c(3000))
#this makes copies of .rav2 files to allow for testing on what is an appropriate value to get stable estimates of rer. sensitivity.runs is the number of copies to make of each rav file:
rav2.filenames <- list.files(pattern = \"*.rav2$\")
sensitivityBuild(rav2.filenames = rav2.filenames, sensitivity.runs = 10)
#list rav files available:
rav2.filenames <- list.files(pattern = \"*.rav2$\")
#run vrap
#a .rds file will be saved for each .rav2 file tested
#the .rds file prefix will match that of the .rav2 file.
lapply(rav2.filenames, FUN=function(x) VRAP2:::Main(x))
###run vrap in parallel:
require(parallel)
no_cores <- detectCores() - 1
# Initiate cluster
cl <- makeCluster(no_cores)
obj.names <- ls()
clusterExport(cl, obj.names)
clusterEvalQ(cl, library(VRAP2))
rav2.filenames <- list.files(pattern = \"*.rav2$\")
parLapply(cl, rav2.filenames, fun=function(x) VRAP2:::Main(x))
stopCluster(cl)
###End of parallel sims
vrap.output1 <- readRDS(\"NRuns_3000sensitivity.run_1.rds\")
vrap.output2 <- readRDS(\"NRuns_3000sensitivity.run_2.rds\")
vrap.output3 <- readRDS(\"NRuns_3000sensitivity.run_3.rds\")
vrap.output <- c(vrap.output1, vrap.output2, vrap.output3)
names(vrap.output)
rer.orig <- calcRERfromsummary(vrap.output = vrap.output, prop.lower.tolerance = 0.05)
cat(vrap.output[[1]]$data.summary$output)
rer.orig[[1]]$rer
rer.modified <- calcRER2(vrap.output = vrap.output, percent.lower = c(.05, 0.1))
rer.modified[[1]]$rer
hist(rer.modified[[1]]$RER.empirical$er.target)
#obtain the rer's for all simulations:
lapply(vrap.output, \"[[\", \"rer\")
")
write(script.str, file="vrap_script.R")
file.edit("vrap_script.R" )
}#writeScriptVRAP
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.