tests/testthat/test_recapr.R

In mbtyers/recapr: Two Event Mark-Recapture Experiment

test_that("funcs",{
  expect_equal(NChapman(n1=100, n2=100, m2=20), 484.7619, tolerance=0.001)
  expect_equal(NPetersen(n1=100, n2=100, m2=20), 500, tolerance=0.001)
  expect_equal(NBailey(n1=100, n2=100, m2=20), 480.9524, tolerance=0.001)
  expect_equal(vChapman(n1=100, n2=100, m2=20), 6729.169, tolerance=0.001)
  expect_equal(vPetersen(n1=100, n2=100, m2=20), 10000, tolerance=0.001)
  expect_equal(vBailey(n1=100, n2=100, m2=20), 8328.18, tolerance=0.001)
  expect_equal(seChapman(n1=100, n2=100, m2=20), 82.03151, tolerance=0.001)
  expect_equal(sePetersen(n1=100, n2=100, m2=20), 100, tolerance=0.001)
  expect_equal(seBailey(n1=100, n2=100, m2=20), 91.25886, tolerance=0.001)

  # should probably do something with lengths
})

draws1 <- rChapman(length=100000, N=500, n1=100, n2=100)
draws2 <- rPetersen(length=100000, N=500, n1=100, n2=100)
draws3 <- rBailey(length=100000, N=500, n1=100, n2=100)
out1 <- pChapman(nullN=500, n1=100, n2=100, m2=28)
out12 <- pChapman(nullN=500, n1=100, n2=100, m2=28,alternative="2-sided")
out2 <- pPetersen(nullN=500, n1=100, n2=100, m2=28)
out22 <- pPetersen(nullN=500, n1=100, n2=100, m2=28,alternative="2-sided")
out3 <- pBailey(nullN=500, n1=100, n2=100, m2=28)
out32 <- pBailey(nullN=500, n1=100, n2=100, m2=28,alternative="2-sided")
test_that("rdraws_pvals",{
  expect_equal(mean(draws1), 499.9744, tolerance=1)
  expect_equal(mean(draws2), 517.7152, tolerance=1)
  expect_equal(mean(draws3), 500.371, tolerance=1)

  expect_error(rChapman(length=10, N=100, n1=1:10, n2=1:4), "n2 must be a single number, or a vector of length equal to the number of draws")
  expect_error(rChapman(length=10, N=100, n1=10), "need to supply either n2 or p2")
  expect_warning(rChapman(length=10, N=100, n1=10, n2=5, p2=.5), "both n2 and p2 specified - only n2 used")
  expect_equal(mean(suppressWarnings(rChapman(length=100000, N=100, n1=10, n2=5, p2=.5))), 50.38573, tolerance=10)

  expect_equal(length(out1),2,tolerance=0.001)
  expect_equal(length(out2),2,tolerance=0.001)
  expect_equal(length(out3),2,tolerance=0.001)
  expect_equal(out1$pval, 0.02036, tolerance=0.05)
  expect_equal(out12$pval, 0.07989, tolerance=0.05)
  expect_equal(out2$pval, 0.02086, tolerance=0.05)
  expect_equal(out22$pval, 0.11992, tolerance=0.05)
  expect_equal(out3$pval, 0.03459, tolerance=0.05)
  expect_equal(out32$pval, 0.1145, tolerance=0.05)

  expect_equal(powChapman(nullN=500, trueN=400, n1=100, n2=100, nsim=5000), 0.3358, 0.05)
  expect_equal(powPetersen(nullN=500, trueN=400, n1=100, n2=100, nsim=5000), 0.3354, 0.05)
  expect_equal(powBailey(nullN=500, trueN=400, n1=100, n2=100, nsim=5000), 0.2774, 0.05)
})

out1 <- ciChapman(n1=100, n2=100, m2=20)
out2 <- ciPetersen(n1=100, n2=100, m2=20)
out3 <- ciBailey(n1=100, n2=100, m2=20)
test_that("ci",{
  expect_equal(length(out1),3,tolerance=0.001)
  expect_equal(length(out2),3,tolerance=0.001)
  expect_equal(length(out3),3,tolerance=0.001)
  expect_equal(sum(out1$ciNorm), 969.5238, tolerance=1)
  expect_equal(sum(out2$ciNorm), 1000, tolerance=1)
  expect_equal(sum(out3$ciNorm), 961.9048, tolerance=1)
  expect_equal(sum(out1$ciBoot), 1078.401, tolerance=1)
  expect_equal(sum(out2$ciBoot), 1190.476, tolerance=1)
  expect_equal(sum(out3$ciBoot), 1125.199, tolerance=1)
})

rr <- n2RR(N=1000, n1=100)
rr2 <- n2RR(N=1000, n1=100, conf=c(0.99,0.95,0.92), acc=c(0.5,0.25,0.22))
test_that("n2",{
  expect_equal(length(rr), 6, tolerance=0.001)
  expect_equal(sum(rr$conf_0.99), 8155.26, tolerance=0.1)
  expect_equal(sum(rr$conf_0.95), 7070.26, tolerance=0.1)
  expect_equal(sum(rr$conf_0.9), 6389.26, tolerance=0.1)
  expect_equal(sum(rr$conf_0.85), 5895.26, tolerance=0.1)
  expect_equal(sum(rr$conf_0.8), 5492.26, tolerance=0.1)
  expect_equal(sum(rr$conf_0.75), 5139.26, tolerance=0.1)

  expect_equal(length(rr2), 2, tolerance=0.001)
  expect_equal(sum(rr2$conf_0.99), 1306.75, tolerance=0.1)

  expect_error(n2RR(N=1000,n1=100,conf=c(.2,.4)),"invalid conf")
  expect_error(n2RR(N=1000,n1=100,acc=c(.3,.4)),"invalid acc")
  expect_error(plotn2sim(N=1000, n1=100,conf=.2),"invalid conf")
})

mat <- matrix(c(30,15,1,0,22,15), nrow=2, ncol=3, byrow=TRUE)
ctest <- consistencytest(n1=c(284,199), n2=c(347,3616,1489), stratamat=mat)
stest <- strattest(n1=c(100,100), n2=c(50,200), m2=c(20,15))
test_that("consistency",{
  expect_equal(length(ctest), 12, tolerance=0.001)
  expect_equal(sum(ctest$test1_tab), 483, tolerance=0.001)
  expect_equal(ctest$test1_Xsqd, 44.43179, tolerance=0.001)
  expect_equal(ctest$test1_df,  3, tolerance=0.001)
  expect_equal(ctest$test1_pval, 0, tolerance=0.001)
  expect_equal(sum(ctest$test2_tab), 5452, tolerance=0.001)
  expect_equal(ctest$test2_Xsqd, 125.4408, tolerance=0.001)
  expect_equal(ctest$test2_df,  2, tolerance=0.001)
  expect_equal(ctest$test2_pval, 0, tolerance=0.001)
  expect_equal(sum(ctest$test3_tab), 483, tolerance=0.001)
  expect_equal(ctest$test3_Xsqd, 0.3185941, tolerance=0.001)
  expect_equal(ctest$test3_df,  1, tolerance=0.001)
  expect_equal(ctest$test3_pval, 0.5724538, tolerance=0.001)

  expect_error(consistencytest(n1=c(284,199), n2=c(347,3616,1489)),"recapture strata must be specified, either with stratamat or m2strata1 and m2strata2 together")
  expect_error(consistencytest(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2),m2strata2=c(1,2,3)),"m2strata1 and m2strata2 must be the same length")# (each element corresponding to an individual)")
  expect_error(consistencytest(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,4),m2strata2=c(1,2,3)),"strata values larger than count vector in event 1")
  expect_error(consistencytest(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,2.2),m2strata2=c(1,2,3)),"m2strata1 values must be positive integers")
  expect_error(consistencytest(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,"a"),m2strata2=c(1,2,3)),"m2strata1 values must be positive integers")
  expect_error(consistencytest(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,-1),m2strata2=c(1,2,3)),"m2strata1 values must be positive integers")

  expect_equal(length(stest), 12, tolerance=0.001)
  expect_equal(sum(stest$event1_table), 250, tolerance=0.001)
  expect_equal(stest$event1_Xsqd, 32.44394, tolerance=0.001)
  expect_equal(stest$event1_df, 1, tolerance=0.001)
  expect_equal(stest$event1_pval, 0, tolerance=0.001)
  expect_equal(sum(stest$event2_table), 200, tolerance=0.001)
  expect_equal(stest$event2_Xsqd, 0.5541126, tolerance=0.001)
  expect_equal(stest$event2_df, 1, tolerance=0.001)
  expect_equal(stest$event2_pval, 0.4566422, tolerance=0.001)
  expect_equal(sum(stest$event1v2_table), 450, tolerance=0.001)
  expect_equal(stest$event1v2_Xsqd, 43.66013, tolerance=0.001)
  expect_equal(stest$event1v2_df, 1, tolerance=0.001)
  expect_equal(stest$event1v2_pval, 3.906508e-11, tolerance=0.001)
})

pctest <- powconsistencytest(n1=c(100,180), n2=c(100,100,100), pmat=matrix(c(1,2,2,5,2,3,1,5),nrow=2,ncol=4,byrow=T), nsim=100000)
pstest <- powstrattest(N=c(1000,2000,1000), n1=c(100,180,100), n2=c(80,160,70), nsim=100000)
test_that("consistency power", {
  expect_equal(length(pctest), 16, tolerance=0.001)
  expect_equal(pctest$pwr1_c, 0.8082572, tolerance=0.001)
  expect_equal(pctest$pwr1_sim, 0.7755, tolerance=0.01)
  expect_equal(pctest$pwr2_c, 0.9953951, tolerance=0.001)
  expect_equal(pctest$pwr2_sim, 0.99689, tolerance=0.01)
  expect_equal(pctest$pwr3_c, 0.1129534, tolerance=0.001)
  expect_equal(pctest$pwr3_sim, 0.11285, tolerance=0.01)
  expect_equal(pctest$ntest1, 280, tolerance=0.001)
  expect_equal(pctest$ntest2, 300, tolerance=0.001)
  expect_equal(pctest$ntest3, 280, tolerance=0.001)
  expect_equal(pctest$p1test1[1,1], 0.03571429, tolerance=0.001)
  expect_equal(pctest$p0test1[1,1], 0.05449907, tolerance=0.001)
  expect_equal(pctest$p1test2[1,1], 0.1424242, tolerance=0.001)
  expect_equal(pctest$p0test2[1,1], 0.1646465, tolerance=0.001)
  expect_equal(pctest$p1test3[1,1], 0.1785714, tolerance=0.001)
  expect_equal(pctest$p0test3[1,1], 0.1890074, tolerance=0.001)
  expect_equal(pctest$alpha, 0.05, tolerance=0.001)

  expect_equal(length(pstest), 3, tolerance=0.001)
  expect_equal(length(pstest$Test1), 7, tolerance=0.001)
  expect_equal(length(pstest$Test2), 7, tolerance=0.001)
  expect_equal(length(pstest$Test3), 7, tolerance=0.001)
  expect_equal(pstest$Test1$prob[1], 0.1, tolerance=0.001)
  expect_equal(pstest$Test2$prob[1], 0.08, tolerance=0.001)
  expect_equal(pstest$Test3$prob[1], 0.4444444, tolerance=0.001)
  expect_equal(pstest$Test1$prob_null[1], 0.095, tolerance=0.001)
  expect_equal(pstest$Test2$prob_null[1], 0.0775, tolerance=0.001)
  expect_equal(pstest$Test3$prob_null[1], 0.4492754, tolerance=0.001)
  expect_equal(pstest$Test1$n, 310, tolerance=0.001)
  expect_equal(pstest$Test2$n, 380, tolerance=0.001)
  expect_equal(pstest$Test3$n, 690, tolerance=0.001)
  expect_equal(pstest$Test1$alpha, 0.05, tolerance=0.001)
  expect_equal(pstest$Test2$alpha, 0.05, tolerance=0.001)
  expect_equal(pstest$Test3$alpha, 0.05, tolerance=0.001)
  expect_equal(pstest$Test1$power, 0.05682567, tolerance=0.001)
  expect_equal(pstest$Test2$power, 0.05755484, tolerance=0.001)
  expect_equal(pstest$Test3$power, 0.1884238, tolerance=0.001)
  expect_equal(pstest$Test1$power_sim, 0.06139, tolerance=0.01)
  expect_equal(pstest$Test2$power_sim, 0.05257, tolerance=0.01)
  expect_equal(pstest$Test3$power_sim, 0.16954, tolerance=0.01)

  expect_error(powconsistencytest(n1=c(100,180,100), n2=c(100,100,100), pmat=matrix(c(1,2,2,5,2,3,1,5),nrow=2,ncol=4,byrow=T)), "pmat must have the same number of rows as first-event strata")
  expect_error(powconsistencytest(n1=c(100,180), n2=c(100,100,100,1), pmat=matrix(c(1,2,2,5,2,3,1,5),nrow=2,ncol=4,byrow=T)), "pmat must have the same number of columns as second-event strata, plus one")
  expect_error(powstrattest(N=c(1000,2000), n1=c(100,180,100), n2=c(80,160,70)), "N, n1, and n2 vectors must be the same length.")
  expect_error(powstrattest(N=c(1000,2000,1000), n1=c(100,180), n2=c(80,160,70)), "N, n1, and n2 vectors must be the same length.")
  expect_error(powstrattest(N=c(1000,2000,1000), n1=c(100,180,100), n2=c(80,160)), "N, n1, and n2 vectors must be the same length.")
})

ci <- cistrat(n1=c(100,200), n2=c(100,500), m2=c(10,10),bootreps=1000000)
test_that("stratified",{
  expect_equal(Nstrat(n1=c(100,200), n2=c(100,500), m2=c(10,10)), 10080, tolerance=0.1)
  expect_error(Nstrat(n1=c(100,200), n2=c(100,500), m2=c(10,10), estimator="cheese"), "invalid estimator")
  expect_error(Nstrat(n1=c(100,200), n2=c(100,500), m2=c(10,10,3)), "n1, n2, and m2 vectors must be of equal length")
  expect_equal(vstrat(n1=c(100,200), n2=c(100,500), m2=c(10,10)), 6513699, tolerance=0.1)
  expect_equal(sestrat(n1=c(100,200), n2=c(100,500), m2=c(10,10)), 2552.195, tolerance=0.1)
  expect_equal(mean(rstrat(length=1000000, N=c(5000,10000), n1=c(500,200), n2=c(500,200))), 14838.2, tolerance=100)
  expect_equal(length(ci), 4, tolerance=0.001)
  expect_equal(ci$Nstrat, 10080, tolerance=0.1)
  expect_equal(sum(ci$Nhat_by_strat), 10080, tolerance=0.1)
  expect_equal(sum(ci$ciNorm_strat), 20160, tolerance=0.1)
  expect_equal(sum(ci$ciBoot_strat), 27573.12, tolerance=10)

  expect_error(rstrat(length=10, N=c(1000,1000), p1=c(.1,.2), p2=.3), "p1, p2, and N vectors must be of equal length")
  expect_error(rstrat(length=10, N=c(1000,1000), p1=c(.1,.2), n2=3), "p1, n2, and N vectors must be of equal length")
  expect_error(rstrat(length=10, N=c(1000,1000), p1=c(.1,.2)), "either n2 or p2 vectors must be specified")
})

mat <- matrix(c(59,30,1,45,280,38,0,42,25), nrow=3, ncol=3, byrow=TRUE)
darr <- NDarroch(n1=c(484,1468,399), n2=c(847,6616,2489), stratamat=mat)
test_that("Darroch",{
  expect_equal(length(darr), 6, tolerance=0.001)
  expect_equal(length(darr$Nhat_strata1), 3, tolerance=0.01)
  expect_equal(length(darr$se_Nhat_strata1), 3, tolerance=0.01)
  expect_equal(length(darr$Nhat_strata2), 3, tolerance=0.01)
  expect_equal(length(darr$se_Nhat_strata1), 3, tolerance=0.01)
  expect_equal(sum(darr$Nhat_strata1), 54053.47, tolerance=0.01)
  expect_equal(sum(darr$se_Nhat_strata1), 17796.6, tolerance=0.01)
  expect_equal(sum(darr$Nhat_strata2), 54053.47, tolerance=0.01)
  expect_equal(sum(darr$se_Nhat_strata2), 12971.88, tolerance=0.01)
  expect_equal(darr$Nhat, 54053.47, tolerance=0.01)
  expect_equal(darr$se_Nhat, 4890.999, tolerance=0.01)
  expect_equal(darr$Nhat_strata1[1], 3584.554, tolerance=0.01)
  expect_equal(darr$se_Nhat_strata1[1], 1914.125, tolerance=0.01)
  expect_equal(darr$Nhat_strata2[1], 6360.409, tolerance=0.01)
  expect_equal(darr$se_Nhat_strata2[1], 848.4842, tolerance=0.01)

  expect_error(NDarroch(n1=c(284,199), n2=c(347,3616,1489)),"recapture strata must be specified, either with stratamat or m2strata1 and m2strata2 together")
  expect_error(NDarroch(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2),m2strata2=c(1,2,3)),"m2strata1 and m2strata2 must be the same length")# (each element corresponding to an individual)")
  expect_error(NDarroch(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,4),m2strata2=c(1,2,3)),"strata values larger than count vector in event 1")
  expect_error(NDarroch(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,2.2),m2strata2=c(1,2,3)),"m2strata1 values must be positive integers")
  expect_error(NDarroch(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,"a"),m2strata2=c(1,2,3)),"m2strata1 values must be positive integers")
  expect_error(NDarroch(n1=c(284,199), n2=c(347,3616,1489), m2strata1=c(1,2,-1),m2strata2=c(1,2,3)),"m2strata1 values must be positive integers")
})