Nothing
R/SSplotComps.R
In r4ss: R code for Stock Synthesis
Defines functions
SSplotComps
Documented in
SSplotComps
SSplotComps <-
function(replist, subplots=1:12,
kind="LEN", sizemethod=1, aalyear=-1, aalbin=-1, plot=TRUE, print=FALSE,
fleets="all", fleetnames="default", sexes="all",
datonly=FALSE, samplesizeplots=TRUE, compresidplots=TRUE, bub=FALSE,
showsampsize=TRUE, showeffN=TRUE, sampsizeline=FALSE,effNline=FALSE,
minnbubble=8, pntscalar=NULL,
scalebubbles=FALSE,cexZ1=1.5,bublegend=TRUE,blue=rgb(0,0,1,0.7),
pwidth=7, pheight=7, punits="in", ptsize=12, res=300,
plotdir="default", cex.main=1, linepos=1, fitbar=FALSE,
do.sqrt=TRUE, smooth=TRUE, cohortlines=c(),
labels = c("Length (cm)", #1
"Age (yr)", #2
"Year", #3
"Observed sample size", #4
"Effective sample size", #5
"Proportion", #6
"cm", #7
"Frequency", #8
"Weight", #9
"Length", #10
"(mt)", #11
"(numbers x1000)", #12
"Stdev (Age) (yr)", #13
"Andre's conditional AAL plot, "), #14
printmkt=TRUE,printsex=TRUE,
maxrows=6,maxcols=6,maxrows2=2,maxcols2=4,rows=1,cols=1,andrerows=3,
fixdims=TRUE,fixdims2=FALSE,maxneff=5000,verbose=TRUE,
scalebins=FALSE,addMeans=TRUE,...)
{
################################################################################
# SSplotComps
################################################################################
if(!exists("make_multifig")) stop("you are missing the function 'make_mulitifig'")
pngfun <- function(file,caption=NA){
png(filename=file,width=pwidth,height=pheight,
units=punits,res=res,pointsize=ptsize)
plotinfo <- rbind(plotinfo,data.frame(file=file,caption=caption))
return(plotinfo)
}
plotinfo <- NULL
SS_versionNumeric <- replist$SS_versionNumeric
lendbase <- replist$lendbase
sizedbase <- replist$sizedbase
agedbase <- replist$agedbase
condbase <- replist$condbase
ghostagedbase <- replist$ghostagedbase
ghostlendbase <- replist$ghostlendbase
ladbase <- replist$ladbase
wadbase <- replist$wadbase
tagdbase1 <- replist$tagdbase1
tagdbase2 <- replist$tagdbase2
nfleets <- replist$nfleets
nseasons <- replist$nseasons
seasfracs <- replist$seasfracs
FleetNames <- replist$FleetNames
nsexes <- replist$nsexes
accuage <- replist$accuage
Age_tuning <- replist$Age_comp_Eff_N_tuning_check
titles <- NULL
titlemkt <- ""
if(plotdir=="default") plotdir <- replist$inputs$dir
if(fleets[1]=="all"){
fleets <- 1:nfleets
}else{
if(length(intersect(fleets,1:nfleets))!=length(fleets)){
stop("Input 'fleets' should be 'all' or a vector of values between 1 and nfleets.")
}
}
if(sexes[1]=="all") sexes <- 1:2
if(fleetnames[1]=="default") fleetnames <- FleetNames
# a few quantities related to data type and plot number
if(kind=="LEN"){
dbase_kind <- lendbase
kindlab=labels[1]
if(datonly){
filenamestart <- "comp_lendat_"
titledata <- "length comp data, "
}else{
filenamestart <- "comp_lenfit_"
titledata <- "length comps, "
}
}
if(kind=="GSTLEN"){
dbase_kind <- ghostlendbase
kindlab=labels[1]
if(datonly){
filenamestart <- "comp_gstlendat_"
titledata <- "ghost length comp data, "
}else{
filenamestart <- "comp_gstlenfit_"
titledata <- "ghost length comps, "
}
}
if(kind=="SIZE"){
dbase_kind <- sizedbase[sizedbase$method==sizemethod,]
sizeunits <- unique(dbase_kind$units)
if(length(sizeunits)>1)
stop("!error with size units in generalized size comp plots:\n",
" more than one unit value per method.\n")
if(sizeunits %in% c("in","cm"))
kindlab <- paste(labels[10]," (",sizeunits,")",sep="")
if(sizeunits %in% c("lb","kg"))
kindlab <- paste(labels[9]," (",sizeunits,")",sep="")
if(datonly){
filenamestart <- "comp_sizedat_"
titledata <- "size comp data, "
}else{
filenamestart <- "comp_sizefit_"
titledata <- "size comps, "
}
}
if(kind=="AGE"){
dbase_kind <- agedbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_agedat_"
titledata <- "age comp data, "
}else{
filenamestart <- "comp_agefit_"
titledata <- "age comps, "
}
}
if(kind=="cond"){
dbase_kind <- condbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_condAALdat_"
titledata <- "conditional age-at-length data, "
}else{
filenamestart <- "comp_condAALfit_"
titledata <- "conditional age-at-length, "
}
}
if(kind=="GSTAGE"){
dbase_kind <- ghostagedbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_gstagedat_"
titledata <- "ghost age comp data, "
}else{
filenamestart <- "comp_gstagefit_"
titledata <- "ghost age comps, "
}
}
if(kind=="GSTcond"){
dbase_kind <- ghostagedbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_gstCAALdat_"
titledata <- "ghost conditional age-at-length data, "
}else{
filenamestart <- "comp_gstCAALfit_"
titledata <- "ghost conditional age-at-length comps, "
}
}
if(kind=="L@A"){
dbase_kind <- ladbase[ladbase$N!=0,] # remove values with 0 sample size
kindlab=labels[2]
filenamestart <- "comp_LAAfit_"
titledata <- "mean length at age, "
dbase_kind$SD <- dbase_kind$Lbin_lo/dbase_kind$N
}
if(kind=="W@A"){
dbase_kind <- wadbase[wadbase$N!=0,] # remove values with 0 sample size
kindlab=labels[2]
filenamestart <- "comp_WAAfit_"
titledata <- "mean weight at age, "
}
if(!(kind%in%c("LEN","SIZE","AGE","cond","GSTAGE","GSTLEN","L@A","W@A"))) stop("Input 'kind' to SSplotComps is not right.")
# add asterix to indicate super periods and then remove rows labeled "skip"
# would be better to somehow show the range of years, but that seems difficult
# at this point
if(any(dbase_kind$SuprPer=="Sup" & dbase_kind$Used=="skip")){
cat("Note: removing super-period composition values labeled 'skip'\n",
" and designating super-period values with a '*'\n")
dbase_kind <- dbase_kind[dbase_kind$SuprPer=="No" | dbase_kind$Used!="skip",]
dbase_kind$YrSeasName <- paste(dbase_kind$YrSeasName,ifelse(dbase_kind$SuprPer=="Sup","*",""),sep="")
}
ageerr_warning <- TRUE
# loop over fleets
for(f in fleets)
{
# check for the presence of data
if(length(dbase_kind$Obs[dbase_kind$Fleet==f])>0)
{
dbasef <- dbase_kind[dbase_kind$Fleet==f,]
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
#get mean sample quantities to show in conditional age-at-length figures
if(kind %in% c("cond","GSTcond") && f %in% Age_tuning$Fleet){
#### these values not to be trusted in the presence of ghost data:
## HarmEffNage <- Age_tuning$"HarMean(effN)"[Age_tuning$Fleet==f]
## MeanNage <- Age_tuning$"mean(inputN*Adj)"[Age_tuning$Fleet==f]
HarmEffNage <- NULL
MeanNage <- NULL
}else{
HarmEffNage <- NULL
MeanNage <- NULL
}
# loop over genders combinations
for(k in (1:3)[testor])
{
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
if(sex %in% sexes){
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part))
{
dbase <- dbase_k[dbase_k$Part==j,]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
## # starting with SSv3.24a, the Yr.S column is already in the output, otherwise fill it in
## if(!"Yr.S" %in% names(dbase)){
## # add fraction of season to distinguish between samples
## dbase$Yr.S <- dbase$Yr + (0.5/nseasons)*dbase$Seas
## }
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(table(dbase$Yr.S,dbase$Ageerr)>0,1,sum))
if(max_n_ageerr > 1){
if(ageerr_warning){
cat("Note: multiple samples with different ageing error types within fleet/year.\n",
" Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
" Bubble plots may be misleading with overlapping bubbles.\n")
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish between types within a year
dbase$Yr.S <- dbase$Yr.S + dbase$Ageerr/(1000*max_n_ageerr)
dbase$YrSeasName <- paste(dbase$YrSeasName,"a",dbase$Ageerr,sep="")
}
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("flt",f,"sex",k,"mkt",j,sep="")
### subplot 1: multi-panel composition plot
if(1 %in% subplots & kind!="cond"){ # for age or length comps, but not conditional AAL
ptitle <- paste(titledata,title_sexmkt, fleetnames[f],sep="") # total title
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun <- function(ipage,...){
# a function to combine a bunch of repeated commands
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
sampsize=dbase$N,effN=dbase$effN,showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(kind=="GSTAGE"){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(kind=="GSTLEN"){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(kind %in% c("L@A","W@A")){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
ptsSD=dbase$SD,
sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
nlegends=1,legtext=list(dbase$YrSeasName),
bars=bars,linepos=(1-datonly)*linepos,
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
} # end tempfun
if(plot) tempfun(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(dbase$Yr.S))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,"/",filenamestart,
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun(ipage=ipage,...)
dev.off()
}
}
} # end subplot 1
# some things related to the next two bubble plots (single or multi-panel)
if(datonly){
z <- dbase$Obs
if(scalebubbles) z <- dbase$N*dbase$Obs # if requested, scale by sample sizes
col <- 1
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
}else{
z <- dbase$Pearson
col <- blue
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
### subplot 2: single panel bubble plot for numbers at length or age
if(2 %in% subplots & bub & kind!="cond"){
# get growth curves if requested
if(length(cohortlines)>0){
growdat <- replist$endgrowth
growdatF <- growdat[growdat$Gender==1 & growdat$Morph==min(growdat$Morph[growdat$Gender==1]),]
if(nsexes > 1){
growdatM <- growdat[growdat$Gender==2 & growdat$Morph==min(growdat$Morph[growdat$Gender==2]),]
}
}
ptitle <- paste(titletype, title_sexmkt, fleetnames[f],sep="")
ptitle <- paste(ptitle," (max=",round(max(z),digits=2),")",sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun2 <- function(){
bubble3(x=dbase$Yr.S, y=dbase$Bin, z=z, xlab=labels[3],
ylab=kindlab,col=col,cexZ1=cexZ1,
legend=bublegend,
las=1,main=ptitle,cex.main=cex.main,maxsize=pntscalar,
allopen=allopen,minnbubble=minnbubble)
# add lines for growth of individual cohorts if requested
if(length(cohortlines)>0){
for(icohort in 1:length(cohortlines)){
cat(" Adding line for",cohortlines[icohort],"cohort\n")
if(kind=="LEN"){
if(k %in% c(1,2))
lines(growdatF$Age+cohortlines[icohort],growdatF$Len_Mid, col="red") #females
if(nsexes>1 & k %in% c(1,3))
lines(growdatM$Age+cohortlines[icohort],growdatM$Len_Mid, col="blue") #males
}
if(kind=="AGE"){
lines(c(cohortlines[icohort],cohortlines[icohort]+accuage),
c(0,accuage),col="red")
}
}
}
}
if(plot) tempfun2()
if(print){ # set up plotting to png file if required
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
file <- paste(plotdir,"/",filenamestart,filetype,
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun2()
dev.off() # close device if png
}
} # end bubble plot
### subplot 3: multi-panel bubble plots for conditional age-at-length
if(3 %in% subplots & kind=="cond"){
ptitle <- paste(titletype, title_sexmkt, fleetnames[f],sep="")
ptitle <- paste(ptitle," (max=",round(max(z),digits=2),")",sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
# calculate scaling of lines showing effect and input sample size
sampsizeline.old <- sampsizeline
effNline.old <- effNline
if(is.logical(sampsizeline) && sampsizeline){
# scaling when displaying only adjusted input sample size
sampsizeline <- max(dbase$Bin)/max(dbase$N,na.rm=TRUE)
if(!datonly && is.logical(effNline) && effNline){
# scaling when displaying both input and effective
sampsizeline <- effNline <- max(dbase$Bin)/max(dbase$N,dbase$effN,na.rm=TRUE)
cat(" Fleet ",f," ",titlesex,"adj. input & effective N in red & green scaled by ",effNline,"\n",sep="")
}else{
cat(" Fleet ",f," ",titlesex,"adj. input N in red scaled by ",sampsizeline,"\n",sep="")
}
}
# function to make plots
tempfun3 <- function(ipage,...){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Lbin_mid,yr=dbase$Yr.S,size=z,
sampsize=dbase$N,showsampsize=showsampsize,effN=dbase$effN,
showeffN=FALSE,
cexZ1=cexZ1,
bublegend=bublegend,
nlegends=1,legtext=list(dbase$YrSeasName),
bars=FALSE,linepos=0,main=ptitle,cex.main=cex.main,
xlab=labels[2],ylab=labels[1],ymin0=FALSE,maxrows=maxrows2,maxcols=maxcols2,
fixdims=fixdims,allopen=allopen,minnbubble=minnbubble,
ptscol=col[1],ptscol2=col[2],ipage=ipage,scalebins=scalebins,
sampsizeline=sampsizeline,effNline=effNline,
sampsizemean=MeanNage,effNmean=HarmEffNage,...)
}
if(plot) tempfun3(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(dbase$Yr.S))/maxrows2/maxcols2)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,"/",filenamestart,filetype,
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun3(ipage=ipage,...)
dev.off() # close device if png
}
}
sampsizeline <- sampsizeline.old
effNline <- effNline.old
} # end conditional bubble plot
### subplots 4 and 5: multi-panel plot of point and line fit to conditional age-at-length
# and Pearson residuals of A-L key for specific years
if((4 %in% subplots | 5 %in% subplots) & aalyear[1] > 0 & kind=="cond"){
for(y in 1:length(aalyear)){
aalyr <- aalyear[y]
if(length(dbase$Obs[dbase$Yr==aalyr])>0){
if(4 %in% subplots){
### subplot 4: multi-panel plot of fit to conditional age-at-length for specific years
ptitle <- paste(aalyr," age-at-length bin, ",title_sexmkt,fleetnames[f],sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
ydbase <- dbase[dbase$Yr==aalyr,]
lenbinlegend <- paste(ydbase$Lbin_lo,labels[7],sep="")
lenbinlegend[ydbase$Lbin_range>0] <- paste(ydbase$Lbin_lo,"-",ydbase$Lbin_hi,labels[7],sep="")
tempfun4 <- function(ipage,...){ # temporary function to aid repeating the big function call
make_multifig(ptsx=ydbase$Bin,ptsy=ydbase$Obs,yr=ydbase$Lbin_lo,
linesx=ydbase$Bin,linesy=ydbase$Exp,
sampsize=ydbase$N,effN=ydbase$effN,showsampsize=showsampsize,showeffN=showeffN,
nlegends=3,legtext=list(lenbinlegend,"sampsize","effN"),
bars=FALSE,linepos=linepos,main=ptitle,cex.main=cex.main,
xlab=labels[2],ylab=labels[6],maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(plot) tempfun4(ipage=0,...)
if(print){
npages <- ceiling(length(unique(ydbase$Yr.S))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
file <- paste(plotdir,"/",filenamestart,filename_fltsexmkt,
"_",aalyr,"_",pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun4(ipage=ipage,...)
dev.off() # close device if print
}
}
} # end if 4 in subplots
if(5 %in% subplots){
### subplot 5: Pearson residuals for A-L key
z <- ydbase$Pearson
ptitle <- paste(aalyr," Pearson residuals for A-L key, ",title_sexmkt,fleetnames[f],sep="")
ptitle <- paste(ptitle," (max=",round(abs(max(z)),digits=2),")",sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun5 <- function(){
bubble3(x=ydbase$Bin,y=ydbase$Lbin_lo,z=z,xlab=labels[2],
ylab=labels[1],col=blue,las=1,main=ptitle,
cex.main=cex.main,maxsize=pntscalar,
cexZ1=cexZ1,
legend=bublegend,
allopen=FALSE,minnbubble=minnbubble)
}
if(plot) tempfun5()
if(print){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
file <- paste(plotdir,"/",filenamestart,"yearresids_",
filename_fltsexmkt,"_",aalyr,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun5()
dev.off() # close device if print
}
} # end if 5 in subplots
}
}
}
### subplot 6: multi-panel plot of point and line fit to conditional age-at-length
# for specific length bins
if(6 %in% subplots & aalbin[1] > 0){
badbins <- setdiff(aalbin, dbase$Lbin_hi)
goodbins <- intersect(aalbin, dbase$Lbin_hi)
if(length(goodbins)>0){
if(length(badbins)>0){
cat("Error! the following inputs for 'aalbin' do not match the Lbin_hi values for the conditional age-at-length data:",badbins,"\n",
" the following inputs for 'aalbin' are fine:",goodbins,"\n")
}
for(ibin in 1:length(goodbins)){ # loop over good bins
ilenbin <- goodbins[ibin]
abindbase <- dbase[dbase$Lbin_hi==ilenbin,]
if(nrow(abindbase)>0){ # check for data associated with this bin
ptitle <- paste("Age-at-length ",ilenbin,labels[7],", ",title_sexmkt,fleetnames[f],sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun6 <- function(ipage,...){ # temporary function to aid repeating the big function call
make_multifig(ptsx=abindbase$Bin,ptsy=abindbase$Obs,yr=abindbase$Yr.S,linesx=abindbase$Bin,linesy=abindbase$Exp,
sampsize=abindbase$N,effN=abindbase$effN,showsampsize=showsampsize,showeffN=showeffN,
nlegends=3,legtext=list(abindbase$YrSeasName,"sampsize","effN"),
bars=bars,linepos=(1-datonly)*linepos,
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(plot) tempfun6(ipage=0,...)
if(print){
npages <- ceiling(length(unique(abindbase$Yr.S))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,
"_length",ilenbin,labels[7],pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun6(ipage=ipage,...)
dev.off() # close device if print
}
} # end print
} # end if data
} # end loop over length bins
} # end if length(goodbins)>0
} # end if plot requested
### subplot 7: sample size plot
if(7 %in% subplots & samplesizeplots & !datonly & !(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
ptitle <- paste("N-EffN comparison, ",titledata,title_sexmkt,fleetnames[f], sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
lfitfunc <- function(){
if(kind=="cond"){
# trap nonrobust effective n's
# should this only be for conditional age-at-length or all plots?
dbasegood <- dbase[dbase$Obs>=0.0001 & dbase$Exp<0.99 & !is.na(dbase$effN) & dbase$effN<maxneff,]
}else{
dbasegood <- dbase
}
if(nrow(dbasegood)>0){
plot(dbasegood$N,dbasegood$effN,xlab=labels[4],main=ptitle,cex.main=cex.main,
ylim=c(0,1.05*max(dbasegood$effN)),xlim=c(0,1.05*max(dbasegood$N)),
col=blue,pch=19,ylab=labels[5],xaxs="i",yaxs="i")
abline(h=0,col="grey")
abline(0,1,col="black")
# add loess smoother if there's at least 6 points with a range greater than 2
if(smooth & length(unique(dbasegood$N)) > 6 & diff(range(dbasegood$N))>2){
old_warn <- options()$warn # previous warnings setting
options(warn=-1) # turn off loess warnings
psmooth <- loess(dbasegood$effN~dbasegood$N,degree=1)
options(warn=old_warn) #returning to old value
lines(psmooth$x[order(psmooth$x)],psmooth$fit[order(psmooth$x)],lwd=1.2,col="red",lty="dashed")
}
if(addMeans){
abline(v=mean(dbasegood$N),lty="22",col='green3')
abline(h=1/mean(1/dbasegood$effN),lty="22",col='green3')
}
}
}
if(plot) lfitfunc()
if(print){ # set up plotting to png file if required
file <- paste(plotdir,filenamestart,"sampsize_",filename_fltsexmkt,".png",sep="")
caption <- ptitle
plotinfo <- pngfun(file=file, caption=caption)
lfitfunc()
dev.off()
}
} # end subplot 7
### subplot 8: Andre's mean age and std. dev. in conditional AAL
if(8 %in% subplots & kind=="cond"){
ptitle <- paste(labels[14], title_sexmkt, fleetnames[f],sep="")
andrefun <- function(ipage=0){
Lens <-sort(unique(dbase$Lbin_lo))
Yrs <- sort(unique(dbase$Yr.S))
ymax <- 1.1*max(dbase$Bin,na.rm=TRUE)
xmax <- max(condbase$Lbin_hi,na.rm=TRUE)
xmin <- min(condbase$Lbin_lo,na.rm=TRUE)
# do some stuff so that figures that span multiple pages can be output as separate PNG files
npanels <- length(Yrs)
npages <- npanels/andrerows
panelrange <- 1:npanels
if(npages > 1 & ipage!=0) panelrange <- intersect(panelrange, 1:andrerows + andrerows*(ipage-1))
Yrs2 <- Yrs[panelrange]
par(mfrow=c(andrerows,2),mar=c(2,4,1,1),oma=c(3,0,3,0))
for (Yr in Yrs2){
y <- dbase[dbase$Yr.S==Yr,]
Size <- NULL; Size2 <- NULL
Obs <- NULL; Obs2 <- NULL
Pred <- NULL; Pred2 <- NULL
Upp <- NULL; Low <- NULL; Upp2 <- NULL; Low2 <- NULL
for (Ilen in Lens){
z <- y[y$Lbin_lo == Ilen,]
if (length(z[,1]) > 0){
weightsPred <- z$Exp/sum(z$Exp)
weightsObs <- z$Obs/sum(z$Obs)
ObsV <- sum(z$Bin*weightsObs)
ObsV2 <- sum(z$Bin*z$Bin*weightsObs)
PredV <- sum(z$Bin*weightsPred)
PredV2 <- sum(z$Bin*z$Bin*weightsPred)
# Overdispersion on N
# NN <- z$N[1]*0.01 # Andre did this for reasons unknown
NN <- z$N[1]
if (max(z$Obs) > 1.0e-4){
Size <- c(Size,Ilen)
Obs <- c(Obs,ObsV)
Pred <- c(Pred,PredV)
varn <-sqrt(PredV2-PredV*PredV)/sqrt(NN)
Pred2 <- c(Pred2,varn)
varn <-sqrt(max(0,ObsV2-ObsV*ObsV))/sqrt(NN)
Obs2 <- c(Obs2,varn)
Low <- c(Low,ObsV-1.64*varn)
Upp <- c(Upp,ObsV+1.64*varn)
if (NN > 1){
Size2 <- c(Size2,Ilen)
Low2 <- c(Low2,varn*sqrt((NN-1)/qchisq(0.95,NN)))
Upp2 <- c(Upp2,varn*sqrt((NN-1)/qchisq(0.05,NN)))
}
}
}
}
if (length(Obs) > 0){
## next line was replaced with setting at the top,
## for consistency across years
#ymax <- max(Pred,Obs,Upp)*1.1
plot(Size,Obs,type='n',xlab="",ylab="Age",xlim=c(xmin,xmax),ylim=c(0,ymax),yaxs="i")
label <- ifelse(nseasons==1, floor(Yr), Yr)
text(x=par("usr")[1],y=.9*ymax,labels=label,adj=c(-.5,0),font=2,cex=1.2)
if(length(Low)>1) polygon(c(Size,rev(Size)),c(Low,rev(Upp)),col='grey95',border=NA)
if(!datonly) lines(Size,Pred,col=4,lwd=3)
points(Size,Obs,pch=16)
lines(Size,Low,lty=3)
lines(Size,Upp,lty=3)
#title(paste("Year = ",Yr,"; Gender = ",Gender))
if(par("mfg")[1]==1){
title(main=ptitle,xlab=labels[1],outer=TRUE,line=1)
}
box()
ymax2 <- max(Obs2,Pred2)*1.1
plot(Size,Obs2,type='n',xlab=labels[1],ylab=labels[13],xlim=c(xmin,xmax),ylim=c(0,ymax2),yaxs="i")
if(length(Low2)>1) polygon(c(Size2,rev(Size2)),c(Low2,rev(Upp2)),col='grey95',border=NA)
if(!datonly) lines(Size,Pred2,col=4,lwd=3)
points(Size,Obs2,pch=16)
lines(Size2,Low2,lty=3)
lines(Size2,Upp2,lty=3)
if(!datonly & par("mfg")[1]==1){
legend('topleft',legend=c("Observed (with 95% interval)","Expected"),
bty='n',col=c(1,4),pch=c(16,NA),lty=c(NA,1),lwd=3)
}
box()
} # end if data exist
} # end loop over years
} # end andrefun
if(plot) andrefun()
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(dbase$Yr.S))/andrerows)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,"/",filenamestart,"Andre_plots",
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
andrefun(ipage=ipage)
dev.off() # close device if png
} # end loop over pages
} # end test for print to PNG option
} # end subplot 8
} # end loop over partitions
} # end test for whether gender in vector of requested sexes
} # end loop over combined/not-combined genders
} # end if data
} # end loop over fleets
### subplot 9: by fleet aggregating across years
if(9 %in% subplots & kind!="cond") # for age or length comps, but not conditional AAL
{
dbasef <- dbase_kind[dbase_kind$Fleet %in% fleets,]
# check for the presence of data
if(nrow(dbasef)>0)
{
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor])
{
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part))
{
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
dbase <- dbase_k[dbase_k$Part==j,]
if(nrow(dbase)>0){
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("sex",k,"mkt",j,sep="")
ptitle <- paste(titledata,title_sexmkt, "aggregated across time by fleet",sep="") # total title
titles <- c(ptitle,titles) # compiling list of all plot titles
Bins <- sort(unique(dbase$Bin))
nbins <- length(Bins)
df <- data.frame(N=dbase$N,
effN=dbase$effN,
obs=dbase$Obs*dbase$N,
exp=dbase$Exp*dbase$N)
agg <- aggregate(x=df, by=list(bin=dbase$Bin,f=dbase$Fleet), FUN=sum)
agg <- agg[agg$f %in% fleets,]
agg$obs <- agg$obs/agg$N
agg$exp <- agg$exp/agg$N
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for(f in unique(agg$f)){
infleet <- agg$f==f
agg$N[infleet] <- max(agg$N[infleet])
agg$effN[infleet] <- max(agg$effN[infleet])
}
namesvec <- fleetnames[agg$f]
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
# group remaining calculations as a function
tempfun7 <- function(ipage,...){
make_multifig(ptsx=agg$bin,ptsy=agg$obs,yr=agg$f,
linesx=agg$bin,linesy=agg$exp,
sampsize=agg$N,effN=agg$effN,
showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,
legtext=list(namesvec,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims2,ipage=ipage,lwd=2,scalebins=scalebins,...)
}
if(plot) tempfun7(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(agg$f))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,
pagetext,"aggregated across time.png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun7(ipage=ipage,...)
dev.off()
}
} # end print function
}else{
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase$YrSeasName),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined genders
} # end if data
} # end subplot 9
### subplot 10: by fleet aggregating across years within each season
if(10 %in% subplots & kind!="cond" & nseasons>1) # for age or length comps, but not conditional AAL
{
dbasef <- dbase_kind[dbase_kind$Fleet %in% fleets,]
# check for the presence of data
if(nrow(dbasef)>0)
{
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor])
{
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part))
{
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
dbase <- dbase_k[dbase_k$Part==j,]
if(nrow(dbase)>0){
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("sex",k,"mkt",j,sep="")
ptitle <- paste(titledata,title_sexmkt, "\naggregated within season by fleet",sep="") # total title
titles <- c(ptitle,titles) # compiling list of all plot titles
Bins <- sort(unique(dbase$Bin))
nbins <- length(Bins)
df <- data.frame(N=dbase$N,
effN=dbase$effN,
obs=dbase$Obs*dbase$N,
exp=dbase$Exp*dbase$N)
agg <- aggregate(x=df, by=list(bin=dbase$Bin,f=dbase$Fleet,s=dbase$Seas), FUN=sum)
agg <- agg[agg$f %in% fleets,]
if(any(agg$s<=0)){
cat("super-periods may not work correctly in plots of aggregated comps\n")
agg <- agg[agg$s > 0,]
}
agg$obs <- agg$obs/agg$N
agg$exp <- agg$exp/agg$N
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for(f in unique(agg$f)){ # loop over fleets
for(s in unique(agg$s[agg$f==f])){ # loop over seasons within fleet
infleetseas <- agg$f==f & agg$s==s
agg$N[infleetseas] <- max(agg$N[infleetseas])
agg$effN[infleetseas] <- max(agg$effN[infleetseas])
}
}
agg$fseas <- agg$f + seasfracs[agg$s]
namesvec <- paste(fleetnames[agg$f]," s",agg$s,sep="")
# group remaining calculations as a function
tempfun8 <- function(ipage,...){
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
make_multifig(ptsx=agg$bin,ptsy=agg$obs,yr=agg$fseas,
linesx=agg$bin,linesy=agg$exp,
sampsize=agg$N,effN=agg$effN,
showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,
legtext=list(namesvec,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims2,ipage=ipage,lwd=2,scalebins=scalebins,...)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase$YrSeasName),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
if(plot) tempfun8(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(agg$fseas))/maxrows/maxcols)
for(ipage in 1:npages)
{
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,pagetext,
"aggregated within season.png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun8(ipage=ipage,...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined genders
} # end if data
} # end subplot 10
### subplot 11: by fleet aggregating across seasons within a year
if(11 %in% subplots & kind!="cond" & nseasons>1){ # for age or length comps, but not conditional AAL
# loop over fleets
for(f in fleets){
dbasef <- dbase_kind[dbase_kind$Fleet==f,]
# check for the presence of data
if(nrow(dbasef)>0){
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor]){
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part)){
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
dbase <- dbase_k[dbase_k$Part==j,]
if(nrow(dbase)>0){
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("flt",f,"sex",k,"mkt",j,sep="")
Bins <- sort(unique(dbase$Bin))
nbins <- length(Bins)
df <- data.frame(N=dbase$N,
effN=dbase$effN,
obs=dbase$Obs*dbase$N,
exp=dbase$Exp*dbase$N)
agg <- aggregate(x=df, by=list(bin=dbase$Bin,f=dbase$Fleet,y=floor(dbase$Yr.S)), FUN=sum)
agg <- agg[agg$f %in% fleets,]
agg$obs <- agg$obs/agg$N
agg$exp <- agg$exp/agg$N
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for(f in unique(agg$f)){
for(y in unique(agg$y[agg$f==f])){
infleetyr <- agg$f==f & agg$y==y
agg$N[infleetyr] <- max(agg$N[infleetyr])
agg$effN[infleetyr] <- max(agg$effN[infleetyr])
}
}
agg$fy <- agg$f + agg$y/10000
# total title
ptitle <- paste(titledata,title_sexmkt,fleetnames[f],
"\naggregated across seasons within year",sep="")
# group remaining calculations as a function
tempfun9 <- function(ipage,...){
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
make_multifig(ptsx=agg$bin,ptsy=agg$obs,yr=agg$fy,
linesx=agg$bin,linesy=agg$exp,
sampsize=agg$N,effN=agg$effN,
showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,
legtext=list(agg$y,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims2,ipage=ipage,lwd=2,scalebins=scalebins,...)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase$YrSeasName),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
} # end tempfun
if(plot) tempfun9(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(agg$fy))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,pagetext,
"aggregated across seasons within year.png",sep="")
pngfun(file=file, caption=caption)
tempfun9(ipage=ipage,...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined genders
} # end if data
} # end loop over fleets
} # end subplot 11
### subplot 12: bubble plot comparison of length or age residuals
# across fleets within gender/partition
if(12 %in% subplots & kind %in% c("LEN","AGE")){
# check for the presence of data
testor <- length(dbase_kind$Gender[dbase_kind$Gender==1 &
dbase_kind$Pick_gender==0 ])>0
testor[2] <- length(dbase_kind$Gender[dbase_kind$Gender==1 &
dbase_kind$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbase_kind$Gender[dbase_kind$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor]){
if(k==1){dbase_k <- dbase_kind[dbase_kind$Gender==1 & dbase_kind$Pick_gender==0,]}
if(k==2){dbase_k <- dbase_kind[dbase_kind$Gender==1 & dbase_kind$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbase_kind[dbase_kind$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
if(sex %in% sexes){
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part)){
dbase_fleets <- dbase_k[dbase_k$Part==j,]
fleetvec <- intersect(fleets,dbase_fleets$Fleet)
npanels <- length(fleetvec)
xlim <- range(dbase_fleets$Yr.S) # set xlim based on range across all fleets
xaxislab <- sort(unique(floor(dbase_fleets$Yr.S))) # label with all years
# get growth curves if requested
if(length(cohortlines)>0){
growdat <- replist$endgrowth
growdatF <- growdat[growdat$Gender==1 & growdat$Morph==min(growdat$Morph[growdat$Gender==1]),]
if(nsexes > 1){
growdatM <- growdat[growdat$Gender==2 & growdat$Morph==min(growdat$Morph[growdat$Gender==2]),]
}
}
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch"
if(j==1) titlemkt <- "discard"
if(j==2) titlemkt <- "retained"
titlemkt <- ifelse(printmkt,titlemkt,"")
title_sexmkt <- paste(titlesex,titlemkt,sep="")
ptitle <- paste(titletype, title_sexmkt, ", comparing across fleets", sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
filename_sexmkt <- paste("sex",k,"mkt",j,sep="")
tempfun11 <- function(ipage=0){
# a function to wrap up multi-fleet bubble plots
# multi-figure plot with as many rows as fleets, or the maxrows value
par(mfrow=c(min(npanels,maxrows),1), mar=c(0.5,0,0,0),oma=c(4,6,3,1))
# set up some stuff for cases where there are more fleets than panels in one plot
panelrange <- 1:npanels
npages <- ceiling(npanels/maxrows) # how many pages of plots
if(npages > 1 & ipage!=0) # range of which panels to print for each page
panelrange <- intersect(panelrange, 1:maxrows + maxrows*(ipage-1))
# loop over fleets
for(f in fleetvec[panelrange]){
dbase <- dbase_fleets[dbase_fleets$Fleet==f,]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), gender, and partition,
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(table(dbase$Yr.S,dbase$Ageerr)>0,1,sum))
if(max_n_ageerr > 1){
if(ageerr_warning){
cat("Note: multiple samples with different ageing error types within fleet/year.\n",
" Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
" Bubble plots may be misleading with overlapping bubbles.\n")
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish between types within a year
dbase$Yr.S <- dbase$Yr.S + dbase$Ageerr/(1000*max_n_ageerr)
dbase$YrSeasName <- paste(dbase$YrSeasName,"a",dbase$Ageerr,sep="")
}
# determine bubble size and colors
if(datonly){
z <- dbase$Obs
if(scalebubbles) z <- dbase$N*dbase$Obs # if requested, scale by sample sizes
col <- 1
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
}else{
z <- dbase$Pearson
col <- blue
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
# make bubbles for a single fleet
# this section is a modified version of tempfun2 above
ylim <- range(dbase$Bin)
ylim[2] <- ylim[2]+0.2*diff(ylim) # add buffer of 10% at the top for fleet name
bubble3(x=dbase$Yr.S, y=dbase$Bin, z=z, col=col, cexZ1=cexZ1,
legend=bublegend,
las=1,main="",cex.main=cex.main,maxsize=pntscalar,allopen=allopen,
xlim=xlim,ylim=ylim,axis1=FALSE)
#legend('top',title=fleetnames[f],legend=NA,bty='n') # old way with label within each panel
mtext(fleetnames[f],side=2,line=4.5,cex=par()$cex)
# add lines for growth of individual cohorts if requested
if(length(cohortlines)>0){
for(icohort in 1:length(cohortlines)){
cat(" Adding line for",cohortlines[icohort],"cohort\n")
if(kind=="LEN"){
if(k %in% c(1,2))
lines(growdatF$Age+cohortlines[icohort],growdatF$Len_Mid, col="red") #females
if(nsexes>1 & k %in% c(1,3))
lines(growdatM$Age+cohortlines[icohort],growdatM$Len_Mid, col="blue") #males
}
if(kind=="AGE"){
lines(c(cohortlines[icohort],cohortlines[icohort]+accuage),
c(0,accuage),col="red")
}
}
}
if(par()$mfg[1]==par()$mfg[3] | f==tail(fleetvec,1)){
# label all years on x-axis of last panel
axis(1,at=xaxislab)
}else{
# or just tick marks for other panels
axis(1,at=xaxislab,labels=rep("",length(xaxislab)))
}
if(par()$mfg[1]==1)
# add title after making first panel
title(main=ptitle, outer=TRUE, xlab=labels[3], ylab=kindlab)
} # end loop over fleets
} # end function wrapping up a single page of the residual comparison plot
# make plots or write to PNG file
if(length(fleetvec)>0){
if(plot) tempfun11(ipage=0)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(fleetvec)/maxrows)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
caption <- paste(caption,
"<br>Note: bubble sizes are scaled to maximum within each panel.",
"<br>Thus, comparisons across panels should focus on patterns, not bubble sizes.")
file <- paste(plotdir,filenamestart,filename_sexmkt,pagetext,
"_multi-fleet_comparison.png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun11(ipage=ipage)
dev.off()
} # end loop over pages within printing PNG
} # end printing to PNG files
} # end test for non-zero number of fleets
} # end loop over partitions
} # end loop over sexes
} # end loop over gender combinations
# restore default single panel settings
par(mfcol=c(rows,cols),mar=c(5,4,4,2)+.1,oma=rep(0,4))
} # end subplot 12
if(!is.null(plotinfo)) plotinfo$category <- "Comp"
return(invisible(plotinfo))
} # end embedded SSplotComps function
###########################
Try the r4ss package in your browser
Any scripts or data that you put into this service are public.
r4ss documentation built on May 2, 2019, 4:56 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 SSplotComps
Documented in SSplotComps
SSplotComps <-
function(replist, subplots=1:12,
kind="LEN", sizemethod=1, aalyear=-1, aalbin=-1, plot=TRUE, print=FALSE,
fleets="all", fleetnames="default", sexes="all",
datonly=FALSE, samplesizeplots=TRUE, compresidplots=TRUE, bub=FALSE,
showsampsize=TRUE, showeffN=TRUE, sampsizeline=FALSE,effNline=FALSE,
minnbubble=8, pntscalar=NULL,
scalebubbles=FALSE,cexZ1=1.5,bublegend=TRUE,blue=rgb(0,0,1,0.7),
pwidth=7, pheight=7, punits="in", ptsize=12, res=300,
plotdir="default", cex.main=1, linepos=1, fitbar=FALSE,
do.sqrt=TRUE, smooth=TRUE, cohortlines=c(),
labels = c("Length (cm)", #1
"Age (yr)", #2
"Year", #3
"Observed sample size", #4
"Effective sample size", #5
"Proportion", #6
"cm", #7
"Frequency", #8
"Weight", #9
"Length", #10
"(mt)", #11
"(numbers x1000)", #12
"Stdev (Age) (yr)", #13
"Andre's conditional AAL plot, "), #14
printmkt=TRUE,printsex=TRUE,
maxrows=6,maxcols=6,maxrows2=2,maxcols2=4,rows=1,cols=1,andrerows=3,
fixdims=TRUE,fixdims2=FALSE,maxneff=5000,verbose=TRUE,
scalebins=FALSE,addMeans=TRUE,...)
{
################################################################################
# SSplotComps
################################################################################
if(!exists("make_multifig")) stop("you are missing the function 'make_mulitifig'")
pngfun <- function(file,caption=NA){
png(filename=file,width=pwidth,height=pheight,
units=punits,res=res,pointsize=ptsize)
plotinfo <- rbind(plotinfo,data.frame(file=file,caption=caption))
return(plotinfo)
}
plotinfo <- NULL
SS_versionNumeric <- replist$SS_versionNumeric
lendbase <- replist$lendbase
sizedbase <- replist$sizedbase
agedbase <- replist$agedbase
condbase <- replist$condbase
ghostagedbase <- replist$ghostagedbase
ghostlendbase <- replist$ghostlendbase
ladbase <- replist$ladbase
wadbase <- replist$wadbase
tagdbase1 <- replist$tagdbase1
tagdbase2 <- replist$tagdbase2
nfleets <- replist$nfleets
nseasons <- replist$nseasons
seasfracs <- replist$seasfracs
FleetNames <- replist$FleetNames
nsexes <- replist$nsexes
accuage <- replist$accuage
Age_tuning <- replist$Age_comp_Eff_N_tuning_check
titles <- NULL
titlemkt <- ""
if(plotdir=="default") plotdir <- replist$inputs$dir
if(fleets[1]=="all"){
fleets <- 1:nfleets
}else{
if(length(intersect(fleets,1:nfleets))!=length(fleets)){
stop("Input 'fleets' should be 'all' or a vector of values between 1 and nfleets.")
}
}
if(sexes[1]=="all") sexes <- 1:2
if(fleetnames[1]=="default") fleetnames <- FleetNames
# a few quantities related to data type and plot number
if(kind=="LEN"){
dbase_kind <- lendbase
kindlab=labels[1]
if(datonly){
filenamestart <- "comp_lendat_"
titledata <- "length comp data, "
}else{
filenamestart <- "comp_lenfit_"
titledata <- "length comps, "
}
}
if(kind=="GSTLEN"){
dbase_kind <- ghostlendbase
kindlab=labels[1]
if(datonly){
filenamestart <- "comp_gstlendat_"
titledata <- "ghost length comp data, "
}else{
filenamestart <- "comp_gstlenfit_"
titledata <- "ghost length comps, "
}
}
if(kind=="SIZE"){
dbase_kind <- sizedbase[sizedbase$method==sizemethod,]
sizeunits <- unique(dbase_kind$units)
if(length(sizeunits)>1)
stop("!error with size units in generalized size comp plots:\n",
" more than one unit value per method.\n")
if(sizeunits %in% c("in","cm"))
kindlab <- paste(labels[10]," (",sizeunits,")",sep="")
if(sizeunits %in% c("lb","kg"))
kindlab <- paste(labels[9]," (",sizeunits,")",sep="")
if(datonly){
filenamestart <- "comp_sizedat_"
titledata <- "size comp data, "
}else{
filenamestart <- "comp_sizefit_"
titledata <- "size comps, "
}
}
if(kind=="AGE"){
dbase_kind <- agedbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_agedat_"
titledata <- "age comp data, "
}else{
filenamestart <- "comp_agefit_"
titledata <- "age comps, "
}
}
if(kind=="cond"){
dbase_kind <- condbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_condAALdat_"
titledata <- "conditional age-at-length data, "
}else{
filenamestart <- "comp_condAALfit_"
titledata <- "conditional age-at-length, "
}
}
if(kind=="GSTAGE"){
dbase_kind <- ghostagedbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_gstagedat_"
titledata <- "ghost age comp data, "
}else{
filenamestart <- "comp_gstagefit_"
titledata <- "ghost age comps, "
}
}
if(kind=="GSTcond"){
dbase_kind <- ghostagedbase
kindlab=labels[2]
if(datonly){
filenamestart <- "comp_gstCAALdat_"
titledata <- "ghost conditional age-at-length data, "
}else{
filenamestart <- "comp_gstCAALfit_"
titledata <- "ghost conditional age-at-length comps, "
}
}
if(kind=="L@A"){
dbase_kind <- ladbase[ladbase$N!=0,] # remove values with 0 sample size
kindlab=labels[2]
filenamestart <- "comp_LAAfit_"
titledata <- "mean length at age, "
dbase_kind$SD <- dbase_kind$Lbin_lo/dbase_kind$N
}
if(kind=="W@A"){
dbase_kind <- wadbase[wadbase$N!=0,] # remove values with 0 sample size
kindlab=labels[2]
filenamestart <- "comp_WAAfit_"
titledata <- "mean weight at age, "
}
if(!(kind%in%c("LEN","SIZE","AGE","cond","GSTAGE","GSTLEN","L@A","W@A"))) stop("Input 'kind' to SSplotComps is not right.")
# add asterix to indicate super periods and then remove rows labeled "skip"
# would be better to somehow show the range of years, but that seems difficult
# at this point
if(any(dbase_kind$SuprPer=="Sup" & dbase_kind$Used=="skip")){
cat("Note: removing super-period composition values labeled 'skip'\n",
" and designating super-period values with a '*'\n")
dbase_kind <- dbase_kind[dbase_kind$SuprPer=="No" | dbase_kind$Used!="skip",]
dbase_kind$YrSeasName <- paste(dbase_kind$YrSeasName,ifelse(dbase_kind$SuprPer=="Sup","*",""),sep="")
}
ageerr_warning <- TRUE
# loop over fleets
for(f in fleets)
{
# check for the presence of data
if(length(dbase_kind$Obs[dbase_kind$Fleet==f])>0)
{
dbasef <- dbase_kind[dbase_kind$Fleet==f,]
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
#get mean sample quantities to show in conditional age-at-length figures
if(kind %in% c("cond","GSTcond") && f %in% Age_tuning$Fleet){
#### these values not to be trusted in the presence of ghost data:
## HarmEffNage <- Age_tuning$"HarMean(effN)"[Age_tuning$Fleet==f]
## MeanNage <- Age_tuning$"mean(inputN*Adj)"[Age_tuning$Fleet==f]
HarmEffNage <- NULL
MeanNage <- NULL
}else{
HarmEffNage <- NULL
MeanNage <- NULL
}
# loop over genders combinations
for(k in (1:3)[testor])
{
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
if(sex %in% sexes){
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part))
{
dbase <- dbase_k[dbase_k$Part==j,]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
## # starting with SSv3.24a, the Yr.S column is already in the output, otherwise fill it in
## if(!"Yr.S" %in% names(dbase)){
## # add fraction of season to distinguish between samples
## dbase$Yr.S <- dbase$Yr + (0.5/nseasons)*dbase$Seas
## }
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(table(dbase$Yr.S,dbase$Ageerr)>0,1,sum))
if(max_n_ageerr > 1){
if(ageerr_warning){
cat("Note: multiple samples with different ageing error types within fleet/year.\n",
" Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
" Bubble plots may be misleading with overlapping bubbles.\n")
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish between types within a year
dbase$Yr.S <- dbase$Yr.S + dbase$Ageerr/(1000*max_n_ageerr)
dbase$YrSeasName <- paste(dbase$YrSeasName,"a",dbase$Ageerr,sep="")
}
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("flt",f,"sex",k,"mkt",j,sep="")
### subplot 1: multi-panel composition plot
if(1 %in% subplots & kind!="cond"){ # for age or length comps, but not conditional AAL
ptitle <- paste(titledata,title_sexmkt, fleetnames[f],sep="") # total title
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun <- function(ipage,...){
# a function to combine a bunch of repeated commands
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
sampsize=dbase$N,effN=dbase$effN,showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(kind=="GSTAGE"){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(kind=="GSTLEN"){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(kind %in% c("L@A","W@A")){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
ptsSD=dbase$SD,
sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
nlegends=1,legtext=list(dbase$YrSeasName),
bars=bars,linepos=(1-datonly)*linepos,
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
} # end tempfun
if(plot) tempfun(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(dbase$Yr.S))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,"/",filenamestart,
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun(ipage=ipage,...)
dev.off()
}
}
} # end subplot 1
# some things related to the next two bubble plots (single or multi-panel)
if(datonly){
z <- dbase$Obs
if(scalebubbles) z <- dbase$N*dbase$Obs # if requested, scale by sample sizes
col <- 1
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
}else{
z <- dbase$Pearson
col <- blue
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
### subplot 2: single panel bubble plot for numbers at length or age
if(2 %in% subplots & bub & kind!="cond"){
# get growth curves if requested
if(length(cohortlines)>0){
growdat <- replist$endgrowth
growdatF <- growdat[growdat$Gender==1 & growdat$Morph==min(growdat$Morph[growdat$Gender==1]),]
if(nsexes > 1){
growdatM <- growdat[growdat$Gender==2 & growdat$Morph==min(growdat$Morph[growdat$Gender==2]),]
}
}
ptitle <- paste(titletype, title_sexmkt, fleetnames[f],sep="")
ptitle <- paste(ptitle," (max=",round(max(z),digits=2),")",sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun2 <- function(){
bubble3(x=dbase$Yr.S, y=dbase$Bin, z=z, xlab=labels[3],
ylab=kindlab,col=col,cexZ1=cexZ1,
legend=bublegend,
las=1,main=ptitle,cex.main=cex.main,maxsize=pntscalar,
allopen=allopen,minnbubble=minnbubble)
# add lines for growth of individual cohorts if requested
if(length(cohortlines)>0){
for(icohort in 1:length(cohortlines)){
cat(" Adding line for",cohortlines[icohort],"cohort\n")
if(kind=="LEN"){
if(k %in% c(1,2))
lines(growdatF$Age+cohortlines[icohort],growdatF$Len_Mid, col="red") #females
if(nsexes>1 & k %in% c(1,3))
lines(growdatM$Age+cohortlines[icohort],growdatM$Len_Mid, col="blue") #males
}
if(kind=="AGE"){
lines(c(cohortlines[icohort],cohortlines[icohort]+accuage),
c(0,accuage),col="red")
}
}
}
}
if(plot) tempfun2()
if(print){ # set up plotting to png file if required
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
file <- paste(plotdir,"/",filenamestart,filetype,
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun2()
dev.off() # close device if png
}
} # end bubble plot
### subplot 3: multi-panel bubble plots for conditional age-at-length
if(3 %in% subplots & kind=="cond"){
ptitle <- paste(titletype, title_sexmkt, fleetnames[f],sep="")
ptitle <- paste(ptitle," (max=",round(max(z),digits=2),")",sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
# calculate scaling of lines showing effect and input sample size
sampsizeline.old <- sampsizeline
effNline.old <- effNline
if(is.logical(sampsizeline) && sampsizeline){
# scaling when displaying only adjusted input sample size
sampsizeline <- max(dbase$Bin)/max(dbase$N,na.rm=TRUE)
if(!datonly && is.logical(effNline) && effNline){
# scaling when displaying both input and effective
sampsizeline <- effNline <- max(dbase$Bin)/max(dbase$N,dbase$effN,na.rm=TRUE)
cat(" Fleet ",f," ",titlesex,"adj. input & effective N in red & green scaled by ",effNline,"\n",sep="")
}else{
cat(" Fleet ",f," ",titlesex,"adj. input N in red scaled by ",sampsizeline,"\n",sep="")
}
}
# function to make plots
tempfun3 <- function(ipage,...){
make_multifig(ptsx=dbase$Bin,ptsy=dbase$Lbin_mid,yr=dbase$Yr.S,size=z,
sampsize=dbase$N,showsampsize=showsampsize,effN=dbase$effN,
showeffN=FALSE,
cexZ1=cexZ1,
bublegend=bublegend,
nlegends=1,legtext=list(dbase$YrSeasName),
bars=FALSE,linepos=0,main=ptitle,cex.main=cex.main,
xlab=labels[2],ylab=labels[1],ymin0=FALSE,maxrows=maxrows2,maxcols=maxcols2,
fixdims=fixdims,allopen=allopen,minnbubble=minnbubble,
ptscol=col[1],ptscol2=col[2],ipage=ipage,scalebins=scalebins,
sampsizeline=sampsizeline,effNline=effNline,
sampsizemean=MeanNage,effNmean=HarmEffNage,...)
}
if(plot) tempfun3(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(dbase$Yr.S))/maxrows2/maxcols2)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,"/",filenamestart,filetype,
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun3(ipage=ipage,...)
dev.off() # close device if png
}
}
sampsizeline <- sampsizeline.old
effNline <- effNline.old
} # end conditional bubble plot
### subplots 4 and 5: multi-panel plot of point and line fit to conditional age-at-length
# and Pearson residuals of A-L key for specific years
if((4 %in% subplots | 5 %in% subplots) & aalyear[1] > 0 & kind=="cond"){
for(y in 1:length(aalyear)){
aalyr <- aalyear[y]
if(length(dbase$Obs[dbase$Yr==aalyr])>0){
if(4 %in% subplots){
### subplot 4: multi-panel plot of fit to conditional age-at-length for specific years
ptitle <- paste(aalyr," age-at-length bin, ",title_sexmkt,fleetnames[f],sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
ydbase <- dbase[dbase$Yr==aalyr,]
lenbinlegend <- paste(ydbase$Lbin_lo,labels[7],sep="")
lenbinlegend[ydbase$Lbin_range>0] <- paste(ydbase$Lbin_lo,"-",ydbase$Lbin_hi,labels[7],sep="")
tempfun4 <- function(ipage,...){ # temporary function to aid repeating the big function call
make_multifig(ptsx=ydbase$Bin,ptsy=ydbase$Obs,yr=ydbase$Lbin_lo,
linesx=ydbase$Bin,linesy=ydbase$Exp,
sampsize=ydbase$N,effN=ydbase$effN,showsampsize=showsampsize,showeffN=showeffN,
nlegends=3,legtext=list(lenbinlegend,"sampsize","effN"),
bars=FALSE,linepos=linepos,main=ptitle,cex.main=cex.main,
xlab=labels[2],ylab=labels[6],maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(plot) tempfun4(ipage=0,...)
if(print){
npages <- ceiling(length(unique(ydbase$Yr.S))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
file <- paste(plotdir,"/",filenamestart,filename_fltsexmkt,
"_",aalyr,"_",pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun4(ipage=ipage,...)
dev.off() # close device if print
}
}
} # end if 4 in subplots
if(5 %in% subplots){
### subplot 5: Pearson residuals for A-L key
z <- ydbase$Pearson
ptitle <- paste(aalyr," Pearson residuals for A-L key, ",title_sexmkt,fleetnames[f],sep="")
ptitle <- paste(ptitle," (max=",round(abs(max(z)),digits=2),")",sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun5 <- function(){
bubble3(x=ydbase$Bin,y=ydbase$Lbin_lo,z=z,xlab=labels[2],
ylab=labels[1],col=blue,las=1,main=ptitle,
cex.main=cex.main,maxsize=pntscalar,
cexZ1=cexZ1,
legend=bublegend,
allopen=FALSE,minnbubble=minnbubble)
}
if(plot) tempfun5()
if(print){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
file <- paste(plotdir,"/",filenamestart,"yearresids_",
filename_fltsexmkt,"_",aalyr,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun5()
dev.off() # close device if print
}
} # end if 5 in subplots
}
}
}
### subplot 6: multi-panel plot of point and line fit to conditional age-at-length
# for specific length bins
if(6 %in% subplots & aalbin[1] > 0){
badbins <- setdiff(aalbin, dbase$Lbin_hi)
goodbins <- intersect(aalbin, dbase$Lbin_hi)
if(length(goodbins)>0){
if(length(badbins)>0){
cat("Error! the following inputs for 'aalbin' do not match the Lbin_hi values for the conditional age-at-length data:",badbins,"\n",
" the following inputs for 'aalbin' are fine:",goodbins,"\n")
}
for(ibin in 1:length(goodbins)){ # loop over good bins
ilenbin <- goodbins[ibin]
abindbase <- dbase[dbase$Lbin_hi==ilenbin,]
if(nrow(abindbase)>0){ # check for data associated with this bin
ptitle <- paste("Age-at-length ",ilenbin,labels[7],", ",title_sexmkt,fleetnames[f],sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
tempfun6 <- function(ipage,...){ # temporary function to aid repeating the big function call
make_multifig(ptsx=abindbase$Bin,ptsy=abindbase$Obs,yr=abindbase$Yr.S,linesx=abindbase$Bin,linesy=abindbase$Exp,
sampsize=abindbase$N,effN=abindbase$effN,showsampsize=showsampsize,showeffN=showeffN,
nlegends=3,legtext=list(abindbase$YrSeasName,"sampsize","effN"),
bars=bars,linepos=(1-datonly)*linepos,
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims,ipage=ipage,scalebins=scalebins,...)
}
if(plot) tempfun6(ipage=0,...)
if(print){
npages <- ceiling(length(unique(abindbase$Yr.S))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,
"_length",ilenbin,labels[7],pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun6(ipage=ipage,...)
dev.off() # close device if print
}
} # end print
} # end if data
} # end loop over length bins
} # end if length(goodbins)>0
} # end if plot requested
### subplot 7: sample size plot
if(7 %in% subplots & samplesizeplots & !datonly & !(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
ptitle <- paste("N-EffN comparison, ",titledata,title_sexmkt,fleetnames[f], sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
lfitfunc <- function(){
if(kind=="cond"){
# trap nonrobust effective n's
# should this only be for conditional age-at-length or all plots?
dbasegood <- dbase[dbase$Obs>=0.0001 & dbase$Exp<0.99 & !is.na(dbase$effN) & dbase$effN<maxneff,]
}else{
dbasegood <- dbase
}
if(nrow(dbasegood)>0){
plot(dbasegood$N,dbasegood$effN,xlab=labels[4],main=ptitle,cex.main=cex.main,
ylim=c(0,1.05*max(dbasegood$effN)),xlim=c(0,1.05*max(dbasegood$N)),
col=blue,pch=19,ylab=labels[5],xaxs="i",yaxs="i")
abline(h=0,col="grey")
abline(0,1,col="black")
# add loess smoother if there's at least 6 points with a range greater than 2
if(smooth & length(unique(dbasegood$N)) > 6 & diff(range(dbasegood$N))>2){
old_warn <- options()$warn # previous warnings setting
options(warn=-1) # turn off loess warnings
psmooth <- loess(dbasegood$effN~dbasegood$N,degree=1)
options(warn=old_warn) #returning to old value
lines(psmooth$x[order(psmooth$x)],psmooth$fit[order(psmooth$x)],lwd=1.2,col="red",lty="dashed")
}
if(addMeans){
abline(v=mean(dbasegood$N),lty="22",col='green3')
abline(h=1/mean(1/dbasegood$effN),lty="22",col='green3')
}
}
}
if(plot) lfitfunc()
if(print){ # set up plotting to png file if required
file <- paste(plotdir,filenamestart,"sampsize_",filename_fltsexmkt,".png",sep="")
caption <- ptitle
plotinfo <- pngfun(file=file, caption=caption)
lfitfunc()
dev.off()
}
} # end subplot 7
### subplot 8: Andre's mean age and std. dev. in conditional AAL
if(8 %in% subplots & kind=="cond"){
ptitle <- paste(labels[14], title_sexmkt, fleetnames[f],sep="")
andrefun <- function(ipage=0){
Lens <-sort(unique(dbase$Lbin_lo))
Yrs <- sort(unique(dbase$Yr.S))
ymax <- 1.1*max(dbase$Bin,na.rm=TRUE)
xmax <- max(condbase$Lbin_hi,na.rm=TRUE)
xmin <- min(condbase$Lbin_lo,na.rm=TRUE)
# do some stuff so that figures that span multiple pages can be output as separate PNG files
npanels <- length(Yrs)
npages <- npanels/andrerows
panelrange <- 1:npanels
if(npages > 1 & ipage!=0) panelrange <- intersect(panelrange, 1:andrerows + andrerows*(ipage-1))
Yrs2 <- Yrs[panelrange]
par(mfrow=c(andrerows,2),mar=c(2,4,1,1),oma=c(3,0,3,0))
for (Yr in Yrs2){
y <- dbase[dbase$Yr.S==Yr,]
Size <- NULL; Size2 <- NULL
Obs <- NULL; Obs2 <- NULL
Pred <- NULL; Pred2 <- NULL
Upp <- NULL; Low <- NULL; Upp2 <- NULL; Low2 <- NULL
for (Ilen in Lens){
z <- y[y$Lbin_lo == Ilen,]
if (length(z[,1]) > 0){
weightsPred <- z$Exp/sum(z$Exp)
weightsObs <- z$Obs/sum(z$Obs)
ObsV <- sum(z$Bin*weightsObs)
ObsV2 <- sum(z$Bin*z$Bin*weightsObs)
PredV <- sum(z$Bin*weightsPred)
PredV2 <- sum(z$Bin*z$Bin*weightsPred)
# Overdispersion on N
# NN <- z$N[1]*0.01 # Andre did this for reasons unknown
NN <- z$N[1]
if (max(z$Obs) > 1.0e-4){
Size <- c(Size,Ilen)
Obs <- c(Obs,ObsV)
Pred <- c(Pred,PredV)
varn <-sqrt(PredV2-PredV*PredV)/sqrt(NN)
Pred2 <- c(Pred2,varn)
varn <-sqrt(max(0,ObsV2-ObsV*ObsV))/sqrt(NN)
Obs2 <- c(Obs2,varn)
Low <- c(Low,ObsV-1.64*varn)
Upp <- c(Upp,ObsV+1.64*varn)
if (NN > 1){
Size2 <- c(Size2,Ilen)
Low2 <- c(Low2,varn*sqrt((NN-1)/qchisq(0.95,NN)))
Upp2 <- c(Upp2,varn*sqrt((NN-1)/qchisq(0.05,NN)))
}
}
}
}
if (length(Obs) > 0){
## next line was replaced with setting at the top,
## for consistency across years
#ymax <- max(Pred,Obs,Upp)*1.1
plot(Size,Obs,type='n',xlab="",ylab="Age",xlim=c(xmin,xmax),ylim=c(0,ymax),yaxs="i")
label <- ifelse(nseasons==1, floor(Yr), Yr)
text(x=par("usr")[1],y=.9*ymax,labels=label,adj=c(-.5,0),font=2,cex=1.2)
if(length(Low)>1) polygon(c(Size,rev(Size)),c(Low,rev(Upp)),col='grey95',border=NA)
if(!datonly) lines(Size,Pred,col=4,lwd=3)
points(Size,Obs,pch=16)
lines(Size,Low,lty=3)
lines(Size,Upp,lty=3)
#title(paste("Year = ",Yr,"; Gender = ",Gender))
if(par("mfg")[1]==1){
title(main=ptitle,xlab=labels[1],outer=TRUE,line=1)
}
box()
ymax2 <- max(Obs2,Pred2)*1.1
plot(Size,Obs2,type='n',xlab=labels[1],ylab=labels[13],xlim=c(xmin,xmax),ylim=c(0,ymax2),yaxs="i")
if(length(Low2)>1) polygon(c(Size2,rev(Size2)),c(Low2,rev(Upp2)),col='grey95',border=NA)
if(!datonly) lines(Size,Pred2,col=4,lwd=3)
points(Size,Obs2,pch=16)
lines(Size2,Low2,lty=3)
lines(Size2,Upp2,lty=3)
if(!datonly & par("mfg")[1]==1){
legend('topleft',legend=c("Observed (with 95% interval)","Expected"),
bty='n',col=c(1,4),pch=c(16,NA),lty=c(NA,1),lwd=3)
}
box()
} # end if data exist
} # end loop over years
} # end andrefun
if(plot) andrefun()
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(dbase$Yr.S))/andrerows)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,"/",filenamestart,"Andre_plots",
filename_fltsexmkt,pagetext,".png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
andrefun(ipage=ipage)
dev.off() # close device if png
} # end loop over pages
} # end test for print to PNG option
} # end subplot 8
} # end loop over partitions
} # end test for whether gender in vector of requested sexes
} # end loop over combined/not-combined genders
} # end if data
} # end loop over fleets
### subplot 9: by fleet aggregating across years
if(9 %in% subplots & kind!="cond") # for age or length comps, but not conditional AAL
{
dbasef <- dbase_kind[dbase_kind$Fleet %in% fleets,]
# check for the presence of data
if(nrow(dbasef)>0)
{
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor])
{
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part))
{
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
dbase <- dbase_k[dbase_k$Part==j,]
if(nrow(dbase)>0){
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("sex",k,"mkt",j,sep="")
ptitle <- paste(titledata,title_sexmkt, "aggregated across time by fleet",sep="") # total title
titles <- c(ptitle,titles) # compiling list of all plot titles
Bins <- sort(unique(dbase$Bin))
nbins <- length(Bins)
df <- data.frame(N=dbase$N,
effN=dbase$effN,
obs=dbase$Obs*dbase$N,
exp=dbase$Exp*dbase$N)
agg <- aggregate(x=df, by=list(bin=dbase$Bin,f=dbase$Fleet), FUN=sum)
agg <- agg[agg$f %in% fleets,]
agg$obs <- agg$obs/agg$N
agg$exp <- agg$exp/agg$N
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for(f in unique(agg$f)){
infleet <- agg$f==f
agg$N[infleet] <- max(agg$N[infleet])
agg$effN[infleet] <- max(agg$effN[infleet])
}
namesvec <- fleetnames[agg$f]
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
# group remaining calculations as a function
tempfun7 <- function(ipage,...){
make_multifig(ptsx=agg$bin,ptsy=agg$obs,yr=agg$f,
linesx=agg$bin,linesy=agg$exp,
sampsize=agg$N,effN=agg$effN,
showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,
legtext=list(namesvec,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims2,ipage=ipage,lwd=2,scalebins=scalebins,...)
}
if(plot) tempfun7(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(agg$f))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,
pagetext,"aggregated across time.png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun7(ipage=ipage,...)
dev.off()
}
} # end print function
}else{
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase$YrSeasName),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined genders
} # end if data
} # end subplot 9
### subplot 10: by fleet aggregating across years within each season
if(10 %in% subplots & kind!="cond" & nseasons>1) # for age or length comps, but not conditional AAL
{
dbasef <- dbase_kind[dbase_kind$Fleet %in% fleets,]
# check for the presence of data
if(nrow(dbasef)>0)
{
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor])
{
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part))
{
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
dbase <- dbase_k[dbase_k$Part==j,]
if(nrow(dbase)>0){
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("sex",k,"mkt",j,sep="")
ptitle <- paste(titledata,title_sexmkt, "\naggregated within season by fleet",sep="") # total title
titles <- c(ptitle,titles) # compiling list of all plot titles
Bins <- sort(unique(dbase$Bin))
nbins <- length(Bins)
df <- data.frame(N=dbase$N,
effN=dbase$effN,
obs=dbase$Obs*dbase$N,
exp=dbase$Exp*dbase$N)
agg <- aggregate(x=df, by=list(bin=dbase$Bin,f=dbase$Fleet,s=dbase$Seas), FUN=sum)
agg <- agg[agg$f %in% fleets,]
if(any(agg$s<=0)){
cat("super-periods may not work correctly in plots of aggregated comps\n")
agg <- agg[agg$s > 0,]
}
agg$obs <- agg$obs/agg$N
agg$exp <- agg$exp/agg$N
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for(f in unique(agg$f)){ # loop over fleets
for(s in unique(agg$s[agg$f==f])){ # loop over seasons within fleet
infleetseas <- agg$f==f & agg$s==s
agg$N[infleetseas] <- max(agg$N[infleetseas])
agg$effN[infleetseas] <- max(agg$effN[infleetseas])
}
}
agg$fseas <- agg$f + seasfracs[agg$s]
namesvec <- paste(fleetnames[agg$f]," s",agg$s,sep="")
# group remaining calculations as a function
tempfun8 <- function(ipage,...){
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
make_multifig(ptsx=agg$bin,ptsy=agg$obs,yr=agg$fseas,
linesx=agg$bin,linesy=agg$exp,
sampsize=agg$N,effN=agg$effN,
showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,
legtext=list(namesvec,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims2,ipage=ipage,lwd=2,scalebins=scalebins,...)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase$YrSeasName),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
}
if(plot) tempfun8(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(agg$fseas))/maxrows/maxcols)
for(ipage in 1:npages)
{
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,pagetext,
"aggregated within season.png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun8(ipage=ipage,...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined genders
} # end if data
} # end subplot 10
### subplot 11: by fleet aggregating across seasons within a year
if(11 %in% subplots & kind!="cond" & nseasons>1){ # for age or length comps, but not conditional AAL
# loop over fleets
for(f in fleets){
dbasef <- dbase_kind[dbase_kind$Fleet==f,]
# check for the presence of data
if(nrow(dbasef)>0){
testor <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender==0 ])>0
testor[2] <- length(dbasef$Gender[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbasef$Gender[dbasef$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor]){
if(k==1){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender==0,]}
if(k==2){dbase_k <- dbasef[dbasef$Gender==1 & dbasef$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbasef[dbasef$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part)){
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), fleet, gender, and partition
dbase <- dbase_k[dbase_k$Part==j,]
if(nrow(dbase)>0){
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch, "
if(j==1) titlemkt <- "discard, "
if(j==2) titlemkt <- "retained, "
titlemkt <- ifelse(printmkt,titlemkt,"")
# plot bars for data only or if input 'fitbar=TRUE'
if(datonly | fitbar) bars <- TRUE else bars <- FALSE
# aggregating identifiers for plot titles and filenames
title_sexmkt <- paste(titlesex,titlemkt,sep="")
filename_fltsexmkt <- paste("flt",f,"sex",k,"mkt",j,sep="")
Bins <- sort(unique(dbase$Bin))
nbins <- length(Bins)
df <- data.frame(N=dbase$N,
effN=dbase$effN,
obs=dbase$Obs*dbase$N,
exp=dbase$Exp*dbase$N)
agg <- aggregate(x=df, by=list(bin=dbase$Bin,f=dbase$Fleet,y=floor(dbase$Yr.S)), FUN=sum)
agg <- agg[agg$f %in% fleets,]
agg$obs <- agg$obs/agg$N
agg$exp <- agg$exp/agg$N
# note: sample sizes will be different for each bin if tail compression is used
# printed sample sizes in plot will be maximum, which may or may not
# represent sum of sample sizes over all years/ages
for(f in unique(agg$f)){
for(y in unique(agg$y[agg$f==f])){
infleetyr <- agg$f==f & agg$y==y
agg$N[infleetyr] <- max(agg$N[infleetyr])
agg$effN[infleetyr] <- max(agg$effN[infleetyr])
}
}
agg$fy <- agg$f + agg$y/10000
# total title
ptitle <- paste(titledata,title_sexmkt,fleetnames[f],
"\naggregated across seasons within year",sep="")
# group remaining calculations as a function
tempfun9 <- function(ipage,...){
if(!(kind %in% c("GSTAGE","GSTLEN","L@A","W@A"))){
make_multifig(ptsx=agg$bin,ptsy=agg$obs,yr=agg$fy,
linesx=agg$bin,linesy=agg$exp,
sampsize=agg$N,effN=agg$effN,
showsampsize=showsampsize,showeffN=showeffN,
bars=bars,linepos=(1-datonly)*linepos,
nlegends=3,
legtext=list(agg$y,"sampsize","effN"),
main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
fixdims=fixdims2,ipage=ipage,lwd=2,scalebins=scalebins,...)
}
# haven't configured this aggregated plot for other types
## if(kind=="GSTAGE"){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=dbase$effN,showsampsize=FALSE,showeffN=FALSE,
## bars=bars,linepos=(1-datonly)*linepos,
## nlegends=3,legtext=list(dbase$YrSeasName,"sampsize","effN"),
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=labels[6],
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
## if(kind %in% c("L@A","W@A")){
## make_multifig(ptsx=dbase$Bin,ptsy=dbase$Obs,yr=dbase$Yr.S,linesx=dbase$Bin,linesy=dbase$Exp,
## sampsize=dbase$N,effN=0,showsampsize=FALSE,showeffN=FALSE,
## nlegends=1,legtext=list(dbase$YrSeasName),
## bars=bars,linepos=(1-datonly)*linepos,
## main=ptitle,cex.main=cex.main,xlab=kindlab,ylab=ifelse(kind=="W@A",labels[9],labels[1]),
## maxrows=maxrows,maxcols=maxcols,rows=rows,cols=cols,
## fixdims=fixdims,ipage=ipage,...)
## }
} # end tempfun
if(plot) tempfun9(ipage=0,...)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(unique(agg$fy))/maxrows/maxcols)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
file <- paste(plotdir,filenamestart,filename_fltsexmkt,pagetext,
"aggregated across seasons within year.png",sep="")
pngfun(file=file, caption=caption)
tempfun9(ipage=ipage,...)
dev.off()
}
} # end print function
} # end test for presence of observations in this partition
} # end loop over partitions
} # end loop over combined/not-combined genders
} # end if data
} # end loop over fleets
} # end subplot 11
### subplot 12: bubble plot comparison of length or age residuals
# across fleets within gender/partition
if(12 %in% subplots & kind %in% c("LEN","AGE")){
# check for the presence of data
testor <- length(dbase_kind$Gender[dbase_kind$Gender==1 &
dbase_kind$Pick_gender==0 ])>0
testor[2] <- length(dbase_kind$Gender[dbase_kind$Gender==1 &
dbase_kind$Pick_gender %in% c(1,3)])>0
testor[3] <- length(dbase_kind$Gender[dbase_kind$Gender==2])>0
# loop over genders combinations
for(k in (1:3)[testor]){
if(k==1){dbase_k <- dbase_kind[dbase_kind$Gender==1 & dbase_kind$Pick_gender==0,]}
if(k==2){dbase_k <- dbase_kind[dbase_kind$Gender==1 & dbase_kind$Pick_gender %in% c(1,3),]}
if(k==3){dbase_k <- dbase_kind[dbase_kind$Gender==2,]}
sex <- ifelse(k==3, 2, 1)
if(sex %in% sexes){
# loop over partitions (discard, retain, total)
for(j in unique(dbase_k$Part)){
dbase_fleets <- dbase_k[dbase_k$Part==j,]
fleetvec <- intersect(fleets,dbase_fleets$Fleet)
npanels <- length(fleetvec)
xlim <- range(dbase_fleets$Yr.S) # set xlim based on range across all fleets
xaxislab <- sort(unique(floor(dbase_fleets$Yr.S))) # label with all years
# get growth curves if requested
if(length(cohortlines)>0){
growdat <- replist$endgrowth
growdatF <- growdat[growdat$Gender==1 & growdat$Morph==min(growdat$Morph[growdat$Gender==1]),]
if(nsexes > 1){
growdatM <- growdat[growdat$Gender==2 & growdat$Morph==min(growdat$Morph[growdat$Gender==2]),]
}
}
## assemble pieces of plot title
# sex
if(k==1) titlesex <- "sexes combined, "
if(k==2) titlesex <- "female, "
if(k==3) titlesex <- "male, "
titlesex <- ifelse(printsex,titlesex,"")
# market category
if(j==0) titlemkt <- "whole catch"
if(j==1) titlemkt <- "discard"
if(j==2) titlemkt <- "retained"
titlemkt <- ifelse(printmkt,titlemkt,"")
title_sexmkt <- paste(titlesex,titlemkt,sep="")
ptitle <- paste(titletype, title_sexmkt, ", comparing across fleets", sep="")
titles <- c(ptitle,titles) # compiling list of all plot titles
filename_sexmkt <- paste("sex",k,"mkt",j,sep="")
tempfun11 <- function(ipage=0){
# a function to wrap up multi-fleet bubble plots
# multi-figure plot with as many rows as fleets, or the maxrows value
par(mfrow=c(min(npanels,maxrows),1), mar=c(0.5,0,0,0),oma=c(4,6,3,1))
# set up some stuff for cases where there are more fleets than panels in one plot
panelrange <- 1:npanels
npages <- ceiling(npanels/maxrows) # how many pages of plots
if(npages > 1 & ipage!=0) # range of which panels to print for each page
panelrange <- intersect(panelrange, 1:maxrows + maxrows*(ipage-1))
# loop over fleets
for(f in fleetvec[panelrange]){
dbase <- dbase_fleets[dbase_fleets$Fleet==f,]
# dbase is the final data.frame used in the individual plots
# it is subset based on the kind (age, len, age-at-len), gender, and partition,
# check for multiple ageing error types within a year to plot separately
max_n_ageerr <- max(apply(table(dbase$Yr.S,dbase$Ageerr)>0,1,sum))
if(max_n_ageerr > 1){
if(ageerr_warning){
cat("Note: multiple samples with different ageing error types within fleet/year.\n",
" Plots label '2005a3' indicates ageing error type 3 for 2005 sample.\n",
" Bubble plots may be misleading with overlapping bubbles.\n")
ageerr_warning <- FALSE
}
# add 1/1000 of a year for each ageing error type to distinguish between types within a year
dbase$Yr.S <- dbase$Yr.S + dbase$Ageerr/(1000*max_n_ageerr)
dbase$YrSeasName <- paste(dbase$YrSeasName,"a",dbase$Ageerr,sep="")
}
# determine bubble size and colors
if(datonly){
z <- dbase$Obs
if(scalebubbles) z <- dbase$N*dbase$Obs # if requested, scale by sample sizes
col <- 1
titletype <- titledata
filetype <- "bub"
allopen <- TRUE
}else{
z <- dbase$Pearson
col <- blue
titletype <- "Pearson residuals, "
filetype <- "resids"
allopen <- FALSE
}
# make bubbles for a single fleet
# this section is a modified version of tempfun2 above
ylim <- range(dbase$Bin)
ylim[2] <- ylim[2]+0.2*diff(ylim) # add buffer of 10% at the top for fleet name
bubble3(x=dbase$Yr.S, y=dbase$Bin, z=z, col=col, cexZ1=cexZ1,
legend=bublegend,
las=1,main="",cex.main=cex.main,maxsize=pntscalar,allopen=allopen,
xlim=xlim,ylim=ylim,axis1=FALSE)
#legend('top',title=fleetnames[f],legend=NA,bty='n') # old way with label within each panel
mtext(fleetnames[f],side=2,line=4.5,cex=par()$cex)
# add lines for growth of individual cohorts if requested
if(length(cohortlines)>0){
for(icohort in 1:length(cohortlines)){
cat(" Adding line for",cohortlines[icohort],"cohort\n")
if(kind=="LEN"){
if(k %in% c(1,2))
lines(growdatF$Age+cohortlines[icohort],growdatF$Len_Mid, col="red") #females
if(nsexes>1 & k %in% c(1,3))
lines(growdatM$Age+cohortlines[icohort],growdatM$Len_Mid, col="blue") #males
}
if(kind=="AGE"){
lines(c(cohortlines[icohort],cohortlines[icohort]+accuage),
c(0,accuage),col="red")
}
}
}
if(par()$mfg[1]==par()$mfg[3] | f==tail(fleetvec,1)){
# label all years on x-axis of last panel
axis(1,at=xaxislab)
}else{
# or just tick marks for other panels
axis(1,at=xaxislab,labels=rep("",length(xaxislab)))
}
if(par()$mfg[1]==1)
# add title after making first panel
title(main=ptitle, outer=TRUE, xlab=labels[3], ylab=kindlab)
} # end loop over fleets
} # end function wrapping up a single page of the residual comparison plot
# make plots or write to PNG file
if(length(fleetvec)>0){
if(plot) tempfun11(ipage=0)
if(print){ # set up plotting to png file if required
npages <- ceiling(length(fleetvec)/maxrows)
for(ipage in 1:npages){
caption <- ptitle
pagetext <- ""
if(npages>1){
pagetext <- paste("_page",ipage,sep="")
caption <- paste(caption, " (plot ",ipage," of ",npages,")",sep="")
}
if(length(grep("Pearson",caption))>0){
caption <- paste(caption,
"<br> \nClosed bubbles are positive residuals and open bubbles are negative residuals.")
}
caption <- paste(caption,
"<br>Note: bubble sizes are scaled to maximum within each panel.",
"<br>Thus, comparisons across panels should focus on patterns, not bubble sizes.")
file <- paste(plotdir,filenamestart,filename_sexmkt,pagetext,
"_multi-fleet_comparison.png",sep="")
plotinfo <- pngfun(file=file, caption=caption)
tempfun11(ipage=ipage)
dev.off()
} # end loop over pages within printing PNG
} # end printing to PNG files
} # end test for non-zero number of fleets
} # end loop over partitions
} # end loop over sexes
} # end loop over gender combinations
# restore default single panel settings
par(mfcol=c(rows,cols),mar=c(5,4,4,2)+.1,oma=rep(0,4))
} # end subplot 12
if(!is.null(plotinfo)) plotinfo$category <- "Comp"
return(invisible(plotinfo))
} # end embedded SSplotComps function
###########################
Try the r4ss package in your browser
Any scripts or data that you put into this service are public.
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.