R/cluster.R

In RDelta: Parse and Query DELTA-Format Taxonomic Descriptions

#' Hierarchical Cluster Analysis
#'
#' Performs hierarchical cluster analysis on a distance matrix generated by
#' the dist() function. This function reads the distance matrix and item names
#' from the DIST output files and creates a dendrogram.
#'
#' @param notu Number of operational taxonomic units (items) in the dataset.
#' @param method The agglomeration method to be used. This should be one of
#'   "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty"
#'   (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC). Default is "average".
#'
#' @return Invisibly returns a list with the following components:
#'   \describe{
#'     \item{hclust}{The \code{hclust} object produced by
#'       \code{\link[stats]{hclust}}, which can be further processed or plotted
#'       by the user.}
#'     \item{cophenetic}{The cophenetic correlation coefficient (a single
#'       numeric value between -1 and 1) measuring how faithfully the
#'       dendrogram represents the original distance matrix. Values closer to
#'       1 indicate a better fit.}
#'     \item{dist}{The \code{dist} object representing the input distance
#'       matrix.}
#'   }
#' @export
#'
#' @examples
#' # Requires the CONFOR and DIST executables on the PATH.
#' \donttest{
#' try({
#'   delta <- load_delta(
#'     system.file("extdata", "chars_viola", package = "RDelta"),
#'     system.file("extdata", "items_viola", package = "RDelta"),
#'     system.file("extdata", "specs_viola", package = "RDelta")
#'   )
#'   dist(delta)
#'   res <- cluster(delta$get_items_nb())
#'   res$cophenetic
#'   res <- cluster(delta$get_items_nb(), method = "complete")
#' })
#' }
cluster <- function(notu, method = "average") {
  # Validate input
  if (!is.numeric(notu) || notu <= 0 || notu != round(notu)) {
    stop("'notu' must be a positive integer")
  }

  # Check if distance matrix files exist
  if (!file.exists("dist.dis")) {
    stop("Distance matrix file 'dist.dis' not found. Run dist() first.")
  }
  if (!file.exists("dist.nam")) {
    stop("Item names file 'dist.nam' not found. Run dist() first.")
  }

  # Read distance matrix
  # The dist.dis file contains the lower triangular values (including diagonal)
  mat <- matrix(0, notu, notu)
  mat[row(mat) >= col(mat)] <- scan("dist.dis")

  # Read item names and assign to matrix
  names <- scan("dist.nam", what = "character")
  colnames(mat) <- names

  # Convert to distance object
  df.dist <- stats::as.dist(mat)

  # Perform hierarchical clustering
  hc <- stats::hclust(df.dist, method = method)
  dc <- stats::cophenetic(hc)
  r  <- stats::cor(df.dist, dc)

  # Plot dendrogram
  oldpar <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(oldpar), add = TRUE)
  graphics::par(mar = c(5, 4, 4, 10) + 0.1)
  graphics::plot(
    stats::as.dendrogram(hc),
    xlab = "Distance",
    ylab = "OTU's",
    main = paste("Dendrogram", "(r = ", format(r, digits = 4), ")"),
    horiz = TRUE,
    edgePar = list(col = "blue", lwd = 3)
  )

  # Return results invisibly so the user can access them programmatically
  invisible(list(
    hclust     = hc,
    cophenetic = r,
    dist       = df.dist
  ))
}