Nothing
SSplotPars <-
function(
dir="c:/path/", repfile="Report.sso",
xlab="Parameter value",ylab="Density",
postfile="posteriors.sso", showpost=TRUE, showprior=TRUE,
showmle=TRUE, showinit=TRUE, showrecdev=TRUE, priorinit=TRUE,
priorfinal=TRUE, showlegend=TRUE, fitrange=FALSE, xaxs="i",
xlim=NULL, ylim=NULL, verbose=TRUE, nrows=3, ncols=3,
new=TRUE, pdf=FALSE, pwidth=7, pheight=7, punits="in",
ptsize=12, returntable=FALSE, strings=c(), exact=FALSE,
newheaders=NULL, burn=0, thin=1,
ctlfile="control.ss_new")
{
################################################################################
#
# SSplotPars
# This function comes with no warranty or guarantee of accuracy
#
# Purpose: To make a multi-figure plot of prior distributions
# from a Stock Synthesis control file
# Written: Ian Taylor, NWFSC/UW. Ian.Taylor-at-noaa.gov
# Returns: Plots of prior distributions used in Stock Synthesis model
# Notes: hosted at http://code.google.com/p/r4ss/
# Required packages: none
#
################################################################################
if(verbose){
cat("Check for new code and report problems at http://code.google.com/p/r4ss/\n")
}
# define subfunction
GetPrior <- function(Ptype,Pmin,Pmax,Pr,Psd,Pval){
# function to calculate prior values is direct translation of code in SSv3
# Ptype used to be a numeric value and is now a character string
# hopefully this function is robust to either option
Ptype2 <- NA
if(is.character(Ptype)){
if(Ptype=="No_prior") Ptype2 <- -1
if(Ptype=="Normal") Ptype2 <- 0
if(Ptype=="Sym_Beta") Ptype2 <- 1
if(Ptype=="Full_Beta") Ptype2 <- 2
if(Ptype=="Log_Norm") Ptype2 <- 3
if(Ptype=="Log_Norm_adjusted") Ptype2 <- 4
}else{
Ptype2 <- Ptype
}
# fix cases where numeric value was read as character (from older SS versions, I think)
if(is.na(Ptype2)) Ptype2 <- as.numeric(Ptype)
if(is.na(Ptype2)) cat("problem with prior type interpretation. Ptype:",Ptype," Ptype2:",Ptype2,"\n")
Pconst <- 0.0001
if(Ptype2==-1){ # no prior
Prior_Like <- rep(0.,length(Pval));
}
if(Ptype2==0){ # normal
Prior_Like <- 0.5*((Pval-Pr)/Psd)^2;
}
if(Ptype2==1){ # symmetric beta value of Psd must be >0.0
mu <- -(Psd*(log( (Pmax+Pmin)*0.5 - Pmin))) - (Psd*(log(0.5)));
Prior_Like <- -(mu+ (Psd*(log(Pval-Pmin+Pconst)))+(Psd*(log(1.-((Pval-Pmin-Pconst)/(Pmax-Pmin))))));
}
if(Ptype2==2){ # CASAL's Beta; check to be sure that Aprior and Bprior are OK before running SS2!
mu <- (Pr-Pmin) / (Pmax-Pmin); # CASAL's v
tau <- (Pr-Pmin)*(Pmax-Pr)/(Psd^2)-1.0;
Bprior <- tau*mu; Aprior <- tau*(1-mu); # CASAL's m and n
if(Bprior<=1.0 | Aprior <=1.0) {cat(" bad Beta prior\n");}
Prior_Like <- (1.0-Bprior)*log(Pconst+Pval-Pmin) + (1.0-Aprior)*log(Pconst+Pmax-Pval)
-(1.0-Bprior)*log(Pconst+Pr-Pmin) - (1.0-Aprior)*log(Pconst+Pmax-Pr);
}
if(Ptype2==3){ # lognormal
Prior_Like <- 0.5*((log(Pval)-Pr)/Psd)^2
}
if(Ptype2==4){ # lognormal with bias correction (from Larry Jacobson)
if(Pmin>0.0){
Prior_Like <- 0.5*((log(Pval)-Pr+0.5*Psd^2)/Psd)^2;
}else{
cat("cannot do prior in log space for parm with min <=0.0\n")
}
}
#### need to work out the following code, including replacing "gammln"
## if(Ptype2==5){ # gamma (from Larry Jacobson)
## warnif <- 1e-15;
## if(Pmin<0.0){
## cat("Lower bound for gamma prior must be >=0. Suggestion ",warnif*10.0,"\n")
## }else{
## # Gamma is defined over [0,+inf) but x=zero causes trouble for some mean/variance combos.
## if(Pval < warnif){
## cat("Pval too close to zero in gamma prior - can not guarantee reliable calculations.\n",
## "Suggest rescaling data (e.g. * 1000)?\n")
## }else{
## scale <- (Psd^2)/Pr; # gamma parameters by method of moments
## shape <- Pr/scale;
## Prior_Like <- -shape*log(scale)-gammln(shape)+(shape-1.0)*log(Pval)-Pval/scale;
## }
## }
## }
return(Prior_Like)
} # end GetPrior
fullpostfile <- paste(dir,postfile,sep="/")
fullrepfile <- paste(dir,repfile,sep="/")
fullctlfile <- paste(dir,ctlfile,sep="/")
postfileinfo <- file.info(fullpostfile)$size
repfileinfo <- file.info(fullrepfile)$size
ctlfileinfo <- file.info(fullctlfile)$size
if(is.na(repfileinfo)) stop("Missing rep file:",fullrepfile)
if(repfileinfo==0) stop("Empty rep file:",fullrepfile)
goodctl <- TRUE
if(is.na(ctlfileinfo)){
cat("Missing control.ss_new file. Assuming recdev limits are -5 & 5.\n")
goodctl <- FALSE
}else{
if(ctlfileinfo==0){
cat("Empty control.ss_new file. Assuming recdev limits are -5 & 5.\n")
goodctl <- FALSE
}
}
if(showpost & is.na(postfileinfo)){
cat("Missing posteriors file: ",postfile,", changing input to 'showpost=FALSE'\n",sep="")
showpost <- FALSE
}
if(showpost & !is.na(postfile) & postfileinfo==0){
cat("Empty posteriors file: ",postfile,", changing input to 'showpost=FALSE'\n",sep="")
showpost <- FALSE
}
## get posteriors
if(showpost & !is.na(postfileinfo) & postfileinfo>0){
test <- readLines(fullpostfile,n=10) # test for presence of file with at least 10 rows
if(length(test)>5){
posts <- read.table(fullpostfile,header=TRUE)
names(posts)[names(posts)=="SR_LN.R0."] <- "SR_LN(R0)"
# remove burn-in and thin the posteriors if requested
posts <- posts[seq(burn+1,nrow(posts),thin), ]
}else{
cat("Posteriors file has too few rows, changing input to 'showpost=FALSE'\n")
showpost <- FALSE
}
}
## get parameter estimates
if(!is.na(repfileinfo) & repfileinfo>0){
replines <- readLines(fullrepfile,n=2000)
parstart <- grep("PARAMETERS",replines)[2]
parend <- grep("DERIVED_QUANTITIES",replines)[2]
nrows2 <- parend - parstart - 3
partable <- read.table(fullrepfile,header=FALSE,nrows=nrows2,skip=parstart,as.is=TRUE,
fill=TRUE,row.names=paste(1:nrows2),col.names=1:60)
partable <- partable[,1:15]
temp <- as.character(partable[1,])
# this command was necessary for some intermediate version of SS (but I forget which)
#temp <- c(temp[1:12],"PR_type_code",temp[13:14])
names(partable) <- temp
partable <- partable[-1,]
rownames(partable) <- 1:nrow(partable)
# check for presence of line which would represent the end of the table
test <- grep("Number_of_active_parameters",partable$Num)
if(length(test)>0) partable <- partable[1:(test-1),]
# clean up contents of the table and make certain columns numeric or character
partable[partable=="_"] <- NA
partable$Active_Cnt <- as.numeric(as.character(partable$Active_Cnt))
partable$Label <- as.character(partable$Label)
for(i in (1:ncol(partable))[!names(partable) %in% c("Label","Status","PR_type")] ){
partable[,i] <- as.numeric(as.character(partable[,i]))
}
}
# subset for only active parameters
allnames <- partable$Label[!is.na(partable$Active_Cnt)]
## get list of subset names if vector "strings" is supplied
if(!is.null(strings)){
goodnames <- NULL
if(exact) goodnames <- allnames[allnames %in% strings]
else for(i in 1:length(strings))
goodnames <- c(goodnames,grep(strings[i],allnames,value=TRUE))
goodnames <- unique(goodnames)
cat("parameters matching input vector 'strings':\n")
print(goodnames)
if(length(goodnames)==0){
cat("No active parameters match input vector 'strings'.\n")
return()
}
}else{
goodnames <- allnames
}
badpars <- grep("Impl_err_",goodnames)
if(length(badpars)>0) goodnames <- goodnames[-badpars]
stds <- partable$Parm_StDev[partable$Label %in% goodnames]
if(showmle & (min(is.na(stds))==1 || min(stds, na.rm=TRUE) <= 0)){
cat("Some parameters have std. dev. values in Report.sso equal to 0.\n",
" Asymptotic uncertainty estimates will not be shown.\n",
" Try re-running the model with the Hessian but no MCMC.\n")
}
# Recruitment Devs
recdevmin <- -5
recdevmin <- 5
recdevlabels <- c("Early_RecrDev_","Early_InitAge_","Main_InitAge_",
"Main_RecrDev_","ForeRecr_","Late_RecrDev_")
if(showrecdev & goodctl){
ctllines <- readLines(fullctlfile)
iline <- grep("#min rec_dev",ctllines)
if(length(iline)==1){
# advanced options stuff
recdevmin <- as.numeric(strsplit(ctllines[iline], " #")[[1]][1])
recdevmax <- as.numeric(strsplit(ctllines[iline+1]," #")[[1]][1])
readrecdev <- as.numeric(strsplit(ctllines[iline+2]," #")[[1]][1])
if(is.na(readrecdev) | readrecdev==1)
cat("This function does not yet display recdev values read from ctl file.\n")
}
}else{
goodnames <- goodnames[!substr(goodnames,1,9) %in% substr(recdevlabels,1,9)]
}
npars <- length(goodnames)
### debugging tools:
## print(goodnames)
## return(partable[partable$Label %in% goodnames,])
# make plot
if(verbose & is.null(xlim)){
if(fitrange){
cat("Plotting range is scaled to fit parameter estimates.\n",
" Change input to 'fitrange=FALSE' to get full parameter range.\n")
}else{
cat("Plotting range is equal to input limits on parameters.\n",
" Range can be scaled to fit estimates by setting input 'fitrange=TRUE'.\n")
}
}
## make plot
if(new & !pdf){
OS <- "Mac"
if(length(grep('linux',version$os)) > 0) OS <- "Linux"
if(length(grep('mingw',version$os)) > 0) OS <- "Windows"
### Note: the following line has been commented out because it was identified
### by Brian Ripley as "against CRAN policies".
#if(exists(".SavedPlots",where=1)) rm(.SavedPlots,pos=1)
if(OS=="Windows") windows(width=pwidth,height=pheight,pointsize=ptsize,record=TRUE)
if(OS=="Linux") X11(width=pwidth,height=pheight,pointsize=ptsize)
if(OS=="Mac") quartz(width=pwidth,height=pheight,pointsize=ptsize)
}
if(pdf){
pdffile <- paste(dir,"/SSplotPars_",format(Sys.time(),'%d-%b-%Y_%H.%M' ),".pdf",sep="")
pdf(file=pdffile,width=pwidth,height=pheight)
if(verbose) cat("PDF file with plots will be: ",pdffile,"\n")
}
if(new) par(mfcol=c(nrows,ncols),mar=c(2,1,2,1),oma=c(2,2,0,0))
if(verbose) cat("Making plots of parameters:\n")
for(ipar in 1:npars){
# grab name and full parameter line
parname <- goodnames[ipar]
if(verbose) cat(" ",parname,"\n")
parline <- partable[partable$Label==parname,]
# grab values
initval <- parline$Init
finalval <- parline$Value
parsd <- parline$Parm_StDev
Pmin <- parline$Min
Pmax <- parline$Max
Ptype <- parline$PR_type
Psd <- parline$Pr_SD
Pr <- parline$Prior
if(substr(parname,1,9) %in% substr(recdevlabels,1,9)){
initval <- 0
Pmin <- recdevmin
Pmax <- recdevmax
Ptype <- 0
Pr <- 0
Psd <- partable$Value[partable$Label=="SR_sigmaR"]
}
x <- seq(Pmin,Pmax,length=5000) # x vector for subsequent calcs
# make empty holders for future information
ymax <- 0 # upper y-limit in plot
xmin <- NULL # lower x-limit in plot
xmax <- NULL # upper x-limit in plot
# get prior
negL_prior <- GetPrior(Ptype=Ptype,Pmin=Pmin,Pmax=Pmax,Pr=Pr,Psd=Psd,Pval=x)
prior <- exp(-1*negL_prior)
# prior likelihood at initial and final values
## priorinit <- exp(-1*GetPrior(Ptype=Ptype,Pmin=Pmin,Pmax=Pmax,Pr=Pr,Psd=Psd,Pval=initval))
## priorfinal <- exp(-1*GetPrior(Ptype=Ptype,Pmin=Pmin,Pmax=Pmax,Pr=Pr,Psd=Psd,Pval=finalval))
if(showprior){
prior <- prior/(sum(prior)*mean(diff(x)))
ymax <- max(ymax,max(prior),na.rm=TRUE) # update ymax
}
# get normal distribution associated with ADMB's estimate
# of the parameter's asymptotic std. dev.
if(showmle){
if(!is.na(parsd) && parsd>0){
mle <- dnorm(x,finalval,parsd)
mlescale <- 1/(sum(mle)*mean(diff(x)))
mle <- mle*mlescale
ymax <- max(ymax,max(mle),na.rm=TRUE) # update ymax
# update x range
xmin <- qnorm(0.001,finalval,parsd)
xmax <- qnorm(0.999,finalval,parsd)
}else{
xmin <- xmax <- finalval
}
}
# get posterior
goodpost <- FALSE
if(showpost){
jpar <- (1:ncol(posts))[names(posts)==parname]
if(length(jpar)==1){
post <- posts[,jpar]
xmin <- min(xmin, quantile(post,0.001)) # update x range
xmax <- max(xmax, quantile(post,0.999)) # update x range
goodpost <- TRUE
}else{
cat("Error! parameter '",parname,"', not found in '",postfile,"'.\n",sep="")
}
}
# get x-range
if(is.null(xlim)){
if(fitrange & ((!is.na(parsd) && parsd!=0) | showpost)){
# if rescaling limits,
# make sure initial value is inside limits
if(showinit){
xmin <- min(initval,xmin)
xmax <- max(initval,xmax)
}
# keep range inside parameter limits
xmin <- max(Pmin,xmin)
xmax <- min(Pmax,xmax)
}else{
# or use parameter limits
xmin <- Pmin
xmax <- Pmax
}
xlim2 <- c(xmin,xmax)
}else{
xlim2 <- xlim
}
# get histogram for posterior based on x-range
if(showpost & goodpost){
jpar <- (1:ncol(posts))[names(posts)==parname]
post <- posts[,jpar]
breakvec <- seq(xmin,xmax,length=50)
if(min(breakvec) > min(post)) breakvec <- c(min(post),breakvec)
if(max(breakvec) < max(post)) breakvec <- c(breakvec,max(post))
posthist <- hist(post,plot=FALSE,breaks=breakvec)
ymax <- max(ymax,max(posthist$density),na.rm=FALSE) # update ymax
}
# make plot
if(is.null(newheaders)) header <- parname else header <- newheaders[ipar]
if(is.null(ylim)) ylim2 <- c(0,1.1*ymax) else ylim2 <- ylim
plot(0,type="n",xlim=xlim2,ylim=ylim2,xaxs=xaxs,yaxs="i",
xlab="",ylab="",main=header,cex.main=1,axes=FALSE)
axis(1)
# axis(2) # don't generally show y-axis values because it's just distracting
# add stuff to plot
colval <- "grey"
if(showpost & goodpost) plot(posthist,add=TRUE,freq=FALSE,col=colval,border=colval)
if(showprior) lines(x,prior,lwd=2,lty=1)
if(showmle){
if(!is.na(parsd) && parsd>0){
lines(x,mle,col="blue",lwd=2)
lines(rep(finalval,2),c(0,dnorm(finalval,finalval,parsd)*mlescale),col="blue")
}else{
abline(v=finalval,col="blue")
}
}
if(showinit){
par(xpd=NA) # stop clipping
points(initval,-0.02*ymax,col="red",pch=17,cex=1.2)
par(xpd=FALSE) # restore original value
}
## if(printlike) mtext(side=3,line=0.2,cex=.8,adj=0,paste("prob@init =",round(priorinit,3)))
## if(printlike) mtext(side=3,line=0.2,cex=.8,adj=1,paste("prob@final =",round(priorfinal,3)))
box()
if(max(par("mfg")[1:2])==1){ # first panel on page
mtext(xlab,side=1,line=0.5,outer=TRUE)
mtext(ylab,side=2,line=0.5,outer=TRUE)
if(showlegend){
showvec <- c(showprior,showmle,showpost,showinit)
legend("topleft",cex=1.2,bty="n",pch=c(NA,NA,15,17)[showvec],
lty=c(1,1,NA,NA)[showvec],lwd=c(2,1,NA,NA)[showvec],
col=c("black","blue","grey","red")[showvec],
pt.cex=c(1,1,2,1)[showvec],
legend=c("prior","max. likelihood","posterior","initial value")[showvec])
} # end legend
} # end first panel stuff
} # end loop over parameters
if(pdf) dev.off() # close PDF file if it was open
if(returntable) return(partable[partable$Label %in% goodnames,])
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.