R/plot.constr.hclust.R

In adespatial: Multivariate Multiscale Spatial Analysis

## **************************************************************************
##
##    (c) 2018-2024 Guillaume Guénard
##        Department de sciences biologiques,
##        Université de Montréal
##        Montreal, QC, Canada
##
##    **Plotting method for the constr.hclust-class**
##
##    This file is part of constr.hclust
##
##    constr.hclust is free software: you can redistribute it and/or modify
##    it under the terms of the GNU General Public License as published by
##    the Free Software Foundation, either version 3 of the License, or
##    (at your option) any later version.
##
##    constr.hclust is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##    GNU General Public License for more details.
##
##    You should have received a copy of the GNU General Public License
##    along with constr.hclust. If not, see <https://www.gnu.org/licenses/>.
##
## /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\
## |                                                            |
## |  CONSTRAINED HIERARCHICAL CLUSTERING                       |
## |  using (user-specified) criterion                          |
## |                                                            |
## |  C implementation of the Lance and Williams (1967) method  |
## |  for hierarchical clustering with or without the spatial   |
## |  contiguity constraint.                                    |
## |                                                            |
## \-----------------------------------------------------------*/
##
##    R source code file
##
## **************************************************************************
##
#' Plotting Method For Space- And Time-Constrained Clustering
#' 
#' Method \code{plot.constr.hclust} displays the results of space-constrained or
#' time-constrained agglomerative cluster analyses obtained from multivariate
#' dissimilarity matrices.
#' 
#' @usage \method{plot}{constr.hclust}(x, k, xlim, ylim, xlab, ylab, bg, col,
#' lty, lwd, col.links, links=FALSE, points=TRUE, pch=21L,
#' hybrids=c("change","single","none"), lty.hyb=1L, lwd.hyb=1, col.hyb="black",
#' plot=TRUE, axes=TRUE, cex=1, lwd.pt=1, invert.axes=FALSE, ...)
#' 
#' @param x A \code{\link{constr.hclust-class}} object.
#' @param k The number of clusters to delineate.
#' @param xlim Optional: limits, in abscissa, of the zone to be plotted.
#' @param ylim Optional: limits, in ordinate, of the zone to be plotted.
#' @param xlab Optional: labels for x axis annotation.
#' @param ylab Optional: labels for y axis annotation.
#' @param bg Optional: a colour or set of colours to be used for the point
#' backgrounds (one for each of the k clusters), see Details.
#' @param col Optional: a colour or set of colours to be used for the point
#' outlines (defaults to \code{"black"} when omitted).
#' @param lty Optional: reference line type (see \link{graphical parameters} for
#' details).
#' @param lwd Optional: reference line width (see \link{graphical parameters}
#' for details).
#' @param col.links Optional: a colour or set of colours to be used for the
#' links within a cluster (defaults to the point background colours when
#' omitted).
#' @param links Should segments be drawn to represent the edges (links)
#' (default: FALSE).
#' @param points Should observation points be drawn (default: TRUE).
#' @param pch Point character to display observations (default: 21, a circle
#' with a background colour).
#' @param hybrids How should hybrid segments be drawn (default: "change").
#' @param lty.hyb Line type to use for hybrid segments (default: lty).
#' @param lwd.hyb Width of hybrid segments with respect to lwd (default: 1).
#' @param col.hyb Colour of hybrid segments, when applicable (default: "black").
#' @param plot Should a new plotting window be opened first (default: TRUE).
#' @param axes Should the axes be displayed (default: TRUE).
#' @param cex Text and symbol magnification (see \link{graphical parameters})
#' (default: 1).
#' @param lwd.pt Line width around points with respect to lwd (default: 1).
#' @param invert.axes Should axes be inverted on the plot (default: FALSE).
#' @param ... Other \link{graphical parameters}.
#' 
#' @details The plotting method uses the coordinates provided by the user of
#' \code{\link{constr.hclust}} to display the observations. It cuts the tree
#' (see \link{cutree}) into \code{k} clusters and displays each cluster using
#' the indices returned by \code{\link{cutree}}. The point background colours
#' can be provided using argument \code{bg}. When they are omitted, default
#' colours are provided automatically by the function as either a set of
#' mutually contrasting colours (when \code{k <= 10}) or rainbow colours (when
#' \code{k > 10}). When \code{links = TRUE}, each edge is displayed as a segment
#' with a colour corresponding to the identity of the clusters located at its
#' ends. A special treatment is done for hybrids edges: those whose ends lie in
#' different clusters; it is controlled by argument \code{hybrids}. When
#' argument \code{hybrids="change"} (the default), hybrid links are represented
#' as segments whose colours change halfway. When \code{hybrids="single"},
#' hybrid edges are shown as single-color lines, whose color is given as
#' argument \code{col.hyb}, whereas \code{hybrids="none"} suppresses the drawing
#' of hybrid edges. Whenever hybrid edges are displayed, their width with
#' respect to the lwd value is controlled by argument \code{lwd.hyb}.
#' 
#' When argument \code{plot=FALSE}, no \code{plot} command is issued and the
#' points (and segments when \code{links = TRUE}) are drawn over an existing
#' plotting window. This functionality is to allow one to plot the result of a
#' constrained clustering over an existing map. In that case, arguments
#' \code{xlim}, \code{ylim}, \code{axes}, and all other
#' \link{graphical parameters} to which the method \link{plot} would responds
#' are ignored.
#' 
#' When disjoint clusters are present (i.e., when the graph provided to
#' \code{\link{constr.hclust}} is not entirely connected), the function does not
#' allow one to plot fewer clusters than the number of disjoint subsets; a
#' warning message is issued to notify the user.
#' 
#' @author Guillaume Guénard \email{guillaume.guenard@umontreal.ca}
#' and Pierre Legendre \email{pierre.legendre@@umontreal.ca}
#' 
#' @examples
#' 
#' ## Artificial map data from Legendre & Legendre (2012, Fig. 13.26)
#' ## n = 16
#' 
#' dat <- c(41,42,25,38,50,30,41,43,43,41,30,50,38,25,42,41)
#' coord.dat <- matrix(c(1,3,5,7,2,4,6,8,1,3,5,7,2,4,6,8,
#'                       4.4,4.4,4.4,4.4,3.3,3.3,3.3,3.3,
#'                       2.2,2.2,2.2,2.2,1.1,1.1,1.1,1.1),16,2)
#' 
#' ## Obtaining a list of neighbours:
#' library(spdep)
#' listW <- nb2listw(tri2nb(coord.dat), style="B")
#' links.mat.dat <- listw2mat(listW)
#' neighbors <- listw2sn(listW)[,1:2]
#' 
#' ## Calculating the (Euclidean) distance between points:
#' D.dat <- dist(dat)
#' 
#' ## Display the points:
#' plot(coord.dat, type='n',asp=1)
#' title("Delaunay triangulation")
#' text(coord.dat, labels=as.character(as.matrix(dat)), pos=3)
#' for(i in 1:nrow(neighbors))
#'     lines(rbind(coord.dat[neighbors[i,1],],
#'           coord.dat[neighbors[i,2],]))
#' 
#' ## Clustering with a contiguity constraint described by a list of
#' ## links:
#' grpWD2cst_constr_hclust <-
#'     constr.hclust(
#'         D.dat, method="ward.D2",
#'         neighbors, coord.dat)
#' 
#' ## Plot the results with k=5 clusters on a map:
#' plot(grpWD2cst_constr_hclust, k=5, links=TRUE, las=1,
#'      xlab="Eastings", ylab="Northings", cex=3, lwd=3)
#' 
#' ## Repeat the plot with other values of k (number of groups)
#' 
#' @importFrom grDevices dev.cur
#' @importFrom graphics par
#' @importFrom utils head
#' 
#' @evalNamespace "S3method(plot,constr.hclust)" ## Waiting for a better way...
plot.constr.hclust <- function(x, k, xlim, ylim, xlab, ylab, bg, col, lty, lwd,
                               col.links, links=FALSE, points=TRUE, pch=21L,
                               hybrids=c("change","single","none"), lty.hyb=1L,
                               lwd.hyb=1, col.hyb="black", plot=TRUE,
                               axes=TRUE, cex=1, lwd.pt=1, invert.axes=FALSE,
                               ...) {
  hybrids <- match.arg(hybrids)
  if(missing(lty)) lty <- par()$lty
  if(missing(lwd)) lwd <- par()$lwd
  
  if(!plot&&(dev.cur()==1L))
    stop("Use 'plot=FALSE' only for drawing over as existing plot!")
  if(is.null(x$coords)) {
    class(x) <- "hclust"
    plot(x, cex=cex, ...)
    return(invisible(NULL))
  }
  if(missing(bg)) {
    if(k > 10) {
      bg <- head(rainbow(1.2*k),k)
    } else
      bg <- head(
        c("blue", "gold", "grey70", "cadetblue2", "red", "orange3",
          "coral2", "green", "blueviolet", "grey30")
        ,k
      )
  } else {
    if(length(bg) < k) {
      if(length(bg) != 1L)
        warning("Argument 'bg' has length < k and had to be recycled")
      bg <- rep(bg, length.out=k)
    }
  }
  if(missing(col)) {
    col <- rep("black", k)
  } else
    if(length(col) < k) {
      if(length(col) != 1L)
        warning("Argument 'col' has length < k and had to be recycled")
      col <- rep(col, length.out=k)
    }
  if(any(nna <- is.na(x$height)))
    if(k < sum(nna) + 1L) {
      warning("Impossible to plot the cluster for k = ", k, " because ",
              "the graph provided involves ", sum(nna) + 1L, " non",
              "-connected clusters that cannot be linked at any ",
              "dissimilarity level. This cluster is plotted for k = ",
              sum(nna) + 1L)
      k <- sum(nna) + 1L
    }
  cl <- cutree(x,k)
  coords <- x$coords
  if(invert.axes) coords <- coords[,2L:1L]
  if(plot) {
    if(missing(xlim)) {
      xlim <- range(coords[,1L])
      if((xlim[2L]-xlim[1L])==0)
        xlim <- xlim + c(-0.5,0.5)
    }
    if(missing(ylim)) {
      ylim <- range(coords[,2L])
      if((ylim[2L]-ylim[1L])==0)
        ylim <- ylim + c(-0.5,0.5)
    }
    if(missing(xlab))
      xlab <- if(diff(range(coords[,1L]))==0) "" else "x"
    if(missing(ylab))
      ylab <- if(diff(range(coords[,2L]))==0) "" else "y"
    plot(NA, asp=1, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, axes=axes,
         type="n", cex=cex, ...)
  }
  if(links) {
    if(missing(col.links)) {
      col.links <- bg
    } else
      if(length(col.links) < k) {
        if(length(col.links) != 1L)
          warning("Argument 'col.links' has length < k and had to be recycled")
        col.links <- rep(col.links, length.out=k)
      }
    if(is.null(x$links))
      x$links <- cbind(1L:(nrow(coords)-1L),2L:nrow(coords))
    dev.hold()
    for(i in 1:nrow(x$links)) {
      ## i=1L
      clij <- cl[x$links[i,1L:2L]]
      if(clij[1L]==clij[2L]) {
        segments(coords[x$links[i,1L],1L],
                 coords[x$links[i,1L],2L],
                 coords[x$links[i,2L],1L],
                 coords[x$links[i,2L],2L],
                 col=col.links[clij[1L]], lty=lty, lwd=lwd)
      } else {
        ## The link is an hybrid
        if(hybrids=="change") {
          mid <- c(mean(coords[x$links[i,],1L]),
                   mean(coords[x$links[i,],2L]))
          segments(coords[x$links[i,1L],1L],
                   coords[x$links[i,1L],2L],
                   mid[1L], mid[2L],
                   col=col.links[clij[1L]],
                   lty=lty.hyb,
                   lwd=lwd*lwd.hyb)
          segments(mid[1L], mid[2L],
                   coords[x$links[i,2L],1L],
                   coords[x$links[i,2L],2L],
                   col=col.links[clij[2L]],
                   lty=lty.hyb,
                   lwd=lwd*lwd.hyb)
        } else if(hybrids=="single")
          segments(coords[x$links[i,1L],1L],
                   coords[x$links[i,1L],2L],
                   coords[x$links[i,2L],1L],
                   coords[x$links[i,2L],2L],
                   col=col.hyb,
                   lty=lty.hyb,
                   lwd=lwd*lwd.hyb)
      }
    }
    dev.flush()
  }
  if(points)
    points(x=coords[,1L], y=coords[,2L], bg=bg[cl], lty=lty,
           lwd=lwd.pt*lwd, cex=cex, pch=pch, col=col[cl], ...)
  return(invisible(NULL))
}
##