R/plot_mds.R

In dbrobust: Robust Distance-Based Visualization and Analysis of Mixed-Type Data

#' Plot MDS Results with Grouped Scatter and Density Plots (Internal)
#'
#' This internal function performs classical or weighted Multidimensional Scaling (MDS)
#' on a given distance matrix and visualizes the resulting coordinates using a
#' pairwise scatterplot matrix with density plots on the diagonal. Grouping information
#' can be provided for colored visual separation.
#'
#' @name plot_mds
#'
#' @param dist_mat A distance matrix or object convertible to a distance matrix.
#' @param k Integer. Number of dimensions to retain in MDS (default is 3).
#' @param weights Optional numeric vector of weights for weighted MDS. If \code{NULL},
#'   classical MDS is performed.
#' @param group Optional factor or vector indicating group membership for observations,
#'   used for coloring plots.
#' @param main_title Optional character string for the main plot title.
#'
#' @return A \code{ggmatrix} object from \code{GGally} representing the
#'   pairs plot with scatterplots and density plots.
#'
#' @details
#' This is an internal helper function. It is not recommended to call \code{plot_mds()} directly.
#' Instead, use \code{\link{visualize_distances}()}, which wraps this function.
#'
#' Weighted MDS is performed with \code{vegan::wcmdscale} if \code{weights} are provided;
#' otherwise, classical MDS (\code{cmdscale}) is used. Diagonal panels show density plots
#' by group, and off-diagonal panels show scatter plots by group.
#'
#' @examples
#' # Load example dataset
#' data("Data_HC_contamination", package = "dbrobust")
#' # Subset of 20 rows
#' Data_small <- Data_HC_contamination[1:20, ]
#'
#' # Define variable types
#' cont_vars <- c("V1", "V2", "V3", "V4")
#' cat_vars  <- c("V5", "V6", "V7")
#' bin_vars  <- c("V8", "V9")
#'
#' # Use column 'w_loop' as weights
#' w <- Data_small$w_loop
#'
#' # Compute robust distances using GGower
#' dist_sq_ggower <- dbrobust:::robust_distances(
#'   data = Data_small,
#'   cont_vars = cont_vars,
#'   bin_vars  = bin_vars,
#'   cat_vars  = cat_vars,
#'   w = w,
#'   alpha = 0.10,
#'   method = "ggower"
#' )
#'
#' # Create factor indicating Normal vs Outlier
#' n_obs <- nrow(dist_sq_ggower)
#' group_vec <- rep("Normal", n_obs)
#' group_vec[attr(dist_sq_ggower, "outlier_idx")] <- "Outlier"
#' group_factor <- factor(group_vec, levels = c("Normal", "Outlier"))
#'
#' # Plot MDS
#' dbrobust:::plot_mds(
#'   dist_mat = dist_sq_ggower,
#'   k = 2,
#'   group = group_factor,
#'   main_title = "MDS of Data_HC_contamination (GGower) with Outliers"
#' )
#'
#' @importFrom ggplot2 ggplot aes geom_point geom_density scale_color_manual scale_fill_manual theme theme_classic element_text element_blank element_rect coord_cartesian ggtitle margin element_line
#' @importFrom GGally ggpairs
#' @importFrom vegan wcmdscale
#' @keywords internal
if (getRversion() >= "2.15.1") utils::globalVariables(c("Group"))
plot_mds <- function(dist_mat, k = 3, weights = NULL, group = NULL, main_title = NULL) {
  # Convert to distance object if needed
  dist_obj <- if (inherits(dist_mat, "dist")) dist_mat else as.dist(dist_mat)

  # Perform classical or weighted MDS
  if (is.null(weights)) {
    mds_res <- cmdscale(dist_obj, k = k, eig = TRUE)
  } else {
    mds_res <- vegan::wcmdscale(dist_obj, w = weights, k = k, eig = TRUE)
  }

  # Convert to data frame and label dimensions
  coords <- as.data.frame(mds_res$points)
  colnames(coords) <- paste0("PCo", 1:k)
  coords$Group <- if (!is.null(group)) factor(group) else factor("All")
  plot_data <- coords[, c(paste0("PCo", 1:k), "Group")]

  # Color palette
  n_groups <- length(levels(coords$Group))
  pal <- get_custom_palette(n_groups)

  # Diagonal density plot (filled + thin outline)
  custom_diag <- function(data, mapping, ...) {
    ggplot(data = data, mapping = mapping) +
      geom_density(aes(fill = Group, color = Group), alpha = 0.4, linewidth = 0.8) +  # filled density
      geom_density(aes(color = Group), linewidth = 0.3, fill = NA) +                   # thin outline
      scale_fill_manual(values = pal) +
      scale_color_manual(values = pal) +
      theme_classic(base_size = 14) +
      theme(
        axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.grid.major = ggplot2::element_line(color = "white", linewidth = 0.2),
        panel.grid.minor = ggplot2::element_line(color = "white", linewidth = 0.1),
        legend.position = "none",
        panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
      )
  }

  # Off-diagonal scatter plot
  custom_offdiag <- function(data, mapping, ...) {
    ggplot(data = data, mapping = mapping) +
      geom_point(aes(color = Group), size = 1.5, alpha = 0.7) +
      scale_color_manual(values = pal) +
      theme_classic(base_size = 14) +
      theme(
        legend.position = "none",
        panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
        panel.grid.major = ggplot2::element_line(color = "grey90", linewidth = 0.3),
        panel.grid.minor = ggplot2::element_line(color = "grey95", linewidth = 0.15),
        axis.title = element_text(face = "bold")
      ) +
      coord_cartesian()
  }

  # Create GGally pairs plot
  ggp <- GGally::ggpairs(
    plot_data,
    columns = 1:k,
    mapping = aes(color = Group, fill = Group),
    diag = list(continuous = custom_diag),
    upper = list(continuous = custom_offdiag),
    lower = list(continuous = custom_offdiag)
  )

  # Theme adjustments
  ggp <- ggp + theme(
    strip.text = element_text(face = "bold", size = 12, color = "darkblue", margin = margin(5, 5, 5, 5)),
    strip.background = element_rect(color = "black", linewidth = 1, fill = "azure2"),
    legend.position = "right",
    legend.title = element_text(face = "bold", size = 10),
    legend.text = element_text(size = 12)
  )

  # Add main title
  if (!is.null(main_title)) {
    ggp <- ggp + ggtitle(main_title) +
      theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))
  }

  return(ggp)
}