R/fastPermAsymptoticFunctions.R

In bdsegal/fastperm: Fast approximation of small p-values in permutation tests by partitioning the permutation space

xiFunRatioMean1 <- function(m, nx, ny, x, y){
  #' Utitlity function for fastPermAsym
  #'
  #' Calculates xi(m)
  #' @export
  
  N <- nx + ny
  # lambda <- nx / N
  
  xbar <- mean(x)
  ybar <- mean(y)
  
  Emean <- g(m, nx, ny, xbar, ybar)
  Evar <- VarR(m, nx, ny, x, y, xbar, ybar)
  
  t <- max(xbar / ybar, ybar / xbar)

  xi <- (t - Emean) / sqrt(Evar)
  
  return(xi)
}

xiFunDiffMean1 <- function(m, nx, ny, x, y){
  #' Utitlity function for fastPermAsym
  #'
  #' Calculates xi(m)
  #' @export
  
  N <- nx + ny
  # lambda <- nx / N
  
  xbar <- mean(x)
  ybar <- mean(y)
  
  Emean <- gLin(m = m, nx, ny, xbar = xbar, ybar = ybar)
    
  Evar <-  VarRLin(m, nx, ny, x, y)
  
  t <- abs(xbar - ybar)

  xi <- (t - Emean) / sqrt(Evar)
  
  return(xi)
}

fastPermAsym <- function(x, y, testStat = ratioMean){
#' Fast approximation of small permutation p-values via asymptotic results
#'
#' This function approximates the p-value (two-sided) for a two
#' sample permutation test by partitioning the permutation space,
#' and using asymptotic results on the p-values within each partition.
#' This function is most useful for small p-values, e.g. p < 10^-6, and
#' large sample sizes.
#' @param x First sample
#' @param y Second sample
#' @param testStat Test statistic, defaults to the ratio of the
#' means (ratioMean). The difference in means (diffMean) is also available.
#' @param B Number of Monte Carlo iterations within each partition. Defaults to 1,000.
#' @keywords fast permtution test two sample
#' @export
#' @examples
#' x <- rexp(100, 5)
#' y <- rexp(100, 2)
#' mStopRatioMean(x, y)
#' fastPerm(x, y)
#'
#' x <- rnorm(110, 0, 1)
#' y <- rnorm(110, 1, 1)
#' mStopDiffMean(x, y)
#' fastPerm(x, y, testStat = diffMean)
  
  comparison <- attributes(testStat)$comparison  
  if (comparison == "ratio" & (min(x,y)<0)) {
    stop("Some sample values <=0. With ratio statistics, all values must be greater than 0.")
  }
  
  if(comparison == "ratio"){
    xiFun <- xiFunRatioMean1
  } else if (comparison == "difference"){
    xiFun <- xiFunDiffMean1
  } else {
    stop("testStat must be either ratioMean or diffMean")
  }

  # swap x and y if xbar < ybar
  # if (mean(x) < mean(y)) {
  #   ytemp <- y
  #   y <- x
  #   x <- ytemp
  # }
  
  nx <- length(x)
  ny <- length(y)
  nMax <- min(nx, ny)
  
  t0 <- testStat(x,y)
  
  pmf <- PiLgamma(nx, ny)
  mMax <- min(which(pmf == max(pmf)))

  if (length(mMax)==1 & mMax[1]==nMax){
    m <- 1:mMax
  } else if (length(mMax)==1 & mMax[1]<nMax){
    m <- c(1:mMax,((mMax-1):0)[1:(nMax-mMax)])
  } else if (length(mMax)==2 & mMax[2]==nMax){
    m <- c(1:mMax[1], mMax[1])
  } else {
    m <- c(1:mMax[1], (mMax[1]:0)[1:(nMax - mMax[1])])
  }
  
  # xi <- xiFun(m=m, nx, ny, x, y)
  xi1 <- xiFun(m=m, nx, ny, x, y)
  xi2 <- xiFun(m=m, ny, nx, y, x)

  # if nx != ny, set p-value in 0 partition to 1 and 
  # estimate p-value in partition min(nx,ny)
  # note: for xi2, when we switch x and y, t become 1/t.
  # Consequently, we need the lower tail for xi2
  if (nx != ny) { 
    # pNorm <- c(1, pnorm(xi, lower.tail=FALSE)) %*% pmf
    xiVec <- c(1, pnorm(xi1, lower.tail=FALSE) + pnorm(xi2, lower.tail = FALSE))
    pNorm <- xiVec %*% pmf
    pT <- xiVec %*% pmf
  # if nx==ny, set both the 0 and nx partition to 1
  } else { 
    # pNorm <- c(1, pnorm(xi[-length(xi)], lower.tail=FALSE), 1) %*% pmf
    xiVec <- c(1, pnorm(xi1[-length(xi1)], lower.tail=FALSE) + 
                  pnorm(xi2[-length(xi2)], lower.tail=FALSE), 1)
    pNorm <- xiVec %*% pmf
    pT <- xiVec %*% pmf
  }

  # if (plot) {
    # par(mar = c(5, 5, 4, 2) + 0.1)
    # plot(x = m, y = xi, pch=19, ylab = expression(paste(xi, "(m)")), 
         # main= bquote(paste(m["stop"]^"asym", " = ", .(mStop), sep = "")),
         # cex.axis = 1.4, cex.lab = 1.4, cex.main = 1.4)
    # abline(a = u, b = 0, col="red")
  # }
  fpAsym <- list(pNorm=as.numeric(pNorm),
                 pT=as.numeric(pT),
                 comparison = attributes(testStat)$comparison,
                 summary = attributes(testStat)$summary,
                 t0=t0)

  class(fpAsym) <- "fastPermAsym"
  return(fpAsym)
}

print.fastPermAsym <- function(fp){
  #' Print function for fastPerm
  #'
  #' This function prints the results of fastPermAsym
  #' @param fp Output from the fastPermAsym function
  #' @keywords fastPerm print
  #' @export
  #' @examples
  #' x <- rexp(100, 5)
  #' y <- rexp(100, 2)
  #' print(fastPermAsym(x, y, testStat = ratioMean))
  
  result <- paste("    fastPerm Two Sample Test   \n\n",
                  fp$comparison, " of ", fp$summary, "s\n", 
                  "\nobserved statistic = ", signif(fp$t0,3),
                  "\n\np-value \t Approximation", 
                  "\n-------------------------------",
                  "\n", signif(fp$pT,3), "\t t-distribution",
                  "\n", signif(fp$pNorm,3), "\t normal",
                  sep = "")
  
  writeLines(result)
}

# testing ---------------------------------------------------------------------
# x <- rexp(n = 100, rate = 2)
# y <- rexp(n = 100, rate = 4)
# mStopRatioMean(x, y, B=1000, plot = TRUE)
# fastPermAsym(x, y, testStat = ratioMean)
# fastPerm(x, y, testStat = ratioMean)