R/SSplotMCMC_ExtraSelex.R

In r4ss: R code for Stock Synthesis

SSplotMCMC_ExtraSelex <- function(post, add=FALSE, nsexes=1,shift=0,
                                  fleetname="default",col="blue"){
  # post is a data.frame containing either derived_posteriors.sso or a good subset of it
  #      it can be an element of the list created by the the SSgetMCMC function
  # add will add to existing plot
  # nsexes should match model values but is only used in the title
  # shift will shift the x values by for overplotting on an existing plot
  # fleetname will make title better
  # col will color points and lines
  
  cols <- grep("Selex_std",names(post))
  if(length(cols)==0){
    stop("no columns in posteriors include text 'Selex_std'")
  }else{
    sel <- post[,cols]
    names <- names(post)[cols]
    splitnames <- strsplit(names,"_")
    namesDF <- as.data.frame(matrix(unlist(strsplit(names,"_")),ncol=6,byrow=T))
    i       <- as.numeric(as.character(namesDF$V3))[1]
    m       <- as.character(namesDF$V4)[1]
    agelen  <- as.character(namesDF$V5)[1]
    bin     <- sort(unique(as.numeric(as.character(namesDF$V6))))+shift
    quants <- apply(sel,2,quantile, probs=c(0.025,0.5,0.975))
    
    xlab <- "Age (years)"
    if(fleetname=="default") fleetname <- paste("Fleet",i)

    if(m=="Fem" & nsexes==1) sextitle3 <- ""
    if(m=="Fem" & nsexes==2) sextitle3 <- "females"
    if(m=="Mal") sextitle3 <- "males"
    main <- paste("Uncertainty in selectivity for",fleetname,sextitle3)
    no0 <- as.numeric(quants[3,]-quants[1,])!=0
    upper=quants[3,]
    lower=quants[1,]
    if(!add) matplot(bin,t(quants),lty=c(3,1,3),lwd=c(1,3,1),xlab=xlab,ylim=c(0,1),main=main,
                     ylab="Selectivity",type="n",xlim=c(0,max(bin)))
    
    lines(bin,quants[2,],lty=1,col=col,lwd=1,type="o")
    arrows(x0=bin[no0], y0=quants[1,no0], x1=bin[no0], y1=quants[3,no0],
           length=0.01, angle=90, code=3, col=col)
    abline(h=0,col="grey")
  }
}