tests/testthat/test_chisq.R

In miheerdew/bmdupdate: This package implements the pvalue calculation for the ChiSq test

source("helper.R")

context("Check Chisq Backend")

set.seed(1234567)
sim1 <- generate_data(200, 30, rho, 100, beta, s2)
sim2 <- generate_data(200, 200, rho, 50, beta, s2)

b1 <- new(BackendChisq, sim1$X, sim1$Y)
b2 <- new(BackendChisq, sim2$X, sim2$Y)

check <- function(sim, A, backend, testX = TRUE, method=Method$Auto, mask=c(), total_mask=mask){
  if (!testX) {
    #Testing Y set
    c <-  precomputations(sim$Y[,A, drop=FALSE], sim$X)
    ny <- sim$dx
    A <- A + sim$dx;
    currentMask <- total_mask[total_mask <= sim$dx]
  } else {
    c <- precomputations(sim$X[,A, drop=FALSE], sim$Y)
    ny <- sim$dy
    currentMask <- total_mask[total_mask > sim$dx] - sim$dx
  }

  ## Compute R pvalues from traces
  if(method == Method$BigX_fast) { #method == Method$BigX_fast) {
   traceR <- t(sapply(1:ny, function(i) traceXR_fast(i,c)))
  } else {
   traceR <- t(sapply(1:ny, function(i) traceXR(i,c)))
  }
  Tstat <- rowSums(c$allr^2)
  as <- traceR[,2] / traceR[,1]
  bs <- traceR[,1]^2 / traceR[, 2]
  pvalsR <- pchisq(c$ndata*Tstat/as, df=bs, lower.tail = FALSE)

  if(length(total_mask) > 0) {
    expect_true(setequal(backend$mask(mask),total_mask))
    pvalsR[currentMask] <- NA
  }

  pvalsC <- as.vector(backend$pvals(A,FALSE, method))
  pvalsCpar <- as.vector(backend$pvals(A, TRUE, method))

  expect_equal(pvalsCpar, pvalsR)
  expect_equal(pvalsC, pvalsR)
}

test_that("Chisq calculates trace correctly for X set",{
  check(sim1, 1:sim1$dx, b1)
  check(sim2, 1:sim2$dx, b2)
  check(sim1, 1:sim1$dx, b1, Method$BigX_fast)
  check(sim2, 1:sim2$dx, b2, Method$BigX_fast)
})

test_that("Chisq calculates trace correctly for X subset",{
  check(sim1, 1:(sim1$dx/2), b1)
  check(sim2, 1:(sim2$dx/2), b2)
  check(sim1, 1:(sim1$dx/2), b1, Method$BigX_fast)
  check(sim2, 1:(sim2$dx/2), b2, Method$BigX_fast)
})

test_that("Chisq calculates trace correctly for Y subset",{
  check(sim1, 1:(sim1$dy/2), b1, testX = FALSE)
  check(sim2, 1:(sim2$dy/2), b2, testX = FALSE)
  check(sim1, 1:(sim1$dy/2), b1, testX = FALSE, Method$BigX_fast)
  check(sim2, 1:(sim2$dy/2), b2, testX = FALSE, Method$BigX_fast)
})

test_that("Chisq calculates trace correctly entire Y subset",{
  check(sim1, 1:(sim1$dy), b1, testX = FALSE)
  check(sim2, 1:(sim2$dy), b2, testX = FALSE)
  check(sim1, 1:(sim1$dy), b1, testX = FALSE, Method$BigX_fast)
  check(sim2, 1:(sim2$dy), b2, testX = FALSE, Method$BigX_fast)
})

test_that("Last singleton X passes", {
  check(sim1, sim1$dx, b1, testX = TRUE)
  check(sim2, sim2$dx, b2, testX = TRUE)
  check(sim1, sim1$dx, b1, testX = TRUE, Method$BigX_fast)
  check(sim2, sim2$dx, b2, testX = TRUE, Method$BigX_fast)
})

test_that("Single Mask/Unmask works", {
  check(sim1, sim1$dx, b1, testX = TRUE, mask = sim1$dx + c(29, 3))
  check(sim2, 1:20, b2, testX = FALSE, mask = 1:100)
  b1$unmask()
  expect_length(b1$mask(numeric(0)), 0)
  b2$unmask()
  expect_length(b2$mask(numeric(0)),0)
})

test_that("Multiple Mask/Unmask work", {
  total <- c()

  mask <- sample(1:7) + sim1$dx
  total <- c(total, mask)
  check(sim1, 1:sim1$dx, b1, testX = TRUE, mask = mask, total=total)

  mask <-  sample(1:sim1$dx, 5)
  total <- c(total, mask)
  check(sim1, 1:sim1$dy, b1, testX = FALSE, mask = mask, total=total)

  mask <-  c(sample(1:sim1$dx, 30), sample(1:sim1$dy, 5) + sim1$dx)
  total <- c(total, mask)
  check(sim1, 1:sim1$dy, b1, testX = FALSE, mask = mask, total=total)

  b1$unmask()
  expect_length(b1$mask(numeric(0)), 0)
  total <- c()

  mask <-  c(sample(1:sim1$dx, 30), sample(1:sim1$dy, 5) + sim1$dx)
  total <- c(total, mask)
  check(sim1, 1:sim1$dy, b1, testX = FALSE, mask = mask, total=total)
})