R/Retrospective/rho/03_Mohns_rho.r
In mkapur/ssdiags:
Defines functions
mohns.rho
##
# DC Edits 2/19/2019 Mohn's Rho Adapted from Felipe Carvalho's Excel file
# C:\000\1004_ICCAT_SFM_2019\01_2019_Meeting\SSv324U_02_Diagnostics\2019_04_Dean_reproduce_Diag_run_1\03_code_Mahns_Rho
# Retrospective error_HurtadoAverage_DC_Edits.xls
# see
# Hurtado-Ferro_et al_ICES_DC_Notes_Diagnostics.pdf
# Carvalho_et_al_2017_Fish_Res_Diagnostics_DC_Notes.pdf
#
# DC Note: for plots below and calculations below:
# Data from 1951-2015 are available for all retrospective model runs
# Estimation of terminal year(s) turned off as indicated in starter in starter file
# Need to manually grab the correct endyrvec=c(2015,2014,2013,2012,2011,2010)
# and set the correct legend lables
# e.g.,
# endyrvec=c(2015,2014,2013,2012,2011,2010)
# legendlabels <- c('2015','2014','2013','2012','2011','2010')
##
library(r4ss)
# working dir set in scource file
#dir.work <- "C:\\000\\1004_ICCAT_SFM_2019\\01_2019_Meeting\\SSv324U_02_Diagnostics\\2019_04_Dean_reproduce_Diag_run_1\\MAKO_ATL_DC_run_1\\Retrospective_try_05"
setwd(dirname.Retrospective)
getwd()
csv.yn <- T
plot.dir <-'plots_3'
Produce.plots <- T
read.outputs.from.R <- T #read outputs from 00_Merged_Retro_analysis_examples.R
#read.outputs.from.R <- F #read outputs from directory structure
if (read.outputs.from.R) {
# read outputs from 00_Merged_Retro_analysis_examples.R
outputs <- retroModels
run_sum <- retroSummary
end.yr.vec <- endyrvec
legendlabels <- as.character(end.yr.vec)
} else {
# read outputs from directory structure
#run dir
rundir <- paste(dirname.Retrospective,"retrospectives",sep="\\")
dirout <- file.path( paste0(rundir) )
#dir paths
dir_paths <- list.dirs(dirout, recursive=F)
mydirs <- dir_paths
mydirs
#[1] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-1"
#[2] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-2"
#[3] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-3"
#[4] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-4"
#[5] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-5"
#[6] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro0"
#order "plots"
mydirs <- dir_paths[c(6,1,2,3,4,5)]
mydirs
#end.yr.vec <- c(2015)
end.yr.vec <- c(2015,2014,2013,2012,2011,2010)
#legendlabels <- c('2015')
legendlabels <- c('2015','2014','2013','2012','2011','2010')
# read outputs
outputs <- SSgetoutput(dirvec=mydirs, getcovar=TRUE, forecast=TRUE)
# outputs <- SSgetoutput(dirvec=mydirs[1], getcovar=TRUE, forecast=TRUE)
# outputs <- SSgetoutput(dirvec=mydirs, getcovar=FALSE, forecast=TRUE)
run_sum <- SSsummarize(outputs)
setwd(dirname.Retrospective)
getwd() # should be back at base dir
#save(run_sum,outputs, grid, dirname.Retrospective, rundir, mydirs, file= file.path(dirname.Retrospective, "/run_sumdat.rdata") )
save(run_sum,outputs, dirname.Retrospective, rundir, file= file.path(dirname.Retrospective, "/run_sumdat.rdata") )
}
names(run_sum)
run_sum$modelnames
run_sum$quants
run_sum$quantsSD
setwd(dirname.Retrospective)
getwd() # should be back at base dir
if (!(plot.dir %in% dir())) {
dir.create(plot.dir)
}
#SSplotComparisons(run_sum, col=rich.colors.short(8) ,subplots=1 , legend=TRUE, endyrvec=end.yr.vec)
#SSplotComparisons(run_sum, col=rich.colors.short(8) ,subplots=2 , legend=TRUE, endyrvec=end.yr.vec)
#SSplotComparisons(run_sum, subplots=1, plot=T, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
#SSplotComparisons(run_sum, subplots=2, plot=T, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
SSplotComparisons(run_sum, plot=T, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
SSplotComparisons(run_sum, plot=T, pdf=T, png=F, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
#graphics.off()
outfile <- "SMA_SSF_vals.csv"
SSB_Virgin <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_Virgin"), ]#
SSB_Initial <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_Initial"), ]#
SSB_1950 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_1950"), ]#
SSB_2010 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2010"), ]#
SSB_2011 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2011"), ]#
SSB_2012 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2012"), ]#
SSB_2013 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2013"), ]#
SSB_2014 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2014"), ]#
SSB_2015 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2015"), ]#
SSB.retro <- run_sum$quants[as.numeric(row.names(SSB_Virgin)):as.numeric(row.names(SSB_2015)),] # note that this is spawning stock fecundity (1,000s of pups produced) and depends on how spawning units set up in control file
outfile_SD <- "SMA_SSF_SD.csv"
SSB_Virgin_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_Virgin"), ]#
SSB_Initial_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_Initial"), ]#
SSB_1950_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_1950"), ]#
SSB_2010_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2010"), ]#
SSB_2011_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2011"), ]#
SSB_2012_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2012"), ]#
SSB_2013_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2013"), ]#
SSB_2014_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2014"), ]#
SSB_2015_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2015"), ]#
SSB.retro.SD <- run_sum$quantsSD[as.numeric(row.names(SSB_Virgin_SD)):as.numeric(row.names(SSB_2015_SD)),] # note that this is spawning stock fecundity (1,000s of pups produced) and depends on how spawning units set up in control file
if (csv.yn) {
setwd(plot.dir)
write.table(SSB.retro, file = outfile, append = F, sep = ",", dec = ".", row.names = T, col.names = T)
write.table(SSB.retro.SD, file = outfile_SD, append = F, sep = ",", dec = ".", row.names = T, col.names = T)
setwd(dirname.Retrospective)
}
#----------------------------------------------------------------------------------------------------
# Begin DC Mohns rho Check
#----------------------------------------------------------------------------------------------------
SSB.retro
end.yr.vec
legendlabels
dim(SSB.retro)
SSB.retro.2015 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2015"), ]#
SSB.retro.2014 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2014"), ]#
SSB.retro.2013 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2013"), ]#
SSB.retro.2012 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2012"), ]#
SSB.retro.2011 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2011"), ]#
SSB.retro.2010 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2010"), ]#
retro.i <- 0
mohns.rho.i <- rep(0,length(end.yr.vec[-1]))
SSB.i <- rep(0,length(end.yr.vec[-1]))
SSB.retro.i <- rep(0,length(end.yr.vec[-1]))
for (retro.yr in end.yr.vec[-1]) {
retro.i = retro.i +1
#retro.i <- 1
#retro.yr <- 2015
SSB.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[1]]
SSB.retro.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[retro.i+1]]
mohns.rho.i[retro.i] <- (SSB.retro.i[retro.i] - SSB.i[retro.i])/SSB.i[retro.i]
}
check.mohns.rho.out <- mean(mohns.rho.i)
retro.i
retro.yr
mohns.rho.i
SSB.i
SSB.retro.i
#----------------------------------------------------------------------------------------------------
# End Dc Mohns rho check
#----------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------
# Begin FC Mohns rho function (with DC Edits to match r4ss function "SSmohnsrho" output "AFSC_Hurtado_SSB")
#----------------------------------------------------------------------------------------------------
mohns.rho<-function(x,y=end.yr.vec){
#x<-SSB.retro
#y<-end.yr.vec
retro.i <- 0
mohns.rho.i <- rep(0,length(y[-1]))
SSB.i <- rep(0,length(y[-1]))
SSB.retro.i <- rep(0,length(y[-1]))
mohns.rho.out <- 0
for (retro.yr in end.yr.vec[-1]) {
retro.i = retro.i +1
#retro.i <- 1
#retro.yr <- 2015
SSB.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[1]]
SSB.retro.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[retro.i+1]]
mohns.rho.i[retro.i] <- (SSB.retro.i[retro.i] - SSB.i[retro.i])/SSB.i[retro.i]
}
mohns.rho.out <- mean(mohns.rho.i)
retro.i
retro.yr
mohns.rho.i
SSB.i
SSB.retro.i
# print(Outs)
# return(Output)
return(list(SSB.i=c(SSB.i),
SSB.retro.i=c(SSB.retro.i),
mohns.rho.i=c(mohns.rho.i),
mohns.rho.out= mohns.rho.out))
}
#----------------------------------------------------------------------------------------------------
# End FC Mohns rho function
#----------------------------------------------------------------------------------------------------
#SSB.retro
get_mohns.rho <- mohns.rho(SSB.retro)
get_mohns.rho
check.mohns.rho.out
if (csv.yn) {
setwd(plot.dir)
write.table(get_mohns.rho, file="01_FC_get_mohns_rho.txt", row.names=F)
write.table(check.mohns.rho.out, file="02_DC_check_mohns_rho.txt", row.names=F)
setwd(dirname.Retrospective)
}
#----------------------------------------------------------------------------------------------------
# Begin Plots 1-4 : SSF > SSF_MSY and SSF/SSF_MSY > 1; Derived Quantity Estiamtes (mode) and 95% CI
#----------------------------------------------------------------------------------------------------
modeldesc <- legendlabels
nmodels <- dim(run_sum$quants)[2]-2
# Plots 1 and 2
# Stock Synthesis model run SSF/SSF_MSY
# and approximate 95% confidence intervals (95% CI = mode ± 1.96*std)
# using the maximum likelihood estimate (mode) and standard deviation (std)
# Plots 3 and 4
# Stock Synthesis model run SSF
# and approximate 95% confidence intervals (95% CI = mode ± 1.96*std)
# using the maximum likelihood estimate (mode) and standard deviation (std)
SSB_MSY <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_MSY"), ][1]#
SSB_VIRG <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_Virgin"), ][1]
msy_equilib <- SSB_MSY/SSB_VIRG
msy_equilib
#> msy_equilib
#replist1
#350 0.5212652
#
# Plot 1 Relative spawning biomass without error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=3 , legend=TRUE, legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
text( 1965, 1.05, "SSF_MSY", col=grey(0.5))
mtext('Derived Quantity Estimates', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSBwLINEmsy.png"), type='png' )
#graphics.off()
} #if plots statement end
#
# Plot 2 Relative spawning biomass with approximate 95% error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=4 , legend=TRUE, legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
text( 1965, 1.05, "SSF_MSY", col=grey(0.5))
mtext('Derived Quantity Estimates and Approximate 95% Confidence Intervals', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSBwLINEmsy_w_95_error.png"), type='png' )
#graphics.off()
} #if plots statement end
#
# Plot 3 Spawning biomass without error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=1 , legend=TRUE,
legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
abline(h=SSB_MSY/1, lty=3,col=grey(0.5), lwd=2 )
text( 1960, 1.05 *SSB_MSY, paste0("SSF_MSY = ",round(SSB_MSY/1000,2)), col=grey(0.5), pos=4)
text( 1960, 1.20 *SSB_MSY, paste0("Mohn's rho = ",round(get_mohns.rho$mohns.rho.out,3)), col=grey(0.5), pos=4)
mtext('Derived Quantity Estimates', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSB_total_wLINEmsy_Q50.png"), type='png' )
#graphics.off()
} #if plots statement end
#
# Plot 4 Spawning biomass with error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=2 , legend=TRUE,
legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
abline(h=SSB_MSY/1, lty=3,col=grey(0.5), lwd=2 )
text( 1960, 1.05 *SSB_MSY, paste0("SSF_MSY = ",round(SSB_MSY/1000,2)), col=grey(0.5), pos=4)
text( 1960, 1.20 *SSB_MSY, paste0("Mohn's rho = ",round(get_mohns.rho$mohns.rho.out,3)), col=grey(0.5), pos=4)
mtext('Derived Quantity Estimates and Approximate 95% Confidence Intervals', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSB_total_wLINEmsy_Q50_w_95_error.png"), type='png' )
#graphics.off()
} #if plots statement end
graphics.off()
#----------------------------------------------------------------------------------------------------
# End Plots 1-4 : SSF > SSF_MSY and SSF/SSF_MSY > 1; Derived Quantity Estiamtes (mode) and 95% CI
#----------------------------------------------------------------------------------------------------
mkapur/ssdiags documentation built on April 8, 2020, 7:17 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions mohns.rho
##
# DC Edits 2/19/2019 Mohn's Rho Adapted from Felipe Carvalho's Excel file
# C:\000\1004_ICCAT_SFM_2019\01_2019_Meeting\SSv324U_02_Diagnostics\2019_04_Dean_reproduce_Diag_run_1\03_code_Mahns_Rho
# Retrospective error_HurtadoAverage_DC_Edits.xls
# see
# Hurtado-Ferro_et al_ICES_DC_Notes_Diagnostics.pdf
# Carvalho_et_al_2017_Fish_Res_Diagnostics_DC_Notes.pdf
#
# DC Note: for plots below and calculations below:
# Data from 1951-2015 are available for all retrospective model runs
# Estimation of terminal year(s) turned off as indicated in starter in starter file
# Need to manually grab the correct endyrvec=c(2015,2014,2013,2012,2011,2010)
# and set the correct legend lables
# e.g.,
# endyrvec=c(2015,2014,2013,2012,2011,2010)
# legendlabels <- c('2015','2014','2013','2012','2011','2010')
##
library(r4ss)
# working dir set in scource file
#dir.work <- "C:\\000\\1004_ICCAT_SFM_2019\\01_2019_Meeting\\SSv324U_02_Diagnostics\\2019_04_Dean_reproduce_Diag_run_1\\MAKO_ATL_DC_run_1\\Retrospective_try_05"
setwd(dirname.Retrospective)
getwd()
csv.yn <- T
plot.dir <-'plots_3'
Produce.plots <- T
read.outputs.from.R <- T #read outputs from 00_Merged_Retro_analysis_examples.R
#read.outputs.from.R <- F #read outputs from directory structure
if (read.outputs.from.R) {
# read outputs from 00_Merged_Retro_analysis_examples.R
outputs <- retroModels
run_sum <- retroSummary
end.yr.vec <- endyrvec
legendlabels <- as.character(end.yr.vec)
} else {
# read outputs from directory structure
#run dir
rundir <- paste(dirname.Retrospective,"retrospectives",sep="\\")
dirout <- file.path( paste0(rundir) )
#dir paths
dir_paths <- list.dirs(dirout, recursive=F)
mydirs <- dir_paths
mydirs
#[1] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-1"
#[2] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-2"
#[3] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-3"
#[4] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-4"
#[5] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro-5"
#[6] "C:\\000\\1005_Diagnostics_pub\\04_Retrospective\\MAKO_ATL_r4ss_Retrospective_merged/Retrospective\\retrospectives/retro0"
#order "plots"
mydirs <- dir_paths[c(6,1,2,3,4,5)]
mydirs
#end.yr.vec <- c(2015)
end.yr.vec <- c(2015,2014,2013,2012,2011,2010)
#legendlabels <- c('2015')
legendlabels <- c('2015','2014','2013','2012','2011','2010')
# read outputs
outputs <- SSgetoutput(dirvec=mydirs, getcovar=TRUE, forecast=TRUE)
# outputs <- SSgetoutput(dirvec=mydirs[1], getcovar=TRUE, forecast=TRUE)
# outputs <- SSgetoutput(dirvec=mydirs, getcovar=FALSE, forecast=TRUE)
run_sum <- SSsummarize(outputs)
setwd(dirname.Retrospective)
getwd() # should be back at base dir
#save(run_sum,outputs, grid, dirname.Retrospective, rundir, mydirs, file= file.path(dirname.Retrospective, "/run_sumdat.rdata") )
save(run_sum,outputs, dirname.Retrospective, rundir, file= file.path(dirname.Retrospective, "/run_sumdat.rdata") )
}
names(run_sum)
run_sum$modelnames
run_sum$quants
run_sum$quantsSD
setwd(dirname.Retrospective)
getwd() # should be back at base dir
if (!(plot.dir %in% dir())) {
dir.create(plot.dir)
}
#SSplotComparisons(run_sum, col=rich.colors.short(8) ,subplots=1 , legend=TRUE, endyrvec=end.yr.vec)
#SSplotComparisons(run_sum, col=rich.colors.short(8) ,subplots=2 , legend=TRUE, endyrvec=end.yr.vec)
#SSplotComparisons(run_sum, subplots=1, plot=T, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
#SSplotComparisons(run_sum, subplots=2, plot=T, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
SSplotComparisons(run_sum, plot=T, pdf=F, png=T, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
SSplotComparisons(run_sum, plot=T, pdf=T, png=F, plotdir=plot.dir, legendloc='bottomright', legendlabels=legendlabels, endyrvec=end.yr.vec)
#graphics.off()
outfile <- "SMA_SSF_vals.csv"
SSB_Virgin <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_Virgin"), ]#
SSB_Initial <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_Initial"), ]#
SSB_1950 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_1950"), ]#
SSB_2010 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2010"), ]#
SSB_2011 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2011"), ]#
SSB_2012 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2012"), ]#
SSB_2013 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2013"), ]#
SSB_2014 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2014"), ]#
SSB_2015 <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_2015"), ]#
SSB.retro <- run_sum$quants[as.numeric(row.names(SSB_Virgin)):as.numeric(row.names(SSB_2015)),] # note that this is spawning stock fecundity (1,000s of pups produced) and depends on how spawning units set up in control file
outfile_SD <- "SMA_SSF_SD.csv"
SSB_Virgin_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_Virgin"), ]#
SSB_Initial_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_Initial"), ]#
SSB_1950_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_1950"), ]#
SSB_2010_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2010"), ]#
SSB_2011_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2011"), ]#
SSB_2012_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2012"), ]#
SSB_2013_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2013"), ]#
SSB_2014_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2014"), ]#
SSB_2015_SD <- run_sum$quantsSD[ which(run_sum$quantsSD[,"Label"]=="SSB_2015"), ]#
SSB.retro.SD <- run_sum$quantsSD[as.numeric(row.names(SSB_Virgin_SD)):as.numeric(row.names(SSB_2015_SD)),] # note that this is spawning stock fecundity (1,000s of pups produced) and depends on how spawning units set up in control file
if (csv.yn) {
setwd(plot.dir)
write.table(SSB.retro, file = outfile, append = F, sep = ",", dec = ".", row.names = T, col.names = T)
write.table(SSB.retro.SD, file = outfile_SD, append = F, sep = ",", dec = ".", row.names = T, col.names = T)
setwd(dirname.Retrospective)
}
#----------------------------------------------------------------------------------------------------
# Begin DC Mohns rho Check
#----------------------------------------------------------------------------------------------------
SSB.retro
end.yr.vec
legendlabels
dim(SSB.retro)
SSB.retro.2015 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2015"), ]#
SSB.retro.2014 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2014"), ]#
SSB.retro.2013 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2013"), ]#
SSB.retro.2012 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2012"), ]#
SSB.retro.2011 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2011"), ]#
SSB.retro.2010 <- SSB.retro[ which(SSB.retro[,"Yr"]=="2010"), ]#
retro.i <- 0
mohns.rho.i <- rep(0,length(end.yr.vec[-1]))
SSB.i <- rep(0,length(end.yr.vec[-1]))
SSB.retro.i <- rep(0,length(end.yr.vec[-1]))
for (retro.yr in end.yr.vec[-1]) {
retro.i = retro.i +1
#retro.i <- 1
#retro.yr <- 2015
SSB.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[1]]
SSB.retro.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[retro.i+1]]
mohns.rho.i[retro.i] <- (SSB.retro.i[retro.i] - SSB.i[retro.i])/SSB.i[retro.i]
}
check.mohns.rho.out <- mean(mohns.rho.i)
retro.i
retro.yr
mohns.rho.i
SSB.i
SSB.retro.i
#----------------------------------------------------------------------------------------------------
# End Dc Mohns rho check
#----------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------
# Begin FC Mohns rho function (with DC Edits to match r4ss function "SSmohnsrho" output "AFSC_Hurtado_SSB")
#----------------------------------------------------------------------------------------------------
mohns.rho<-function(x,y=end.yr.vec){
#x<-SSB.retro
#y<-end.yr.vec
retro.i <- 0
mohns.rho.i <- rep(0,length(y[-1]))
SSB.i <- rep(0,length(y[-1]))
SSB.retro.i <- rep(0,length(y[-1]))
mohns.rho.out <- 0
for (retro.yr in end.yr.vec[-1]) {
retro.i = retro.i +1
#retro.i <- 1
#retro.yr <- 2015
SSB.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[1]]
SSB.retro.i[retro.i] <- SSB.retro[ which(SSB.retro[,"Yr"]==retro.yr), ][[retro.i+1]]
mohns.rho.i[retro.i] <- (SSB.retro.i[retro.i] - SSB.i[retro.i])/SSB.i[retro.i]
}
mohns.rho.out <- mean(mohns.rho.i)
retro.i
retro.yr
mohns.rho.i
SSB.i
SSB.retro.i
# print(Outs)
# return(Output)
return(list(SSB.i=c(SSB.i),
SSB.retro.i=c(SSB.retro.i),
mohns.rho.i=c(mohns.rho.i),
mohns.rho.out= mohns.rho.out))
}
#----------------------------------------------------------------------------------------------------
# End FC Mohns rho function
#----------------------------------------------------------------------------------------------------
#SSB.retro
get_mohns.rho <- mohns.rho(SSB.retro)
get_mohns.rho
check.mohns.rho.out
if (csv.yn) {
setwd(plot.dir)
write.table(get_mohns.rho, file="01_FC_get_mohns_rho.txt", row.names=F)
write.table(check.mohns.rho.out, file="02_DC_check_mohns_rho.txt", row.names=F)
setwd(dirname.Retrospective)
}
#----------------------------------------------------------------------------------------------------
# Begin Plots 1-4 : SSF > SSF_MSY and SSF/SSF_MSY > 1; Derived Quantity Estiamtes (mode) and 95% CI
#----------------------------------------------------------------------------------------------------
modeldesc <- legendlabels
nmodels <- dim(run_sum$quants)[2]-2
# Plots 1 and 2
# Stock Synthesis model run SSF/SSF_MSY
# and approximate 95% confidence intervals (95% CI = mode ± 1.96*std)
# using the maximum likelihood estimate (mode) and standard deviation (std)
# Plots 3 and 4
# Stock Synthesis model run SSF
# and approximate 95% confidence intervals (95% CI = mode ± 1.96*std)
# using the maximum likelihood estimate (mode) and standard deviation (std)
SSB_MSY <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_MSY"), ][1]#
SSB_VIRG <- run_sum$quants[ which(run_sum$quants[,"Label"]=="SSB_Virgin"), ][1]
msy_equilib <- SSB_MSY/SSB_VIRG
msy_equilib
#> msy_equilib
#replist1
#350 0.5212652
#
# Plot 1 Relative spawning biomass without error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=3 , legend=TRUE, legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
text( 1965, 1.05, "SSF_MSY", col=grey(0.5))
mtext('Derived Quantity Estimates', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSBwLINEmsy.png"), type='png' )
#graphics.off()
} #if plots statement end
#
# Plot 2 Relative spawning biomass with approximate 95% error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=4 , legend=TRUE, legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
text( 1965, 1.05, "SSF_MSY", col=grey(0.5))
mtext('Derived Quantity Estimates and Approximate 95% Confidence Intervals', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSBwLINEmsy_w_95_error.png"), type='png' )
#graphics.off()
} #if plots statement end
#
# Plot 3 Spawning biomass without error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=1 , legend=TRUE,
legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
abline(h=SSB_MSY/1, lty=3,col=grey(0.5), lwd=2 )
text( 1960, 1.05 *SSB_MSY, paste0("SSF_MSY = ",round(SSB_MSY/1000,2)), col=grey(0.5), pos=4)
text( 1960, 1.20 *SSB_MSY, paste0("Mohn's rho = ",round(get_mohns.rho$mohns.rho.out,3)), col=grey(0.5), pos=4)
mtext('Derived Quantity Estimates', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSB_total_wLINEmsy_Q50.png"), type='png' )
#graphics.off()
} #if plots statement end
#
# Plot 4 Spawning biomass with error
#
# Plots saved to the current plot directory
# Produce.plots=T
# ptype="png"
# ptype="jpeg"
if (Produce.plots) { #if plots statement begin
windows(width=8,height=6,record=T)
par(mfrow=c(1,1),mar=c(4.5,5,1.5,1),oma=c(0,0,1,0), las=1)
SSplotComparisons(run_sum, col=rich.colors.short(nmodels), subplots=2 , legend=TRUE,
legendloc="topright", legendlabels =modeldesc,
shadeForecast=TRUE , endyrvec=end.yr.vec,minbthresh=0, btarg=0 )
abline(h=SSB_MSY/1, lty=3,col=grey(0.5), lwd=2 )
text( 1960, 1.05 *SSB_MSY, paste0("SSF_MSY = ",round(SSB_MSY/1000,2)), col=grey(0.5), pos=4)
text( 1960, 1.20 *SSB_MSY, paste0("Mohn's rho = ",round(get_mohns.rho$mohns.rho.out,3)), col=grey(0.5), pos=4)
mtext('Derived Quantity Estimates and Approximate 95% Confidence Intervals', side=3, outer=T, line=-1)
savePlot( file=paste0(plot.dir,"/SSB_total_wLINEmsy_Q50_w_95_error.png"), type='png' )
#graphics.off()
} #if plots statement end
graphics.off()
#----------------------------------------------------------------------------------------------------
# End Plots 1-4 : SSF > SSF_MSY and SSF/SSF_MSY > 1; Derived Quantity Estiamtes (mode) and 95% CI
#----------------------------------------------------------------------------------------------------
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.