diff --git a/R/get_Z.R b/R/get_Z.R index b9ec4ce..d8e6b3f 100644 --- a/R/get_Z.R +++ b/R/get_Z.R @@ -416,6 +416,10 @@ forward_sim <- function(init, colo, ex, sample = FALSE){ assertthat::assert_that(is.numeric(init)) assertthat::assert_that(length(init) == 1) assertthat::assert_that(init >= 0 & init <= 1) + assertthat::assert_that(all(colo >= 0, na.rm = TRUE)) + assertthat::assert_that(all(colo <= 1, na.rm = TRUE)) + assertthat::assert_that(all(ex >= 0, na.rm = TRUE)) + assertthat::assert_that(all(ex <= 1, na.rm = TRUE)) out <- rep(NA, length_out) if(sample){ diff --git a/R/single_C-prediction-epred.R b/R/single_C-prediction-epred.R index 8033149..d9c3fcd 100644 --- a/R/single_C-prediction-epred.R +++ b/R/single_C-prediction-epred.R @@ -4,6 +4,9 @@ #' @param prep Output of \code{brms::prepare_predictions}. See brms custom #' families vignette at #' https://cran.r-project.org/web/packages/brms/vignettes/brms_customfamilies.html +#' @param ... unused additional arguments. See brms custom +#' families vignette at +#' https://cran.r-project.org/web/packages/brms/vignettes/brms_customfamilies.html #' @return Posterior predictions posterior_predict_occupancy_single_C <- function(i, prep, ...) { mu <- brms::get_dpar(prep, "mu", i = i) diff --git a/man/get_Z.Rd b/man/get_Z.Rd index 1809b5f..c652c42 100644 --- a/man/get_Z.Rd +++ b/man/get_Z.Rd @@ -33,7 +33,8 @@ or bernoulli samples from those probabilities (TRUE)} \item{new_data}{Optional new data at which to predict the Z matrix. Can be the output of `make_flocker_data` or the `unit_covs` input to `make_flocker_data` provided that `history_condition` is `FALSE` and the -occupancy model is for a single season.} +occupancy model is a single-season model (rep-constant, rep-varying, or +data-augmented).} \item{allow_new_levels}{allow new levels for random effect terms in `new_data`? Will error if set to `FALSE` and new levels are provided in `new_data`.} diff --git a/man/log_lik_occupancy_single_C.Rd b/man/log_lik_occupancy_single_C.Rd index 1f25b33..1e1edd2 100644 --- a/man/log_lik_occupancy_single_C.Rd +++ b/man/log_lik_occupancy_single_C.Rd @@ -3,12 +3,12 @@ \name{log_lik_occupancy_single_C} \alias{log_lik_occupancy_single_C} \title{A log-likelihood function for the rep-constant occupancy model, sufficient for -\code{brms::loo(vc_fit)} to work.} +\code{brms::loo()} to work.} \usage{ log_lik_occupancy_single_C(i, prep) } \arguments{ -\item{i}{Posterior iteration} +\item{i}{Observation id} \item{prep}{Output of \code{brms::prepare_predictions}. See brms custom families vignette at @@ -19,5 +19,5 @@ The log-likelihood for observation i } \description{ A log-likelihood function for the rep-constant occupancy model, sufficient for -\code{brms::loo(vc_fit)} to work. +\code{brms::loo()} to work. } diff --git a/man/posterior_epred_occupancy_single_C.Rd b/man/posterior_epred_occupancy_single_C.Rd new file mode 100644 index 0000000..164b5e8 --- /dev/null +++ b/man/posterior_epred_occupancy_single_C.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/single_C-prediction-epred.R +\name{posterior_epred_occupancy_single_C} +\alias{posterior_epred_occupancy_single_C} +\title{A posterior epred function for the rep-constant occupancy model, sufficient +for \code{brms::posterior_epred()} to work.} +\usage{ +posterior_epred_occupancy_single_C(prep) +} +\arguments{ +\item{prep}{Output of \code{brms::prepare_predictions}. See brms custom +families vignette at +https://cran.r-project.org/web/packages/brms/vignettes/brms_customfamilies.html} +} +\value{ +Posterior epreds +} +\description{ +A posterior epred function for the rep-constant occupancy model, sufficient +for \code{brms::posterior_epred()} to work. +} diff --git a/man/posterior_predict_occupancy_single_C.Rd b/man/posterior_predict_occupancy_single_C.Rd new file mode 100644 index 0000000..9fe1c5b --- /dev/null +++ b/man/posterior_predict_occupancy_single_C.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/single_C-prediction-epred.R +\name{posterior_predict_occupancy_single_C} +\alias{posterior_predict_occupancy_single_C} +\title{A posterior predict function for the rep-constant occupancy model, sufficient +for \code{brms::predict.brmsfit()} to work.} +\usage{ +posterior_predict_occupancy_single_C(i, prep, ...) +} +\arguments{ +\item{i}{Observation id} + +\item{prep}{Output of \code{brms::prepare_predictions}. See brms custom +families vignette at +https://cran.r-project.org/web/packages/brms/vignettes/brms_customfamilies.html} + +\item{...}{unused additional arguments. See brms custom +families vignette at +https://cran.r-project.org/web/packages/brms/vignettes/brms_customfamilies.html} +} +\value{ +Posterior predictions +} +\description{ +A posterior predict function for the rep-constant occupancy model, sufficient +for \code{brms::predict.brmsfit()} to work. +} diff --git a/man/predict_flocker.Rd b/man/predict_flocker.Rd index 1e9fa45..7c42d94 100644 --- a/man/predict_flocker.Rd +++ b/man/predict_flocker.Rd @@ -29,7 +29,11 @@ exclusively by the posterior distribution for the occupancy probability psi.} \item{new_data}{Optional new data at which to predict. If `NULL`, predictions -are given at the data points used for model fitting ("retrodictions")} +are given at the data points used for model fitting. +Can be the output of 'make_flocker_data' or the 'unit_covs' input to +'make_flocker_data' provided that 'history_condition' is 'FALSE' and the +occupancy model is a single-season model (rep-constant, rep-varying, or +data-augmented).} \item{mixed}{When `new_data` is not provided, should random effect levels be drawn from their posteriors (`FALSE`, the default) or re-sampled from @@ -46,11 +50,9 @@ levels not present in the original data, how should predictions be handled? Passed directly to brms::prepare_predictions, which see.} } \value{ -An array of posterior predictions. For a single-season rep-varying - model, a 3-dimensional array where the first dimension is the closure-unit, - the second dimension is the rep, the third dimension is the iteration, - and the value is 1, 0, or NA indicating detection, non-detection, or - non-existence of the sampling event. +An array of posterior predictions in the same shape as the + observations passed to `make_flocker_data()` with posterior iterations + stacked along the final dimension. } \description{ Get posterior predictions from a flocker model diff --git a/tests/testthat/test-get_Z.R b/tests/testthat/test-get_Z.R new file mode 100644 index 0000000..5a4cfce --- /dev/null +++ b/tests/testthat/test-get_Z.R @@ -0,0 +1,185 @@ +test_that("get_Z gives valid returns", { + Z_single <- get_Z(example_flocker_model_single) + expect_true(all(Z_single >= 0)) + expect_true(all(Z_single <= 1)) + expect_true(any(Z_single == 1)) + + Z_single_C <- get_Z(example_flocker_model_single_C) + expect_true(all(Z_single_C >= 0)) + expect_true(all(Z_single_C <= 1)) + expect_true(any(Z_single_C == 1)) + + Z_augmented <- get_Z(example_flocker_model_aug) + expect_true(all(Z_augmented >= 0)) + expect_true(all(Z_augmented <= 1)) + expect_true(any(Z_augmented == 1)) + + Z_multi_colex_ex <- get_Z(example_flocker_model_multi_colex_ex) + expect_true(all(Z_multi_colex_ex >= 0, na.rm = T)) + expect_true(all(Z_multi_colex_ex <= 1, na.rm = T)) + expect_true(any(Z_multi_colex_ex == 1, na.rm = T)) + + Z_multi_colex_eq <- get_Z(example_flocker_model_multi_colex_eq) + expect_true(all(Z_multi_colex_eq >= 0, na.rm = T)) + expect_true(all(Z_multi_colex_eq <= 1, na.rm = T)) + expect_true(any(Z_multi_colex_eq == 1, na.rm = T)) + + Z_multi_auto_ex <- get_Z(example_flocker_model_multi_auto_ex) + expect_true(all(Z_multi_auto_ex >= 0, na.rm = T)) + expect_true(all(Z_multi_auto_ex <= 1, na.rm = T)) + expect_true(any(Z_multi_auto_ex == 1, na.rm = T)) + + Z_multi_auto_eq <- get_Z(example_flocker_model_multi_auto_eq) + expect_true(all(Z_multi_auto_eq >= 0, na.rm = T)) + expect_true(all(Z_multi_auto_eq <= 1, na.rm = T)) + expect_true(any(Z_multi_auto_eq == 1, na.rm = T)) + + + + Z_single_nohist <- get_Z(example_flocker_model_single, history_condition = FALSE) + expect_true(all(Z_single_nohist >= 0)) + expect_true(all(Z_single_nohist <= 1)) + expect_true(sum(Z_single == 1) > sum(Z_single_nohist == 1)) + + Z_single_C_nohist <- get_Z(example_flocker_model_single_C, history_condition = FALSE) + expect_true(all(Z_single_C_nohist >= 0)) + expect_true(all(Z_single_C_nohist <= 1)) + expect_true(sum(Z_single_C == 1) > sum(Z_single_C_nohist == 1)) + + Z_augmented_nohist <- get_Z(example_flocker_model_aug, history_condition = FALSE) + expect_true(all(Z_augmented_nohist >= 0)) + expect_true(all(Z_augmented_nohist <= 1)) + expect_true(sum(Z_augmented == 1) > sum(Z_augmented_nohist == 1)) + + Z_multi_colex_ex_nohist <- get_Z(example_flocker_model_multi_colex_ex, history_condition = FALSE) + expect_true(all(Z_multi_colex_ex_nohist >= 0, na.rm = T)) + expect_true(all(Z_multi_colex_ex_nohist <= 1, na.rm = T)) + expect_true(sum(Z_multi_colex_ex == 1, na.rm = TRUE) > sum(Z_multi_colex_ex_nohist == 1, na.rm = TRUE)) + + Z_multi_colex_eq_nohist <- get_Z(example_flocker_model_multi_colex_eq, history_condition = FALSE) + expect_true(all(Z_multi_colex_eq_nohist >= 0, na.rm = T)) + expect_true(all(Z_multi_colex_eq_nohist <= 1, na.rm = T)) + expect_true(sum(Z_multi_colex_eq == 1, na.rm = TRUE) > sum(Z_multi_colex_eq_nohist == 1, na.rm = TRUE)) + + Z_multi_auto_ex_nohist <- get_Z(example_flocker_model_multi_auto_ex, history_condition = FALSE) + expect_true(all(Z_multi_auto_ex_nohist >= 0, na.rm = T)) + expect_true(all(Z_multi_auto_ex_nohist <= 1, na.rm = T)) + expect_true(sum(Z_multi_auto_ex == 1, na.rm = TRUE) > sum(Z_multi_auto_ex_nohist == 1, na.rm = TRUE)) + + Z_multi_auto_eq_nohist <- get_Z(example_flocker_model_multi_auto_eq, history_condition = FALSE) + expect_true(all(Z_multi_auto_eq_nohist >= 0, na.rm = T)) + expect_true(all(Z_multi_auto_eq_nohist <= 1, na.rm = T)) + expect_true(sum(Z_multi_auto_eq == 1, na.rm = TRUE) > sum(Z_multi_auto_eq_nohist == 1, na.rm = TRUE)) + +}) + + +test_that("get_Z with sampling gives valid returns", { + Z_single <- get_Z(example_flocker_model_single, sample = TRUE) + expect_true(all(Z_single %in% c(0, 1, NA))) + + Z_single_C <- get_Z(example_flocker_model_single_C, sample = TRUE) + expect_true(all(Z_single_C %in% c(0, 1, NA))) + + Z_augmented <- get_Z(example_flocker_model_aug, sample = TRUE) + expect_true(all(Z_augmented %in% c(0, 1, NA))) + + Z_multi_colex_ex <- get_Z(example_flocker_model_multi_colex_ex, sample = TRUE) + expect_true(all(Z_multi_colex_ex %in% c(0, 1, NA))) + + Z_multi_colex_eq <- get_Z(example_flocker_model_multi_colex_eq, sample = TRUE) + expect_true(all(Z_multi_colex_eq %in% c(0, 1, NA))) + + Z_multi_auto_ex <- get_Z(example_flocker_model_multi_auto_ex, sample = TRUE) + expect_true(all(Z_multi_auto_ex %in% c(0, 1, NA))) + + Z_multi_auto_eq <- get_Z(example_flocker_model_multi_auto_eq, sample = TRUE) + expect_true(all(Z_multi_auto_eq %in% c(0, 1, NA))) + +}) + +test_that("new_data works as expected", { + fd1 <- simulate_flocker_data() + mfd1 <- make_flocker_data(fd1$obs, fd1$unit_covs, fd1$event_covs) + expect_silent(get_Z(example_flocker_model_single, new_data = mfd1)) +# expect_silent(get_Z(example_flocker_model_single, history_condition = FALSE, new_data = fd1$unit_covs)) + expect_error(get_Z(example_flocker_model_single, history_condition = TRUE, new_data = fd1$unit_covs)) +}) + + +test_that("incorrect types cause an error", { + expect_error(forward_sim("init", c(0.5, 0.5), c(0.5, 0.5))) + expect_error(forward_sim(0.5, "colo", c(0.5, 0.5))) + expect_error(forward_sim(0.5, c(0.5, 0.5), "ex")) +}) + +test_that("negative probabilities cause an error", { + expect_silent(forward_sim(.4, c(.5, .5), c(.5, .5))) + expect_error(forward_sim(.4, c(.5, NA), c(.5, .5))) + expect_error(forward_sim(-0.1, c(0.5, .5), c(0.5, .5))) + expect_error(forward_sim(0.5, c(.5, -0.5), c(.5, 0.5))) + expect_error(forward_sim(0.5, c(NA, 0.5), c(NA, -0.5))) +}) + +test_that("correct inputs produce correct outputs", { + # Simple case with no NA values and no sampling + expect_equal(forward_sim(0.5, c(0.5, 0.5), c(0.5, 0.5)), c(0.5, 0.5)) + # Simple case with sampling, the result is stochastic, so we can't test exact values, just types and lengths + result <- forward_sim(0.5, c(0.5, 0.5), c(0.5, 0.5), sample = TRUE) + expect_type(result, "integer") + expect_length(result, 2) +}) + +test_that("edge cases are handled", { + # Case with all NA values should return NA + expect_error(forward_sim(NA, c(NA, NA), c(NA, NA))) + # Case with empty vectors + expect_error(forward_sim(0.5, numeric(0), numeric(0))) +}) + +# Test for forward_backward_algorithm function +test_that("forward_backward_algorithm returns correct results", { + el0 <- c(0.1, 0.2, 0.3) + el1 <- c(0.9, 0.8, 0.7) + init <- 0.5 + colo <- c(0.6, 0.7, 0.8) + ex <- c(0.4, 0.3, 0.2) + result <- forward_backward_algorithm(el0, el1, init, colo, ex) + expect_type(result, "double") + expect_length(result, length(el0)) +}) + +# Test for forward_backward_sampling function +test_that("forward_backward_sampling returns correct results", { + el0 <- c(0.1, 0.2, 0.3) + el1 <- c(0.9, 0.8, 0.7) + init <- 0.5 + colo <- c(0.6, 0.7, 0.8) + ex <- c(0.4, 0.3, 0.2) + result <- forward_backward_sampling(el0, el1, init, colo, ex) + expect_type(result, "double") + expect_length(result, length(el0)) +}) + +# Test for forward_algorithm function +test_that("forward_algorithm returns correct results", { + el0 <- c(0.1, 0.2, 0.3) + el1 <- c(0.9, 0.8, 0.7) + init <- 0.5 + colo <- c(0.6, 0.7, 0.8) + ex <- c(0.4, 0.3, 0.2) + result <- forward_algorithm(el0, el1, init, colo, ex) + expect_type(result, "double") + expect_equal(dim(result), c(length(el0), 2)) +}) + +# Test for backward_algorithm function +test_that("backward_algorithm returns correct results", { + el0 <- c(0.1, 0.2, 0.3) + el1 <- c(0.9, 0.8, 0.7) + colo <- c(0.6, 0.7, 0.8) + ex <- c(0.4, 0.3, 0.2) + result <- backward_algorithm(el0, el1, colo, ex) + expect_type(result, "double") + expect_equal(dim(result), c(length(el0), 2)) +})