R/plot_sensitivity_results.R

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

#
# plot_sensitivity_results.R
#
#' Create a plot of results from the Morris method sensitivity analysis.
#'
#' Reads processed results generated by the function e2e_run_sens(), or the function e2e_process_sens_mc(), and plots the Elementary Effect mean (x-axis; magnitude of sensitivity)
#' against Elementary Effect standard deviation (y-axis; strength of interactions between parameters).
#'
#' Processed results from the function e2e_run_sens(), or the function e2e_process_sens_mc(), are stored in the csv file
#' /results/Modelname/Variantname/sorted_parameter_elementary_effects-*.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 example data for one of the two North Sea model variants supplied with the package.
#'
#' Each symbol in the plot represents a single parameter in the model. The parameters are colour-coded to indicate 6 different types - fitted and fixed
#' parameters of the ecology model, fishing fleet model parameters, fishery harvest ratios, environmental drivers, and physical configuration parameters.
#'
#' The plot also shows a wedge formed by the two dashed lines. These correspond to +/-2 standard errors of the mean, so for points falling outside of the
#' wedge there is a significant expectation that the distribution of elementary effects is non-zero. For points falling within the wedge the distribution of
#' elementary effects is not significantly different from zero.
#' 
#' For details of how the Elementary Effect values are derived for each parameter see help(e2e_run_sens)
#'
#' @references Morris, M.D. (1991). Factorial sampling plans for preliminary computational experiments. Technometrics, 33, 161-174.
#'
#' @param model R-list object defining the model configuration compiled by the e2e_read() function
#' @param use.example Logical. If TRUE 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 text
#'
#' @noRd
#
# ------------------------------------------------------------------------------

plot_sensitivity_results <- 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")

if(use.example==TRUE){
	Sorted_SENS_results <- get.example.results(model.name, model.variant, "sorted_parameter_elementary_effects")
}

if(use.example==FALSE){
	sensfile	<- csvname(resultsdir, "sorted_parameter_elementary_effects", model.ident)
	message("Reading sensitivity analysis data from '", sensfile, "'")

	check.exists(sensfile)

	Sorted_SENS_results <- readcsv(sensfile)
}

	NTRAJ <- Sorted_SENS_results$Ntrajectories[1]

	senscritname <- Sorted_SENS_results$criterion[1]

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

#	xamin<-1.05*min(Sorted_SENS_results$EEmean,na.rm=TRUE)
#	xamax<-max(1, (1.05*max(Sorted_SENS_results$EEmean,na.rm=TRUE)) )

	xamin<-min(Sorted_SENS_results$EEmean,na.rm=TRUE) - 0.7 *( max(Sorted_SENS_results$EEmean,na.rm=TRUE) - min(Sorted_SENS_results$EEmean,na.rm=TRUE) )
	xamaxt<-max((1.05*max(Sorted_SENS_results$EEmean,na.rm=TRUE)) )
        xaxspan<-xamaxt-xamin
        xamax<-xamaxt + (xaxspan*0.05)

	yamax<-1.05*max(Sorted_SENS_results$EEsd,na.rm=TRUE)

	#black = fitted (0)
	#red = fixed (1)
	#green = fishing parametars (2)
	#blue = harvest rates (3)
	#orange = drivers (4)
	#yellow = physical setup (5)
        leglabels<-c("fitted ecology parameters","fixed ecology parameters","fishing fleet parameters","harvest ratios","environmental drivers","physical setup")
        legcolours<-c("black","red","green","blue","orange","yellow")
        legsymbols<-rep(16,6)
        legsize<-rep(1,6)

	par(mfrow=c(1,1))
	par(mar=c(4.5,5,2.5,0.5))

	plot(Sorted_SENS_results$EEmean[which(Sorted_SENS_results$fixfit==0)],Sorted_SENS_results$EEsd[which(Sorted_SENS_results$fixfit==0)],type="p",pch=legsymbols[1],cex=legsize[1],col=legcolours[1],xlim=c(xamin,xamax),ylim=c(0,yamax),xaxt="n",yaxt="n",ann=FALSE)
	points(Sorted_SENS_results$EEmean[which(Sorted_SENS_results$fixfit==1)],Sorted_SENS_results$EEsd[which(Sorted_SENS_results$fixfit==1)],type="p",pch=legsymbols[2],cex=legsize[2],col=legcolours[2])
	points(Sorted_SENS_results$EEmean[which(Sorted_SENS_results$fixfit==2)],Sorted_SENS_results$EEsd[which(Sorted_SENS_results$fixfit==2)],type="p",pch=legsymbols[3],cex=legsize[3],col=legcolours[3])
	points(Sorted_SENS_results$EEmean[which(Sorted_SENS_results$fixfit==3)],Sorted_SENS_results$EEsd[which(Sorted_SENS_results$fixfit==3)],type="p",pch=legsymbols[4],cex=legsize[4],col=legcolours[4])
	points(Sorted_SENS_results$EEmean[which(Sorted_SENS_results$fixfit==4)],Sorted_SENS_results$EEsd[which(Sorted_SENS_results$fixfit==4)],type="p",pch=legsymbols[5],cex=legsize[5],col=legcolours[5])
	points(Sorted_SENS_results$EEmean[which(Sorted_SENS_results$fixfit==5)],Sorted_SENS_results$EEsd[which(Sorted_SENS_results$fixfit==5)],type="p",pch=legsymbols[6],cex=legsize[6],col=legcolours[6])


	abline(0,(sqrt(NTRAJ)/2),lty="dashed")
	abline(0,-(sqrt(NTRAJ)/2),lty="dashed")
#	axis(side=1,at=c(-5,-4,-3,-2,-1,-0.5,0,0.5,1,2,3,4,5,6),las=1,cex.axis=1)
#	axis(side=2,at=c(0,0.25,0.5,1,1.5,2,3,4,5,6),las=1,cex.axis=1)
	axis(side=1,las=1,cex.axis=1)
	axis(side=2,las=1,cex.axis=1)
	mtext("Elementary Effect mean",cex=1.3,side=1,line=2.25)
	mtext("Elementary Effect standard deviation",cex=1.3,side=2,line=2.8)
	legend("bottomleft",  bg="transparent", leglabels, pch=legsymbols, cex=legsize*0.75, col=legcolours, title= "Parameter types")
	mtext(paste("Sensitivity analysis criterion: ",senscritname,sep="") ,cex=1.0,side=3,line=0.5)


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

}