visualize() --> plot.grasps()

Author: Carol-seven

NAMESPACE

@@ -1,9 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(plot,grasps)
 export(gen_prec_sbm)
 export(grasps)
 export(sparsify_block_banded)
-export(visualize)
 import(ggplot2)
 importFrom(Rcpp,evalCpp)
 importFrom(Rdpack,reprompt)

R/gen_prec_sbm.R

@@ -119,7 +119,7 @@
 #' @importFrom igraph as_adjacency_matrix sample_sbm
 #'
 #' @return
-#' A list with the following components:
+#' An object with S3 class "grasps" containing the following components:
 #' \describe{
 #' \item{Omega}{The precision matrix with SBM block structure.}
 #' \item{Sigma}{The covariance matrix, i.e., the inverse of \code{Omega}.}
@@ -141,6 +141,7 @@
 #'                     weight.dists = list(my_gamma, "unif"),
 #'                     weight.paras = list(NULL, c(min = 0, max = 5)),
 #'                     cond.target = 100)
+#' plot(sim)
 #'
 #' @export
 
@@ -238,7 +239,9 @@ gen_prec_sbm <- function(d,
   ## covariance matrix
   Sigma <- solve(Omega)
 
-  return(list(Omega = Omega, Sigma = Sigma,
-              sparsity = sum(Omega == 0) / length(Omega),
-              membership = membership))
+  result <- list(Omega = Omega, Sigma = Sigma,
+                 sparsity = sum(Omega == 0) / length(Omega),
+                 membership = membership)
+  class(result) <- "grasps"
+  return(result)
 }

R/grasps.R

@@ -115,7 +115,7 @@
 #' @importFrom Rdpack reprompt
 #'
 #' @return
-#' A list containing the following components:
+#' An object with S3 class "grasps" containing the following components:
 #' \describe{
 #' \item{hatOmega}{The estimated precision matrix.}
 #' \item{lambda}{The optimal regularization parameter.}
@@ -137,11 +137,26 @@
 #' \item{IC.score}{The information criterion score for each parameter
 #' combination when \code{crit} is set to \code{"AIC"}, \code{"BIC"},
 #' \code{"EBIC"}, or \code{"HBIC"}.}
+#' \item{membership}{The group membership.}
 #' }
 #'
 #' @references
 #' \insertAllCited{}
 #'
+#' @examples
+#' ## reproducibility for everything
+#' set.seed(1234)
+#'
+#' ## n-by-d data matrix
+#' X <- matrix(rnorm(200), 10, 20)
+#'
+#' ## group membership
+#' membership <- c(rep(1,5), rep(2,5), rep(3,4), rep(4,6))
+#'
+#' ## lasso, BIC
+#' res <- grasps(X = X, membership = membership, penalty = "lasso", crit = "BIC")
+#' plot(res)
+#'
 #' @export
 
 grasps <- function(X, n = nrow(X), membership, penalty,
@@ -156,34 +171,37 @@ grasps <- function(X, n = nrow(X), membership, penalty,
   d <- ncol(X)
 
   if (length(membership) != d) {
-    stop("The length of 'membership' must equal the column dimension of X!")
+    stop('The length of `membership` must equal the column dimension of `X`!')
   }
   if (!penalty %in% c("lasso", "adapt", "mcp", "scad")) {
-    stop("Error in penalty! Available options: 'lasso', 'adapt', 'mcp', 'scad'.")
+    stop('Error in `penalty`!\nAvailable options: "lasso", "adapt", "mcp", "scad".')
+  }
+  if (!crit %in% c("AIC", "BIC", "EBIC", "HBIC", "CV")) {
+    stop('Error in `crit`!\nAvailable options: "AIC", "BIC", "EBIC", "HBIC", "CV".')
   }
   if (!all(lambda > 0)) {
-    stop("The parameter 'lambda' must be positive!")
+    stop('The parameter `lambda` must be positive!')
   }
   if (!all(alpha >= 0 & alpha <= 1)) {
-    stop("The parameter 'alpha' must be in [0,1]!")
+    stop('The parameter `alpha` must be in [0,1]!')
   }
   if (rho <= 0) {
-    stop("The parameter 'rho' must be positive!")
+    stop('The parameter `rho` must be positive!')
   }
   if (tau.incr <= 1) {
-    stop("The parameter 'tau.incr' must be greater than 1!");
+    stop('The parameter `tau.incr` must be greater than 1!')
   }
   if (tau.decr <= 1) {
-    stop("The parameter 'tau.decr' must be greater than 1!");
+    stop('The parameter `tau.decr` must be greater than 1!')
   }
   if (nu <= 1) {
-    stop("The parameter 'nu' must be greater than 1!");
+    stop('The parameter `nu` must be greater than 1!')
   }
   if (tol.abs <= 0) {
-    stop("The parameter 'tol.abs' must be positive!");
+    stop('The parameter `tol.abs` must be positive!')
   }
   if (tol.rel <= 0) {
-    stop("The parameter 'tol.rel' must be positive!");
+    stop('The parameter `tol.rel` must be positive!')
   }
 
   if (isSymmetric(X)) {
@@ -262,11 +280,11 @@ grasps <- function(X, n = nrow(X), membership, penalty,
     if (crit == "CV") {
 
       if(is.null(X)) {
-        stop("CV requires the n-by-p data matrix!")
+        stop('CV requires the n-by-d data matrix!')
       }
 
       if (kfold < 2 | kfold > n) {
-        stop("'kfold' must be between 2 and the row dimension of X!")
+        stop('`kfold` must be between 2 and the row dimension of `X`!')
       }
 
       CV <- ADMMsggm_CV(X = X, group_idx = group.idx, penalty = penalty,
@@ -323,6 +341,8 @@ grasps <- function(X, n = nrow(X), membership, penalty,
 
   }
 
+  result$membership <- membership
+  class(result) <- "grasps"
   return(result)
 
 }

R/visualize.R R/plot.grasps.RR/visualize.R renamed to R/plot.grasps.R

@@ -1,18 +1,10 @@
-#' Visualize a Matrix with Group Boundaries
+#' Plot Function for "grasps" (Visualize a Matrix with Group Boundaries)
 #'
 #' @description
-#' Visualize a matrix (e.g., a precision, covariance, or adjacency matrix) as a
-#' heatmap with optional group-based reordering and dashed boundary lines
+#' Visualize a precision matrix as a heatmap with dashed boundary lines
 #' separating group blocks.
 #'
-#' @param mat A p-by-p matrix.
-#'
-#' @param membership An integer vector specifying the group membership.
-#' The length of \code{membership} must be consistent with the dimension p.
-#'
-#' @param reorder_by_group A boolean (default = \code{FALSE}) specifying whether
-#' to reorder both rows and columns of \code{mat} according to \code{membership},
-#' such that variables from the same group are contiguous.
+#' @param x An object with S3 class "grasps".
 #'
 #' @param colors A vector of colors specifying an n-color gradient scale for
 #' the fill aesthetics.
@@ -35,49 +27,47 @@
 #' }
 #'
 #' ## block-structured precision matrix based on SBM
-#' sim <- gen_prec_sbm(p = 100, K = 5,
+#' sim <- gen_prec_sbm(d = 100, K = 5,
 #'                     within.prob = 0.5, between.prob = 0.05,
 #'                     weight.dists = list(my_gamma, "unif"),
 #'                     weight.paras = list(NULL, c(min = 0, max = 1)),
-#'                     min.eig = 0.1)
+#'                     cond.target = 100)
 #'
 #' ## visualization
-#' visualize(sim$Omega, sim$membership)
+#' plot(sim)
 #'
 #' @export
+#' @noRd
 
-visualize <- function(mat, membership, reorder_by_group = FALSE,
-                      colors = colorRampPalette(
-                        c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F",
-                          "yellow", "#FF7F00", "red", "#7F0000"))(256)) {
+plot.grasps <- function(x, colors = NULL, ...) {
 
-  p <- ncol(mat)
-  if (length(membership) != p) {
-    stop(sprintf("Length of 'membership' (%d) must equal the matrix dimension p (%d).",
-                 length(membership), p))
+  if(is.null(x$hatOmega)) {
+    mat <- x$Omega
+  } else {
+    mat <- x$hatOmega
   }
 
-  ## optionally reorder rows/cols by group membership
-  if (reorder_by_group) {
-    o <- order(membership)
-    mat <- mat[o, o]
-    membership <- membership[o]
+  d <- ncol(mat)
+  if (is.null(colors)) {
+    colors <- colorRampPalette(
+      c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F",
+        "yellow", "#FF7F00", "red", "#7F0000"))(256)
   }
 
   ## compute group sizes and boundary positions for dashed lines
-  grp_sizes <- table(membership)
+  grp_sizes <- table(x$membership)
   cuts <- cumsum(grp_sizes)
   bnds <- cuts[-length(cuts)] + 0.5 ## boundaries between groups
-  y_bnds <- p - bnds + 1  ## flipped y, note: scale_y_discrete(limits = rev)
+  y_bnds <- d - bnds + 1  ## flipped y, note: scale_y_discrete(limits = rev)
 
   ## declare
   Col <- Row <- value <- NULL
 
   ## plot data
-  labs <- paste0("V", seq_len(p))
+  labs <- paste0("V", seq_len(d))
   plotData <- data.frame(
-    Row = factor(rep(labs, times = p),  levels = labs),
-    Col = factor(rep(labs, each = p), levels = labs),
+    Row = factor(rep(labs, times = d),  levels = labs),
+    Col = factor(rep(labs, each = d), levels = labs),
     value = as.vector(mat),
     check.names = FALSE
   )
@@ -87,9 +77,6 @@ visualize <- function(mat, membership, reorder_by_group = FALSE,
   #   pivot_longer(-Row, names_to = "Col", values_to = "value") %>%
   #   mutate(Col = factor(Col, levels = levels(Row)))
 
-  ## sparsity (proportion of zero entries)
-  sparsity <- round(sum(mat == 0) / length(mat), 4)
-
   ## zero -> NA for better white rendering
   plotData$value[plotData$value == 0] <- NA
 
@@ -101,15 +88,17 @@ visualize <- function(mat, membership, reorder_by_group = FALSE,
     coord_fixed() +
     geom_tile() +
     guides(fill = guide_colourbar(title = NULL, barwidth = 0.5, barheight = 5)) +
-    scale_y_discrete(limits = rev) +
+    scale_x_discrete(limits = labs, expand = c(0, 0)) +
+    scale_y_discrete(limits = rev(labs), expand = c(0, 0)) +
     scale_fill_gradientn(colours = colors,
                          values = rescale(c(vmin, 0, vmax)),
                          limits = c(vmin, vmax),
                          na.value = "white") +
     geom_vline(xintercept = bnds, linetype = "dashed") +
     geom_hline(yintercept = y_bnds, linetype = "dashed") +
     labs(x = NULL, y = NULL,
-         title = sprintf("p = %d, Sparsity = %s", ncol(mat), sparsity)) +
+         title = sprintf("Dimension = %d, Sparsity = %s", d,
+                         round(sum(mat == 0) / length(mat), 4))) +
     theme_bw() +
     theme(axis.text = element_blank(),
           axis.ticks = element_blank(),

R/sparsify_block_banded.R

@@ -14,7 +14,7 @@
 #' where groups whose labels differ by at most \code{neighbor.range} are
 #' considered neighbors and kept in the mask.
 #'
-#' @return A list containing:
+#' @return An object with S3 class "grasps" containing the following components:
 #' \describe{
 #' \item{Omega}{The masked precision matrix.}
 #' \item{Sigma}{The covariance matrix, i.e., the inverse of \code{Omega}.}
@@ -33,18 +33,15 @@
 #'
 #' ## default: keep blocks within ±1 of each group
 #' res1 <- sparsify_block_banded(est$hatOmega, membership, neighbor.range = 1)
-#' ## visualization
-#' visualize(res1$Omega, res1$membership)
+#' plot(res1)
 #'
 #' ## wider band: keep blocks within ±2 of each group
 #' res2 <- sparsify_block_banded(est$hatOmega, membership, neighbor.range = 2)
-#' ## visualization
-#' visualize(res2$Omega, res2$membership)
+#' plot(res2)
 #'
 #' ## special case: block-diagonal matrix
 #' res3 <- sparsify_block_banded(est$hatOmega, membership, neighbor.range = 0)
-#' ## visualization
-#' visualize(res3$Omega, res3$membership)
+#' plot(res3)
 #'
 #' @export
 
@@ -64,7 +61,9 @@ sparsify_block_banded <- function(mat, membership, neighbor.range = 1) {
   ## compute covariance
   Sigma <- solve(Omega)
 
-  return(list(Omega = Omega, Sigma = Sigma,
-              sparsity = sum(Omega == 0) / length(Omega),
-              membership = membership))
+  result <- list(Omega = Omega, Sigma = Sigma,
+                 sparsity = sum(Omega == 0) / length(Omega),
+                 membership = membership)
+  class(result) <- "grasps"
+  return(result)
 }

README.Rmd

@@ -68,15 +68,26 @@ devtools::install_github("Carol-seven/grasps")
 ## Example
 
 
-```
+```{r}
 library(grasps)
 
+## reproducibility for everything
 set.seed(1234)
 
-X <- matrix(rnorm(200), 10, 20)
-membership <- c(rep(1,5), rep(2,5), rep(3,4), rep(4,6))
+## block-structured precision matrix based on SBM
+sim <- gen_prec_sbm(d = 100, K = 5,
+                    within.prob = 0.5, between.prob = 0.05,
+                    weight.dists = list("gamma", "unif"),
+                    weight.paras = list(c(shape = 20, scale = 5), c(min = 0, max = 1)),
+                    cond.target = 100)
+
+## synthetic data
+library(MASS)
+X <- MASS::mvrnorm(n = 50, mu = rep(0, ncol(sim$Sigma)), Sigma = sim$Sigma)
 
-res <- grasps(X, membership = membership, penalty = "lasso", crit = "BIC")
+## solution
+res <- grasps(X = X, membership = sim$membership, penalty = "lasso", crit = "BIC")
+plot(res)
 ```
 
 

README.md

@@ -40,12 +40,29 @@ You can install the development version of **grasps** from
 
 ## Example
 
-    library(grasps)
+``` r
+library(grasps)
 
-    X <- matrix(rnorm(200), 10, 20)
-    membership <- c(rep(1,5), rep(2,5), rep(3,4), rep(4,6))
+## reproducibility for everything
+set.seed(1234)
 
-    res <- grasps(X, membership = membership, penalty = "lasso", crit = "BIC")
+## block-structured precision matrix based on SBM
+sim <- gen_prec_sbm(d = 100, K = 5,
+                    within.prob = 0.5, between.prob = 0.05,
+                    weight.dists = list("gamma", "unif"),
+                    weight.paras = list(c(shape = 20, scale = 5), c(min = 0, max = 1)),
+                    cond.target = 100)
+
+## synthetic data
+library(MASS)
+X <- MASS::mvrnorm(n = 50, mu = rep(0, ncol(sim$Sigma)), Sigma = sim$Sigma)
+
+## solution
+res <- grasps(X = X, membership = sim$membership, penalty = "lasso", crit = "BIC")
+plot(res)
+```
+
+<img src="man/figures/README-unnamed-chunk-2-1.png" width="100%" />
 
 ## Reference