Merge pull request #49 from suzannejin/master

change 'tails' argument to take 'right' and 'both'

R/3-shared-getAdjacency.R

@@ -1,18 +1,25 @@
-#' Get Adjacency Matrix as indicated by permutation tests
+#' Get Adjacency Matrix as indicated by permutation tests.
+#' 
+#' This function gets the significant pairs according to the permutation tests. Then it fills 
+#' the adjacency matrix with 1 if pair is significant, otherwise 0. The significance is determined
+#' by the cutoff value for which the false discovery rate (FDR) is less or equal than the given 
+#' value 'fdr'. The significant pairs are those that have a value greater/less or equal than the 
+#' cutoff, depending on the metric.
 #'
-#' This function gets the significant pairs, according to the permutation tests. 
-#' Then it fills the adjacency matrix with 1 if pair is significant, otherwise 0.
 #' @param object A \code{propd} or \code{propr} object.
 #' @param fdr A float value for the false discovery rate. Default is 0.05.
 #' @param window_size An integer. Default is 1. When it is greater than 1, the FDR
 #' values would be smoothed out by a moving average of the given window size.
-#' @param tails 1 for one-sided on the right, and 2 for two-sided. When NULL, use default
-#' test for the given metric. This is only relevant for \code{propr} objects.
+#' @param tails 'right' or 'both'. 'right' is for one-sided on the right. 'both' is to 
+#' combine one-sided on the right (positive values) and left (negative values). This is only 
+#' relevant for \code{propr} objects, as \code{propd} objects are always one-sided and only 
+#' have positive values.
 #' @return An adjacency matrix.
 #'
 #' @export
-getAdjacencyFDR <- function(object, fdr = 0.05, window_size = 1, tails = NULL) {
-
+getAdjacencyFDR <- 
+  function(object, fdr = 0.05, window_size = 1, tails = c('right', 'both')) {
+
   if (inherits(object, "propr")){
     adj <- getAdjacencyFDR.propr(object, fdr=fdr, window_size=window_size, tails=tails)
 
@@ -60,23 +67,18 @@ getAdjacencyFDR.propd <-
 #' pairs from a \code{propd} object.
 #' @export
 getAdjacencyFDR.propr <-
-  function(object, fdr = 0.05, window_size = 1, tails = NULL) {
+  function(object, fdr = 0.05, window_size = 1, tails = c('right', 'both')) {
 
-    if (is.null(tails)) {
-      tails <- ifelse(object@has_meaningful_negative_values, 2, 1)
-    }
-    if (tails != 1 & tails != 2) {
-        stop("Please provide a valid value for tails: 1 or 2.")
-    }
-    if (tails == 1 & object@has_meaningful_negative_values) {
-      warning("Significant pairs are chosen based on one-sided FDR test.")
-    }
-    if (tails == 2 & !object@direct) {
-      stop("Two-sided FDR is not available for this metric.")
+    # handle tails argument
+    tails <- match.arg(tails)
+    if (tails == 'both' & !object@direct) {
+      warning("Running tails='right' instead")  
+      tails <- 'right'   # set to right, when non negative values can be expected for a given metric.
     }
 
-    # correct fdr when two-sided
-    if (tails == 2) fdr <- fdr / 2
+    # correct fdr when considering both tails
+    # so that the total fdr is the given value
+    if (tails == 'both') fdr <- fdr / 2
 
     # create empty matrix
     adj <- matrix(0, nrow = ncol(object@matrix), ncol = ncol(object@matrix))
@@ -90,12 +92,12 @@ getAdjacencyFDR.propr <-
       if (object@direct) {
         adj[object@matrix >= 0 & object@matrix >= cutoff] <- 1
       } else {
-        adj[object@matrix <= cutoff] <- 1
+        adj[object@matrix <= cutoff] <- 1  # don't have negative values
       }
     }
 
     # for negative tail
-    if (tails == 2){
+    if (tails == 'both'){
       cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size, positive=F)
       if (cutoff) adj[object@matrix < 0 & object@matrix <= cutoff] <- 1
     }

R/3-shared-getCutoff.R

@@ -61,7 +61,7 @@ getCutoffFDR.propr <- function(object, fdr = 0.05, window_size = 1, positive = T
   # get index of FDR values below the threshold
   index <- (df$FDR <= fdr) & (is.finite(df$FDR))
   if (!any(index)) {
-    warning("No significant cutoff found for the given FDR, when positive = ", positive)
+    warning("No significant cutoff found for the given FDR, when positive tail = ", positive)
     return(FALSE)
   }
 
@@ -130,7 +130,7 @@ getCutoffFstat <- function(object, pval = 0.05, fdr_adjusted = FALSE) {
 
   if (fdr_adjusted) {
     message("Alert: Returning an empiric cutoff based on the $FDR slot.")
-    index <- (object@results$FDR < pval) & (is.finite(object@results$FDR))
+    index <- (object@results$FDR <= pval) & (is.finite(object@results$FDR))
     if (any(index)) {
       cutoff <- max(object@results$theta[index])
     } else{

R/3-shared-getResults.R

@@ -19,16 +19,17 @@ getResults <-
 
 #' Get Significant Results from Object based on the permutation tests.
 #'
-#' This function provides a unified wrapper to retrieve results
-#'  from a \code{propr} or \code{propd} object keeping only the
-#'  statistically significant pairs.
+#' This function retrieves results from a \code{propr} or \code{propd} object keeping only the
+#' statistically significant pairs. The significance is determined by the cutoff value for which
+#' the false discovery rate (FDR) is less or equal than the given value 'fdr'. The significant 
+#' pairs are those that have a value greater/less or equal than the cutoff, depending on the metric.
 #'
 #' @inheritParams getAdjacencyFDR
 #' @return A \code{data.frame} of results.
 #'
 #' @export
 getSignificantResultsFDR <-
-  function(object, fdr = 0.05, window_size = 1, tails = NULL) {
+  function(object, fdr = 0.05, window_size = 1, tails = c('right', 'both')) {
 
     if (inherits(object, "propr")) {
       results <- getSignificantResultsFDR.propr(object, fdr=fdr, window_size=window_size, tails=tails)
@@ -68,19 +69,13 @@ getSignificantResultsFDR.propd <-
 #' only the statistically significant pairs.
 #' @export
 getSignificantResultsFDR.propr <- 
-  function(object, fdr = 0.05, window_size = 1, tails = NULL) {
+  function(object, fdr = 0.05, window_size = 1, tails = c('right', 'both')) {
 
-    if (is.null(tails)) {
-      tails <- ifelse(object@has_meaningful_negative_values, 2, 1)
-    }
-    if (tails != 1 & tails != 2) {
-        stop("Please provide a valid value for tails: 1 or 2.")
-    }
-    if (tails == 1 & object@has_meaningful_negative_values) {
-      warning("Significant pairs are chosen based on one-sided FDR test.")
-    }
-    if (tails == 2 & !object@direct) {
-      stop("Two-sided FDR is not available for this metric.")
+    # handle tails argument
+    tails <- match.arg(tails)
+    if (tails == 'both' & !object@direct) {
+      warning("Running tails='right' instead")  
+      tails <- 'right'   # set to right, when non negative values can be expected for a given metric.
     }
 
     # function to subset the results data frame based on the cutoff
@@ -93,8 +88,9 @@ getSignificantResultsFDR.propr <-
       }
     }
 
-    # correct fdr when two-sided
-    if (tails == 2) fdr <- fdr / 2
+    # correct fdr when considering both tails
+    # so that the total fdr is the given value
+    if (tails == 'both') fdr <- fdr / 2
 
     # define results data frame
     df <- getResults(object)
@@ -107,7 +103,7 @@ getSignificantResultsFDR.propr <-
     results <- rbind(results, subsetBeyondCutoff(part, cutoff))
 
     # get the significant negative values
-    if (tails == 2) {
+    if (tails == 'both') {
       cutoff <- getCutoffFDR(object, fdr=fdr, window_size=window_size, positive=F)
       part <- df[which(df$propr < 0),]
       results <- rbind(results, subsetBeyondCutoff(part, cutoff))

R/3-shared-updateCutoffs.R

@@ -1,50 +1,58 @@
 #' Update FDR by Permutation
 #'
-#' This function updates FDR for a set of cutoffs.
-#'
-#' This function wraps \code{updateCutoffs.propr} and
-#'  \code{updateCutoffs.propd}.
+#' This function updates the FDR for each cutoff. By default, the set of cutoffs are determined
+#' based on the quantile of the data, so that the cutoffs are evenly spaced across the data.
+#' The FDR is calculated as the ratio between the number of permuted values beyond the cutoff
+#' and the number of true values beyond the the cutoff. 
+#' When tails is set to 'right', the FDR is calculated only on the positive side of the data.
+#' When tails is set to 'both', the FDR is calculated on both sides of the data.
 #'
 #' @param object A \code{propr} or \code{propd} object.
-#' @param number_of_cutoffs An integer. The number of cutoffs to test.
-#' Given this number, the cutoffs will be determined based on the quantile of the data.
-#' In this way, the cutoffs will be evenly spaced across the data.
-#' @param custom_cutoffs A numeric vector. When provided, this vector is used
-#' as the FDR cutoffs to test, and number_of_cutoffs is ignored.
-#' @param tails 1 for one-sided on the right, and 2 for two-sided. When two-sided,
-#' the FDR is calculated for both positive and negative values. When NULL, use default
-#' option for the given metric. This is only relevant for \code{propr} objects.
+#' @param number_of_cutoffs An integer. The number of cutoffs to test. Given this number, 
+#' the cutoffs will be determined based on the quantile of the data. In this way, the 
+#' cutoffs will be evenly spaced across the data.
+#' @param custom_cutoffs A numeric vector. When provided, this vector is used as the set of 
+#' cutoffs to test, and 'number_of_cutoffs' is ignored.
+#' @param tails 'right' or 'both'. 'right' is for one-sided on the right. 'both' is to 
+#' combine one-sided on the right (positive values) and left (negative values). This 
+#' is only relevant for \code{propr} objects, as \code{propd} objects are always one-sided
+#' and only have positive values.
 #' @param ncores An integer. The number of parallel cores to use.
 #' @return A \code{propr} or \code{propd} object with the FDR slot updated.
+#' 
 #' @export
 updateCutoffs <-
   function(object,
            number_of_cutoffs = 100,
            custom_cutoffs = NULL,
-           tails = NULL,
+           tails = c('right','both'),
            ncores = 1) {
+    tails <- match.arg(tails)
 
-    if (!is.null(tails)){
-      if (tails != 1 & tails != 2) stop("Tails only accept values 1 or 2, if not NULL.")
-    }
-
-    get_cutoffs <- function(values, number_of_cutoffs, custom_cutoffs, tails=1) {
+    # function to get the set of cutoffs
+    get_cutoffs <- function(values, number_of_cutoffs, custom_cutoffs, tails='right') {
       if (!is.null(custom_cutoffs)) return(custom_cutoffs)
-      if (tails == 1) values <- values[values >= 0]
+      if (tails == 'right') {
+        values <- values[values >= 0]
+        if (length(values) == 0) stop("No positive values found.")
+      }
       return(as.numeric(quantile(values, probs = seq(0, 1, length.out = number_of_cutoffs))))
     }
 
+    # update FDR values for propr object
     if (inherits(object, "propr")) {
-      if (is.null(tails)) tails <- ifelse(object@has_meaningful_negative_values, 2, 1)
-      if (tails == 1 & object@has_meaningful_negative_values) warning("One-sided FDR test is performed.")
-      if (tails == 2 & !object@direct) stop("Two-sided FDR is not available for this metric.")
+      if (tails == 'both' & !object@direct) {
+        warning("Running tails='right' instead")  
+        tails <- 'right'   # set to right, when non negative values can be expected for a given metric.
+      }
       values <- object@results$propr
       cutoffs <- get_cutoffs(values, number_of_cutoffs, custom_cutoffs, tails)
       updateCutoffs.propr(object, cutoffs, ncores)
 
+    # update FDR values for propd object
     } else if (inherits(object, "propd")) {
       values <- object@results$theta
-      cutoffs <- get_cutoffs(values, number_of_cutoffs, custom_cutoffs, tails=1)
+      cutoffs <- get_cutoffs(values, number_of_cutoffs, custom_cutoffs, tails='right')
       updateCutoffs.propd(object, cutoffs, ncores)
 
     } else{
@@ -64,7 +72,6 @@ updateCutoffs <-
 #' @export
 updateCutoffs.propr <-
   function(object, cutoffs, ncores) {
-
     if (identical(object@permutes, list(NULL))) {
       stop("Permutation testing is disabled.")
     }

README.Rmd

@@ -54,7 +54,7 @@ pr <- updateCutoffs(
         pr,
         number_of_cutoffs = 100,  # number of cutoffs to estimate FDR
         custom_cutoffs = NULL,  # or specify custom cutoffs
-        tails = 1,  # 1 or 2-tailed test
+        tails = 'right',  # consider only the positive values ('right') or both sides ('both')
         ncores = 4  # parallelize here
       ) 
 ```
@@ -81,7 +81,7 @@ pr <- updateCutoffs(
         pr,
         number_of_cutoffs = 100,  # number of cutoffs to estimate FDR
         custom_cutoffs = NULL,  # or specify custom cutoffs
-        tails = 1,  # 1 or 2-tailed test
+        tails = 'right',  # consider only the positive values ('right') or both sides ('both')
         ncores = 4  # parallelize here
       ) 
 ```

README.md

@@ -101,7 +101,7 @@ pr <- updateCutoffs(
         pr,
         number_of_cutoffs = 100,  # number of cutoffs to estimate FDR
         custom_cutoffs = NULL,  # or specify custom cutoffs
-        tails = 1,  # 1 or 2-tailed test
+        tails = 'right',  # consider only the positive values ('right') or both sides ('both')
         ncores = 4  # parallelize here
       ) 
 ```
@@ -131,7 +131,7 @@ pr <- updateCutoffs(
         pr,
         number_of_cutoffs = 100,  # number of cutoffs to estimate FDR
         custom_cutoffs = NULL,  # or specify custom cutoffs
-        tails = 1,  # 1 or 2-tailed test
+        tails = 'right',  # consider only the positive values ('right') or both sides ('both')
         ncores = 4  # parallelize here
       ) 
 ```