R/addSSsummarize.R

In r4ss: R code for Stock Synthesis

addSSsummarize <- function(origModels,newModels) {
########################################################################
# This function adds the specified components to a model summary object
# If there are no original models, it will create an object that can be used with r4ss model comparisons
#
# Written by Allan Hicks, 1/12/2011, allan.hicks@noaa.gov
#
# origModels is the list of original models that are being added to
#   it is the result of the SSsumarize function
#   if origModels is NULL, then a new list will be returned that can be used in r4ss model comparisons
# 
# newModels is a list of models to be added to the origModels
#   each new model is an element of the list, and is a list itself
#   each new model is a list with possible components named:
#       npars:          the number of parameters in the model
#       maxgrad:        the maximum gradient component (if used)
#       nsexes:         the number of sexes
#       likelihoods:    likelihoods from the model. A data.frame with the 2nd column as names, which matches on names from origModels
#                           SS uses the following names
#                            TOTAL
#                            Equil_catch
#                            Survey
#                            Length_comp
#                            Age_comp
#                            Recruitment
#                            recast_Recruitment
#                            Parm_priors
#                            Parm_softbounds
#                            Parm_devs
#                            Crash_Pen
#                            Size_at_age
#       likelambdas     NOT IMPLEMENTED
#       #pars#:         NOT IMPLEMENTED YET FOR DIFFICULTY IN MATCHING PARAMETERS
#       #parsSD#:       NOT IMPLEMENTED YET FOR DIFFICULTY IN MATCHING PARAMETERS
#       #parsphases#:   NOT IMPLEMENTED YET FOR DIFFICULTY IN MATCHING PARAMETERS
#       SpawnBio:       Spawning biomass matrix 
#                           1st column is year
#                           2nd column is spawning biomass in same units as original models (SS reports female spawning biomass)
#                           3rd column is the standard deviation of estimated spawning biomass
#                           4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
#                           5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
#       Bratio:         Depletion matrix
#                           1st column is year
#                           2nd column is depletion
#                           3rd column is the standard deviation of depletion (optional)
#                           4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
#                           5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
#       SPRratio:         SPR ratio matrix
#                           1st column is year
#                           2nd column is depletion
#                           3rd column is the standard deviation (optional)
#                           4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
#                           5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
#       recruits:       Recruitment matrix
#                           1st column is year
#                           2nd column is recruitment as in original models (SS reports age-0 recruits)
#                           3rd column is the standard deviation (optional)
#                           4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
#                           5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
#       recdevs:        Recruitment deviate matrix
#                           1st column is year
#                           2nd column is deviate (matched with original models)
#       growth:         #NOT IMPLEMENTED
#       indices:        Matrix of fits to indices
#                           1st column is year
#                           2nd column is observed index (data)
#                           3rd column is expected index (prediction)
#                           4th column is catchability coefficient (q)
#                           5th column is standard error of index (total used in fitting)
#                           6th column is a likelihood for this point, or enter any value to make sure it plots, or enter NA not to plot the estimate
#                           
#       InitAgeYrs      #NOT IMPLEMENTED

    if(!is.list(newModels[[1]])) stop("The new models needs to be a list of lists. A list of new models, each with a list of quantities to enter.")

    cat("Adding",length(newModels),"new models to",origModels$n,"original models\n")
    
    if(is.null(origModels)) {
        stop("Not implemented yet for no original models\n")
        #origModels <- list()
    }
    
    addColumn <- function(origDF,newDF,matchOrig,matchNew,nOrig,valNew=1) {
        #adds a column for the new model to the dataframe and matches using the original column labeled with matchOrig and the new 
        #nOrig is the number of original models
        if(is.null(newDF)) {
            newDF <- data.frame(rep(NA,nrow(origDF)),origDF[,matchOrig])
            valNew <- 1  
            matchNew <- 2
        }
        x.add <- newDF[!(newDF[,matchNew] %in% origDF[,matchOrig]),]
        if(nrow(x.add)>0) { #add additional rows for missing names
            origDFnrow <- nrow(origDF)
            for(i in 1:nrow(x.add)) {origDF <- rbind(origDF,NA)} #add neessary number of rows of NA
            origDF[,matchOrig] <- c(origDF[1:origDFnrow,matchOrig],as.character(x.add[,matchNew]))
        }
        x.ind <- match(origDF[,matchOrig],newDF[,matchNew])   #should have all rows of newDF and origDF
        out <- cbind(origDF[,1:nOrig],newDF[x.ind,valNew],origDF[,(nOrig+1):ncol(origDF)])
        colnames(out) <- c(1:(nOrig+1),colnames(origDF)[(nOrig+1):ncol(origDF)])
        out
    }
    
    addNewModel.fn <- function(x,models) {
        #adds a single model to the origModels
        #models is the set of models formatted as in SSsummarize
        #written as a function to use lapply, avoid loops and make it easier when origModels=NULL
        n <- models$n
        models$n <- n+1
        if(!is.null(x$npars))   { models$npars[n+1] <- x$npars }else{ models$npars[n+1] <- NA }
        if(!is.null(x$listnames))   { models$npars[n+1] <- x$listnames }else{ models$npars[n+1] <- NA }
        if(!is.null(x$keyvec))   { models$npars[n+1] <- x$keyvec }else{ models$npars[n+1] <- NA }
        if(!is.null(x$maxgrad)) { models$maxgrad[n+1] <- x$maxgrad }else{ models$maxgrad[n+1] <- NA }
        if(!is.null(x$nsexes))  { models$nsexes[n+1] <- x$nsexes }else{ cat("nsexes must be defined! Inserting a value of 3.\n"); models$nsexes[n+1] <- 3}
        
        models$pars <- addColumn(models$pars,NULL,"Label",2,n)   ###NEED TO IMPLEMENT
        models$parsSD <- addColumn(models$parsSD,NULL,"Label",2,n)   ###NEED TO IMPLEMENT
        models$parphases <- addColumn(models$parphases,NULL,"Label",2,n)   ###NEED TO IMPLEMENT
        models$quants <- addColumn(models$quants,NULL,"Label",2,n)   ###NEED TO IMPLEMENT
        models$quantsSD <- addColumn(models$quantsSD,NULL,"Label",2,n)   ###NEED TO IMPLEMENT        
        models$likelihoods <- addColumn(models$likelihoods,x$likelihoods,"Label",2,n)
        models$likelambdas <- addColumn(models$likelambdas,NULL,"Label",2,n)   ###NEED TO IMPLEMENT
        #add Spawning Biomasses, note that it is only the column that changes
        if(!is.null(x$SpawnBio)) print("adding SpawnBio")
        if(!is.null(x$SpawnBio)){if(ncol(x$SpawnBio)<5) {stop("There must be at least 5 columns in SpawnBio: Yr, SB, SD, Lower, Upper\n")}}
        models$SpawnBio <- addColumn(models$SpawnBio,x$SpawnBio,"Yr",1,n,2)
        models$SpawnBio$Yr <- as.numeric(models$SpawnBio$Yr)
        models$SpawnBio <- models$SpawnBio[order(models$SpawnBio[,"Yr"]),]
        models$SpawnBioSD <- addColumn(models$SpawnBioSD,x$SpawnBio,"Yr",1,n,3)
        models$SpawnBioSD$Yr <- as.numeric(models$SpawnBioSD$Yr)
        models$SpawnBioSD <- models$SpawnBioSD[order(models$SpawnBioSD[,"Yr"]),]
        models$SpawnBioLower <- addColumn(models$SpawnBioLower,x$SpawnBio,"Yr",1,n,4)
        models$SpawnBioLower$Yr <- as.numeric(models$SpawnBioLower$Yr)
        models$SpawnBioLower <- models$SpawnBioLower[order(models$SpawnBioLower[,"Yr"]),]
        models$SpawnBioUpper <- addColumn(models$SpawnBioUpper,x$SpawnBio,"Yr",1,n,5)
        models$SpawnBioUpper$Yr <- as.numeric(models$SpawnBioUpper$Yr)
        models$SpawnBioUpper <- models$SpawnBioUpper[order(models$SpawnBioUpper[,"Yr"]),]
        #add Bratio (depletion), note that it is only the column that changes
        if(!is.null(x$Bratio)) {
            if(ncol(x$Bratio)<5) {stop("There must be at least 5 columns in Bratio: Yr, SB, SD, Lower, Upper\n")}
            print("adding Bratio")
        }
        models$Bratio <- addColumn(models$Bratio,x$Bratio,"Yr",1,n,2)
        models$Bratio$Yr <- as.numeric(models$Bratio$Yr)
        models$Bratio <- models$Bratio[order(models$Bratio[,"Yr"]),]
        models$BratioSD <- addColumn(models$BratioSD,x$Bratio,"Yr",1,n,3)
        models$BratioSD$Yr <- as.numeric(models$BratioSD$Yr)
        models$BratioSD <- models$BratioSD[order(models$BratioSD[,"Yr"]),]
        models$BratioLower <- addColumn(models$BratioLower,x$Bratio,"Yr",1,n,4)
        models$BratioLower$Yr <- as.numeric(models$BratioLower$Yr)
        models$BratioLower <- models$BratioLower[order(models$BratioLower[,"Yr"]),]
        models$BratioUpper <- addColumn(models$BratioUpper,x$Bratio,"Yr",1,n,5)
        models$BratioUpper$Yr <- as.numeric(models$BratioUpper$Yr)
        models$BratioUpper <- models$BratioUpper[order(models$BratioUpper[,"Yr"]),]
        #add recruits        
        models$recruits <- addColumn(models$recruits,x$recruits,"Yr",1,n,2)
        models$recruitsSD <- addColumn(models$recruitsSD,x$recruits,"Yr",1,n,3)
        models$recruitsLower <- addColumn(models$recruitsLower,x$recruits,"Yr",1,n,4)
        models$recruitsUpper <- addColumn(models$recruitsUpper,x$recruits,"Yr",1,n,5)
        models$recdevs <- addColumn(models$recdevs,x$recdevs,"Yr",1,n,2)
        #add SPR
        models$SPRratio <- addColumn(models$SPRratio,x$SPRratio,"Yr",1,n,2)
        models$SPRratioSD <- addColumn(models$SPRratioSD,x$SPRratio,"Yr",1,n,3)
        models$SPRratioLower <- addColumn(models$SPRratioLower,x$SPRratio,"Yr",1,n,4)
        models$SPRratioUpper <- addColumn(models$SPRratioUpper,x$SPRratio,"Yr",1,n,5)


        models$growth <- cbind(models$growth,NA)  #NEED TO IMPLEMENT
        
        #Indices
        if(!is.null(x$indices)) {
            newIndices <- as.data.frame(cbind(NA,x$indices[,1],NA,NA,x$indices[,2],x$indices[,3],x$indices[,4],NA,x$indices[,5],NA,x$indices[,6],NA,NA,NA,NA,n+1,n+1))
            names(newIndices) <- names(models$indices)
            models$indices <- rbind(models$indices,newIndices)
        }
        
        #models$InitAgeYrs <- cbind(models$InitAgeYrs,NA)  #NEED TO IMPLEMENT
        
        return(models)
    }
    
    for(ii in 1:length(newModels)) {  #lapply doesn't work because it returns a separate list for each newModel
        origModels <- addNewModel.fn(newModels[[ii]],models=origModels)
    }
        
    return(origModels)
}

if(F) {

    likes <- c(4242,5555,9933,3377)
    likes <- data.frame(likes=likes,label=c("Survey","Other","TOTAL","Other2"))
    SB <- data.frame(c(1971,1972,1973),c(3333,5555,7777),c(100,200,300),c(1111,3333,5555),c(5555,7777,9999))
    BR <- data.frame(c(1971,1972,1973),c(0.3333,0.5555,0.7777),c(0.100,.200,.300),c(.1111,.3333,.5555),c(.5555,.7777,.9999))
    newModel1 <- list(npars=42,likelihoods=likes,SpawnBio=SB,Bratio=BR)  #data.frame(1,"NatM_p_1_Fem_GP_1"))
    likes <- c(4233,5544,9988,3344)
    likes <- data.frame(likes=likes,label=c("Recruitment","Other","TOTAL","Other2"))
    newModel2 <- list(npars=42,likelihoods=likes)
    newModels <- list(newModel1,newModel2)
    tmp <- addSSsummarize(mysummary,newModels)

        x.ind <- match(models$likelihoods[,"Label"],names(x$likelihoods))
        models$likelihoods <- cbind(models$likelihoods[,1:n],x$likelihoods[x.ind],models$likelihoods[,n+1])
        colnames(models$likelihoods) <- c(1:(n+1),"Label")
        x.add <- x$likelihoods[!(names(x$likelihoods) %in% models$likelihoods[,"Label"])]
        if(length(x.add)>0) { #add additional rows for missing names
            tmp <- as.data.frame(matrix(NA,nrow=length(x.add),ncol=ncol(models$likelihoods),dimnames=list(NULL,colnames(models$likelihoods))))
            tmp[,n+1] <- x.add
            tmp[,"Label"] <- names(x.add)
            models$likelihoods <- rbind(models$likelihoods,tmp)
        }

    addNewModel.fn <- function(x,models) {
        #adds a single model to the origModels
        #models is the set of models formatted as in SSsummarize
        #written as a function to use lapply, avoid loops and make it easier when origModels=NULL
        n <- models$n
        models$n <- n+1
        if(!is.null(x$npars))   { models$npars[n+1] <- x$npars }else{ models$npars[n+1] <- NA }
        if(!is.null(x$maxgrad)) { models$maxgrad[n+1] <- x$maxgrad }else{ models$maxgrad[n+1] <- NA }
        if(!is.null(x$nsexes))  { models$nsexes[n+1] <- x$nsexes }else{ models$nsexes[n+1] <- NA }
        
        models$likelihoods <- addColumn(models$likelihoods,x$likelihoods,"Label",2,n)
        
        return(models)
    }

}