R/plotKS.R

In lfabreti/convenience: Convergence assessment for phylogenetic inference

#' Histogram of the KS values
#' 
#' Plots the histogram of the KS values for the combination of all runs. The MCMC must have at least 2 runs
#' 
#' @importFrom grDevices dev.off pdf
#' @importFrom graphics abline hist layout legend lines par points polygon rect title
#' 
#' @param x A list of convenience.diag type
#' @param precision The precision of the mean estimates. Default is 0.01
#' @param fill_color The color to fill the histogram bars
#' @param filename The name of the file to save the plot
#' @param xlab Defines the label for the x-axis. To remove the label set to NA. Default = NULL
#' @param ylab Defines the label for the y-axis. To remove the label set to NA. Default = NULL
#' @param title Defines the title of the plot. To remove the title set to NA. Default = NULL
#' @param ... (various) Additional arguments passed to plot().
#' 
#' @return Histogram
#' 
#' @export


plotKS <- function(x, precision = 0.01, fill_color = NULL, filename = NULL, xlab = NULL, ylab = NULL, title = NULL, ...){
  
  col_threshold <- "gray69"
  
  if( is.null(fill_color) ){
    fill_color <- "cadetblue4"
  }
  
  if( !(is.null(filename)) ){
    pdf(file = filename, width = 4.5, height = 4.5)
  }
  
  minimumESS <- minESS(precision)
  minimumKS <- ksThreshold(0.01,minimumESS)
  
  KS_values <- vector()
  
  for (i in 1:ncol(x$continuous_parameters$compare_runs)) {
    KS_values <- c(KS_values, x$continuous_parameters$compare_runs[[i]])
  }
  y_topLim <- max(hist(KS_values, plot = FALSE)$counts)
  
  par(mar = c(4.1, 3.9, 2.1, 1.0))
  if( is.null(xlab)) xlab <- "Kolmogorov-Smirnov score" else xlab <- xlab
  if( is.null(ylab)) ylab <- "Counts" else ylab <- ylab
  if( is.null(title)) title <- "Kolmogorov-Smirnov test" else title <- title
  plot <- plot(NA,
               xlab = xlab,
               ylab = ylab,
               main = title,
               cex.main = 0.9,
               xlim = c( (min(minimumKS, KS_values) - 0.01), (max(minimumKS, KS_values) + 0.01)  ),
               ylim = c(0,y_topLim + 1),
               las=1,
               bty="l")
  plot <- rect(xleft = min(minimumKS, KS_values) - 0.01, ybottom = 0, xright = minimumKS, ytop = y_topLim+1, border = NA, col = "gray89")
  plot <- lines(x = c(minimumKS,minimumKS),y=c(0,y_topLim+1), col = col_threshold, lwd= 2, lty = 2)
  
  
  plot <- hist(KS_values, 
               xlim = c( (min(minimumKS, KS_values) - 0.01), (max(minimumKS, KS_values) + 0.01)  ),
               ylim = c(0,y_topLim + 1),
               breaks = seq(0, max(minimumKS, KS_values)+0.01, 0.0023),
               col = fill_color,
               border=F,
               las = 1,
               add=T,
               ...)

  
  if( !(is.null(filename)) ){
    dev.off()
  } 
  
}