R/utils.R

Defines functions AT2_one_perm col_t_test_fun pirwin.hall

Documented in pirwin.hall

#' Distribution function of the Irwin-Hall distribution
#' @description Computes the distribution function of Irwin-Hall with parameter n. This function is NOT vectorized.
#' @param x A number between 0 and n.
#' @param n The number of i.i.d. U(0,1) to sum. It should be an integer.
#'
#' @return The CDF of the Irwin-Hall(n) at x.
#' @export
#' @references
#' Irwin-Hall Distribution. (n.d.). Randomservices.org. Retrieved May 8, 2023, from https://www.randomservices.org/random/special/IrwinHall.html
#' @examples
#' pirwin.hall(1, 2)
#' pirwin.hall(2.5, 4)

pirwin.hall <- function(x, n) {
  if (any(!is.numeric(x) | !is.numeric(n) | (length(x) != 1) | (length(n) != 1) | (x < 0) | (x > n) | (n%%1 != 0)))
    stop("Invalid inputs.")

  k <- 0:n
  return ( 1/2 + 1/(2*factorial(n)) * sum( (-1)^k * choose(n, k) * sign(x-k) * (x-k)^n ) )
}


col_t_test_fun <- function(col_data, label) {
  # conducts a two sample t-test on a data vector [col_data] with assigned label (2 factors). Return the t-statistic.
  t.test(formula = col_data ~ as.factor(label),
         alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE)$statistic
}


AT2_one_perm <- function(perm_idx, edge_t_stat_selected, edge_dist_mat) {
  # Compute AT2 with one random permutation of the label.
  perm_label <- sample(edge_dist_mat[,1])  # recall that the first column is sample label

  t_stat_perm <- apply(edge_dist_mat[,-1], 2, col_t_test_fun, label = perm_label)  # t-statistic for ALL edges with permuted label
  t_perc_perm <- ecdf(t_stat_perm)(t_stat_perm)  # percentile of each t-statistic among ALL t-statistic

  names(t_stat_perm) <- names(t_perc_perm) <- colnames(edge_dist_mat[,-1])
  t_perc_perm_selected <- t_perc_perm[paste(edge_t_stat_selected$V1, edge_t_stat_selected$V2, sep = "-")]

  AT2 <- mean(abs(t_perc_perm_selected - 0.5))  # Compute AT2 associated with the permuted label

  return (AT2)
}

Try the ExprNet package in your browser

Any scripts or data that you put into this service are public.

ExprNet documentation built on Aug. 15, 2023, 9:08 a.m.