R/cutoffplot.R

In matrixcut: Determines Clustering Threshold Based on Similarity Values

#' @name cutoffplot
#' @title Draws a Plot Showing the Number of Similarity Values as a Function
#' of the Cutoff Value.
#' @description The function takes a square sequence similarity matrix and creates a plot
#' showing the number of similarity values above a given cutoff value, drawn on the x-axis,
#' ranging from 0 to 1. The function also depicts the inflection point with a vertical
#' line. Upper and lower bounds can be provided between which the inflection point will
#' be found (if it exists).
#'
#' Version 0.0.1.
#' Author: Dr. Matthew Cserhati
#' Email: csmatyi@protonmail.com
#' March 19, 2023
#'
#' @importFrom grDevices colorRampPalette dev.off jpeg
#' @importFrom graphics box legend
#' @importFrom inflection stats
#' @importFrom graphics abline lines
#' @importFrom loess predict
#' @importFrom bese graph_from_adjacency_matrix
#'
#' @param mx a square sequence similarity matrix
#' @param lower_bound lower bounds for calculating the inflection point in, default = 0
#' @param upper_bound upper bounds for calculating the inflection point in, default = 1
#' @return A plot showing the number of similarity values as a function of the cutoff threshold.
#'
#' @references Mullner. <ArXiv:1109.2378>
#' @references Cserhati, Carter. (2020, Journal of Creation 34(3):41-50), <https://dl0.creation.com/articles/p137/c13759/j34-3_64-73.pdf>
#'
#' @examples
#' cutoffplot(liliales,0.7,0.85)
#' cutoffplot(liliales)
#'
#'@export
utils::globalVariables(c("exit","abline","loess","lines","predict","bese","graph_from_adjacency_matrix"))
cutoffplot <- function(mx,lower_bound=0,upper_bound=1) {
  if (!is.numeric(mx)) {
    print("Data frame contains non-numeric values!")
    exit()
  }
  if (lower_bound < -1 | lower_bound > 1 | upper_bound < -1 | upper_bound > 1) {
    print("One of the boundaries is out of range (0-1)!")
    exit()
  }

  mxv <- mx[upper.tri(mx)]
  minv <- floor(10*min(mxv))/10
  xs <- c()
  xc <- seq(minv,1,0.001)
  for (i in xc) {
    nt <- sum(as.integer(mxv >= i))
    xs <- c(xs,nt)
  }

  xc2 <- xc[xc >= lower_bound & xc <= upper_bound]
  xs2 <- xs[xc >= lower_bound & xc <= upper_bound]
  b <- inflection::bese(xc2,xs2,0)

  lx <- loess(xs~xc, span=0.1)
  xs_pred <- predict(lx,xc)

  plot(xc,xs,pch=20,xlab="cutoff",ylab="no. pairs",cex=0.5)
  lines(xc,xs_pred,col="blue")
  if (!is.nan(b$iplast)) {
    if (lower_bound > minv) {abline(v=lower_bound,col="green")}
    if (upper_bound < 1) {abline(v=upper_bound,col="green")}
    abline(v=b$iplast,col="red")
  }
}