aboutsummaryrefslogtreecommitdiff
path: root/tests/testthat/test_occuCOP.R
diff options
context:
space:
mode:
Diffstat (limited to 'tests/testthat/test_occuCOP.R')
-rw-r--r--tests/testthat/test_occuCOP.R485
1 files changed, 485 insertions, 0 deletions
diff --git a/tests/testthat/test_occuCOP.R b/tests/testthat/test_occuCOP.R
new file mode 100644
index 0000000..ca880bd
--- /dev/null
+++ b/tests/testthat/test_occuCOP.R
@@ -0,0 +1,485 @@
+context("occuCOP fitting function")
+skip_on_cran()
+
+COPsimul <- function(psi = 0.5,
+ lambda = 1,
+ M = 100,
+ J = 5) {
+
+ z_i <- sample(
+ x = c(0, 1),
+ size = M,
+ prob = c(1 - psi, psi),
+ replace = T
+ )
+
+ y = matrix(rpois(n = M * J, lambda = lambda), nrow = M, ncol = J) * z_i
+
+ return(y)
+}
+
+
+test_that("unmarkedFrameOccuCOP is constructed correctly", {
+ set.seed(123)
+
+ # Simulate data
+ M = 100
+ J = 5
+ y = COPsimul(psi = .5,
+ lambda = 1,
+ M = M,
+ J = J)
+ L = y * 0 + 1
+
+ psiCovs = data.frame(
+ "psiNum" = rnorm(M),
+ "psiInt" = as.integer(rpois(n = M, lambda = 5)),
+ "psiBool" = sample(c(T, F), size = M, replace = T),
+ "psiChar" = sample(letters[1:5], size = M, replace = T),
+ "psiFactUnord" = factor(sample(
+ letters[5:10], size = M, replace = T
+ )),
+ "psiFactOrd" = sample(factor(c("A", "B", "C"), ordered = T), size =
+ M, replace = T)
+ )
+
+ lambdaCovs = list(
+ "lambdaNum" = matrix(
+ rnorm(M * J),
+ nrow = M, ncol = J
+ ),
+ "lambdaInt" = matrix(
+ as.integer(rpois(n = M * J, lambda = 1e5)),
+ nrow = M, ncol = J
+ ),
+ "lambdaBool" = matrix(
+ sample(c(T, F), size = M * J, replace = T),
+ nrow = M, ncol = J
+ ),
+ "lambdaChar" = matrix(
+ sample(letters[1:5], size = M * J, replace = T),
+ nrow = M, ncol = J
+ ),
+ "lambdaFactUnord" = matrix(
+ factor(sample(letters[5:10], size = M * J, replace = T)),
+ nrow = M, ncol = J
+ ),
+ "lambdaFactOrd" = matrix(
+ sample(factor(c("A", "B", "C"), ordered = T), size = M * J, replace = T),
+ nrow = M, ncol = J
+ )
+ )
+
+
+ # Creating a unmarkedFrameOccuCOP object
+ expect_warning(umf <- unmarkedFrameOccuCOP(y = y))
+ expect_no_error(umf <- unmarkedFrameOccuCOP(y = y, L = L))
+
+
+ # Create subsets
+ expect_no_error(umf_sub_i <- umf[1:3, ])
+ expect_no_error(umf_sub_j <- umf[, 1:2])
+ expect_no_error(umf_sub_ij <- umf[1:3, 1:2])
+
+ # unmarkedFrameOccuCOP organisation ----------------------------------------------
+ expect_true(inherits(umf, "unmarkedFrameOccuCOP"))
+ expect_equivalent(numSites(umf_sub_i), 3)
+ expect_equivalent(obsNum(umf_sub_j), 2)
+ expect_equivalent(numSites(umf_sub_ij), 3)
+ expect_equivalent(obsNum(umf_sub_ij), 2)
+
+ # unmarkedFrameOccuCOP display ---------------------------------------------------
+
+ # print method
+ expect_output(print(umf), "Data frame representation of unmarkedFrame object")
+ expect_output(print(umf_sub_i), "Data frame representation of unmarkedFrame object")
+ expect_output(print(umf[1,]), "Data frame representation of unmarkedFrame object")
+ expect_output(print(umf[,1]), "Data frame representation of unmarkedFrame object")
+ expect_output(print(umf[1,1]), "Data frame representation of unmarkedFrame object")
+
+ # summary method for unmarkedFrameOccuCOP
+ expect_output(summary(umf), "unmarkedFrameOccuCOP Object")
+ expect_output(summary(umf_sub_ij), "unmarkedFrameOccuCOP Object")
+
+ # plot method for unmarkedFrameOccuCOP
+ expect_no_error(plot(umf))
+ expect_no_error(plot(umf_sub_ij))
+
+
+ # Input handling: covariates -------------------------------------------------
+ expect_no_error(umfCovs <- unmarkedFrameOccuCOP(
+ y = y,
+ L = L,
+ siteCovs = psiCovs,
+ obsCovs = lambdaCovs
+ ))
+ expect_output(print(umfCovs), "Data frame representation of unmarkedFrame object")
+ expect_output(summary(umfCovs), "unmarkedFrameOccuCOP Object")
+ expect_no_error(plot(umfCovs))
+
+ # Input handling: NA ---------------------------------------------------------
+
+ # NA should not pose problems when creating the unmarkedFrameOccuCOP object.
+ # The warnings and potential errors will be displayed when fitting the model.
+ # Except when y only contains NA: then there's an error.
+
+ ## NA in y
+ yNA <- y
+ yNA[1:2,] <- NA
+ yNA[3:5, 3:4] <- NA
+ yNA[,3] <- NA
+ expect_no_error(umfNA <- unmarkedFrameOccuCOP(y = yNA, L = L))
+ expect_output(print(umfNA), "Data frame representation of unmarkedFrame object")
+ expect_output(summary(umfNA), "unmarkedFrameOccuCOP Object")
+ expect_no_error(plot(umfNA))
+
+ ## NA in L
+ obsLengthNA <- L
+ obsLengthNA[1:2,] <- NA
+ obsLengthNA[3:5, 3:4] <- NA
+ obsLengthNA[,5] <- NA
+ expect_no_error(umfNA <- unmarkedFrameOccuCOP(y = y, L = obsLengthNA))
+ expect_output(print(umfNA), "Data frame representation of unmarkedFrame object")
+ expect_output(summary(umfNA), "unmarkedFrameOccuCOP Object")
+
+ expect_no_error(plot(umfNA))
+
+ ## NA also in covariates
+ psiCovsNA <- psiCovs
+ psiCovsNA[1:5,] <- NA
+ psiCovsNA[c(8,10,12), 3] <- NA
+ psiCovsNA[,1] <- NA
+ lambdaCovsNA <- lambdaCovs
+ lambdaCovsNA[[1]][1:5,] <- NA
+ lambdaCovsNA[[1]][,3] <- NA
+ lambdaCovsNA[[3]][,5] <- NA
+ expect_no_error(umfCovsNA <- unmarkedFrameOccuCOP(
+ y = yNA,
+ L = obsLengthNA,
+ siteCovs = psiCovsNA,
+ obsCovs = lambdaCovsNA
+ ))
+ expect_output(print(umfCovsNA), "Data frame representation of unmarkedFrame object")
+ expect_output(summary(umfCovsNA), "unmarkedFrameOccuCOP Object")
+ expect_no_error(plot(umfCovsNA))
+
+ ## all NA in y
+ yallNA <- y
+ yallNA[1:M, 1:J] <- NA
+ expect_error(unmarkedFrameOccuCOP(y = yallNA, L = L))
+
+ # Input handling: error and warnings -----------------------------------------
+ # Error when there are decimals in y
+ y_with_decimals = y
+ y_with_decimals[1, 1] = .5
+ expect_error(unmarkedFrameOccuCOP(y = y_with_decimals, L = L))
+
+ # Warning if y is detection/non-detection instead of count
+ y_detec_nodetec = (y > 0) * 1
+ expect_warning(unmarkedFrameOccuCOP(y = y_detec_nodetec, L = L))
+
+ # Error if the dimension of L is different than that of y
+ expect_error(unmarkedFrameOccuCOP(y = y, L = L[1:2, 1:2]))
+})
+
+
+test_that("occuCOP can fit simple models", {
+ # Simulate data
+ set.seed(123)
+ M = 100
+ J = 5
+ y = COPsimul(psi = .5,
+ lambda = 1,
+ M = M,
+ J = J)
+ L = y * 0 + 1
+
+ # Create umf
+ umf <- unmarkedFrameOccuCOP(y = y, L = L)
+
+ # Fitting options ----
+
+ ## With default parameters ----
+ expect_no_error(fit_default <- occuCOP(data = umf, L1 = TRUE))
+ expect_warning(occuCOP(data = umf, psiformula = ~ 1, lambdaformula = ~ 1, psistarts = 0, lambdastarts = 0))
+
+ ## With chosen starting points ----
+ expect_no_error(occuCOP(data = umf,
+ psiformula = ~ 1, lambdaformula = ~ 1,
+ psistarts = qlogis(.7),
+ lambdastarts = log(0.1), L1=T))
+ expect_no_error(occuCOP(data = umf,
+ psiformula = ~ 1, lambdaformula = ~ 1,
+ starts = c(qlogis(.7), log(0.1)), L1 = T))
+ # warning if all starts and psistarts and lambdastarts were furnished
+ # and starts != c(psistarts, lambdastarts)
+ expect_no_error(occuCOP(data = umf, starts = c(0, 0),
+ psistarts = c(0), lambdastarts = c(0), L1 = T))
+ expect_warning(occuCOP(data = umf, starts = c(0, 1),
+ psistarts = c(0), lambdastarts = c(0), L1 = T))
+ # errors if starting vectors of the wrong length
+ expect_error(occuCOP(data = umf, starts = c(0), L1 = T))
+ expect_error(occuCOP(data = umf, psistarts = c(0, 0), lambdastarts = 0, L1 = T))
+ expect_error(occuCOP(data = umf, lambdastarts = c(0, 0), L1 = T))
+
+ # With different options
+ expect_no_error(occuCOP(data = umf, method = "Nelder-Mead", psistarts = 0, lambdastarts = 0, L1=T))
+ expect_error(occuCOP(data = umf, method = "ABC", psistarts = 0, lambdastarts = 0, L1=T))
+
+ expect_no_error(occuCOP(data = umf, se = F, psistarts = 0, lambdastarts = 0, L1=T))
+ expect_error(occuCOP(data = umf, se = "ABC"))
+
+ expect_no_error(occuCOP(data = umf, engine = "R", psistarts = 0, lambdastarts = 0, L1=T))
+ expect_error(occuCOP(data = umf, engine = "julia", psistarts = 0, lambdastarts = 0, L1=T))
+
+ expect_no_error(occuCOP(data = umf, na.rm = F, psistarts = 0, lambdastarts = 0, L1=T))
+ expect_error(occuCOP(data = umf, na.rm = "no", psistarts = 0, lambdastarts = 0, L1=T))
+
+ # Looking at at COP model outputs ----
+ expect_is(fit_default, "unmarkedFitOccuCOP")
+ expect_equivalent(coef(fit_default), c(0.13067954, 0.06077929), tol = 1e-5)
+
+ ## backTransform
+ expect_no_error(backTransform(fit_default, type = "psi"))
+ expect_no_error(backTransform(fit_default, type = "lambda"))
+ expect_error(backTransform(fit_default, type = "state"))
+ expect_error(backTransform(fit_default, type = "det"))
+ expect_is(backTransform(fit_default, type = "psi"), "unmarkedBackTrans")
+
+ ## predict with newdata = fit@data
+ expect_no_error(umpredpsi <- predict(object = fit_default, type = "psi"))
+ expect_equal(umpredpsi$Predicted[1], 0.5326235, tol = 1e-5)
+ expect_no_error(umpredlambda <- predict(object = fit_default, type = "lambda"))
+ expect_equal(umpredlambda$Predicted[1], 1.062664, tol = 1e-5)
+ expect_error(predict(object = fit_default, type = "state"))
+
+ ## predict with newdata = 1
+ expect_no_error(predict(
+ object = fit_default,
+ newdata = data.frame(1),
+ type = "psi"
+ ))
+ expect_no_error(predict(
+ object = fit_default,
+ newdata = data.frame(1),
+ type = "lambda"
+ ))
+ expect_no_error(predict(
+ object = fit_default,
+ newdata = data.frame("a"=1:5,"b"=10:14),
+ type = "psi"
+ ))
+
+ # Fitting accurately ----
+ ## psi = 0.50, lambda = 1 ----
+ psi_test = .5
+ lambda_test = 1
+ fit_accur <- occuCOP(data = unmarkedFrameOccuCOP(
+ y = COPsimul(
+ psi = psi_test,
+ lambda = lambda_test,
+ M = 1000,
+ J = 10
+ ),
+ L = matrix(1, nrow = 1000, ncol = 10)
+ ), psistarts = 0, lambdastarts = 0, L1=T)
+ psi_estimate = backTransform(fit_accur, type = "psi")@estimate
+ lambda_estimate = backTransform(fit_accur, type = "lambda")@estimate
+ expect_equivalent(
+ psi_estimate,
+ psi_test,
+ tol = 0.05
+ )
+ expect_equivalent(
+ lambda_estimate,
+ lambda_test,
+ tol = 0.05
+ )
+
+ ## psi = 0.25, lambda = 5 ----
+ psi_test = 0.25
+ lambda_test = 5
+ fit_accur <- occuCOP(data = unmarkedFrameOccuCOP(
+ y = COPsimul(
+ psi = psi_test,
+ lambda = lambda_test,
+ M = 1000,
+ J = 10
+ ),
+ L = matrix(1, nrow = 1000, ncol = 10)
+ ), psistarts = 0, lambdastarts = 0, L1=T)
+ psi_estimate = backTransform(fit_accur, type = "psi")@estimate
+ lambda_estimate = backTransform(fit_accur, type = "lambda")@estimate
+ expect_equivalent(
+ psi_estimate,
+ psi_test,
+ tol = 0.05
+ )
+ expect_equivalent(
+ lambda_estimate,
+ lambda_test,
+ tol = 0.05
+ )
+
+ ## psi = 0.75, lambda = 0.5 ----
+ psi_test = 0.75
+ lambda_test = 0.5
+ fit_accur <- occuCOP(data = unmarkedFrameOccuCOP(
+ y = COPsimul(
+ psi = psi_test,
+ lambda = lambda_test,
+ M = 1000,
+ J = 10
+ ),
+ L = matrix(1, nrow = 1000, ncol = 10)
+ ), psistarts = 0, lambdastarts = 0, L1=T)
+ psi_estimate = backTransform(fit_accur, type = "psi")@estimate
+ lambda_estimate = backTransform(fit_accur, type = "lambda")@estimate
+ expect_equivalent(
+ psi_estimate,
+ psi_test,
+ tol = 0.05
+ )
+ expect_equivalent(
+ lambda_estimate,
+ lambda_test,
+ tol = 0.05
+ )
+
+ # With NAs ----
+ yNA <- y
+ yNA[1,] <- NA
+ yNA[3, 1] <- NA
+ yNA[4, 3] <- NA
+ yNA[, 5] <- NA
+ expect_no_error(umfNA <- unmarkedFrameOccuCOP(y = yNA, L = L))
+
+ expect_warning(fit_NA <- occuCOP(data = umfNA, psistarts = 0, lambdastarts = 0, L1=T))
+ expect_error(occuCOP(data = umfNA, psistarts = 0, lambdastarts = 0, na.rm = F))
+})
+
+test_that("We can simulate COP data", {
+
+ # From scratch ----
+
+ # With no covariates
+ expect_no_error(simulate(
+ "occuCOP",
+ formulas = list(psi = ~ 1, lambda = ~ 1),
+ coefs = list(
+ psi = c(intercept = 0),
+ lambda = c(intercept = 0)
+ ),
+ design = list(M = 100, J = 100)
+ ))
+
+ # With quantitative covariates
+ expect_no_error(simulate(
+ "occuCOP",
+ formulas = list(psi = ~ elev, lambda = ~ rain),
+ coefs = list(
+ psi = c(intercept = qlogis(.5), elev = -0.5),
+ lambda = c(intercept = log(3), rain = -1)
+ ),
+ design = list(M = 100, J = 5)
+ ))
+
+ # With guides
+ expect_no_error(simulate(
+ "occuCOP",
+ formulas = list(psi = ~ elev, lambda = ~ rain),
+ coefs = list(
+ psi = c(intercept = qlogis(.5), elev = -0.5),
+ lambda = c(intercept = log(3), rain = -1)
+ ),
+ design = list(M = 100, J = 5),
+ guide = list(elev=list(dist=rnorm, mean=12, sd=0.5))
+ ))
+
+ # With qualitative covariates
+ expect_no_error(umf <- simulate(
+ "occuCOP",
+ formulas = list(psi = ~ elev + habitat, lambda = ~ 1),
+ coefs = list(
+ psi = c(
+ intercept = qlogis(.2),
+ elev = -0.5,
+ habitatB = .5,
+ habitatC = .8
+ ),
+ lambda = c(intercept = log(3))
+ ),
+ design = list(M = 100, J = 5),
+ guide = list(habitat = factor(levels = c("A", "B", "C")))
+ ))
+
+ # From unmarkedFitOccuCOP ----
+ expect_no_error(umfit <- occuCOP(
+ umf,
+ psiformula = ~ habitat,
+ L1 = T,
+ psistarts = c(0,0,0),
+ lambdastarts = 0
+ ))
+ expect_no_error(simulate(umfit))
+})
+
+test_that("occuCOP can fit and predict models with covariates", {
+ # Simulate data with covariates ----
+ set.seed(123)
+ expect_no_error(umf <- simulate(
+ "occuCOP",
+ formulas = list(psi = ~ elev + habitat, lambda = ~ rain),
+ coefs = list(
+ psi = c(
+ intercept = qlogis(.2),
+ elev = -0.5,
+ habitatB = .5,
+ habitatC = .8
+ ),
+ lambda = c(intercept = log(3), rain = -1)
+ ),
+ design = list(M = 100, J = 5),
+ guide = list(habitat = factor(levels = c("A", "B", "C")))
+ ))
+
+ # Fit ----
+ expect_no_error(umfit <- occuCOP(
+ umf,
+ psiformula = ~ habitat + elev,
+ lambdaformula = ~ rain,
+ L1 = T,
+ psistarts = c(0,0,0,0),
+ lambdastarts = c(0,0)
+ ))
+
+ expect_error(occuCOP(
+ umf,
+ psiformula = ~ habitat+elev,
+ lambdaformula = ~ rain,
+ L1 = T,
+ psistarts = c(0),
+ lambdastarts = c(0,0)
+ ))
+
+ expect_equivalent(
+ coef(umfit),
+ c(-1.5350679, 0.4229763, 0.7398768, -1.0456397, 1.2333424, -0.8344109),
+ tol = 1e-5
+ )
+
+ # Predict ----
+ expect_no_error(predict(umfit, type = "psi"))
+ expect_no_error(umpredpsi <- predict(
+ umfit,
+ type = "psi",
+ newdata = data.frame("habitat" = c("A", "B", "C"), "elev" = c(0, 0, 0)),
+ appendData = TRUE
+ ))
+ expect_equivalent(umpredpsi$Predicted, c(0.1772534, 0.2474811, 0.3110551), tol = 1e-5)
+
+ expect_no_error(umpredlambda <- predict(umfit, type = "lambda", appendData = TRUE))
+ expect_no_error(predict(umfit, type = "lambda", level = 0.5))
+ expect_equal(umpredlambda$Predicted[1], 1.092008, tol = 1e-5)
+})
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-02 15:35:54 +0000