tests/testthat/test-alpha.R

context("alpha diversity")

test_that("Alpha diversity values are consistent with QIIME2 tests", {
  x_qiime2 <- c(0, 1, 1, 4, 2, 5, 2, 4, 1, 2)
  expect_equal(berger_parker_d(x_qiime2), 5 / 22)
  expect_equal(brillouin_d(c(1, 2, 0, 0, 3, 1)), 0.863, tol=0.001)
  expect_equal(chao1(x_qiime2), 9.75)
  expect_equal(chao1(x_qiime2, bias_corrected=FALSE), 10.5)
  x_no_singles <- c(0, 2, 2, 4, 5, 0, 0, 0, 0, 0)
  expect_equal(chao1(x_no_singles), 4)
  expect_equal(chao1(x_no_singles, bias_corrected=FALSE), 4)
  x_no_doubles <- c(0, 1, 1, 4, 5, 0, 0, 0, 0, 0)
  expect_equal(chao1(x_no_doubles), 5)
  expect_equal(chao1(x_no_doubles, bias_corrected=FALSE), 5)
  expect_equal(dominance(c(1, 0, 2, 5, 2)), 0.340, tol=0.001)
  expect_equal(doubles(x_qiime2), 3)
  expect_equal(doubles(c(0, 3, 4)), 0)
  expect_equal(doubles(2), 1)
  expect_equal(enspie(c(1, 1, 1, 1, 1, 1)), 6)
  expect_equal(enspie(c(13, 13, 13, 13)), 4)
  x_enspie1 <- c(1, 41, 0, 0, 12, 13)
  expect_equal(enspie(x_enspie1), 2.250, tol=0.001)
  expect_equal(enspie(x_enspie1), 1 / dominance(x_enspie1))
  x_enspie2 <- c(1, 0, 2, 5, 2)
  expect_equal(enspie(x_enspie2), 1 / dominance(x_enspie2))
  # TODO: esty_ci
  #expect_equal(fisher_alpha(c(4, 3, 4, 0, 1, 0, 2)), 2.782, tol=0.001)
  #expect_equal(fisher_alpha(c(1, 6, 1, 0, 1, 0, 5)), 2.782, tol=0.001)
  #expect_equal(fisher_alpha(c(61, 0, 0, 1)), 0.395, tol=0.001)
  #expect_equal(fisher_alpha(c(999, 0, 10)), 0.240, tol=0.001)
  # TODO: gini_index
  arr_good <- c(rep(1, 75), 2, 2, 2, 2, 2, 2, 3, 4, 4)
  expect_equal(goods_coverage(arr_good), 0.235, tol=0.001)
  x_heip <- c(1, 2, 3, 1)
  expect_equal(heip_e(x_heip), (exp(shannon(x_heip)) - 1) / 3)
  expect_equal(heip_e(c(500, 300, 200)), 0.90, tol=0.01)
  expect_equal(heip_e(c(500, 299, 200, 1)), 0.61, tol=0.01)
  x_kempton <- c(
    2, 3, 3, 3, 3, 3, 4, 4, 4, 6, 6, 7, 7, 9, 9, 11, 14,
    15, 15, 20, 29, 33, 34, 36, 37, 53, 57, 138, 146, 170)
  expect_equal(kempton_taylor_q(x_kempton), 14 / log(34 / 4))
  expect_equal(kempton_taylor_q(sample(x_kempton)), 14 / log(34 / 4))
  expect_equal(margalef(x_qiime2), 8 / log(22))
  expect_equal(mcintosh_d(c(1, 2, 3)), 0.636, tol=0.001)
  expect_equal(mcintosh_e(c(1, 2, 3, 1)), sqrt(15 / 19))
  expect_equal(menhinick(x_qiime2), 9 / sqrt(22))
  # TODO: michaelis_menten_fit
  expect_equal(richness(c(4, 3, 4, 0, 1, 0, 2)), 5)
  expect_equal(richness(x_qiime2), 9)
  x_pielou <- c(1, 2, 3, 1)
  expect_equal(pielou_e(x_pielou), shannon(x_pielou) / log(4))
  expect_equal(pielou_e(x_qiime2), 0.925, tol=0.001)
  expect_equal(pielou_e(c(1, 1)), 1)
  expect_equal(pielou_e(c(1, 1, 196, 1, 1)), 0.078, tol=0.001)
  expect_equal(pielou_e(c(0, 0, 200, 0, 0)), NaN)
  expect_equal(robbins(c(1, 2, 3, 0, 1)), 2 / 7)
  expect_equal(shannon(5, base = 2), 0)
  expect_equal(shannon(c(5, 5), base = 2), 1)
  expect_equal(shannon(c(1, 1, 1, 1, 0), base = 2), 2)
  expect_equal(simpson(c(1, 0, 2, 5, 2)), 0.66)
  expect_equal(simpson(5), 0)
  expect_equal(simpson_e(c(1, 1, 1, 1, 1, 1, 1)), 1)
  expect_equal(simpson_e(c(500, 400, 600, 500)), 0.980, tol=0.001)
  expect_equal(singles(x_qiime2), 3)
  expect_equal(singles(c(0, 3, 4)), 0)
  expect_equal(singles(1), 1)
  expect_equal(strong(c(1, 2, 3, 1)), 0.214, tol=0.001)
})

test_that("Invalid count vectors are detected", {
  expect_true(check_positive(c(0, 0, 1, 2))) # counts pass
  expect_true(check_positive(c(0, 0.1))) # proportions pass
  expect_error(check_positive(c(0, 3, -1)))
  expect_error(richness(c(0, 0, -0.01, 0.99)))
})

test_that("Empty count vectors produce warnings or valid result", {
  x_empty <- c(0, 0, 0, 0)
  expect_warning(berger_parker_d(x_empty), "not defined")
  expect_warning(brillouin_d(x_empty), "not defined")
  expect_equal(chao1(x_empty), 0)
  expect_warning(dominance(x_empty), "not defined")
  expect_equal(doubles(x_empty), 0)
  expect_warning(enspie(x_empty), "not defined")
  expect_warning(etsy_ci(x_empty), "not defined")
  #expect_equal(fisher_alpha(x_empty), 1)
  expect_warning(goods_coverage(x_empty), "not defined")
  expect_warning(heip_e(x_empty), "not defined")
  expect_warning(invsimpson(x_empty), "not defined")
  expect_warning(kempton_taylor_q(x_empty), "not defined")
  expect_warning(margalef(x_empty), "not defined")
  expect_warning(mcintosh_d(x_empty), "not defined")
  expect_warning(mcintosh_e(x_empty), "not defined")
  expect_warning(menhinick(x_empty), "not defined")
  # TODO: michaelis_menten_fit
  expect_warning(pielou_e(x_empty), "not defined")
  expect_equal(richness(x_empty), 0)
  expect_warning(robbins(x_empty), "not defined")
  expect_warning(shannon(x_empty), "not defined")
  expect_warning(simpson(x_empty), "not defined")
  expect_warning(simpson_e(x_empty), "not defined")
  expect_equal(singles(x_empty), 0)
  expect_warning(strong(x_empty), "not defined")
})

test_that("Empty count vectors give correct values", {
  x_empty <- c(0, 0, 0, 0)
  expect_equal(suppressWarnings(berger_parker_d(x_empty)), NA)
  expect_equal(suppressWarnings(brillouin_d(x_empty)), NA)
  expect_equal(suppressWarnings(dominance(x_empty)), NA)
  expect_equal(doubles(x_empty), 0)
  expect_equal(suppressWarnings(enspie(x_empty)), NA)
  expect_equal(suppressWarnings(etsy_ci(x_empty)), NA)
  #expect_equal(fisher_alpha(x_empty), 1)
  expect_equal(suppressWarnings(goods_coverage(x_empty)), NA)
  expect_equal(suppressWarnings(heip_e(x_empty)), NA)
  expect_equal(suppressWarnings(invsimpson(x_empty)), NA)
  expect_equal(suppressWarnings(kempton_taylor_q(x_empty)), NA)
  expect_equal(suppressWarnings(margalef(x_empty)), NA)
  expect_equal(suppressWarnings(mcintosh_d(x_empty)), NA)
  expect_equal(suppressWarnings(mcintosh_e(x_empty)), NA)
  expect_equal(suppressWarnings(menhinick(x_empty)), NA)
  # TODO: michaelis_menten_fit
  expect_equal(suppressWarnings(pielou_e(x_empty)), NA)
  expect_equal(richness(x_empty), 0)
  expect_equal(suppressWarnings(robbins(x_empty)), NA)
  expect_equal(suppressWarnings(shannon(x_empty)), NA)
  expect_equal(suppressWarnings(simpson(x_empty)), NA)
  expect_equal(suppressWarnings(simpson_e(x_empty)), NA)
  expect_equal(singles(x_empty), 0)
  expect_equal(suppressWarnings(strong(x_empty)), NA)
})

test_that("Alpha diversity values are correct for simple example", {
  x_simple <- c(0, 1, 1, 4, 3, 2)
  expect_equal(berger_parker_d(x_simple), 4 / 11)
  expect_equal(brillouin_d(x_simple), 1.076, tol=0.001)
  expect_equal(dominance(x_simple), 0.256, tol=0.001)
  expect_equal(doubles(x_simple), 1)
  expect_equal(enspie(x_simple), invsimpson(x_simple))
  expect_equal(etsy_ci(x_simple), c(-0.158, 0.522), tol=0.001)
  #expect_equal(fisher_alpha(x_simple), 3.538, tol=0.001)
  expect_equal(goods_coverage(x_simple), 0.818, tol=0.001)
  expect_equal(heip_e(x_simple), 0.835, tol=0.001)
  expect_equal(invsimpson(x_simple), 3.903, tol=0.001)
  expect_equal(kempton_taylor_q(x_simple), 1.820, tol=0.001)
  expect_equal(margalef(x_simple), 1.668, tol=0.001)
  expect_equal(mcintosh_d(x_simple), 0.707, tol=0.01)
  expect_equal(mcintosh_e(x_simple), 0.765, tol=0.001)
  expect_equal(menhinick(x_simple), 1.508, tol=0.001)
  # TODO: michaelis_menten_fit
  expect_equal(pielou_e(x_simple), 0.912, tol=0.001)
  expect_equal(richness(x_simple), 5)
  expect_equal(robbins(x_simple), 2 /11)
  expect_equal(shannon(x_simple), 1.468, tol=0.001)
  expect_equal(simpson(x_simple), 0.744, tol=0.001)
  expect_equal(simpson_e(x_simple), 0.781, tol=0.001)
  expect_equal(singles(x_simple), 2)
  expect_equal(strong(x_simple), 0.236, tol=0.001)
})
kylebittinger/abdiv documentation built on Dec. 15, 2019, 12:32 a.m.