Add parameters rank and mode_orth to enable using rank-R TGCCA

* Important change in this version: tau and sparsity are converted to matrices directly in select_analysis

NAMESPACE

@@ -48,6 +48,7 @@ export(rgcca_stability)
 export(rgcca_transform)
 importFrom(Deriv,Deriv)
 importFrom(MASS,ginv)
+importFrom(abind,abind)
 importFrom(caret,confusionMatrix)
 importFrom(caret,multiClassSummary)
 importFrom(caret,postResample)
@@ -66,7 +67,6 @@ importFrom(matrixStats,rowSds)
 importFrom(matrixStats,rowSums2)
 importFrom(methods,is)
 importFrom(rlang,.data)
-importFrom(stats,complete.cases)
 importFrom(stats,cor)
 importFrom(stats,median)
 importFrom(stats,model.matrix)

R/block.R

@@ -45,19 +45,30 @@ new_sparse_block <- function(x, j, sparsity, tol = 1e-08, ...) {
   )
 }
 
-new_tensor_block <- function(x, j, ..., class = character()) {
+new_tensor_block <- function(x, j, rank, mode_orth, ..., class = character()) {
   new_block(
-    x, j, factors = NULL, weights = NULL, ..., class = c(class, "tensor_block")
+    x, j, rank = rank, mode_orth = mode_orth, factors = NULL,
+    weights = NULL, ..., class = c(class, "tensor_block")
+  )
+}
+
+new_regularized_tensor_block <- function(x, j, rank, mode_orth, tau, ...) {
+  new_tensor_block(
+    x, j, rank = rank, mode_orth = mode_orth, tau = tau,
+    M = NULL, ..., class = "tensor_regularized_block"
   )
 }
 
 ### Utility method to choose the adequate class
-create_block <- function(x, j, bias, na.rm, tau, sparsity, tol) {
+create_block <- function(x, j, bias, na.rm, tau, sparsity,
+                         tol, rank, mode_orth) {
   if (length(dim(x)) > 2) {         # TGCCA
     if (tau < 1) {
-      res <- new_regularized_tensor_block(x, j, tau, bias = bias, na.rm = na.rm)
+      res <- new_regularized_tensor_block(
+        x, j, rank, mode_orth, tau, bias = bias, na.rm = na.rm
+      )
     } else {
-      res <- new_tensor_block(x, j, bias = bias, na.rm = na.rm)
+      res <- new_tensor_block(x, j, rank, mode_orth, bias = bias, na.rm = na.rm)
     }
   } else {
     if (sparsity < 1) {             # SGCCA

R/block_init.R

@@ -44,14 +44,20 @@ block_init.dual_regularized_block <- function(x, init = "svd") {
 block_init.tensor_block <- function(x, init = "svd") {
   if (init == "svd") {
     x$factors <- lapply(seq_along(dim(x$x))[-1], function(m) {
-      svd(apply(x$x, m, c), nu = 0, nv = 1)$v
+      svd(apply(x$x, m, c), nu = 0, nv = x$rank)$v
     })
   } else {
     x$factors <- lapply(seq_along(dim(x$x))[-1], function(m) {
-      svd(matrix(rnorm(dim(x$x)[m]), dim(x$x)[m]), nu = 1, nv = 0)$u
+      if (m == x$mode_orth) {
+        svd(matrix(
+          rnorm(dim(x$x)[m] * x$rank), dim(x$x)[m]
+        ), nu = x$rank, nv = 0)$u
+      } else {
+        matrix(rnorm(dim(x$x)[m] * x$rank), dim(x$x)[m])
+      }
     })
   }
-  x$weights <- rep(1 / sqrt(1), 1)
+  x$weights <- rep(1 / sqrt(x$rank), x$rank)
 
   return(block_project(x))
 }

R/block_update.R

@@ -29,8 +29,17 @@ block_update.tensor_block <- function(x, grad) {
     grad_m <- grad_m %*% khatri_rao(
       Reduce(khatri_rao, rev(x$factors[-seq_len(m)])), other_factors
     )
-    x$factors[[m]] <- grad_m %*% diag(x$weights, nrow = length(x$weights))
-    x$factors[[m]] <- x$factors[[m]] / norm(x$factors[[m]], type = "2")
+    if (m == x$mode_orth) {
+      SVD <- svd(
+        grad_m %*% diag(x$weights, nrow = x$rank), nu = x$rank, nv = x$rank
+      )
+      x$factors[[m]] <- SVD$u %*% t(SVD$v)
+    } else {
+      x$factors[[m]] <- grad_m %*% diag(x$weights, nrow = x$rank)
+      x$factors[[m]] <- apply(
+        x$factors[[m]], 2, function(y) y / norm(y, type = "2")
+      )
+    }
 
     other_factors <- khatri_rao(x$factors[[m]], other_factors)
   }

R/check_blocks.r

@@ -33,6 +33,7 @@
 #'   dimnames in the same order.
 #' }
 #' @inheritParams rgcca
+#' @importFrom abind abind
 #' @noRd
 check_blocks <- function(blocks, quiet = FALSE, response = NULL) {
   blocks <- check_blocks_is_list(blocks)
@@ -258,7 +259,7 @@ check_blocks_align <- function(blocks, m = 1) {
     y <- array(NA, dim = extra_dim)
     dimnames(y)[[m]] <- missing_names
     rownames(y) <- missing_names
-    return(abind::abind(x, y, along = m))
+    return(abind(x, y, along = m))
   })
 
   # Align blocks using names on dimension m

R/checks.R

@@ -219,8 +219,12 @@ check_ncomp <- function(ncomp, blocks, min = 1, superblock = FALSE,
         "only one number of components must be specified (superblock)."
       )
     }
+    is_superblock_present <- "superblock" %in% names(blocks)
+    ncomp <- rep(max(ncomp), ifelse(
+      is_superblock_present, length(blocks), length(blocks) + 1
+    ))
     max_ncomp <- ifelse(
-      "superblock" %in% names(blocks),
+      is_superblock_present,
       NCOL(blocks[[length(blocks)]]),
       sum(vapply(blocks, NCOL, FUN.VALUE = integer(1)))
     )
@@ -240,6 +244,11 @@ check_ncomp <- function(ncomp, blocks, min = 1, superblock = FALSE,
 
   ncomp <- elongate_arg(ncomp, blocks)
   check_size_blocks(blocks, "ncomp", ncomp)
+
+  if (!is.null(response)) {
+    ncomp[response] <- max(ncomp[-response])
+  }
+
   ncomp <- vapply(
     seq_along(ncomp),
     function(x) {
@@ -319,7 +328,7 @@ check_size_blocks <- function(blocks, x, y = x, n_row = NULL,
 }
 
 check_penalty <- function(penalty, blocks, method = "rgcca", superblock = FALSE,
-                          ncomp = NULL) {
+                          ncomp = 1) {
   penalty <- elongate_arg(penalty, blocks)
   is_matrix <- is.matrix(penalty)
   DIM <- dim(penalty)
@@ -348,7 +357,10 @@ check_penalty <- function(penalty, blocks, method = "rgcca", superblock = FALSE,
     )
   }
 
-  if (is_matrix) penalty <- matrix(penalty, DIM[1], DIM[2])
+  penalty <- matrix(penalty, ncol = size)
+  if (nrow(penalty) < ncomp) {
+    penalty <- matrix(penalty, nrow = ncomp, ncol = size, byrow = TRUE)
+  }
 
   return(penalty)
 }
@@ -378,6 +390,49 @@ check_tau <- function(tau) {
   invisible(tau)
 }
 
+check_mode_orth <- function(mode_orth, blocks) {
+  mode_orth <- elongate_arg(mode_orth, blocks)
+  mode_orth <- check_integer("mode_orth", mode_orth, min = 1, type = "vector")
+  check_size_blocks(blocks, "mode_orth", mode_orth)
+  lapply(seq_along(blocks), function(j) {
+    if (length(dim(blocks[[j]])) > 2 &&
+        mode_orth[j] > length(dim(blocks[[j]])) - 1) {
+      stop_rgcca(
+        "mode_orth should be comprise between ", 1 ,
+        " and ", length(dim(blocks[[j]])) - 1, " (that is the number of
+          modes of block ", j, " minus 1)."
+      )
+    }
+  })
+  return(mode_orth)
+}
+
+check_rank <- function(rank, blocks, mode_orth, ncomp) {
+  rank <- elongate_arg(rank, blocks)
+  p <- ifelse(is.vector(rank), 1, nrow(rank))
+  J <- length(blocks)
+  rank <- vapply(seq_along(rank), function(x) {
+    y <- check_integer("rank", rank[x], min = 1)
+    j <- (x - 1) %% J + 1
+    if (length(dim(blocks[[j]])) > 2 &&
+        y > dim(blocks[[j]])[mode_orth[j] + 1]) {
+      stop_rgcca(
+        "rank[", x, "] should be comprise between ", 1 ,
+        " and ", dim(blocks[[j]])[mode_orth[j] + 1], " (that is the number of",
+        " variables of the mode bearing the orthogonality for block ", j, ")."
+      )
+    }
+    return(y)
+  }, FUN.VALUE = 1L)
+
+  rank <- matrix(rank, nrow = p)
+  check_size_blocks(blocks, "rank", rank, n_row = ncomp)
+  if (p < ncomp) {
+    rank <- matrix(rank, nrow = ncomp, ncol = J, byrow = TRUE)
+  }
+  return(rank)
+}
+
 check_scheme <- function(scheme) {
   if (mode(scheme) != "function") {
     scheme <- tolower(scheme)

R/format_output.R

@@ -43,11 +43,11 @@ format_output <- function(func_out, rgcca_args, opt, blocks) {
   names(func_out$AVE$AVE_X_cor) <- names_AVE
 
   ### Set names and shave
-  array_idx <- which(
-    vapply(blocks, function(x) length(dim(x)) > 2, FUN.VALUE = logical(1L))
+  array_idx <- vapply(
+    blocks, function(x) length(dim(x)) > 2, FUN.VALUE = logical(1L)
   )
 
-  for (j in array_idx) {
+  for (j in seq_along(blocks)[array_idx]) {
     func_out$factors[[j]] <- lapply(
       seq_along(func_out$factors[[j]]),
       function(m) {
@@ -60,7 +60,7 @@ format_output <- function(func_out, rgcca_args, opt, blocks) {
     rownames(func_out$a[[j]]) <- do.call(paste, c(grid, sep = "_"))
   }
 
-  for (j in seq_along(blocks)[-array_idx]) {
+  for (j in seq_along(blocks)[!array_idx]) {
     rownames(func_out$a[[j]]) <- colnames(blocks[[j]])
   }
 

R/get_rgcca_args.R

@@ -23,6 +23,7 @@ get_rgcca_args <- function(object, default_args = list()) {
       tol = default_args$tol,
       init = tolower(default_args$init),
       bias = default_args$bias,
+      rank = default_args$rank,
       quiet = default_args$quiet,
       scale = default_args$scale,
       ncomp = default_args$ncomp,
@@ -34,6 +35,7 @@ get_rgcca_args <- function(object, default_args = list()) {
       response = default_args$response,
       NA_method = tolower(default_args$NA_method),
       comp_orth = default_args$comp_orth,
+      mode_orth = default_args$mode_orth,
       n_iter_max = default_args$n_iter_max,
       connection = default_args$connection,
       superblock = default_args$superblock,

R/plot.rgcca_cv.R

@@ -30,7 +30,7 @@ plot.rgcca_cv <- function(x, type = c("sd", "quantile"),
   }
 
   best <- which(apply(
-    x$params, 1, function(z) identical(z, x$best_params)
+    x$params, 1, function(z) identical(z, drop(x$best_params))
   ))
 
   idx_order <- seq_len(nrow(df))

R/plot.rgcca_permutation.R

@@ -32,7 +32,7 @@ plot.rgcca_permutation <- function(x,
 
   # Mark the best parameter set
   best <- which(apply(
-    x$params[idx_order, ], 1, function(z) identical(z, x$best_params)
+    x$params[idx_order, ], 1, function(z) identical(z, drop(x$best_params))
   ))
   df$label[best] <- "Best parameter set"
 

R/rgcca.R

@@ -174,6 +174,18 @@
 #' orthogonal block components or orthogonal block weight vectors.
 #' @param A Deprecated argument, please use blocks instead.
 #' @param C Deprecated argument, please use connection instead.
+#' @param rank Either an integer, an integer vector of
+#' size \eqn{J} or an integer matrix
+#' of dimension \eqn{\textrm{max}(\textrm{ncomp}) \times J} giving the rank
+#' of the decomposition sought for the canonical vectors in TGCCA.
+#' If block \eqn{j} is an array with at least three dimensions, rank must be
+#' comprised between 1 and the number of variables on the mode bearing the
+#' orthogonality constraint. See \textrm{mode_orth}.
+#' @param mode_orth Either an integer or an integer vector of size \eqn{J}
+#' designating the mode which associated set of factors will be orthogonal
+#' in the decomposition sought by TGCCA. If block \eqn{j} is an array with
+#' \eqn{d > 2} dimensions, \textrm{mode_orth} must be comprised between
+#' 1 and \eqn{d - 1}.
 #' @return A fitted rgcca object.
 #' @return \item{Y}{A list of \eqn{J} elements. The jth element
 #' of the list \eqn{Y}
@@ -433,6 +445,7 @@ rgcca <- function(blocks, connection = NULL, tau = 1, ncomp = 1,
                   superblock = FALSE,
                   NA_method = "na.ignore", quiet = TRUE,
                   n_iter_max = 1000, comp_orth = TRUE,
+                  rank = 1, mode_orth = 1,
                   A = NULL, C = NULL) {
   # Check for deprecated arguments
   if (!missing(A)) {
@@ -478,14 +491,12 @@ rgcca <- function(blocks, connection = NULL, tau = 1, ncomp = 1,
   }
 
   ### Call the gcca function
-  gcca_args <- rgcca_args[c(
-    "connection", "ncomp", "scheme", "init", "bias", "tol",
-    "verbose", "superblock", "response", "n_iter_max", "comp_orth"
-  )]
+  gcca_args <- rgcca_args[-which(names(rgcca_args) %in% c(
+    "quiet", "scale", "scale_block", "method", "NA_method"
+  ))]
   gcca_args[["na.rm"]] <- na.rm
   gcca_args[["blocks"]] <- blocks
   gcca_args[["disjunction"]] <- opt$disjunction
-  gcca_args[[opt$param]] <- rgcca_args[[opt$param]]
   func_out <- do.call(rgcca_outer_loop, gcca_args)
 
   ### Format the output

R/rgcca_bootstrap.R

@@ -105,8 +105,7 @@ rgcca_bootstrap <- function(rgcca_res, n_boot = 100,
     rgcca_res$call$blocks <- Map(
       function(x, y) x[, y, drop = FALSE], rgcca_res$call$blocks, keep_var
     )
-    rgcca_res$call$tau <-
-      rgcca_res$call$sparsity <- rep(1, length(rgcca_res$blocks))
+    rgcca_res$call$tau[] <- rgcca_res$call$sparsity[] <- 1
 
     rgcca_res <- rgcca(rgcca_res)
   }

R/rgcca_cv.r

@@ -186,6 +186,8 @@ rgcca_cv <- function(blocks,
                      verbose = TRUE,
                      n_iter_max = 1000,
                      comp_orth = TRUE,
+                     rank = 1,
+                     mode_orth = 1,
                      ...) {
   ### Try to retrieve parameters from a rgcca object
   rgcca_args <- as.list(environment())
@@ -338,7 +340,7 @@ rgcca_cv <- function(blocks,
     par_type = param$par_type,
     params = param$par_value,
     validation = validation,
-    best_params = param$par_value[best_param_idx, ],
+    best_params = param$par_value[best_param_idx, , drop = FALSE],
     classification = model$classification,
     prediction_model = model$model_name
   )

R/rgcca_inner_loop.R

@@ -1,7 +1,9 @@
 rgcca_inner_loop <- function(A, C, g, dg, tau = rep(1, length(A)),
                              sparsity = rep(1, length(A)),
                              verbose = FALSE, init = "svd", bias = TRUE,
-                             tol = 1e-08, na.rm = TRUE, n_iter_max = 1000) {
+                             tol = 1e-08, na.rm = TRUE, n_iter_max = 1000,
+                             rank = rep(1, length(A)),
+                             mode_orth = rep(1, length(A))) {
   if (!is.numeric(tau)) {
     # From Schafer and Strimmer, 2005
     tau <- vapply(A, tau.estimate, na.rm = na.rm, FUN.VALUE = 1.0)
@@ -15,7 +17,9 @@ rgcca_inner_loop <- function(A, C, g, dg, tau = rep(1, length(A)),
 
   ### Initialization
   block_objects <- lapply(seq_along(A), function(j) {
-    create_block(A[[j]], j, bias, na.rm, tau[j], sparsity[j], tol)
+    create_block(
+      A[[j]], j, bias, na.rm, tau[j], sparsity[j], tol, rank[j], mode_orth[j]
+    )
   })
 
   block_objects <- lapply(block_objects, block_init, init = init)