update function VS() to allow a matrix/array input

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: BayesSurvive
 Title: Bayesian Survival Models for High-Dimensional Data
-Version: 0.0.4
-Date: 2024-07-10
+Version: 0.0.5
+Date: 2024-09-27
 Authors@R: c(person("Zhi", "Zhao", role=c("aut","cre"), email = "[email protected]"),
              person("Waldir", "Leoncio", role=c("aut")),
              person("Katrin", "Madjar", role=c("aut")),
@@ -22,7 +22,7 @@ Imports: Rcpp, ggplot2, GGally, mvtnorm, survival, riskRegression,
 Suggests: knitr, testthat
 LazyData: true
 NeedsCompilation: yes
-Packaged: 2024-07-10 08:34:27 UTC; zhiz
+Packaged: 2024-09-27 15:26:16 UTC; zhiz
 Author: Zhi Zhao [aut, cre],
   Waldir Leoncio [aut],
   Katrin Madjar [aut],

NEWS.md

@@ -6,6 +6,7 @@
 * Rename the output of MPM coefficients in function `coef.BayesSurvive()`
 * Added `cpp` argument to `BayesSurvive()` to allow for faster computation using `Rcpp`
 * Added validation to some `BayesSurvive()` arguments
+* Allow function `VS()` for an input with a matrix or array or a list consisting of matrices/arrays
 
 # BayesSurvive 0.0.3
 

R/BayesSurvive-package.R

@@ -14,7 +14,6 @@ NULL
 .onLoad <- function(libname, pkgname) {
   # CRAN OMP THREAD LIMIT
   Sys.setenv("OMP_THREAD_LIMIT" = 1)
-
 }
 
-## nocov end
+## nocov end

R/RcppExports.R

@@ -9,14 +9,6 @@ settingInterval_cpp <- function(y, delta_, s_, J_) {
     .Call(`_BayesSurvive_settingInterval_cpp`, y, delta_, s_, J_)
 }
 
-matProdVec <- function(x, y) {
-    .Call(`_BayesSurvive_matProdVec`, x, y)
-}
-
-sumMatProdVec <- function(x, y) {
-    .Call(`_BayesSurvive_sumMatProdVec`, x, y)
-}
-
 updateBH_cpp <- function(x_, beta_, J_, ind_r_d_, hPriorSh_, d_, c0_) {
     .Call(`_BayesSurvive_updateBH_cpp`, x_, beta_, J_, ind_r_d_, hPriorSh_, d_, c0_)
 }

R/VS.R

@@ -8,7 +8,8 @@
 #'
 #' @name VS
 #'
-#' @param x fitted object obtained with \code{BayesSurvive}
+#' @param x fitted object obtained with \code{BayesSurvive}, or a matrix/array,
+#' or a list consisting of matrices and arrays
 #' @param method variable selection method to choose from
 #' \code{c("CI", "SNC", "MPM", "FDR")}. Default is "FDR"
 #' @param threshold SNC threshold value (default 0.5) or the Bayesian expected false
@@ -66,65 +67,202 @@
 #'
 #' @export
 VS <- function(x, method = "FDR", threshold = NA, subgroup = 1) {
-  if (!inherits(x, "BayesSurvive")) {
-    stop("Use only with 'BayesSurvive' object!")
+  if (!(inherits(x, "BayesSurvive") || is.array(x) || is.list(x))) {
+    stop("Use only with 'BayesSurvive' object or a matrix or an array or a list!")
   }
-  if (!method %in% c("CI", "SNC", "SNC-BIC", "MPM", "FDR")) {
-    stop("'method' should be one of c('CI', 'SNC', 'SNC-BIC', 'MPM', 'FDR')!")
+  if (!inherits(x, "BayesSurvive") && is.list(x)) {
+    if (any(!sapply(x, function(xx) is.array(xx)))) {
+      stop("The list input has to consist of matrices and/or arrays!")
+    }
+    # If the input is a matrix or array, reformat it to be an list
+    x <- list(x)
+  }
+  if (!method %in% c("CI", "SNC", "MPM", "FDR")) { # "SNC-BIC",
+    stop("'method' should be one of c('CI', 'SNC', 'MPM', 'FDR')!")
   }
 
   if (method == "CI") {
-    betas <- coef.BayesSurvive(x,
-      type = "mean", CI = 95,
-      subgroup = subgroup
-    )
-    ret <- (betas$CI.lower > 0) | (betas$CI.upper < 0)
+    # for an input from "BayesSurvive"
+    if (inherits(x, "BayesSurvive")) {
+      betas <- coef.BayesSurvive(x,
+        type = "mean", CI = 95,
+        subgroup = subgroup
+      )
+      ret <- (betas$CI.lower > 0) | (betas$CI.upper < 0)
+    } else {
+      # for an output from a matrix or array or list
+      if (!is.list(x)) {
+        betas <- list()
+        # if (is.matrix(x)) {
+        #   betas$CI.lower <- apply(x, 2, quantile, 0.025)
+        #   betas$CI.upper <- apply(x, 2, quantile, 0.975)
+        # }
+        # if (is.array(x)) {
+        #   betas$CI.lower <- apply(x, c(2,3), quantile, 0.025)
+        #   betas$CI.upper <- apply(x, c(2,3), quantile, 0.975)
+        # }
+        # ret <- (betas$CI.lower > 0) | (betas$CI.upper < 0)
+      } else {
+        betas <- rep(list(NULL), length(x))
+        ret <- rep(list(NULL), length(x))
+        for (l in seq_len(length(x))) {
+          # if (length(dim(x[[l]])) == 2) {
+          #   betas[[l]]$CI.lower <- apply(x[[l]], 2, quantile, 0.025)
+          #   betas[[l]]$CI.upper <- apply(x[[l]], 2, quantile, 0.975)
+          # }
+          # if (length(dim(x[[l]])) == 3) {
+          #   betas[[l]]$CI.lower <- apply(x[[l]], c(2,3), quantile, 0.025)
+          #   betas[[l]]$CI.upper <- apply(x[[l]], c(2,3), quantile, 0.975)
+          # }
+          betas[[l]]$CI.lower <- apply(x[[l]], seq_len(length(dim(x[[l]])))[-1], quantile, 0.025)
+          betas[[l]]$CI.upper <- apply(x[[l]], seq_len(length(dim(x[[l]])))[-1], quantile, 0.975)
+          ret[[l]] <- (betas[[l]]$CI.lower > 0) | (betas[[l]]$CI.upper < 0)
+        }
+      }
+    }
   }
 
   if (method == "SNC") {
     if (is.na(threshold)) {
       threshold <- 0.5
     }
-    if (x$input$S > 1 || !x$input$MRF.G) {
-      x$output$beta.p <- x$output$beta.p[[subgroup]]
-    } 
-    beta_p <- x$output$beta.p[-(1:(x$input$burnin / x$input$thin + 1)), ]
 
-    ret <- rep(FALSE, NCOL(beta_p))
+    if (inherits(x, "BayesSurvive")) {
+      if (x$input$S > 1 || !x$input$MRF.G) {
+        x$output$beta.p <- x$output$beta.p[[subgroup]]
+      }
+      beta_p <- x$output$beta.p[-(1:(x$input$burnin / x$input$thin + 1)), ]
+
+      ret <- rep(FALSE, NCOL(beta_p))
 
-    for (j in seq_len(NCOL(beta_p))) {
-      if (sum(abs(beta_p[, j]) > sd(beta_p[, j])) / NROW(beta_p) > threshold) {
-        ret[j] <- TRUE
+      for (j in seq_len(NCOL(beta_p))) {
+        if (sum(abs(beta_p[, j]) > sd(beta_p[, j])) / NROW(beta_p) > threshold) {
+          ret[j] <- TRUE
+        }
+      }
+    } else {
+      # count the total number of parameters in the list
+      total_num <- 0
+      for (l in seq_len(length(x))) {
+        total_num <- total_num + prod(dim(x[[l]])[-1])
+      }
+
+      ret <- rep(list(NULL), length(x))
+      for (l in seq_len(length(x))) {
+        # define the size of each component in the list
+        ret[[l]] <- array(FALSE, dim = dim(x[[l]])[-1])
+
+        if (is.matrix(x[[l]])) { # for an matrix
+          for (j in seq_len(NCOL(x[[l]]))) {
+            if (sum(abs(x[[l]][, j]) > sd(x[[l]][, j])) / total_num > threshold) {
+              ret[[l]][j] <- TRUE
+            }
+          }
+        } else { # for an array
+          for (j in seq_len(dim(x[[l]])[2])) {
+            for (k in seq_len(dim(x[[l]])[3])) {
+              if (sum(abs(x[[l]][, j, k]) > sd(x[[l]][, j, k])) / total_num > threshold) {
+                ret[[l]][j, k] <- TRUE
+              }
+            }
+          }
+        }
       }
     }
   }
 
   if (method == "MPM") {
-    if (x$input$S > 1 || !x$input$MRF.G) {
-      x$output$gamma.margin <- x$output$gamma.margin[[subgroup]]
+    if (inherits(x, "BayesSurvive")) {
+      if (x$input$S > 1 || !x$input$MRF.G) {
+        x$output$gamma.margin <- x$output$gamma.margin[[subgroup]]
+      }
+      ret <- (x$output$gamma.margin >= 0.5)
+    } else {
+      ret <- rep(list(NULL), length(x))
+
+      for (l in seq_len(length(x))) {
+        # define the size of each component in the list
+        ret[[l]] <- array(FALSE, dim = dim(x[[l]])[-1])
+
+        if (is.matrix(x[[l]])) { # for an matrix
+          for (j in seq_len(NCOL(x[[l]]))) {
+            ret[[l]][j] <- (mean(x[[l]][, j] != 0) >= 0.5)
+          }
+        } else { # for an array
+          for (j in seq_len(dim(x[[l]])[2])) {
+            for (k in seq_len(dim(x[[l]])[3])) {
+              ret[[l]][j, k] <- (mean(x[[l]][, j, k] != 0) >= 0.5)
+            }
+          }
+        }
+      }
     }
-    ret <- (x$output$gamma.margin >= 0.5)
   }
 
   if (method == "FDR") {
     if (is.na(threshold)) {
       threshold <- 0.05
     }
-    if (x$input$S > 1 || !x$input$MRF.G) {
-      x$output$gamma.margin <- x$output$gamma.margin[[subgroup]]
-    }
-    gammas <- x$output$gamma.margin
-
-    sorted_gammas <- sort(gammas, decreasing = TRUE)
-    # computing the fdr
-    fdr <- cumsum((1 - sorted_gammas)) / seq_len(length(sorted_gammas))
-    # determine index of the largest fdr less than threshold
-    if (min(fdr) >= threshold) {
-      ret <- rep(FALSE, length(gammas))
+    if (inherits(x, "BayesSurvive")) {
+      if (x$input$S > 1 || !x$input$MRF.G) {
+        x$output$gamma.margin <- x$output$gamma.margin[[subgroup]]
+      }
+      gammas <- x$output$gamma.margin
+
+      sorted_gammas <- sort(gammas, decreasing = TRUE)
+      # computing the fdr
+      fdr <- cumsum((1 - sorted_gammas)) / seq_len(length(sorted_gammas))
+      # determine index of the largest fdr less than threshold
+      if (min(fdr) >= threshold) {
+        ret <- rep(FALSE, length(gammas))
+      } else {
+        thecut.index <- max(which(fdr < threshold))
+        gammas_threshold <- sorted_gammas[thecut.index]
+        ret <- (gammas > gammas_threshold)
+      }
     } else {
-      thecut.index <- max(which(fdr < threshold))
-      gammas_threshold <- sorted_gammas[thecut.index]
-      ret <- (gammas > gammas_threshold)
+      ret <- gammas <- rep(list(NULL), length(x))
+      gammas_vec <- NULL
+
+      # compute mPIPs
+      for (l in seq_len(length(x))) {
+        # define the size of each component in the list
+        gammas[[l]] <- array(FALSE, dim = dim(x[[l]])[-1])
+
+        if (is.matrix(x[[l]])) { # for an matrix
+          for (j in seq_len(NCOL(x[[l]]))) {
+            gammas[[l]][j] <- mean(x[[l]][, j] != 0)
+          }
+        } else { # for an array
+          for (j in seq_len(dim(x[[l]])[2])) {
+            for (k in seq_len(dim(x[[l]])[3])) {
+              gammas[[l]][j, k] <- mean(x[[l]][, j, k] != 0)
+            }
+          }
+        }
+        # save all mPIPs into a vector
+        gammas_vec <- c(gammas_vec, as.vector(gammas[[l]]))
+      }
+
+      sorted_gammas <- sort(gammas_vec, decreasing = TRUE)
+      # computing the fdr
+      fdr <- cumsum((1 - sorted_gammas)) / seq_len(length(sorted_gammas))
+      # determine index of the largest fdr less than threshold
+      if (min(fdr) >= threshold) {
+        ret <- rep(FALSE, length(gammas_vec))
+      } else {
+        thecut.index <- max(which(fdr < threshold))
+        gammas_threshold <- sorted_gammas[thecut.index]
+        ret_vec <- (gammas_vec > gammas_threshold)
+      }
+
+      # reformat the results into a list
+      for (l in seq_len(length(x))) {
+        ret[[l]] <- array(FALSE, dim = dim(x[[l]])[-1])
+        len <- prod(dim(x[[l]])[-1])
+        ret[[l]] <- array(ret_vec[1:len], dim = dim(x[[l]])[-1])
+        ret_vec <- ret_vec[-c(1:len)]
+      }
     }
   }
 

R/coef.BayesSurvive.R

@@ -105,7 +105,7 @@ coef.BayesSurvive <- function(object, MPM = FALSE, type = "mean", CI = 95,
     )
     names(tbl)[2] <- type
     if (SD) tbl <- data.frame(tbl, SD = apply(beta_p, 2, sd))
-    #tbl$term <- factor(tbl$term, levels = tbl$term)
+    # tbl$term <- factor(tbl$term, levels = tbl$term)
   } else {
     # MPM coefficients
     if (object$input$S > 1 || !object$input$MRF.G) {
@@ -121,11 +121,10 @@ coef.BayesSurvive <- function(object, MPM = FALSE, type = "mean", CI = 95,
       term = x_names,
       estimate = object$output$beta.margin
     )
-    tbl$estimate <- tbl$estimate / object$output$gamma.margin * 
+    tbl$estimate <- tbl$estimate / object$output$gamma.margin *
       (object$output$gamma.margin >= 0.5)
     tbl$estimate[is.na(tbl$estimate)] <- 0
   }
 
   tbl
-
 }

R/plot.BayesSurvive.R

@@ -107,9 +107,8 @@ plot.BayesSurvive <- function(x, type = "mean", interval = TRUE,
   # pdf("psbcBeta.pdf", height = 5, width = 3.5)
 
   # Sys.setenv(`_R_S3_METHOD_REGISTRATION_NOTE_OVERWRITES_` = "false")
-  pCoef <- ggcoef(tbl, conf.int = interval, ...) + 
+  pCoef <- ggcoef(tbl, conf.int = interval, ...) +
     xlab(expression(Posterior ~ ~beta)) + ylab("")
   pCoef
   # dev.off()
-
 }

R/predict.BayesSurvive.R

@@ -28,10 +28,10 @@
 #' @param \dots not used
 #'
 #' @return A list object with 5 components if \code{type="brier"} including
-#'"model", "times", "Brier", "IBS" and "IPA" (Index of Prediction Accuracy), 
+#' "model", "times", "Brier", "IBS" and "IPA" (Index of Prediction Accuracy),
 #' otherwise  a list of 7 components with the first component as
-#' the specified argument \code{type} and "se", "band", "type", "diag", 
-#' "baseline" and "times", see function \code{riskRegression::predictCox} for 
+#' the specified argument \code{type} and "se", "band", "type", "diag",
+#' "baseline" and "times", see function \code{riskRegression::predictCox} for
 #' details
 #'
 #' @examples
@@ -215,7 +215,7 @@ predict.BayesSurvive <- function(object, survObj.new, type = "brier",
           metrics = "brier", summary = c("ibs", "ipa"),
           null.model = TRUE, times = times
         )$Brier$score
-        
+
         # Brier <- BrierScore[BrierScore$model != "Null model", ]
         # extract scores for Null model and the Bayesian Cox model
         ibs[1, ] <- Brier$Brier[c(length(times), length(times) * 2)]
@@ -251,7 +251,7 @@ predict.BayesSurvive <- function(object, survObj.new, type = "brier",
       # Brier scores of NULL.model do not change
       Brier.null <- BrierScore[1:(nrow(BrierScore) / 2), -c(1:2)]
       Brier <- BrierScore[-c(1:(nrow(BrierScore) / 2)), -c(1:2)]
-      
+
       # calculate Brier scores based on other MCMC estimates
       for (i in 2:nrow(betas)) {
         data_train$lp <- survObj$lp.all[, i] # lp_all_train[, i]
@@ -271,18 +271,20 @@ predict.BayesSurvive <- function(object, survObj.new, type = "brier",
         Brier <- Brier + BrierScore[-c(1:(nrow(BrierScore) / 2)), -c(1:2)]
       }
       Brier <- rbind(Brier.null, Brier / nrow(betas))
-      
+
       # extract IBS for Null model and the Bayesian Cox model
       ibs[1, ] <- Brier$Brier[c(length(times), length(times) * 2)]
       ibs[2, ] <- Brier$IBS[c(length(times), length(times) * 2)]
       ibs[3, ] <- Brier$IPA[c(length(times), length(times) * 2)]
     }
-    
-    #ibs <- rbind(rep(max(times), 2), ibs)
+
+    # ibs <- rbind(rep(max(times), 2), ibs)
     t0 <- round(max(times), digits = 1)
-    rownames(ibs) <- c(paste0("Brier(t=", t0, ")"), 
-                       paste0("IBS(t:0~", t0, ")"), 
-                       paste0("IPA(t=", t0, ")"))
+    rownames(ibs) <- c(
+      paste0("Brier(t=", t0, ")"),
+      paste0("IBS(t:0~", t0, ")"),
+      paste0("IPA(t=", t0, ")")
+    )
     if (verbose) {
       # cat("                      IBS\n",
       #   "  Null model          ", ibs[1, 1],
@@ -294,5 +296,4 @@ predict.BayesSurvive <- function(object, survObj.new, type = "brier",
     invisible(Brier)
     # return(Brier)
   }
-
 }