R/exploit.plt.R

In LobsterScience/bio.lobster: Lobster Stock Assessment Tools for DFO Maritimes

#' @export
exploit.plt <- function(output, years, plt=c('f','m'),graphic="R", type="mcmc",path=''){

	
# Fit plots
if(missing(years))years<-output$data$iyr
if(graphic=='pdf')pdf(file.path(path,"exploit.pdf"), width = 6.5, height = 3.5)
if(graphic=="R")windows(6.5,3.5)

par(mar = c(1, 4, 1, 4), omi = c(0.2, 0.1, 0.2, 0.1))

if("m"%in%plt)ymax=1
if(!"m"%in%plt)ymax=ceiling(max(apply(output$sims.list$mu, 2, quantile, 0.975))*10)/10

if("f"%in%plt){
	if(type=="mcmc"){
		plot(years, output$median$mu, type = 'o', pch = 16, ylab = "", las = 1, mgp = c(0.5, 0.5, 0), xlab = "", tcl = -0.3, ylim = c(0, ymax),xlim=c(min(years),max(years)+1), asp = 'xy')
		#abline(h=0.25, lty=2)
		lines(years, apply(output$sims.list$mu, 2, quantile, 0.025), lty = 2)
		lines(years, apply(output$sims.list$mu, 2, quantile, 0.975), lty = 2)
	}
	else if(type=="mle"){
		plot(years, output$MLE$mu, type = 'o', pch = 16, ylab = "", las = 1, mgp = c(0.5, 0.5, 0), xlab = "", tcl = -0.3, ylim = c(0, ymax),xlim=c(min(years),max(years)+1), asp = 'xy')
	}
	mtext("Exploitation Rate", 2, 2.5, cex = 1.25)
}
if("mR"%in%plt){
	if(type=="mcmc"){
		points(years, exp(-output$median$mR), type = 'o', pch = 18,col='green')
		lines(years, apply(exp(-output$sims.list$mR), 2, quantile, 0.025), lty = 2,col='green')
		lines(years, apply(exp(-output$sims.list$mR), 2, quantile, 0.975), lty = 2,col='green')
	}
	else if(type=="mle"){
		points(years, exp(-output$MLE$mR), type = 'o', pch = 18,col='green')
	}
}
if("m"%in%plt){
	if(type=="mcmc"){
		points(years, exp(-output$median$m), type = 'o', pch = 18,col='grey50')
		lines(years, apply(exp(-output$sims.list$m), 2, quantile, 0.025), lty = 2,col='grey50')
		lines(years, apply(exp(-output$sims.list$m), 2, quantile, 0.975), lty = 2,col='grey50')
	}
	else if(type=="mle"){
		points(years, exp(-output$MLE$m), type = 'o', pch = 18,col='grey50')
	}
	mtext("Natural suvival fraction", 4, 2.5, cex = 1.25)
	axis(4, tcl = -0.3, las=1)
}
if(!"m"%in%plt){
	axis(4, at=1-exp(-pretty(c(0,-log(1-ymax)))),lab=pretty(c(0,-log(1-ymax))),tcl = -0.3, las=1, mgp = c(0.5, 0.5, 0))
	mtext("F", 4, 2.5, cex = 1.25)

}

if(graphic=='pdf')dev.off()
}