Update von Louis Files für CRAN. Danach merge mit der Version von Lukas vorgesehen.

.gitignore

@@ -4,3 +4,4 @@
 *.DS_Store*
 *.Rproj
 docs
+scripts

DESCRIPTION

@@ -1,12 +1,16 @@
 Package: CVN
 Title: Covariate-varying Networks
 Version: 1.1
-Author: Louis Dijkstra [aut, cre]
-Maintainer: Louis Dijkstra <[email protected]>
-Description: A package for inferring high-dimensional Gaussian graphical
-    networks that change with multiple discrete covariates
+Author: Louis Dijkstra [aut],
+        Lukas Burk [ctb],
+        DFG [fnd],
+        Ronja Foraita [cre]
+Maintainer: Ronja Foraita <[email protected]>
+Description: Inferring high-dimensional Gaussian graphical networks that 
+   change with multiple discrete covariates. 
+   Louis Dijkstra, Arne Godt, Ronja Foraita (2024) <arXiv:2407.19978>.
 License: GPL (>= 3)
-URL: https://github.com/bips-hb/CVN
+URL: https://github.com/bips-hb/CVN, https://arxiv.org/abs/2407.19978
 BugReports: https://github.com/bips-hb/CVN/issues
 Depends:
     R (>= 4.0.2),
@@ -33,4 +37,7 @@ Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.3.2
 Suggests: 
-    igraph
+    igraph,
+    knitr,
+    rmarkdown   
+VignetteBuilder: knitr

NAMESPACE

@@ -4,20 +4,13 @@ S3method(plot,cvn)
 S3method(print,cvn)
 export("%>%")
 export(CVN)
-export(check_correctness_input)
 export(create_edges_visnetwork)
-export(create_matrix_D)
 export(create_nodes_visnetwork)
 export(create_weight_matrix)
-export(determine_information_criterion)
-export(determine_information_criterion_cvn)
-export(estimate)
 export(find_core_graph)
 export(genlasso_wrapper)
-export(glasso)
 export(hamming_distance)
 export(hamming_distance_adj_matrices)
-export(hits_end_lambda_intervals)
 export(interpolate)
 export(matrix_A_inner_ADMM)
 export(plot_hamming_distances)
@@ -27,16 +20,20 @@ export(plot_weight_matrix)
 export(relative_difference_precision_matrices)
 export(set_attributes_to_edges_visnetwork)
 export(strip_cvn)
-export(updateTheta)
-export(updateY)
-export(updateZ_wrapper)
 export(visnetwork)
 export(visnetwork_cvn)
 import(Rcpp)
 import(ggplot2)
 import(visNetwork)
 importFrom(Matrix,Matrix)
+importFrom(crayon,green)
+importFrom(crayon,red)
 importFrom(doSNOW,registerDoSNOW)
+importFrom(dplyr,"%>%")
+importFrom(dplyr,arrange)
+importFrom(dplyr,filter)
+importFrom(dplyr,rename)
+importFrom(dplyr,select)
 importFrom(foreach,"%dopar%")
 importFrom(foreach,foreach)
 importFrom(magrittr,"%>%")
@@ -47,4 +44,7 @@ importFrom(stats,optim)
 importFrom(stats,runif)
 importFrom(stats,var)
 importFrom(utils,combn)
+importFrom(visNetwork,visIgraphLayout)
+importFrom(visNetwork,visNetwork)
+importFrom(visNetwork,visOptions)
 useDynLib(CVN)

R/CVN-package.R

@@ -1,3 +1,4 @@
+#' @aliases CVN-package NULL
 #' @keywords internal
 "_PACKAGE"
 
@@ -16,10 +17,13 @@
 #' @importFrom stats runif
 #' @importFrom stats var
 #' @importFrom utils combn
+#' @importFrom dplyr %>% filter
+#' @importFrom crayon green red
 #' @useDynLib CVN
 ## usethis namespace: end
 NULL
 
+
 # Silence global variable warning
 # Ideally this would be solved more cleanly, but this suffices
 utils::globalVariables(c(
@@ -31,3 +35,28 @@ utils::globalVariables(c(
   "value",
   "aic", "bic"
 ))
+
+
+#' Covariate-varying Networks
+#'
+#' Inferring high-dimensional Gaussian graphical networks that 
+#' change with multiple discrete covariates. 
+#'
+#' @name CVN-package
+#' @docType package
+#' @aliases CVN-package
+#' @author Louis Dijkstra\cr Maintainer and contributors:
+#' Lukas Burk, Ronja Foraita <foraita@@leibniz-bips.de>
+#' @references Dijkstra L, Godt A, Foraita R 
+#' \emph{xxx (2024), Arxiv},
+#' \url{https://arxiv.org/abs/2407.19978}.
+#' @keywords graphical models 
+#' @rdname CVN-package
+#' @examples
+#'
+#' data(grid)
+#' W <- create_weight_matrix(type = "grid")
+#' 
+#' cvn <- CVN(grid, W, lambda1 = 1, lambda2 = 1:2, 
+#'            eps = 1e-3, maxiter = 1000, verbose = TRUE)
+NULL

R/CVN.R

@@ -14,9 +14,9 @@
 #'             Number of observations can differ
 #' @param W The \eqn{(m \times m)}-dimensional symmetric
 #'          weight matrix \eqn{W}
-#' @param lambda1 Vector with different \eqn{\lambda_1} LASSO penalty terms
+#' @param lambda1 Vector with different \eqn{\lambda_1}. LASSO penalty terms 
 #'                (Default: \code{1:2})
-#' @param lambda2 Vector with different \eqn{\lambda_2} global smoothing parameter values
+#' @param lambda2 <- Vector with different \eqn{\lambda_2}. The global smoothing parameter values 
 #'                (Default: \code{1:2})
 #' @param gamma1 A vector of \eqn{\gamma_1}'s LASSO penalty terms, where
 #'              \eqn{\gamma_1 = \frac{2 \lambda_1}{m p (1 - p)}}. If \code{gamma1}
@@ -97,7 +97,10 @@
 #'   \item{\code{minimal}}{If \code{TRUE}, \code{data}, \code{Theta} and \code{Sigma} are not added}
 #'   \item{\code{hits_border_aic}}{If \code{TRUE}, the optimal model based on the AIC hits the border of \eqn{(\lambda_1, \lambda_2)}}
 #'   \item{\code{hits_border_bic}}{If \code{TRUE}, the optimal model based on the BIC hits the border of \eqn{(\lambda_1, \lambda_2)}}
-#' @examples
+#'   
+#' @name CVN function
+#' @aliases CVN
+#' @examples 
 #' data(grid)
 #' m <- 9 # must be 9 for this example
 #'
@@ -109,10 +112,10 @@
 #' diag(W) <- 0
 #'
 #' # lambdas:
-#' lambda1 = 1:2
-#' lambda2 = 1:2
-#'
-#' (cvn <- CVN::CVN(grid, W, lambda1 = lambda1, lambda2 = lambda2,
+#' lambda1 = 1  # can also be lambda1 = 1:2 
+#' lambda2 = 1
+#' 
+#' (cvn <- CVN::CVN(grid, W, lambda1 = lambda1, lambda2 = lambda2, 
 #'                  eps = 1e-3, maxiter = 1000, verbose = TRUE))
 #' @export
 CVN <- function(data, W, lambda1 = 1:2, lambda2 = 1:2,
@@ -190,6 +193,7 @@ CVN <- function(data, W, lambda1 = 1:2, lambda2 = 1:2,
   # Compute the empirical covariance matrices --------------
   Sigma <- lapply(1:m, function(i) cov(data[[i]])*(n_obs[i] - 1) / n_obs[i])
 
+
   # initialize results list ------------
   global_res <- list(
     Theta    = list(),
@@ -248,13 +252,13 @@ CVN <- function(data, W, lambda1 = 1:2, lambda2 = 1:2,
       a <- (res$lambda1[i] / rho)^2 + 1
     } else {
       # Determine the value of the diagonal matrix A such that A - D'D > 0 (positive definite)
-      a <- CVN::matrix_A_inner_ADMM(W, res$lambda1[i] / rho, res$lambda2[i] / rho) + 1
+      a <- matrix_A_inner_ADMM(W, res$lambda1[i] / rho, res$lambda2[i] / rho) + 1
     }
 
     # Estimate the graphs -------------------------------------
     eta1 <- res$lambda1[i] / rho
     eta2 <- res$lambda2[i] / rho
-    CVN::estimate(m, p, W, Theta, Z, Y, a, eta1, eta2, Sigma, n_obs,
+         estimate(m, p, W, Theta, Z, Y, a, eta1, eta2, Sigma, n_obs,
                   rho, rho_genlasso,
                   eps, eps_genlasso,
                   maxiter, maxiter_genlasso, truncate = truncate,
@@ -265,7 +269,8 @@ CVN <- function(data, W, lambda1 = 1:2, lambda2 = 1:2,
   # go over each pair of penalty terms
   if (n_cores > 1) { # parallel
     est <- foreach(i = 1:(length(lambda1) * length(lambda2)),
-                   .options.snow = opts) %dopar% {
+                   .options.snow = opts,
+                   .export = c("estimate")) %dopar% {
                      estimate_lambda_values(i)
                    }
   } else { # sequential
@@ -284,14 +289,14 @@ CVN <- function(data, W, lambda1 = 1:2, lambda2 = 1:2,
     res$n_iterations[i] <- est[[i]]$n_iterations
 
     # determine the AIC
-    res$aic[i] <- CVN::determine_information_criterion(Theta = est[[i]]$Z,
+    res$aic[i] <- determine_information_criterion(Theta = est[[i]]$Z,
                                                      adj_matrices = est[[i]]$adj_matrices,
                                                      Sigma = Sigma,
                                                      n_obs = n_obs,
                                                      type = "AIC")
 
     # determine the BIC
-    res$bic[i] <- CVN::determine_information_criterion(Theta = est[[i]]$Z,
+    res$bic[i] <- determine_information_criterion(Theta = est[[i]]$Z,
                                                        adj_matrices = est[[i]]$adj_matrices,
                                                        Sigma = Sigma,
                                                        n_obs = n_obs,
@@ -310,7 +315,7 @@ CVN <- function(data, W, lambda1 = 1:2, lambda2 = 1:2,
 
   # determine whether the optimal model based on the AIC or BIC hits the border
   # of lambda1 and/or lambda2
-  hit <- CVN::hits_end_lambda_intervals(res)
+  hit <- hits_end_lambda_intervals(res)
   global_res$hits_border_aic <- hit$hits_border_aic
   global_res$hits_border_bic <- hit$hits_border_bic
 
@@ -326,70 +331,4 @@ CVN <- function(data, W, lambda1 = 1:2, lambda2 = 1:2,
   return(global_res)
 }
 
-#' Print Function for the CVN Object Class
-#'
-#' @export
-print.cvn <- function(x, ...) {  # TODO
-  cat(sprintf("Covariate-varying Network (CVN)\n\n"))
-
-  if (all(x$results$converged)) {
-    cat(crayon::green(sprintf("\u2713 all converged\n\n")))
-  } else {
-    cat(crayon::red(sprintf("\u2717 did not converge (maxiter of %d not sufficient)\n\n", x$maxiter)))
-  }
-
-  # print following variables
-  cat(sprintf("Number of graphs (m)    : %d\n", x$m))
-  cat(sprintf("Number of variables (p) : %d\n", x$p))
-  cat(sprintf("Number of lambda pairs  : %d\n\n", x$n_lambda_values))
-
-  # If x is an object returned by CVN::glasso() it doesn't have a $W element and this errors
-  # Causes example in glasso.R to error during R CMD check
-  # TODO: Check if glasso() behaves correctly, possibly adjust print method for CVN:glasso subclass?
-  if (!is.null(x$W)) {
-    cat(sprintf("Weight matrix (W):\n"))
-    print(Matrix::Matrix(x$W, sparse = T))
-  }
-
-  cat(sprintf("\n"))
-  print(x$results)
-}
-
-#' Strip CVN
-#'
-#' Function that removes most of the items to make the CVN object
-#' more memory sufficient. This is especially important when the
-#' graphs are rather larger
-#'
-#' @param cvn \code{cvn} object
-#'
-#' @return Reduced cvn where \code{Theta}, \code{data} and \code{Sigma}
-#'         are removed
-#' @export
-strip_cvn <- function(cvn) {
-
-  if ('Theta' %in% names(cvn)) {
-    cvn <- within(cvn, rm(Theta))
-  }
-
-  if ('data' %in% names(cvn)) {
-    cvn <- within(cvn, rm(data))
-  }
-
-  if ('Sigma' %in% names(cvn)) {
-    cvn <- within(cvn, rm(Sigma))
-  }
-
-  # set the variable keeping track of whether the cvn is
-  # striped to TRUE
-  cvn$minimal <- TRUE
-
-  return(cvn)
-}
 
-#' Plot Function for CVN Object Class
-#'
-#' @export
-plot.cvn <- function(x, ...) {
-  CVN::visnetwork_cvn(x)
-}

R/RcppExports.R

@@ -2,7 +2,6 @@
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #' Solving Generalized LASSO with fixed \eqn{\lambda = 1}
-#'
 #' Solves efficiently the generalized LASSO problem of the form
 #' \deqn{
 #'   \hat{\beta} = \text{argmin } \frac{1}{2} || y - \beta ||_2^2 + ||D\beta||_1
@@ -12,7 +11,7 @@
 #' We solve this optimization problem using an adaption of the ADMM
 #' algorithm presented in Zhu (2017).
 #'
-#' @param y The \eqn{y} vector of length \eqn{m}
+#' @param Y The \eqn{y} vector of length \eqn{m}
 #' @param W The weight matrix \eqn{W} of dimensions \eqn{m x m}
 #' @param m The number of graphs
 #' @param eta1 Equals \eqn{\lambda_1 / rho}
@@ -41,7 +40,6 @@ genlassoRcpp <- function(Y, W, m, eta1, eta2, a, rho, max_iter, eps, truncate) {
 }
 
 #' The \eqn{Z}-update Step
-#'
 #' A \code{C} implementation of the \eqn{Z}-update step. We
 #' solve a generalized LASSO problem repeatedly for each of the
 #' individual edges
@@ -50,6 +48,7 @@ genlassoRcpp <- function(Y, W, m, eta1, eta2, a, rho, max_iter, eps, truncate) {
 #' @param p The number of variables
 #' @param Theta A list of matrices with the \eqn{\Theta}-matrices
 #' @param Y A list of matrices with the \eqn{Y}-matrices
+#' @param W The weight matrix \eqn{W} of dimensions \eqn{m x m}
 #' @param eta1 Equals \eqn{\lambda_1 / rho}
 #' @param eta2 Equals \eqn{\lambda_2 / rho}
 #' @param a Value added to the diagonal of \eqn{-D'D} so that

R/check-correctness-input.R

@@ -15,9 +15,11 @@
 #' @param gamma2 A vector of \eqn{\gamma_2}'s global smoothing parameters. Note 
 #'               that \eqn{\gamma_2 = \frac{4 \lambda_2}{m(m-1)p(p-1)}}
 #' @param rho The \eqn{\rho} ADMM's penalty parameter 
+#' 
+#' @importFrom stats var
 #'        
 #' @seealso \code{\link{CVN}}
-#' @export
+#' @keywords internal
 check_correctness_input <- function(raw_data, W, 
                                     lambda1, lambda2, 
                                     gamma1, gamma2, 
@@ -36,7 +38,7 @@ check_correctness_input <- function(raw_data, W,
 
   ncols_per_dataset <- sapply(raw_data, function(X) ncol(X))
 
-  if (m > 1 && var(ncols_per_dataset) != 0) { 
+  if (m > 1 && stats::var(ncols_per_dataset) != 0) { 
     stop("The number of columns (variables) should 
               be the same for all datasets in raw_data") 
   }

R/create-matrix-D.R

@@ -13,6 +13,8 @@
 #'                        (Default: \code{TRUE})
 #' 
 #' @return A \eqn{((m \cdot (m+1)/2) \times m)}-dimensional matrix 
+#' 
+#' @importFrom utils combn
 #'         
 #' @references Tibshirani, R. J., & Taylor, J. (2011). 
 #'             The solution path of the generalized lasso. 
@@ -27,19 +29,8 @@
 #' lambda2 <- .4
 #' rho <- 1
 #' 
-#' CVN::create_matrix_D(W, lambda1, lambda2, rho) 
-#' #      [,1] [,2] [,3] [,4]
-#' # [1,]  0.2  0.0  0.0  0.0
-#' # [2,]  0.0  0.2  0.0  0.0
-#' # [3,]  0.0  0.0  0.2  0.0
-#' # [4,]  0.0  0.0  0.0  0.2
-#' # [5,]  0.4 -0.4  0.0  0.0
-#' # [6,]  0.4  0.0 -0.4  0.0
-#' # [7,]  0.4  0.0  0.0 -0.4
-#' # [8,]  0.0  0.4 -0.4  0.0
-#' # [9,]  0.0  0.4  0.0 -0.4
-#' # [10,]  0.0  0.0  0.4 -0.4
-#' @export
+#' CVN:::create_matrix_D(W, lambda1, lambda2, rho) 
+#' @keywords internal
 create_matrix_D <- function(W, lambda1, lambda2, rho = 1, remove_zero_row = TRUE) { 
 
   eta1 <- lambda1 / rho