prec_to_adj

Author: Carol-seven

NAMESPACE

@@ -8,6 +8,7 @@ export(compute_penalty)
 export(gen_prec_sbm)
 export(grasps)
 export(performance)
+export(prec_to_adj)
 export(sparsify_block_banded)
 import(ggforce)
 import(ggplot2)

R/performance.R

@@ -91,7 +91,7 @@
 #' }
 #'
 #' @example
-#' inst/example/ex-performance.R
+#' inst/example/ex-grasps.R
 #'
 #' @export
 

R/prec_to_adj.R

@@ -0,0 +1,90 @@
+#' Adjacency Matrix from Precision Matrix
+#'
+#' @description
+#' Convert a precision matrix to a partial-correlation-based adjacency matrix.
+#'
+#' @param prec.mat A numeric precision matrix.
+#'
+#' @param diag.zero A logical value (default = TRUE) specifying whether to
+#' set the diagonal entries of the adjacency matrix to 0.
+#' If \code{diag.zero = FALSE}, the diagonal entries are set to 1 for a weighted
+#' network. For an unweighted network (\code{weighted = FALSE}), the diagonal is
+#' always forced to 0 to avoid self-loops.
+#'
+#' @param absolute A logical value (default = FALSE) specifying whether to
+#' take the absolute values of the partial correlations.
+#'
+#' @param threshold A nonnegative numeric value (default = \code{NULL})
+#' specifying the threshold for edge filtering.
+#' Entries with absolute values smaller than the threshold are set to 0.
+#'
+#' @param weighted A logical value (default = TRUE) specifying whether to
+#' return a weighted adjacency matrix.
+#' If \code{weighted = FALSE}, the matrix is a binary adjacency matrix with
+#' entries equal to 0 or 1.
+#'
+#' @return
+#' A numeric adjacency matrix.
+#'
+#' @details
+#' For a precision matrix \eqn{\Omega}, the partial correlation between nodes
+#' \eqn{i} and \eqn{j} is computed as
+#' \deqn{\rho_{ij} = - \Omega_{ij} / \sqrt{\Omega_{ii}\Omega_{jj}}.}
+#'
+#' @example
+#' inst/example/ex-grasps.R
+#'
+#' @export
+
+prec_to_adj <- function(prec.mat, diag.zero = TRUE, absolute = FALSE,
+                        threshold = NULL, weighted = TRUE) {
+
+  if (!is.matrix(prec.mat)) {
+    prec.mat <- as.matrix(prec.mat)
+  }
+  if (any(!is.finite(prec.mat))) {
+    stop("The precision matrix contains non-finite values!")
+  }
+  if (!isSymmetric(prec.mat)) {
+    stop("The precision matrix must be symmetric!")
+  }
+  d <- diag(prec.mat)
+  if (any(d <= 0)) {
+    stop("All diagonal entries of the precision matrix must be positive!")
+  }
+
+  ## partial-correlation-based adjacency matrix
+  norm_mat <- outer(sqrt(d), sqrt(d), "*")
+  adj_mat <- -prec.mat / norm_mat
+
+  ## optional: diagonal entries
+  if (diag.zero) {
+    diag(adj_mat) <- 0
+  } else {
+    diag(adj_mat) <- 1
+  }
+
+  ## optional: absolute value
+  if (absolute) {
+    adj_mat <- abs(adj_mat)
+  }
+
+  ## optional: thresholding
+  if (!is.null(threshold)) {
+    if (!is.numeric(threshold) || length(threshold) != 1 || threshold < 0) {
+      stop("The threshold must be a nonnegative scalar!")
+    }
+    adj_mat[abs(adj_mat) < threshold] <- 0
+  }
+
+  ## optional: binary adjacency matrix
+  if (!weighted) {
+    adj_mat <- 1 * (adj_mat != 0)
+    diag(adj_mat) <- 0
+    if (!diag.zero) {
+      message("For an unweighted network, the diagonal is forced to 0 to avoid self-loops.")
+    }
+  }
+
+  return(adj_mat)
+}

inst/example/ex-grasps.R

@@ -13,15 +13,21 @@ sim <- gen_prec_sbm(p = 30, K = 3,
 ## visualization
 plot(sim)
 
-## n-by-d data matrix
+## n-by-p data matrix
 library(MASS)
 X <- mvrnorm(n = 20, mu = rep(0, 30), Sigma = sim$Sigma)
 
-## adapt, HBIC
-res <- grasps(X = X, membership = sim$membership, penalty = "adapt", crit = "HBIC")
+## precision matrix: adaptive lasso; BIC
+prec <- grasps(X = X, membership = sim$membership, penalty = "adapt", crit = "BIC")
 
-## visualization
-plot(res)
+## precision matrix visualization
+plot(prec)
 
 ## performance
-performance(hatOmega = res$hatOmega, Omega = sim$Omega)
+performance(hatOmega = prec$hatOmega, Omega = sim$Omega)
+
+## adjacency matrix: diagonal = 0; raw partial correlations;
+##                   no thresholding; weighted network
+adj <- prec_to_adj(prec$hatOmega,
+                   diag.zero = TRUE, absolute = FALSE,
+                   threshold = NULL, weighted = TRUE)