R/plot_final_year_trophic_data_with_ci.R

In StrathE2E2: End-to-End Marine Food Web Model

#
# plot_final_year_trophic_data_with_ci.R
#
#' Plots of credible value distributions mean trophic level and omnivory index generated by the Monte Carlo e2e_run_mc() function
#'
#' Create a two-panel plot showing the credible intervals of whole domain mean trophic level and omnivory index of each food web component over the final year of a simulation
#' generated by the Monte Carlo function. The trophic indices originate from the NetIndices R- package.
#'
#' Post-processed data from the e2e_run_mc() function are stored in the file
#' /results/Modelname/Variantname/CredInt/CredInt_processed_networkresults-*.csv, where * represents the model run identifier
#' (model.ident) text embedded in the R-list object created by the e2e_read() function.
#'
#' Optionally the function can read an example data set for one of the two North Sea model variants supplied with the package.
#'
#' The plots shows a set of box and whisker plots of a) mean trophic level, and b) omnivory index derived by the NetIndices package
#' from the annual flux-matrix for the final year of each model run. In each case the maximum likelihood model is shown by a red bar. The median of
#' the credible values distribution
#' is shown by a black bar. The box indicates the 50% credible interval (spanning quartiles of the cumulative likelihood
#' of simulated values). Whiskers span the 99% credible interval.
#'
#' Within each panel the box and whisker plots for each food web component are ordered by the maximum likelihood mean trophic level.
#'
#' NOTE that the NetIndices package assigns detritis and dissolved nutrients as trophic level 1, so that the phytoplankton and macrophytes are assigned trophic level 2.
#'
#' For details of how the distribution of credible output values from StrathE2E are calculated see the help information for the e2e_run_mc() function.
#'
#' @param model R-list object defining the model configuration compiled by the e2e_read() function
#' @param use.example Logical. If TRUE then use pre-computed example data from the internal North Sea model rather than user-generated data (default=FALSE)
#'
#' @return Graphical display in a new graphics window
#'
#' @importFrom graphics polygon
#' @importFrom graphics text
#'
#' @noRd
#
# ------------------------------------------------------------------------------

plot_final_year_trophic_data_with_ci <- function(model, use.example=FALSE) {

start_par = par()$mfrow
on.exit(par(mfrow = start_par))

	resultsdir	<- elt(model, "setup", "resultsdir")
	model.ident	<- elt(model, "setup", "model.ident")
	model.name 	<- elt(model, "setup", "model.name")
	model.variant 	<- elt(model, "setup", "model.variant")

corefilename<-"CredInt_processed_networkresults"

if(use.example==TRUE){
	networkresults <- get.example.results(model.name, model.variant, corefilename, CREDINT_DIR)
}

if(use.example==FALSE){
	credpath	<- makepath(resultsdir, CREDINT_DIR)
	credfile	<- csvname(credpath, corefilename, model.ident)
	message("Reading credible interval processed data from '", credfile, "'")

	networkresults<- readcsv(credfile, row.names=1)
}


#Columns 11-26 = mean trophic level
#Columns 37-52 = ominvory indices

#Code to assemble a dataframe of trophic index data, ordered by baseline trophic level

baseline_trophiclevel<-networkresults[1,11:26]

baseline_omnivory<-networkresults[1,37:52]

asc_trophiclevel<-order(baseline_trophiclevel)
trlevel_names<- c("Macrophytes",
                  "Phytoplankton",
                  "Omniv.zooplankton",
                  "Carn.zooplankton",
                  "Planktiv.fish larvae",
                  "Dem.fish larvae",
                  "Planktiv.fish",
                  "Mig.fish",
                  "Dem.fish",
                  "Susp/dep.benth larvae",
                  "Carn/scav.benth larvae",
                  "Susp/dep.benth",
                  "Carn/scav.benth",
                  "Birds",
                  "Pinnipeds",
                  "Cetaceans")
#trlevel_names[asc_trophiclevel]
#baseline_trophiclevel[asc_trophiclevel]

#...............
#function to get trophic level data 
get_trophicdata<-function(zz){

gotdata<-data.frame(tl=rep(NA,6),oz=rep(NA,6))

gotdata$tl<-networkresults[,(10+zz)]
gotdata$oz<-networkresults[,(36+zz)]

return(gotdata)
}
#...............


trophiclevel_table <-  data.frame(rep(NA,6),rep(NA,6),rep(NA,6),rep(NA,6),rep(NA,6),rep(NA,6),
                                  rep(NA,6),rep(NA,6),rep(NA,6),rep(NA,6),rep(NA,6),rep(NA,6),
                                  rep(NA,6),rep(NA,6),rep(NA,6),rep(NA,6))
names(trophiclevel_table) <- trlevel_names[asc_trophiclevel]
rownames(trophiclevel_table)<-rownames(networkresults)

omnivory_table<-trophiclevel_table

baseline_trophiclevel_table<-trophiclevel_table
baseline_omnivory_table<-omnivory_table
for(jjk in 1:16){
baseline_trophiclevel_table[,jjk]<-baseline_trophiclevel[asc_trophiclevel[jjk]]
baseline_omnivory_table[,jjk]<-baseline_omnivory[asc_trophiclevel[jjk]]
}

for(jjk in 1:16){

tloxdata<-get_trophicdata(asc_trophiclevel[jjk])
trophiclevel_table[,jjk]<-tloxdata$tl
omnivory_table[,jjk]<-tloxdata$oz

}






trlevel_main<-list(stats=trophiclevel_table[2:6,],n=rep(100,16),conf=NULL,out=numeric(length=0),names=names(trophiclevel_table))
trlevel_base<-list(stats=baseline_trophiclevel_table[2:6,],n=rep(100,16),conf=NULL,out=numeric(length=0),names=names(trophiclevel_table))

omniv_main<-list(stats=omnivory_table[2:6,],n=rep(100,16),conf=NULL,out=numeric(length=0),names=names(omnivory_table))
omniv_base<-list(stats=baseline_omnivory_table[2:6,],n=rep(100,16),conf=NULL,out=numeric(length=0),names=names(omnivory_table))


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


par(mfrow=c(2,1))
par(mar=c(4,10,.5,1.75))
   bxp(trlevel_main,boxwex=0.75,at=seq(1,16),ylim=c(2,max(trophiclevel_table[2:6,],na.rm=T)),cex.axis=0.75,show.names=TRUE,horizontal=TRUE,las=1,boxcol="black",whiskcol="black",whisklty="solid",medcol="black",staplecol="black")
   boxplot(baseline_trophiclevel_table[2:6,],add=TRUE,boxwex=0.6,boxcol="red",whiskcol="red",whisklty="solid",medlty=0,staplecol="red",horizontal=TRUE,yaxt="n",xaxt="n")
   mtext("Mean trophic level",cex=1,side=1,line=2)

   bxp(omniv_main,boxwex=0.75,at=seq(1,16),ylim=c(0,max(omnivory_table[2:6,],na.rm=T)),cex.axis=0.75,show.names=TRUE,horizontal=TRUE,las=1,boxcol="black",whiskcol="black",whisklty="solid",medcol="black",staplecol="black")
   boxplot(baseline_omnivory_table[2:6,],add=TRUE,boxwex=0.6,boxcol="red",whiskcol="red",whisklty="solid",medlty=0,staplecol="red",horizontal=TRUE,yaxt="n",xaxt="n")
   mtext("Omnivory index",cex=1,side=1,line=2)


}