In wStockhausen/rTCSAM2013: Functions to run and plot output from TCSAM2013 models.

knitr::opts_chunk$set(echo=FALSE, error=FALSE, message=FALSE, warning=FALSE,
                      dpi=300,fig.width=6.5,out.width=6.5)

    require(pander);
    require(reshape2);
    require(rTCSAM2013);
    panderOptions("table.caption.prefix","Table ")
    panderOptions("table.split.table",Inf)
    paths<-params$paths;
    obj<-params$obj;
    out.dir<-params$out.dir;
    tblno<-1;

Input model cases

    if (is.null(obj)){
        if (!is.character(paths)){
            cat("'paths' parameter is not a character vector!\n");
            stop();
        }
        obj<-list();
        for (case in names(paths)){
            cat(case,": '",paths[case],"'\n",sep='');
            res<-getResLst(inp.dir=paths[case]);
            if (is.null(res)) {
                #remove path from paths
                paths<-paths[names(paths)!=case];
            } else {
                #add res to obj list
                obj[[case]]<-res;
            }
            rm(res);
        }
    } else {
        obj<-convertToListOfResults(obj);
    }
    cases<-names(obj);
    ncases<-length(cases);

    if (!is.null(paths)){ 
        t<-as.data.frame(list(case=cases, path=getSubPaths(paths,last=3)));
    } else {
        t<-as.data.frame(list(case=cases));
    }
    row.names(t)<-NULL;
    pander(t,caption=paste0(tblno,". Model cases for comparison."));
    tblno<-tblno+1;

Objective function components

    dfr<-getMDFR.ObjFunComponents(obj);
    if (ncases==1){
        #single model case
        cols<-c("category","description","weight","likelihood","objFun");
        fn<-file.path(out.dir,"ModelComparison.OFCs.csv");
        write.csv(dfr[,cols],file=fn);
    } else {
        dfro<-reshape2::dcast(dfr,formula="category+description~case",fun.aggregate=sum,value.var="objFun");
        fn<-file.path(out.dir,"ModelComparison.OFCs.csv");
        write.csv(dfro,row.names=FALSE,file=fn);
        rm(fn,dfro);
    }
    panderOptions('knitr.auto.asis',FALSE);##need to do this for loops below
    cats<-as.character(unique(dfr$category))
    for (i in 1:length(cats)){
        dfrp<-dfr[dfr$category==cats[i],];
        row.names(dfrp)<-NULL;
        if (ncases==1){
            #single model case
            cols<-c("description","weight","likelihood","objFun");
            pandoc.table(dfrp[,cols],digits=2,round=2,keep.trailing.zeros=TRUE,justify='lccc',
                   caption=paste0(tblno,". Objective function ",cats[i]," components."))
            tblno<-tblno+1;
        } else {
            #multiple models
            ##objective function values
            dfro<-reshape2::dcast(dfrp,formula="description~case",fun.aggregate=sum,value.var="objFun");
            just<-'l';for (j in 1:ncases) just<-paste0(just,'c');
            pandoc.table(dfro,digits=2,round=2,keep.trailing.zeros=TRUE,justify=just,
                         caption=paste0(tblno,". Objective function ",cats[i]," components."))
            tblno<-tblno+1;
        }
    }
    if (ncases>1){
        ##multiple models, differences tables
        dfrd<-reshape2::dcast(dfr,formula="category+description~case",fun.aggregate=sum,value.var="diff")[,-3];
        names(dfrd)[1+2:ncases]<-paste0(cases[2:ncases],"-",cases[1]);
        write.csv(dfrd,row.names=FALSE,file=file.path(out.dir,"ModelComparison.OFCs.diffs.csv"));
        for (i in 1:length(cats)){
            dfrp<-dfr[dfr$category==cats[i],];
            row.names(dfrp)<-NULL;
            dfrd<-reshape2::dcast(dfrp,formula="description~case",fun.aggregate=sum,value.var="diff")[,-2];
            names(dfrd)[2:ncases]<-paste0(cases[2:ncases],"-",cases[1]);
            just<-'l';for (j in 2:ncases) just<-paste0(just,'c');
            pandoc.table(dfrd,digits=2,round=2,keep.trailing.zeros=TRUE,justify=just,
                         caption=paste0(tblno,". Objective function ",cats[i]," component differences."))
            tblno<-tblno+1;
        }
    }
    panderOptions('knitr.auto.asis',TRUE)

Parameter estimates

    dfr<-getMDFR.ParamsPlusStdDevs(obj);
    dfr$param<-gsub("\\[..\\]","",dfr$param,perl=TRUE);#remove [XX] from parameter vectors
    dfr$index<-as.character(dfr$index);
    dfr$phase<-as.character(dfr$phase);
    if (ncases==1){
        cols<-c("category","description","param","index","phase",
                "min","max","init","value","stdv","check");
        fn<-file.path(out.dir,"ModelComparison.Params.csv");
        write.csv(dfr[,cols],row.names=FALSE,file=fn);
    } else {
        mdfr<-melt(dfr,id.vars=c("category","process","case","description","param","index"),measure.vars=c("value","stdv"));
        ddfr<-dcast(mdfr,category+process+description+param+index~variable+case,value.var='value');
        names(ddfr)<-gsub("_","\n",names(ddfr),fixed=TRUE);
        fn<-file.path(out.dir,"ModelComparison.Params.csv");
        write.csv(ddfr,row.names=FALSE,file=fn);
    }
    panderOptions('knitr.auto.asis',FALSE);##need to do this for loops below
    cats<-as.character(unique(dfr$category))
    for (i in 1:length(cats)){
        dfrp<-dfr[dfr$category==cats[i],];
        prcs<-as.character(unique(dfrp$process));
        for (j in 1:length(prcs)){
            dfrpp<-dfrp[dfrp$process==prcs[j],];
            row.names(dfrpp)<-NULL;
            if (ncases==1){
                cols<-c("description","param","index","phase",
                        "min","max","init","value","stdv","check");
                pandoc.table(dfrpp[,cols],keep.trailing.zeros=TRUE,justify='llcccccccc',
                       caption=paste0(tblno,". Parameter estimates for ",cats[i]," ",prcs[j]," ."));
                tblno<-tblno+1;
            } else {
                mdfr<-melt(dfrpp,id.vars=c("case","description","param","index"),measure.vars=c("value","stdv"));
                ddfr<-dcast(mdfr,description+param+index~variable+case,value.var='value');
                just<-'llc';
                for (k in 1:ncases) just<-paste0(just,'rr');
                names(ddfr)<-gsub("_","\n",names(ddfr),fixed=TRUE);
                pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
                       caption=paste0(tblno,". Parameter estimates for ",cats[i]," ",prcs[j]," ."));
                tblno<-tblno+1;
            }
        }
    }
    panderOptions('knitr.auto.asis',TRUE);

Mature biomass-at-mating

    dfr<-getMDFR.PopQuantities(obj,type='MB_yx');
    ddfr<-dcast(dfr,y~x+case,value.var='val');
    names(ddfr)[1]<-'year';
    fn<-file.path(out.dir,paste0("ModelComparison.MatureBiomassAtMating.csv"));
    write.csv(ddfr,row.names=FALSE,file=fn);
    panderOptions('knitr.auto.asis',FALSE);##need to do this for loops below
    #males
    ddfr<-dcast(dfr[dfr$x=='male',],y~case,value.var='val');
    names(ddfr)[1]<-'year';
    just<-'l';
    for (i in 1:ncases) just<-paste0(just,'r');
    pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
           caption=paste0(tblno,". Estimated MMB-at-mating time (1000's t)."));
    tblno<-tblno+1;
    #females
    ddfr<-dcast(dfr[dfr$x=='female',],y~case,value.var='val');
    names(ddfr)[1]<-'year';
    just<-'l';
    for (i in 1:ncases) just<-paste0(just,'r');
    pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
           caption=paste0(tblno,". Estimated MFB-at-mating time (1000's t)."));
    tblno<-tblno+1;
    panderOptions('knitr.auto.asis',TRUE);

Recruitment

    dfr<-getMDFR.PopQuantities(obj,type='R_y');
    panderOptions('knitr.auto.asis',FALSE);##need to do this for loops below
    ddfr<-dcast(dfr,y~case,value.var='val');
    names(ddfr)[1]<-'year';
    fn<-file.path(out.dir,paste0("ModelComparison.Recruitment.csv"));
    write.csv(ddfr,row.names=FALSE,file=fn);
    just<-'l';
    for (i in 1:ncases) just<-paste0(just,'r');
    pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
           caption=paste0(tblno,". Estimated recruitment (millions)."));
    tblno<-tblno+1;
    panderOptions('knitr.auto.asis',TRUE);

Mature survey biomass

    dfr<-getMDFR.SurveyQuantities(obj,type='MB_yx');
    dfrp<-dcast(dfr,x+y~category+case,value.var='val')
    dfrp<-dfrp[,c(1:3,2+ncases+(1:ncases))];
    names(dfrp)[2:3]<-c('year','observed');
    names(dfrp)[3+(1:ncases)]<-gsub("predicted_","",names(dfrp)[3+(1:ncases)],fixed=TRUE);
    fn<-file.path(out.dir,paste0("ModelComparison.MatureSurveyBiomass.csv"));
    write.csv(dfrp,row.names=FALSE,file=fn);
    panderOptions('knitr.auto.asis',FALSE);##need to do this for loops below
    just<-'lr'; for (i in 1:ncases) just<-paste0(just,'r');
    #males
    ddfr<-dfrp[dfrp$x=='male',2:ncol(dfrp)]; rownames(ddfr)<-NULL;
    pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
           caption=paste0(tblno,". Observed and estimated mature male survey biomass (1000's t)."));
    tblno<-tblno+1;
    #females
    ddfr<-dfrp[dfrp$x=='female',2:ncol(dfrp)]; rownames(ddfr)<-NULL;
    pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
           caption=paste0(tblno,". Observed and estimated mature female survey biomass (1000's t)."));
    tblno<-tblno+1;
    panderOptions('knitr.auto.asis',TRUE);

Retained catch

    dfr<-getMDFR.FisheryQuantities(obj,type='bio.retm');
    dfrp<-dcast(dfr,x+y~category+case,value.var='val')
    dfrp<-dfrp[,c(1:3,2+ncases+(1:ncases))];
    names(dfrp)[2:3]<-c('year','observed');
    names(dfrp)[3+(1:ncases)]<-gsub("predicted_","",names(dfrp)[3+(1:ncases)],fixed=TRUE);
    panderOptions('knitr.auto.asis',FALSE);##need to do this for loops below
    just<-'lr'; for (i in 1:ncases) just<-paste0(just,'r');
    #males
    ddfr<-dfrp[dfrp$x=='male',2:ncol(dfrp)]; rownames(ddfr)<-NULL;
    fn<-file.path(out.dir,paste0("ModelComparison.RetainedCatch.csv"));
    write.csv(ddfr,row.names=FALSE,file=fn);
    pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
           caption=paste0(tblno,". Observed and estimated retained catch (1000's t)."));
    tblno<-tblno+1;
    panderOptions('knitr.auto.asis',TRUE);

Total catch mortality

    panderOptions('knitr.auto.asis',FALSE);##need to do this for loops below
    dfr<-getMDFR.FisheryQuantities(obj,type='bio.totm');
    for (f in c("TCF","SCF","RKF","GTF")){
        dfrp<-dcast(dfr[dfr$fishery==f,],x+y~category+case,value.var='val')
        dfrp<-dfrp[,c(1:3,2+ncases+(1:ncases))];
        names(dfrp)[2:3]<-c('year','observed');
        names(dfrp)[3+(1:ncases)]<-gsub("predicted_","",names(dfrp)[3+(1:ncases)],fixed=TRUE);
        just<-'lr'; for (i in 1:ncases) just<-paste0(just,'r');
        for (x in c('male','female')){ 
            ddfr<-dfrp[dfrp$x==x,2:ncol(dfrp)]; rownames(ddfr)<-NULL;
            fn<-file.path(out.dir,paste("ModelComparison.TotalCatch",f,x,"csv",sep='.'));
            cat(fn,"\n")
            write.csv(ddfr,file=fn);
            pandoc.table(ddfr,keep.trailing.zeros=FALSE,justify=just,
                   caption=paste0(tblno,". Observed and estimated total ",x," catch mortality biomass (1000's t) in ",f,"."));
            tblno<-tblno+1;
        }#x
    }#f
    panderOptions('knitr.auto.asis',TRUE);

    dfr.ft<-getMDFR.FisheryQuantities(obj,type="effSS.tot");
    dfr.ft$type<-'total catch';
    dfr.fr<-getMDFR.FisheryQuantities(obj,type="effSS.ret");
    dfr.fr$type<-'retained catch';
    dfr.f<-rbind(dfr.fr,dfr.ft);
    write.csv(dfr.f,file.path(out.dir,"ModelComparison.FisherySampleSizes.csv"));
    dfr.s <-getMDFR.SurveyQuantities(obj,type="effSS_y");
    write.csv(dfr.s,file.path(out.dir,"ModelComparison.SurveySampleSizes.csv"));