R/SSplotRetroDevs.R

In r4ss/r4ss: R Code for Stock Synthesis

# this function replaced by SSplotRetroRecruits

SSplotRetroDevs <- function(
  retroSummary,
  endyrvec,
  cohorts,
  ylim = c(-3, 3),
  uncertainty = FALSE,
  labels = c(
    "Recruitment deviation relative to recent estimate",
    "Recruitment deviation",
    "Age"
  ),
  main = "Retrospective analysis of recruitment deviations",
  mcmcVec = FALSE,
  relative = FALSE,
  labelyears = TRUE,
  legend = FALSE,
  leg.ncols = 4
) {
  addpoly <- function(
    yrvec,
    lower,
    upper,
    shadecol = rgb(0, 0, 0, .1),
    col = 1
  ) {
    # add shaded uncertainty intervals behind line
    # modified from SSplotComparisons in r4ss package
    polygon(
      x = c(yrvec, rev(yrvec)),
      y = c(lower, rev(upper)),
      border = NA,
      col = shadecol
    )
    lines(yrvec, lower, lty = 3, col = col)
    lines(yrvec, upper, lty = 3, col = col)
  }

  n <- retroSummary[["n"]]
  recdevs <- retroSummary[["recdevs"]]
  recdevsLower <- retroSummary[["recdevsLower"]]
  recdevsUpper <- retroSummary[["recdevsUpper"]]
  lowerCI <- retroSummary[["lowerCI"]]
  upperCI <- retroSummary[["upperCI"]]

  colvec <- rich.colors.short(length(cohorts), alpha = .7)
  shadecolvec <- rich.colors.short(length(cohorts), alpha = .1)
  colvec <- rainbow(length(cohorts), alpha = .7)
  shadecolvec <- rainbow(length(cohorts), alpha = .1)
  colvec.txt <- colvec
  # make text darker
  for (i in seq_along(colvec)) {
    tmp <- col2rgb(colvec[i]) / 255
    colvec.txt[i] <- rgb(tmp[1] / 2, tmp[2] / 2, tmp[3] / 2, alpha = .7)
  }
  print(cbind(colvec, colvec.txt))

  ylab <- ifelse(relative, labels[1], labels[2])
  maxage <- max(endyrvec) - min(cohorts)
  xlim <- c(0, maxage)
  if (labelyears) {
    xlim <- xlim + c(-.8, .8)
  } # expand x-axis to make room for labels

  plot(
    0,
    type = "n",
    xlim = xlim,
    ylim = ylim,
    xlab = labels[3],
    ylab = ylab,
    main = main,
    axes = FALSE
  )
  axis(1, at = 0:maxage)
  axis(2, at = ylim[1]:ylim[2], las = 1)
  abline(h = 0, col = "grey")
  box()

  if (legend) {
    ylim <- ylim + c(0, 1)
  }

  if (length(mcmcVec) == 1) {
    mcmcVec <- rep(mcmcVec, n)
  }
  if (any(mcmcVec)) {
    mcmc <- retroSummary[["mcmc"]]
  }
  for (imodel in (1:n)[mcmcVec]) {
    tmp <- unique(c(
      grep("_RecrDev_", names(mcmc[[imodel]])),
      grep("_InitAge_", names(mcmc[[imodel]])),
      grep("ForeRecr_", names(mcmc[[imodel]]))
    ))
    if (length(tmp) > 0) {
      # there are some mcmc values to use
      mcmc.tmp <- mcmc[[imodel]][, tmp] # subset of columns from MCMC for this model
      mcmclabs <- names(mcmc.tmp)
      lower <- apply(mcmc.tmp, 2, quantile, prob = lowerCI) # hard-wired probability
      med <- apply(mcmc.tmp, 2, quantile, prob = 0.5) # hard-wired probability
      upper <- apply(mcmc.tmp, 2, quantile, prob = upperCI) # hard-wired probability
      recdevs[, imodel] <- med[match(recdevs[["Label"]], mcmclabs)]
      recdevsLower[, imodel] <- lower[match(recdevsLower[["Label"]], mcmclabs)]
      recdevsUpper[, imodel] <- upper[match(recdevsUpper[["Label"]], mcmclabs)]
    }
  }

  for (iy in seq_along(cohorts)) {
    y <- cohorts[iy]
    cohortdevs <- recdevs[recdevs[["Yr"]] == y, 1:n]
    cohortdevsLower <- recdevsLower[recdevsLower[["Yr"]] == y, 1:n]
    cohortdevsUpper <- recdevsUpper[recdevsUpper[["Yr"]] == y, 1:n]
    # combine rows where the parameter labels may differ
    if (nrow(cohortdevs) > 1) {
      cohortdevs2 <- rep(NA, n)
      cohortdevsLower2 <- rep(NA, n)
      cohortdevsUpper2 <- rep(NA, n)
      for (icol in 1:n) {
        cohortdevs2[icol] <- cohortdevs[!is.na(cohortdevs[, icol]), icol]
        cohortdevsLower2[icol] <- cohortdevsLower[
          !is.na(cohortdevsLower[, icol]),
          icol
        ]
        cohortdevsUpper2[icol] <- cohortdevsUpper[
          !is.na(cohortdevsUpper[, icol]),
          icol
        ]
      }
      cohortdevs <- cohortdevs2
      cohortdevsLower <- cohortdevsLower2
      cohortdevsUpper <- cohortdevsUpper2
    }
    cohortdevs <- as.numeric(cohortdevs)
    cohortdevsLower <- as.numeric(cohortdevsLower)
    cohortdevsUpper <- as.numeric(cohortdevsUpper)

    goodmodels <- (1:n)[endyrvec - y >= 0]
    if (relative) {
      # relative to final estimate
      if (uncertainty) {
        addpoly(
          yrvec = endyrvec[goodmodels] - y,
          lower = cohortdevsLower[goodmodels] - cohortdevs[max(goodmodels)],
          upper = cohortdevsUpper[goodmodels] - cohortdevs[max(goodmodels)],
          shadecol = shadecolvec[iy],
          col = colvec[iy]
        )
      }
      lines(
        endyrvec[goodmodels] - y,
        cohortdevs[goodmodels] - cohortdevs[max(goodmodels)],
        type = "o",
        col = colvec[iy],
        lwd = 3,
        pch = 16
      )
      if (labelyears) {
        text(
          x = (endyrvec[goodmodels] - y)[1] - 0.5,
          y = (cohortdevs[goodmodels] - cohortdevs[max(goodmodels)])[1],
          labels = y,
          col = colvec.txt[iy],
          cex = .7
        )
      }
    } else {
      # true value
      if (uncertainty) {
        addpoly(
          yrvec = endyrvec[goodmodels] - y,
          lower = cohortdevsLower[goodmodels],
          upper = cohortdevsUpper[goodmodels],
          shadecol = shadecolvec[iy],
          col = colvec[iy]
        )
      }
      lines(
        endyrvec[goodmodels] - y,
        cohortdevs[goodmodels],
        type = "o",
        col = colvec[iy],
        lwd = 3,
        pch = 16
      )
      if (labelyears) {
        text(
          x = rev(endyrvec[goodmodels] - y)[1] + 0.5,
          y = rev(cohortdevs[goodmodels])[1],
          labels = y,
          col = colvec.txt[iy],
          cex = .7
        )
      }
    }
  }
  if (legend) {
    legend(
      "topright",
      lwd = 3,
      lty = 1,
      pch = 16,
      col = colvec,
      legend = cohorts,
      title = "Cohort birth year",
      ncol = leg.ncols,
      bg = rgb(1, 1, 1, .3),
      box.col = NA
    )
  }
}

## if(FALSE){
##   #### example use
##   # source this file
##   source('c:/SS/hake/Hake_2012/retro/retro_script.R')

##   # move to directory one level above existing model run
##   setwd('C:/ss/hake/Hake_2013/runs/')

##   # run the function above
##   retro(olddir='2013hake_12',years=0:-10)
##   # read in output
##   retroModels <- SSgetoutput(dirvec=paste('retrospectives/retro',-10:0,sep=''))
##   # summarize output
##   retroSummary <- SSsummarize(retroModels)

##   # set the ending year of each model in the set
##   endyrvec <- retroModels[[1]][["endyr"]]-10:0
##   # make comparison plot
##   pdf('retrospectives/retrospective_comparison_plots.pdf')
##   SSplotComparisons(retroSummary,endyrvec=endyrvec,new=FALSE)
##   dev.off()

##   # make Ianelli-style plot of recdev retrospectives
##   pdf('retrospectives/retrospective_dev_plots.pdf',width=7,height=10)
##   par(mfrow=c(2,1))
##   # first scaled relative to most recent estimate
##   SSplotRetroDevs(retroSummary, endyrvec=endyrvec, cohorts=1999:2012, relative=TRUE, legend=FALSE)
##   # second without scaling
##   SSplotRetroDevs(retroSummary, endyrvec=endyrvec, cohorts=1999:2012, relative=FALSE, legend=FALSE)
##   dev.off()

## }