cgrandin / iscam-gui

0 stars 5 forks source link

Add plots of leading parameter estimates on the ComparisonFigs tab #14

Open cgrandin opened 6 years ago

cgrandin commented 6 years ago

Attached is some standalone code I wrote last week (with example plots). In my code you have to make the tables for each scenario first and, because they all have different names, I am just selecting the first file in the folder (L24), so the code will break if it isn’t the first file. I was just doing a quick and dirty because I wanted to look at the results quickly. In hindsight it would have been easier and preferable to make box plots based on the MCMC output (iscam_mcmc.csv). Then you get a nicer plot and the filename is always the same.

#stopgap code to plot parameters
#Robyn Forrest July 6 2018

graphics.off()
rm(list=ls(all=TRUE))
library(tidyverse)
library(readr)

setwd("C:/GitHub/iscam-gui/Scenarios/")
outDir <- "C:/GitHub/iscam-gui/ExtraPlots/"

Area="5AB"

#1. read scenario folder names
scNames <- list.dirs(path = "C:/GitHub/iscam-gui/Scenarios", full.names = TRUE, recursive = FALSE)
nScenarios <- length(scNames)

#Make a list to put in the tables of parameter estimates
parTables <- list()

#Put parameter estimates into list
for(i in seq_along(scNames)) {
  fileName <- list.files(path=paste0(scNames[i],"/Tables"))
  fileName<-fileName[1]
  print(fileName)
  tmp <- as.data.frame(read.csv(paste0(scNames[i],"/Tables/",fileName), header=T))
  colnames(tmp) <- c("2.5%","50%","97.5%","MPD")
  parTables[[i]] <- tmp
}

#Now get them out and plot them    seq_along(parNames)
parNames <- rownames(parTables[[1]])

for(j in seq_along(parNames)){
  parTable <- as.data.frame(matrix(ncol=4, nrow=nScenarios))
  colnames(tmp) <- c("2.5%","50%","97.5%","MPD")
  Parameter <- parNames[j]

  for(i in seq_along(scNames)){
    parEsts      <- parTables[[i]]
    parTable[i,] <-  parEsts[j,]

  }#end i

  png(paste0(outDir,Area,"_MCMC_",Parameter,".png"))
  plot(1:nScenarios, parTable[,2], pch=19, col=4, las=1, 
       main=paste("MCMC:",Parameter), ylim=c(0,1.1*max(parTable[,3])), 
       xlab="Scenario", ylab=Parameter, cex.lab=1.1, cex=2)
  arrows(1:nScenarios,parTable[,2],1:nScenarios,parTable[,3],col=4, angle=90, lwd=2)
  arrows(1:nScenarios,parTable[,2],1:nScenarios,parTable[,1],col=4, angle=90, lwd=2)
  dev.off()

  png(paste0(outDir,Area,"_MPD_",Parameter,".png"))
  barplot(parTable[,4], names.arg=1:nScenarios, col=4, las=1, 
       main=paste("MPD:",Parameter), ylim=c(0,1.1*max(parTable[,4])), 
       xlab="Scenario", ylab=Parameter, cex.lab=1.1)

  dev.off()
}#end j
cgrandin commented 6 years ago

Add MPD to each box in the boxplot. Use all samples for boxplot.