tests/testthat/test-round-robin.R

In poppr: Genetic Analysis of Populations with Mixed Reproduction

context("Round Robin Tests")

x <- structure(
  c("B", "C", "B", "C", "B", "C", "A", "C", "A", "A", "C", "B", "C", "A", "A"),
  .Dim = c(5L, 3L)
)
suppressWarnings(x <- df2genind(x, ploidy = 1))
rrx_m <- rrmlg(x)
mlg_truth <- structure(
  c(3L, 2L, 3L, 1L, 1L, 2L, 4L, 2L, 3L, 1L, 2L, 3L, 2L, 3L, 1L),
  .Dim = c(5L, 3L),
  .Dimnames = list(NULL, c("L1", "L2", "L3"))
)
rrx_f <- rraf(x)
freq_truth <- structure(
  list(
    L1 = structure(
      c(0.333333333333333, 0.666666666666667),
      .Names = c("L1.B", "L1.C")
    ),
    L2 = structure(c(0.25, 0.75), .Names = c("L2.C", "L2.A")),
    L3 = structure(
      c(0.333333333333333, 0.333333333333333, 0.333333333333333),
      .Names = c("L3.C", "L3.B", "L3.A")
    )
  ),
  .Names = c("L1", "L2", "L3")
)

test_that("rrmlg produces the correct mlgs", {
  expect_equivalent(rrx_m, mlg_truth)
})

test_that("rrmlg will not work on genlight objects", {
  skip_on_cran()
  expect_error(rrmlg(glSim(10, 10, 10, parallel = FALSE)))
})

test_that("rrmlg can take loci objects", {
  skip_on_cran()
  expect_equivalent(rrmlg(pegas::as.loci(x)), rrx_m)
})

test_that("rraf produces correct allele frequencies", {
  skip_on_cran()
  expect_equivalent(rrx_f, freq_truth)
  expect_equivalent(rraf(x, res = "vector"), unlist(freq_truth))
})

test_that("rare_allele_correction gets angry if called from global environment", {
  skip_on_cran()
  expect_warning(do.call(
    poppr:::rare_allele_correction,
    args = list(rrx_f, rrx_m),
    envir = .GlobalEnv
  ))
})

test_that("correction is properly applied in rraf", {
  skip_on_cran()
  data(monpop)
  mll(monpop) <- "original"
  monrrmlg <- 1 / colSums(!apply(rrmlg(monpop), 2, duplicated))

  monc <- rraf(monpop) # default
  monnc <- rraf(monpop, correction = FALSE) # No correction
  monc_sto <- rraf(monpop, sum_to_one = TRUE) # summed to one
  monc_mlg <- rraf(monpop, d = "mlg") # by multilocus genotype
  monc_rrmlg <- rraf(monpop, d = "rrmlg") # by round-robin multilocus genotype
  monc_m <- rraf(monpop, m = 0.5) # assume diploid (not true here, though)
  monc_e <- rraf(monpop, e = pi) # A ridiculous correction
  monc_e2 <- rraf(monpop, e = pi, d = "rrmlg", m = 0.5) # e overrides

  SER <- function(x) x$SER["SER.147"]
  SED <- function(x) x$SED["SED.187"]

  # Correction can be turned off and on
  monc_vec <- vapply(monc, sum, numeric(1))
  monnc_vec <- vapply(monnc, sum, numeric(1))
  expect_gt(sum(monc_vec), sum(monnc_vec))
  expect_equivalent(sum(monnc_vec), nLoc(monpop))

  # sum_to_one argument augments does what it says
  monc_sum2one <- vapply(monc_sto, sum, numeric(1))
  # If the sum to one works, then we should have leftover tolerance that is
  # greater than zero but less than one.
  tol <- sum(monc_vec - monc_sum2one)
  expect_lt(tol, 1)
  expect_gt(tol, 0)
  expect_equal(sum(monc_sum2one), nLoc(monpop))
  # The default is 1/n
  expect_equivalent(SER(monc), 1 / nInd(monpop))
  expect_true(identical(SER(monc)[[1]], SED(monc)[[1]]))

  # The multiplier works
  expect_equivalent(SER(monc) / 2, SER(monc_m))

  # Works by MLG
  expect_equivalent(SER(monc_mlg), 1 / nmll(monpop))

  # Works by rrMLG
  expect_equivalent(SER(monc_rrmlg), monrrmlg["SER"])
  # We would not expect cross-loci corrections to be equal
  expect_false(identical(SER(monc_rrmlg)[[1]], SED(monc_rrmlg)[[1]]))

  # Works by custom value
  expect_equivalent(SER(monc_e), pi)

  # Custom value overrides anything else
  expect_equivalent(SED(monc_e), SED(monc_e2))
})

test_that("rraf produces a data frame", {
  skip_on_cran()
  rrx_df <- rraf(x, res = "data.frame")
  expect_is(rrx_df, "data.frame")
  expect_equivalent(names(rrx_df), c("frequency", "locus", "allele"))
  # The length of the data should be equql
  expect_equal(nrow(rrx_df), sum(lengths(rrx_f)))
  # The content should be equal
  expect_equal(rrx_df$frequency, unlist(rrx_f, use.names = FALSE))
})

test_that("rraf calculates per population", {
  skip_on_cran()
  data(Pram)
  rrx_matrix <- rraf(Pram, by_pop = TRUE)
  nout <- numeric(ncol(tab(Pram)))
  exp_matrix <- t(vapply(seppop(Pram), rraf, nout, res = "vector"))
  expect_is(rrx_matrix, "matrix")
  expect_equivalent(dim(rrx_matrix), c(nPop(Pram), ncol(tab(Pram))))
  expect_equivalent(rownames(rrx_matrix), popNames(Pram))
  # The correction is 1/N per pop. PistolRSF_OR has a pop size of 4
  expect_equivalent(rrx_matrix[nrow(rrx_matrix), ncol(rrx_matrix)], 1 / 4)
  expect_equal(rrx_matrix, exp_matrix)
})

test_that("rraf calculates per population when supplied with a population factor", {
  skip_on_cran()
  data(Pram)
  rrx_matrix <- rraf(Pram, by_pop = TRUE, correction = FALSE)
  rr_pop_factor <- rraf(Pram, pop = as.character(pop(Pram)), correction = FALSE)
  rr_formula <- rraf(Pram, pop = ~ SOURCE / STATE, correction = FALSE)
  expect_equivalent(rrx_matrix, rr_pop_factor)
  expect_equivalent(rrx_matrix, rr_formula)
})

test_that("psex and pgen internals produce expected results", {
  skip_on_cran()

  # setup -------------------------------------------------------------------

  # From Parks and Werth, 1993
  x <- "
   Hk Lap Mdh2 Pgm1 Pgm2 X6Pgd2
  54 12 12 12 23 22 11
  36 22 22 11 22 33 11
  10 23 22 11 33 13 13"

  xtab <- read.table(text = x, header = TRUE, row.names = 1)
  xgid <- df2genind(
    xtab[rep(rownames(xtab), as.integer(rownames(xtab))), ],
    ncode = 1
  )
  afreq <- c(
    Hk.1 = 0.167,
    Hk.2 = 0.795,
    Hk.3 = 0.038,
    Lap.1 = 0.190,
    Lap.2 = 0.798,
    Lap.3 = 0.012,
    Mdh2.0 = 0.011,
    Mdh2.1 = 0.967,
    Mdh2.2 = 0.022,
    Pgm1.2 = 0.279,
    Pgm1.3 = 0.529,
    Pgm1.4 = 0.162,
    Pgm1.5 = 0.029,
    Pgm2.1 = 0.128,
    Pgm2.2 = 0.385,
    Pgm2.3 = 0.487,
    X6Pgd2.1 = 0.526,
    X6Pgd2.2 = 0.051,
    X6Pgd2.3 = 0.423
  )
  freqs <- afreq[colnames(tab(xgid))]
  mfreq <- matrix(freqs, nrow = 1, dimnames = list(NULL, names(freqs)))
  pNotGen <- psex(xgid, by_pop = FALSE, freq = freqs, G = 45, method = "single")
  pGen <- exp(rowSums(pgen(xgid, by_pop = FALSE, freq = freqs)))
  pGenm <- exp(rowSums(pgen(xgid, by_pop = FALSE, freq = mfreq)))
  res <- matrix(c(unique(pGen), unique(pNotGen)), ncol = 2)
  expected_result <- structure(
    c(
      4.14726753733799e-05,
      0.00192234266749449,
      0.000558567714432373,
      0.00186456862059092,
      0.0829457730256321,
      0.024829127718338
    ),
    .Dim = c(3L, 2L)
  )

  # tests -------------------------------------------------------------------
  # Testing that pgen can take both a matrix and vector of af
  expect_equal(pGen, pGenm)
  expect_error(
    pgen(xgid, by_pop = FALSE, freq = mfreq[, -1, drop = FALSE]),
    "frequency matrix"
  )
  expect_error(
    pgen(xgid, by_pop = FALSE, freq = mfreq[, -1, drop = TRUE]),
    "frequencies"
  )
  # Testing single encounter
  expect_equivalent(res, expected_result)
})


test_that("psex produces a vector", {
  skip_on_cran()
  data(Pram)
  psexpram <- psex(Pram)
  expect_is(psexpram, "numeric")
  expect_equal(length(psexpram), nInd(Pram))
})

test_that("psex can use the multiple method (not testing for accuracy)", {
  skip_on_cran()
  data(Pram)
  psexpram <- psex(Pram, method = "multiple")
  expect_is(psexpram, "numeric")
  expect_equal(length(psexpram), nInd(Pram))
})
test_that("psex can take population factors", {
  skip_on_cran()
  data(Pram)
  psexpram <- psex(Pram)
  psexpop <- psex(Pram, pop = as.character(pop(Pram)))
  psexform <- psex(Pram, pop = ~ SOURCE / STATE)
  expect_equivalent(psexpram, psexpop)
  expect_equivalent(psexpram, psexform)
})

test_that("pgen produces a matrix", {
  skip_on_cran()
  data(Pram)
  pgenlog <- pgen(Pram, by_pop = FALSE)
  expect_is(pgenlog, "matrix")
  expect_equivalent(exp(pgenlog), pgen(Pram, by_pop = FALSE, log = FALSE))
})

test_that("psex can take population factors", {
  skip_on_cran()
  data(Pram)
  pgenpram <- pgen(Pram)
  pgenpop <- pgen(Pram, pop = as.character(pop(Pram)))
  pgenform <- pgen(Pram, pop = ~ SOURCE / STATE)
  expect_equivalent(pgenpram, pgenpop)
  expect_equivalent(pgenpram, pgenform)
})

test_that("pgen and psex work for haploids", {
  skip_on_cran()
  data(monpop)
  expect_is(psex(monpop, by_pop = FALSE), "numeric")
})

test_that("pgen can't work with polyploids", {
  skip_on_cran()
  data(Pinf)
  expect_error(pgen(Pinf))
})

test_that("psex is accurate", {
  skip_on_cran()
  skip_if_not_installed("RClone")
  # Setup
  #
  data("zostera", package = "RClone")
  rzos <- RClone::convert_GC(zostera[, -c(1:3)], num = 3)
  pzos <- as.genclone(df2genind(
    zostera[, -c(1:3)],
    ncode = 3,
    pop = zostera[, 1],
    strata = zostera[, 1, drop = FALSE]
  ))
  extract_psex <- function(x) {
    as.numeric(x$psex)
  }
  #
  # Testing basic equality, no population
  rzpsex <- RClone::psex(rzos, RR = TRUE) %>%
    extract_psex %>%
    setNames(indNames(pzos))
  pzpsex <- poppr::psex(pzos, by_pop = FALSE, method = "multiple")
  nas <- is.na(rzpsex)

  expect_equal(rzpsex[!nas], pzpsex[!nas])
  #
  # Testing equality with population
  rzpsex <- RClone::psex(rzos, RR = TRUE, vecpop = pop(pzos)) %>%
    lapply(extract_psex) %>%
    unlist() %>%
    setNames(indNames(pzos))
  gtab <- table(pop(pzos))
  pzpsex <- poppr::psex(pzos, by_pop = TRUE, method = "multiple")
  pzpsex2 <- poppr::psex(pzos, by_pop = TRUE, G = gtab, method = "multiple")
  expect_error(poppr::psex(
    pzos,
    by_pop = TRUE,
    G = as.vector(gtab),
    method = "multiple"
  ))
  nas <- is.na(rzpsex)
  expect_equal(rzpsex[!nas], pzpsex[!nas])
  expect_equal(rzpsex[!nas], pzpsex2[!nas])
})