diff options
author | Ken Kellner <kenkellner@users.noreply.github.com> | 2022-10-18 12:57:36 -0400 |
---|---|---|
committer | GitHub <noreply@github.com> | 2022-10-18 12:57:36 -0400 |
commit | 1cc3e845c8d610512da7402edceb8d103154216b (patch) | |
tree | f460cd37331cfb0969b731c998fc05cd03b5cf98 | |
parent | 5063b596ad880b1cba772090c1e48eb13a7c7bdd (diff) | |
parent | bb3198f1a214a58b42ecb4ad4222f376088a641f (diff) |
Merge pull request #241 from kenkellner/parboot_update
Add error handling to parboot
-rw-r--r-- | DESCRIPTION | 4 | ||||
-rw-r--r-- | R/boot.R | 135 | ||||
-rw-r--r-- | tests/testthat/test_parboot.R | 44 |
3 files changed, 101 insertions, 82 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 1bbd46b..a878727 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: unmarked -Version: 1.2.5.9006 -Date: 2022-08-24 +Version: 1.2.5.9007 +Date: 2022-10-12 Type: Package Title: Models for Data from Unmarked Animals Authors@R: c( @@ -30,81 +30,66 @@ setMethod("replaceY", "unmarkedFrameOccuMulti", object }) -setMethod("parboot", "unmarkedFit", - function(object, statistic=SSE, nsim=10, report, seed = NULL, parallel = TRUE, ncores, ...) -{ - dots <- list(...) - statistic <- match.fun(statistic) - call <- match.call(call = sys.call(-1)) - formula <- object@formula - umf <- getData(object) - y <- getY(object) - ests <- as.numeric(coef(object)) - starts <- ests - if(methods::.hasSlot(object, "TMB") && !is.null(object@TMB)) starts <- NULL - t0 <- statistic(object, ...) - lt0 <- length(t0) - t.star <- matrix(NA, nsim, lt0) - if(!missing(report)) - cat("t0 =", t0, "\n") - simdata <- umf - if (!is.null(seed)) set.seed(seed) - simList <- simulate(object, nsim = nsim, na.rm = FALSE) - availcores <- detectCores() - if(missing(ncores)) ncores <- availcores - 1 - if(ncores > availcores) ncores <- availcores - - no_par <- ncores < 2 || nsim < 100 || !parallel - - if (no_par) { - if (!missing(report)) { - for(i in 1:nsim) { - simdata <- replaceY(simdata, simList[[i]]) - fit <- update(object, data=simdata, starts=starts, se=FALSE) - t.star[i,] <- statistic(fit, ...) - if(!missing(report)) { - if (nsim > report && i %in% seq(report, nsim, by=report)) - cat("iter", i, ": ", t.star[i, ], "\n") - } - } - } else { - t.star <- pbsapply(1:nsim, function(i) { - simdata <- replaceY(simdata, simList[[i]]) - fit <- update(object, data=simdata, starts=starts, se=FALSE) - t.star.tmp <- statistic(fit, ...) - }) - if (lt0 > 1) - t.star <- t(t.star) - else - t.star <- matrix(t.star, ncol = lt0) - } - } else { - message("Running parametric bootstrap in parallel on ", ncores, " cores.") - if (!missing(report)) message("Bootstrapped statistics not reported during parallel processing.") - cl <- makeCluster(ncores) - if (!is.null(seed)) parallel::clusterSetRNGStream(cl, iseed = seed) - on.exit(stopCluster(cl)) - varList <- c("simList", "y", "object", "simdata", "starts", "statistic", "dots") - # If call formula is an object, include it too - fm.nms <- all.names(object@call) - if (!any(grepl("~", fm.nms))) varList <- c(varList, fm.nms[2]) - ## Hack to get piFun for unmarkedFitGMM and unmarkedFitMPois - if(.hasSlot(umf, "piFun")) varList <- c(varList, umf@piFun) - clusterExport(cl, varList, envir = environment()) - clusterEvalQ(cl, library(unmarked)) - clusterEvalQ(cl, list2env(dots)) - t.star.parallel <- pblapply(1:nsim, function(i) { - simdata <- replaceY(simdata, simList[[i]]) - fit <- update(object, data = simdata, starts = starts, se = FALSE) - t.star <- statistic(fit, ...) - }, cl = cl) - t.star <- matrix(unlist(t.star.parallel), nrow = length(t.star.parallel), byrow = TRUE) - } - if (!is.null(names(t0))) - colnames(t.star) <- names(t0) - else colnames(t.star) <- paste("t*", 1:lt0, sep="") - out <- new("parboot", call = call, t0 = t0, t.star = t.star) - return(out) + +setMethod("parboot", "unmarkedFit", function(object, statistic=SSE, nsim=10, + report, seed = NULL, parallel = FALSE, ncores, ...){ + + if(!missing(report)){ + warning("report argument is non-functional and will be deprecated in the next version", call.=FALSE) + } + + dots <- list(...) + call <- match.call(call = sys.call(-1)) + stopifnot(is.function(statistic)) + starts <- as.numeric(coef(object)) + # Get rid of starting values if model was fit with TMB + if(methods::.hasSlot(object, "TMB") && !is.null(object@TMB)) starts <- NULL + + t0 <- statistic(object, ...) + + simList <- simulate(object, nsim = nsim, na.rm = FALSE) + + availcores <- parallel::detectCores() - 1 + if(missing(ncores) || ncores > availcores) ncores <- availcores + + cl <- NULL + if(parallel){ + cl <- parallel::makeCluster(ncores) + on.exit(parallel::stopCluster(cl)) + parallel::clusterEvalQ(cl, library(unmarked)) + env_vars <- c("dots", "replaceY") + fm.nms <- all.names(object@call) + if (!any(grepl("~", fm.nms))) env_vars <- c(env_vars, fm.nms[2]) + if(.hasSlot(object@data, "piFun")) env_vars <- c(env_vars, object@data@piFun) + parallel::clusterExport(cl, env_vars, envir = environment()) + parallel::clusterEvalQ(cl, list2env(dots)) + } + + run_sim <- function(x, object, statistic, starts, t0, ...){ + simdata <- replaceY(object@data, x) + tryCatch({ + #if(runif(1,0,1) < 0.5) stop("fail") # for testing error trapping + fit <- update(object, data=simdata, starts=starts, se=FALSE) + statistic(fit, ...) + }, error=function(e){ + t0[] <- NA + t0 + }) + } + + t.star <- t(pbapply::pbsapply(simList, run_sim, object=object, + statistic=statistic, starts=starts, t0=t0, + cl=cl, ...)) + if(length(t0) == 1) t.star <- matrix(t.star, ncol=1) + + failed <- apply(t.star, 1, function(x) any(is.na(x))) + if(sum(failed) > 0){ + warning(paste0("Model fitting failed in ",sum(failed), " sims."), call.=FALSE) + t.star <- t.star[!failed,,drop=FALSE] + } + + new("parboot", call = call, t0 = t0, t.star = t.star) + }) diff --git a/tests/testthat/test_parboot.R b/tests/testthat/test_parboot.R index c3ece25..af62658 100644 --- a/tests/testthat/test_parboot.R +++ b/tests/testthat/test_parboot.R @@ -32,15 +32,49 @@ test_that("parboot works", { dev.off() expect_equal(pl, NULL) - # check that report works - rep_output <- capture.output(parboot(fm, fitstats, nsim=3, report=TRUE)) - expect_equal(substr(rep_output[1], 1,2), "t0") + # check that report argument gives warning + expect_warning(parboot(fm, fitstats, nsim=3, report=TRUE)) }) test_that("parboot works in parallel",{ skip_on_cran() skip_on_ci() # check parallel - pb <- parboot(fm, nsim=101, parallel=TRUE, ncores=2) - expect_equal(length(pb@t.star), 101) + pb <- parboot(fm, nsim=10, parallel=TRUE, ncores=2) + expect_equal(length(pb@t.star), 10) +}) + +test_that("parboot handles failing model fits", { + + fail_func <- function(x){ + rand <- rnorm(1) + if(rand > 0.5){ + stop("fail") + } + return(rand) + } + + set.seed(123) + expect_warning(pb <- parboot(fm, nsim=20, statistic=fail_func)) + expect_equal(nrow(pb@t.star), 13) + + expect_warning(pb <- parboot(fm, nsim=20, statistic=fail_func, parallel=TRUE)) + expect_true(nrow(pb@t.star) < 20) + +}) + +test_that("parboot handles statistic functions with additional arguments", { + + opt_func <- function(x, y){ + res <- mean(residuals(x), na.rm=TRUE) + c(res=res, y=y) + } + + pb <- parboot(fm, nsim=10, statistic=opt_func, y=0.1) + expect_equal(colnames(pb@t.star), c("res", "y")) + expect_true(all(pb@t.star[,"y"]==0.1)) + + pb <- parboot(fm, nsim=10, statistic=opt_func, y=0.1, parallel=TRUE) + expect_equal(colnames(pb@t.star), c("res", "y")) + expect_true(all(pb@t.star[,"y"]==0.1)) }) |