tests/testthat/test-nullsims.R

In segtest: Tests for Segregation Distortion in Polyploids

test_that("null simulations produce uniform p-values", {
  skip("not for unit testing")
  p1_ploidy <- 6
  p1 <- 0
  p2_ploidy <- 6
  p2 <- 2
  q <- gf_freq(
    p1_g = p1,
    p1_ploidy = p1_ploidy,
    p1_gamma = 1,
    p1_beta = NULL,
    p1_alpha = NULL,
    p1_type = "mix",
    p2_g = p2,
    p2_ploidy = p2_ploidy,
    p2_gamma= c(0.5, 0.5),
    p2_beta = NULL,
    p2_alpha = NULL,
    p2_type = "mix",
    pi = 0.03)
  niter <- 200
  nsamp <- 200

  pval_vec <- rep(NA_real_, times = niter)
  df_vec <- rep(NA_real_, times = niter)
  stat_vec <- rep(NA_real_, times = niter)
  df1_vec <- rep(NA_real_, times = niter)
  df0_vec <- rep(NA_real_, times = niter)
  pi_vec <- rep(NA_real_, times = niter)

  for (i in seq_len(niter)) {
    if (i %% 10 == 0) {
      cat(i, "out of", niter, "\n")
    }
    nvec <- c(stats::rmultinom(n = 1, size = nsamp, prob = q))
    system.time({
    sout <- seg_lrt(
      x = nvec,
      p1_ploidy = p1_ploidy,
      p2_ploidy = p2_ploidy,
      p1 = p1,
      p2 = p2,
      model = "seg",
      outlier = TRUE,
      opt = "bobyqa",
      optg = "NLOPT_GN_MLSL_LDS",
      chisq = FALSE)
    })
    # if (sout$p_value < 0.01) stop("here")
    pval_vec[[i]] <- sout$p_value
    df_vec[[i]] <- sout$df
    stat_vec[[i]] <- sout$stat
    df1_vec[[i]] <- sout$alt$df1
    df0_vec[[i]] <- sout$null$df0
    ## pi_vec[[i]] <- sout$null$gam[[3]]$pi
  }

  plot(
    x = ppoints(n = sum(!is.na(pval_vec))),
    y = sort(pval_vec[!is.na(pval_vec)]),
    xlab = "Theoretical",
    ylab = "Sample")
  abline(a = 0, b = 1, lty = 2, col = 2)

  ## Ideal circumstances is we know q
  pvec_chisq <- replicate(n = 1000, expr = {
    nvec <- c(stats::rmultinom(n = 1, size = nsamp, prob = q))
    suppressWarnings(
      stats::chisq.test(x = nvec, p = q)$p.value
    )
  })
  plot(
    x = ppoints(n = length(pvec_chisq)),
    y = sort(pvec_chisq),
    xlab = "Theoretical",
    ylab = "Sample")
  abline(a = 0, b = 1, lty = 2, col = 2)

  ## Train my eye
  plot(
    x = ppoints(n = niter),
    y = sort(stats::runif(niter)),
    xlab = "Theoretical",
    ylab = "Sample")
  abline(a = 0, b = 1, lty = 2, col = 2)
})

test_that("corner cases work", {
  sout <- seg_lrt(x = c(1, 0, 0, 0, 0), p1_ploidy = 4, p2_ploidy = 4, p1 = 0, p2 = 0, outlier = FALSE)
  expect_equal(sout$p_value, 1, tolerance = 1e-3)

  sout <- seg_lrt(x = c(0, 0, 1, 0, 0), p1_ploidy = 4, p2_ploidy = 4, p1 = 0, p2 = 4, outlier = FALSE)
  expect_equal(sout$p_value, 1, tolerance = 1e-3)

  sout <- seg_lrt(x = c(0, 0, 0, 0, 1), p1_ploidy = 4, p2_ploidy = 4, p1 = 4, p2 = 4, outlier = FALSE)
  expect_equal(sout$p_value, 1, tolerance = 1e-3)

  sout <- seg_lrt(x = c(0, 1, 0, 0, 0), p1_ploidy = 4, p2_ploidy = 4, p1 = 0, p2 = 0, outlier = FALSE)
  expect_equal(sout$p_value, 0, tolerance = 1e-3)

  sout <- seg_lrt(x = c(0, 0, 1, 0, 0), p1_ploidy = 4, p2_ploidy = 4, p1 = 0, p2 = 0, outlier = FALSE)
  expect_equal(sout$p_value, 0, tolerance = 1e-3)

  sout <- seg_lrt(x = c(0, 0, 0, 1, 0), p1_ploidy = 4, p2_ploidy = 4, p1 = 0, p2 = 0, outlier = FALSE)
  expect_equal(sout$p_value, 0, tolerance = 1e-3)

  sout <- seg_lrt(x = c(0, 0, 0, 0, 1), p1_ploidy = 4, p2_ploidy = 4, p1 = 0, p2 = 0, outlier = FALSE)
  expect_equal(sout$p_value, 0, tolerance = 1e-3)
})

test_that("skip", {
  skip("testing out frequencies")
  # library(numDeriv)
  # fn <- function(par) {
  #   gf_freq(
  #     p1_g = 3,
  #     p1_ploidy = 4,
  #     p1_gamma = 1,
  #     p1_beta = par[[1]],
  #     p1_alpha = NULL,
  #     p1_type = "mix",
  #     p2_g = 1,
  #     p2_ploidy = 4,
  #     p2_gamma= 1,
  #     p2_beta = par[[2]],
  #     p2_alpha = NULL,
  #     p2_type = "mix",
  #     pi = par[[3]])
  # }
  # q <- fn(par = c(0, 0, 0.015))
  # par <- c(beta_bounds(4), beta_bounds(4), 0.03)
  # svd(numDeriv::jacobian(func = fn, x = par/2))$d
  #
  # obj <- function(par) {
  #   sum((fn(par) - q)^2)
  # }
  # oout <- optim(par = par, fn = obj, lower = c(0, 0, 0), upper = c(0.1, 0.1, 0.01), method = "L-BFGS-B")
  #
  # oout$par
  # fn(oout$par) - q
  # q
})

test_that("seg at simplex and auto_dr are same", {
  q1 <- gamfreq(g = 1, ploidy = 4, gamma = 1, alpha = NULL, beta = 0.01, type = "mix", add_dr = TRUE)
  q2 <- gamfreq(g = 1, ploidy = 4, gamma = NULL, alpha = 0.04, beta = NULL, type = "polysomic", add_dr = TRUE)
  expect_equal(q1, q2)

  q1 <- gamfreq(g = 3, ploidy = 4, gamma = 1, alpha = NULL, beta = 0.01, type = "mix", add_dr = TRUE)
  q2 <- gamfreq(g = 3, ploidy = 4, gamma = NULL, alpha = 0.04, beta = NULL, type = "polysomic", add_dr = TRUE)
  expect_equal(q1, q2)

  q1 <- gf_freq(p1_g = 1, p1_ploidy = 4, p1_gamma = 1, p1_alpha = NULL, p1_beta = 0.01, p1_type = "mix", p1_add_dr = TRUE, p2_g = 3, p2_ploidy = 4, p2_gamma = 1, p2_alpha = NULL, p2_beta = 0.01, p2_type = "mix", p2_add_dr = TRUE, pi = 0.01)
  q2 <- gf_freq(p1_g = 1, p1_ploidy = 4, p1_gamma = NULL, p1_alpha = 0.04, p1_beta = NULL, p1_type = "polysomic", p1_add_dr = TRUE, p2_g = 3, p2_ploidy = 4, p2_gamma = NULL, p2_alpha = 0.04, p2_beta = NULL, p2_type = "polysomic", p2_add_dr = TRUE, pi = 0.01)
  expect_equal(q1, q2)

  q1 <- gf_freq(p1_g = 3, p1_ploidy = 4, p1_gamma = 1, p1_alpha = NULL, p1_beta = 0.01, p1_type = "mix", p1_add_dr = TRUE, p2_g = 1, p2_ploidy = 4, p2_gamma = 1, p2_alpha = NULL, p2_beta = 0.01, p2_type = "mix", p2_add_dr = TRUE, pi = 0.01)
  q2 <- gf_freq(p1_g = 3, p1_ploidy = 4, p1_gamma = NULL, p1_alpha = 0.04, p1_beta = NULL, p1_type = "polysomic", p1_add_dr = TRUE, p2_g = 1, p2_ploidy = 4, p2_gamma = NULL, p2_alpha = 0.04, p2_beta = NULL, p2_type = "polysomic", p2_add_dr = TRUE, pi = 0.01)
  expect_equal(q1, q2)
})

test_that("auto_dr and seg are same", {
  set.seed(1)
  nvec <- c(32L, 2431L, 4979L, 2520L, 38L)
  sout1 <- seg_lrt(
    x = nvec,
    p1_ploidy = 4,
    p2_ploidy = 4,
    p1 = 3,
    p2 = 1,
    model = "auto_dr",
    outlier = TRUE)
  sout2 <- seg_lrt(
    x = nvec,
    p1_ploidy = 4,
    p2_ploidy = 4,
    p1 = 3,
    p2 = 1,
    model = "seg",
    outlier = TRUE)
  expect_equal(sout1$null$l0, sout2$null$l0)
})