diff options
| author | Ken Kellner <kenkellner@users.noreply.github.com> | 2022-11-02 13:43:45 -0400 |
|---|---|---|
| committer | GitHub <noreply@github.com> | 2022-11-02 13:43:45 -0400 |
| commit | 38c0806179d15388443555578b83bc367ebabb89 (patch) | |
| tree | 518761b944a3a126230946f802c88a6d45c7908b | |
| parent | 06f15cef6579b931307fb5bd6fb386d63abfb6b5 (diff) | |
| parent | 971eb1127559d48f78d3a292b6321926b97e5bca (diff) | |
Merge pull request #65 from kenkellner/kfold
Add K-fold cross-validation method
45 files changed, 981 insertions, 76 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 105f6d4..b272a4b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ubms -Version: 1.1.9006 -Date: 2022-09-20 +Version: 1.1.9007 +Date: 2022-11-02 Title: Bayesian Models for Data from Unmarked Animals using 'Stan' Authors@R: person("Ken", "Kellner", email="contact@kenkellner.com", role=c("cre","aut")) @@ -12,7 +12,7 @@ Imports: gridExtra, lme4, loo, - Matrix (>= 1.5-0), + Matrix, methods, pbapply, Rcpp (>= 0.12.0), @@ -34,7 +34,7 @@ License: GPL (>=3) URL: https://kenkellner.com/ubms/ BugReports: https://github.com/kenkellner/ubms/issues Encoding: UTF-8 -RoxygenNote: 7.1.2 +RoxygenNote: 7.2.0 Biarch: true LinkingTo: BH (>= 1.66.0), @@ -60,11 +60,13 @@ Collate: 'missing.R' 'distsamp.R' 'fitlist.R' - 'occuRN.R' - 'mb_chisq.R' - 'multinomPois.R' + 'kfold.R' 'occuTTD.R' + 'multinomPois.R' + 'occuRN.R' 'pcount.R' + 'loglik.R' + 'mb_chisq.R' 'plot_effects.R' 'plot_posteriors.R' 'predict.R' @@ -2,6 +2,7 @@ export(RSR) export(cauchy) +export(extract_log_lik) export(gamma) export(gof) export(laplace) @@ -29,6 +30,7 @@ exportMethods(fitted) exportMethods(getY) exportMethods(get_elapsed_time) exportMethods(get_stancode) +exportMethods(kfold) exportMethods(loo) exportMethods(modSel) exportMethods(names) @@ -81,6 +83,7 @@ importFrom(graphics,plot) importFrom(grid,gpar) importFrom(grid,textGrob) importFrom(gridExtra,grid.arrange) +importFrom(loo,kfold) importFrom(loo,loo) importFrom(loo,loo_compare) importFrom(loo,loo_model_weights) diff --git a/R/RcppExports.R b/R/RcppExports.R index a0cd5ce..ae97b81 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,6 +9,30 @@ exp_counts_occuRN <- function(obs, Kmin, lam, r) { .Call(`_ubms_exp_counts_occuRN`, obs, Kmin, lam, r) } +get_loglik_occu <- function(y, M, si, psimat, pmat, Kmin) { + .Call(`_ubms_get_loglik_occu`, y, M, si, psimat, pmat, Kmin) +} + +get_loglik_pcount <- function(y, M, si, lammat, pmat, K, Kmin) { + .Call(`_ubms_get_loglik_pcount`, y, M, si, lammat, pmat, K, Kmin) +} + +get_loglik_occuRN <- function(y, M, si, lammat, rmat, K, Kmin) { + .Call(`_ubms_get_loglik_occuRN`, y, M, si, lammat, rmat, K, Kmin) +} + +get_loglik_multinomPois <- function(y, M, si, lammat, pmat, pi_type) { + .Call(`_ubms_get_loglik_multinomPois`, y, M, si, lammat, pmat, pi_type) +} + +get_loglik_colext <- function(y, M, Tsamp, J, si, psicube, phicube, pmat, nd) { + .Call(`_ubms_get_loglik_colext`, y, M, Tsamp, J, si, psicube, phicube, pmat, nd) +} + +get_loglik_occuTTD <- function(y, M, si, psimat, lammat, k, delta, ydist) { + .Call(`_ubms_get_loglik_occuTTD`, y, M, si, psimat, lammat, k, delta, ydist) +} + simz_pcount <- function(y, lam_post, p_post, K, Kmin, kvals) { .Call(`_ubms_simz_pcount`, y, lam_post, p_post, K, Kmin, kvals) } @@ -26,6 +26,9 @@ #' @param prior_coef_det Prior distribution for the regression coefficients of #' the detection model #' @param prior_sigma Prior distribution on random effect standard deviations +#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values +#' in the output. This can greatly increase the size of the model. If +#' \code{FALSE}, the values are calculated post-hoc from the posteriors #' @param ... Arguments passed to the \code{\link{stan}} call, such as #' number of chains \code{chains} or iterations \code{iter} #' @@ -60,6 +63,7 @@ stan_colext <- function(psiformula = ~1, prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ...){ umf <- process_umf(data) @@ -80,7 +84,7 @@ stan_colext <- function(psiformula = ~1, prior_intercept_det, prior_coef_det, prior_sigma) submodels <- ubmsSubmodelList(state, col, ext, det) - ubmsFit("colext", match.call(), data, response, submodels, ...) + ubmsFit("colext", match.call(), data, response, submodels, log_lik, ...) } @@ -19,20 +19,24 @@ setClass("ubmsFitOccu", contains = "ubmsFit") # Child class for abundance/N-mixture type models setClass("ubmsFitAbun", contains = "ubmsFit") -ubmsFit <- function(model, call, data, response, submodels, ...){ +ubmsFit <- function(model, call, data, response, submodels, log_lik=TRUE, ...){ #Find missing response <- update_missing(response, submodels) submodels <- update_missing(submodels, response) #Fit model - fit <- fit_model(model, response, submodels, ...) + fit <- fit_model(model, response, submodels, log_lik, ...) + # Exit early if just returning Stan inputs + if(check_return_inputs(...)) return(fit) #Remove placeholder submodels submodels <- remove_placeholders(submodels) #Construct output - new(fit_class(model), call=call, data=data, stanfit=fit, - response=response, submodels=submodels, loo=get_loo(fit)) + out <- new(fit_class(model), call=call, data=data, stanfit=fit, + response=response, submodels=submodels, loo=empty_loo()) + out@loo <- get_loo(out) + out } remove_placeholders <- function(submodels){ @@ -41,6 +45,29 @@ remove_placeholders <- function(submodels){ submodels } +empty_loo <- function(){ + out <- list() + class(out) <- "psis_loo" + out +} + +get_loo <- function(object, cores=getOption("mc.cores", 1)){ + loglik <- extract_log_lik(object, merge_chains=FALSE) + r_eff <- loo::relative_eff(exp(loglik), cores=cores) + loo::loo(loglik, r_eff=r_eff, cores=cores) +} + +# Function to check if just Stan inputs should be returned +# Used by kfold method +# Pass return_inputs=TRUE in ... when calling stan_* +check_return_inputs <- function(...){ + args <- list(...) + if("return_inputs" %in% names(args)){ + if(args$return_inputs) return(TRUE) + } + FALSE +} + fit_class <- function(mod){ cap <- paste0(toupper(substr(mod,1,1)), substr(mod,2,nchar(mod))) paste0("ubmsFit",cap) @@ -48,9 +75,14 @@ fit_class <- function(mod){ #Fit stan model #' @include inputs.R -fit_model <- function(name, response, submodels, ...){ +fit_model <- function(name, response, submodels, log_lik, ...){ model <- name_to_stanmodel(name, submodels) - inp <- build_stan_inputs(name, response, submodels) + inp <- build_stan_inputs(name, response, submodels, log_lik) + # Should just Stan inputs be returned? + if(check_return_inputs(...)){ + inp$submodels <- submodels + return(inp) + } mod <- stanmodels[[model]] mod@model_name <- name fit <- sampling(mod, data=inp$stan_data, pars=inp$pars, ...) @@ -92,9 +124,3 @@ setMethod("stanfit_names", "ubmsSubmodelList", function(object, ...){ out <- unname(out) out }) - -get_loo <- function(stanfit, cores=getOption("mc.cores", 1)){ - loglik <- loo::extract_log_lik(stanfit, merge_chains=FALSE) - r_eff <- loo::relative_eff(exp(loglik), cores=cores) - loo::loo(loglik, r_eff=r_eff, cores=cores) -} diff --git a/R/fitlist.R b/R/fitlist.R index 5da3889..ce8ae27 100644 --- a/R/fitlist.R +++ b/R/fitlist.R @@ -49,7 +49,8 @@ setMethod("fitList", "list", function(...){ #' @importFrom unmarked modSel #' @export setMethod("modSel", "ubmsFitList", function(object, ...){ - loos <- lapply(object@models, loo, ...) + #loos <- lapply(object@models, loo, ...) + loos <- lapply(object@models, function(x) x@loo) elpd <- sapply(loos, function(x) x$estimates[1]) p_loo <- sapply(loos, function(x) x$estimates[2]) compare <- loo::loo_compare(loos)[names(elpd),] @@ -1,9 +1,9 @@ - build_stan_inputs <- function(name, response, submodels, ...){ + build_stan_inputs <- function(name, response, submodels, log_lik, ...){ model_code <- name_to_modelcode(name) y_data <- get_stan_data(response) - pars <- get_pars(submodels) + pars <- get_pars(submodels, name, log_lik) submodels <- unname(submodels@submodels) types <- sapply(submodels, function(x) x@type) submodel_data <- lapply(submodels, get_stan_data) @@ -39,13 +39,18 @@ add_placeholder_priors <- function(submodel_data, types){ setGeneric("get_pars", function(object, ...) standardGeneric("get_pars")) -setMethod("get_pars", "ubmsSubmodelList", function(object, ...){ +setMethod("get_pars", "ubmsSubmodelList", function(object, model_name, log_lik, ...){ #Remove placeholder submodels - we don't want to save those parameters submodels <- object@submodels submodels <- submodels[!sapply(submodels, is_placeholder)] submodels <- unname(submodels) out <- unlist(lapply(submodels, get_pars)) - c(out, "log_lik") + + if(model_name == "distsamp") log_lik <- TRUE + if(any(sapply(submodels, has_spatial))) log_lik <- TRUE + + if(log_lik) out <- c(out, "log_lik") + out }) setMethod("get_pars", "ubmsSubmodel", function(object, ...){ diff --git a/R/kfold.R b/R/kfold.R new file mode 100644 index 0000000..3dc0500 --- /dev/null +++ b/R/kfold.R @@ -0,0 +1,68 @@ +#' K-fold Cross-validation of a ubmsFit Model +#' +#' Randomly partition data into K subsets of equal size (by site). Re-fit the model +#' K times, each time leaving out one of the subsets. Calculate the log-likelihood +#' for each of the sites that was left out. This function is an alternative +#' to \code{loo} (leave-one-out cross validation). +#' +#' @param x A \code{ubmsFit} model +#' @param K Number of folds into which the data will be partitioned +#' @param folds An optional vector with length equal to the number of sites in the data and containing integers from 1 to K, to manually assign sites to folds. You should use this if you plan to compare multiple models, since the folds for each model should be identical. You can use \code{loo::kfold_split_random} to generate this vector +#' @param quiet If \code{TRUE}, suppress progress bar +#' @param ... Currently ignored +#' +#' @return An object of class \code{elpd_generic} that is compatible with \code{loo::loo_compare} +#' +#' @include fit.R +#' @importFrom loo kfold +#' @export +setMethod("kfold", "ubmsFit", function(x, K=10, folds=NULL, quiet=FALSE, ...){ + if(is.null(folds)){ + folds <- loo::kfold_split_random(K, unmarked::numSites(x@data)) + } else { + stopifnot(length(folds) == unmarked::numSites(x@data)) + stopifnot(max(folds) == K) + } + if(has_spatial(x)){ + stop("kfold does not work with spatial models", call.=FALSE) + } + if(inherits(x, "ubmsFitDistsamp")){ + stop("kfold method not yet supported for stan_distsamp models", call.=FALSE) + } + + op <- pbapply::pboptions() + if(quiet) pbapply::pboptions(type = "none") + ll_unsort <- pbapply::pblapply(1:K, ll_fold, x, folds) + pbapply::pboptions(op) + ll_unsort <- do.call("cbind", ll_unsort) + + # Sort sites back to original order + sort_ind <- unlist(lapply(1:K, function(i) which(folds==i))) + ll <- matrix(NA, nrow=nrow(ll_unsort), ncol=ncol(ll_unsort)) + ll[,sort_ind] <- ll_unsort + ll <- ll[,apply(ll, 2, function(x) !any(is.na(x))),drop=FALSE] + + loo::elpd(ll) + +}) + +ll_fold <- function(i, object, folds){ + cl <- object@call + train_data <- object@data[which(!folds==i),] + cl$data <- train_data + cl$refresh <- 0 + refit <- suppressWarnings(eval(cl)) + + test_data <- object@data[which(folds==i),] + cl$data <- test_data + cl$return_inputs <- TRUE + inps <- eval(cl) + refit@submodels <- inps$submodels + ll <- get_loglik(refit, inps) + + # Handle missing sites + out <- matrix(NA, nrow=nrow(ll), ncol=sum(folds==i)) + not_missing <- which(folds==i) %in% which(!removed_sites(object)) + out[,not_missing] <- ll + out +} diff --git a/R/loglik.R b/R/loglik.R new file mode 100644 index 0000000..9d7644f --- /dev/null +++ b/R/loglik.R @@ -0,0 +1,147 @@ +# Calculate parameter from X, Z, offset +# Can't use posterior_linpred for this because it doesn't drop NAs +# Inps are output from get_stan_inputs() +# Right now the permute argument isn't used +calculate_par <- function(object, inps, submodel, permute=FALSE){ + X <- inps$stan_data[[paste0("X_",submodel)]] + beta <- extract_posterior(object, paste0("beta_", submodel)) + lp <- X %*% t(beta) + if(inps$stan_data[[paste0("has_random_",submodel)]]){ + Z <- Z_matrix(inps$submodels[submodel], na.rm=TRUE) + b <- extract_posterior(object, paste0("b_", submodel)) + lp <- lp + Z %*% t(b) + } + lp <- lp + inps$stan_data[[paste0("offset_",submodel)]] + do.call(inps$submodels[submodel]@link, list(lp)) +} + +extract_posterior <- function(object, par, permute=FALSE){ + if(permute) return(extract(object@stanfit, par)[[1]]) + ex <- extract(object@stanfit, par, permuted=FALSE) + smp <- lapply(1:ncol(ex), function(x){ + sapply(1:dim(ex)[3], function(i) ex[,x,i,drop=FALSE]) + }) + do.call("rbind", smp) +} + +setGeneric("get_loglik", function(object, ...) standardGeneric("get_loglik")) + +#' @include fit.R +setMethod("get_loglik", "ubmsFit", function(object, inps, ...){ + stop("Method not available for this fit type", call.=FALSE) +}) + +#' @include occu.R +setMethod("get_loglik", "ubmsFitOccu", function(object, inps, ...){ + psi <- calculate_par(object, inps, submodel="state") + p <- calculate_par(object, inps, submodel="det") + get_loglik_occu(inps$stan_data$y, inps$stan_data$M, inps$stan_data$si-1, + psi, p, inps$stan_data$Kmin) +}) + +#' @include pcount.R +setMethod("get_loglik", "ubmsFitPcount", function(object, inps, ...){ + lam <- calculate_par(object, inps, submodel="state") + p <- calculate_par(object, inps, submodel="det") + get_loglik_pcount(inps$stan_data$y, inps$stan_data$M, inps$stan_data$si-1, + lam, p, inps$stan_data$K, inps$stan_data$Kmin) +}) + +#' @include occuRN.R +setMethod("get_loglik", "ubmsFitOccuRN", function(object, inps, ...){ + lam <- calculate_par(object, inps, submodel="state") + p <- calculate_par(object, inps, submodel="det") + get_loglik_occuRN(inps$stan_data$y, inps$stan_data$M, inps$stan_data$si-1, + lam, p, inps$stan_data$K, inps$stan_data$Kmin) +}) + +#' @include multinomPois.R +setMethod("get_loglik", "ubmsFitMultinomPois", function(object, inps, ...){ + lam <- calculate_par(object, inps, submodel="state") + p <- calculate_par(object, inps, submodel="det") + get_loglik_multinomPois(inps$stan_data$y, inps$stan_data$M, inps$stan_data$si-1, + lam, p, inps$stan_data$y_dist) +}) + +#' @include colext.R +setMethod("get_loglik", "ubmsFitColext", function(object, inps, ...){ + psi <- calculate_par(object, inps, submodel="state") + psicube <- array(NA, c(dim(psi), 2)) + psicube[,,1] <- 1 - psi + psicube[,,2] <- psi + psicube <- aperm(psicube, c(1,3,2)) + + col <- calculate_par(object, inps, submodel="col") + ext <- calculate_par(object, inps, submodel="ext") + phicube <- array(NA, c(dim(col), 4)) + phicube[,,1] <- 1 - col + phicube[,,2] <- col + phicube[,,3] <- ext + phicube[,,4] <- 1 - ext + phicube <- aperm(phicube, c(1,3,2)) + pmat <- calculate_par(object, inps, submodel="det") + + get_loglik_colext(inps$stan_data$y, inps$stan_data$M, inps$stan_data$Tsamp-1, + inps$stan_data$J, inps$stan_data$si-1, psicube, + phicube, pmat, 1 - inps$stan_data$Kmin) +}) + +#' @include occuTTD.R +setMethod("get_loglik", "ubmsFitOccuTTD", function(object, inps, ...){ + psi <- calculate_par(object, inps, submodel="state") + lam <- calculate_par(object, inps, submodel="det") + k <- rep(0, ncol(psi)) + if(inps$stan_data$y_dist == 3){ + k <- extract(object, "beta_shape", permute=FALSE) + k <- exp(as.vector(k)) + } + get_loglik_occuTTD(inps$stan_data$aux2, inps$stan_data$M, inps$stan_data$si-1, + psi, lam, k, inps$stan_data$aux1, inps$stan_data$y_dist) +}) + +#' Extract Pointwise Log-likelihood From Model +#' +#' Extracts the pointwise log-likelihood matrix or array from a model. +#' This is useful as an input for functions in the \code{loo} package. +#' If called on a \code{ubmsFit} object, the log-likelihood matrix or array +#' is calculated using the posterior distributions of the parameters; the +#' log-likelihood values are not actually saved inside the model object. +#' If called on a \code{stanfit} object, \code{loo::extract_log_lik} is used. +#' In this case, the log-likelihood must be saved as one of the output +#' parameters in Stan. +#' +#' @param object A \code{ubmsFit} or \code{stanfit} object +#' @param parameter_name The name of the log-likelihood parameter in the +#' Stan output; ignored when \code{object} is a \code{ubmsFit} +#' @param merge_chains If \code{TRUE} (the default), all Markov chains are +#' merged together (i.e., stacked) and a matrix is returned. If ‘FALSE’ +#' they are kept separate and an array is returned. +#' +#' @return A matrix (samples x sites) or array (samples x chains x sites) +#' +#' @aliases extract_log_lik,ubmsFit-method extract_log_lik,ubmsFitDistsamp-method +#' @export +setGeneric("extract_log_lik", + function(object, parameter_name="log_lik", merge_chains=TRUE){ + loo::extract_log_lik(object, parameter_name, merge_chains) + }) + +setMethod("extract_log_lik", "ubmsFit", + function(object, parameter_name = "log_lik", merge_chains=TRUE){ + #if log_lik was saved as an output parameter, just return that + if(loglik_saved(object)){ + return(loo::extract_log_lik(object@stanfit, parameter_name, merge_chains)) + } + inps <- rebuild_inputs(object) + ll <- get_loglik(object, inps, merge_chains) + if(!merge_chains){ + nchain <- length(rstan::get_inits(object@stanfit)) + niter <- nrow(ll) / nchain + ll <- array(ll, c(niter, nchain, ncol(ll))) + } + ll +}) + +loglik_saved <- function(fit){ + "log_lik" %in% fit@stanfit@sim$pars_oi +} diff --git a/R/missing.R b/R/missing.R index b5615af..3f4c739 100644 --- a/R/missing.R +++ b/R/missing.R @@ -88,3 +88,8 @@ get_row_reps <- function(submodel, response){ if(yl %% n != 0) stop("Invalid model matrix dimensions", call.=FALSE) yl / n } + +#' @include fit.R +removed_sites <- function(object){ + object['state']@missing +} diff --git a/R/multinomPois.R b/R/multinomPois.R index 9b05ef6..9b99f39 100644 --- a/R/multinomPois.R +++ b/R/multinomPois.R @@ -15,6 +15,9 @@ #' @param prior_coef_det Prior distribution for the regression coefficients of #' the detection model #' @param prior_sigma Prior distribution on random effect standard deviations +#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values +#' in the output. This can greatly increase the size of the model. If +#' \code{FALSE}, the values are calculated post-hoc from the posteriors #' @param ... Arguments passed to the \code{\link{stan}} call, such as #' number of chains \code{chains} or iterations \code{iter} #' @@ -39,6 +42,7 @@ stan_multinomPois <- function(formula, prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ...){ forms <- split_formula(formula) @@ -62,7 +66,7 @@ stan_multinomPois <- function(formula, prior_intercept_det, prior_coef_det, prior_sigma) submodels <- ubmsSubmodelList(state, det) - ubmsFit("multinomPois", match.call(), data, response, submodels, ...) + ubmsFit("multinomPois", match.call(), data, response, submodels, log_lik, ...) } get_pifun_type <- function(umf){ @@ -15,6 +15,9 @@ #' @param prior_coef_det Prior distribution for the regression coefficients of #' the detection model #' @param prior_sigma Prior distribution on random effect standard deviations +#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values +#' in the output. This can greatly increase the size of the model. If +#' \code{FALSE}, the values are calculated post-hoc from the posteriors #' @param ... Arguments passed to the \code{\link{stan}} call, such as #' number of chains \code{chains} or iterations \code{iter} #' @@ -46,6 +49,7 @@ stan_occu <- function(formula, prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ...){ forms <- split_formula(formula) @@ -68,7 +72,7 @@ stan_occu <- function(formula, prior_intercept_det, prior_coef_det, prior_sigma) submodels <- ubmsSubmodelList(state, det) - ubmsFit("occu", match.call(), data, response, submodels, ...) + ubmsFit("occu", match.call(), data, response, submodels, log_lik, ...) } @@ -19,6 +19,9 @@ #' @param prior_coef_det Prior distribution for the regression coefficients of #' the detection model #' @param prior_sigma Prior distribution on random effect standard deviations +#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values +#' in the output. This can greatly increase the size of the model. If +#' \code{FALSE}, the values are calculated post-hoc from the posteriors #' @param ... Arguments passed to the \code{\link{stan}} call, such as #' number of chains \code{chains} or iterations \code{iter} #' @@ -47,6 +50,7 @@ stan_occuRN <- function(formula, prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ...){ forms <- split_formula(formula) @@ -69,7 +73,7 @@ stan_occuRN <- function(formula, prior_intercept_det, prior_coef_det, prior_sigma) submodels <- ubmsSubmodelList(state, det) - ubmsFit("occuRN", match.call(), data, response, submodels, ...) + ubmsFit("occuRN", match.call(), data, response, submodels, log_lik, ...) } diff --git a/R/occuTTD.R b/R/occuTTD.R index 605bec3..9eff184 100644 --- a/R/occuTTD.R +++ b/R/occuTTD.R @@ -30,6 +30,9 @@ #' @param prior_intercept_shape Prior distribution for the intercept of the #' shape parameter (i.e., log(shape)) for Weibull TTD models #' @param prior_sigma Prior distribution on random effect standard deviations +#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values +#' in the output. This can greatly increase the size of the model. If +#' \code{FALSE}, the values are calculated post-hoc from the posteriors #' @param ... Arguments passed to the \code{\link{stan}} call, such as #' number of chains \code{chains} or iterations \code{iter} #' @@ -88,6 +91,7 @@ stan_occuTTD <- function(psiformula=~1, prior_coef_det = normal(0, 2.5), prior_intercept_shape = normal(0,2.5), prior_sigma = gamma(1, 1), + log_lik = TRUE, ...){ if(data@numPrimary > 1) stop("Dynamic models not yet supported", call.=FALSE) @@ -119,7 +123,7 @@ stan_occuTTD <- function(psiformula=~1, submodels <- ubmsSubmodelList(state, det, shape) - ubmsFit("occuTTD", match.call(), data, response, submodels, ...) + ubmsFit("occuTTD", match.call(), data, response, submodels, log_lik, ...) } @@ -19,6 +19,9 @@ #' @param prior_coef_det Prior distribution for the regression coefficients of #' the detection model #' @param prior_sigma Prior distribution on random effect standard deviations +#' @param log_lik If \code{TRUE}, Stan will save pointwise log-likelihood values +#' in the output. This can greatly increase the size of the model. If +#' \code{FALSE}, the values are calculated post-hoc from the posteriors #' @param ... Arguments passed to the \code{\link{stan}} call, such as #' number of chains \code{chains} or iterations \code{iter} #' @@ -47,6 +50,7 @@ stan_pcount <- function(formula, prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ...){ forms <- split_formula(formula) @@ -69,7 +73,7 @@ stan_pcount <- function(formula, prior_intercept_det, prior_coef_det, prior_sigma) submodels <- ubmsSubmodelList(state, det) - ubmsFit("pcount", match.call(), data, response, submodels, ...) + ubmsFit("pcount", match.call(), data, response, submodels, log_lik, ...) } #Output object----------------------------------------------------------------- diff --git a/R/spatial.R b/R/spatial.R index a77000c..5c76f06 100644 --- a/R/spatial.R +++ b/R/spatial.R @@ -140,6 +140,10 @@ setMethod("has_spatial", "ubmsSubmodel", function(object, ...){ methods::.hasSlot(object, "spatial") }) +setMethod("has_spatial", "ubmsFit", function(object, ...){ + any(sapply(object@submodels@submodels, has_spatial)) +}) + setClass("ubmsSubmodelSpatial", contains = "ubmsSubmodel", slots=c(data_aug="data.frame", sites_aug="logical", spatial="formula")) diff --git a/R/submodel.R b/R/submodel.R index f7b6c8d..8090113 100644 --- a/R/submodel.R +++ b/R/submodel.R @@ -94,7 +94,7 @@ setMethod("model.matrix", "ubmsSubmodel", formula <- lme4::nobars(object@formula) out <- model.matrix(formula, mf) if(na.rm) out <- out[!object@missing,,drop=FALSE] - if(warn){ + if(warn & !all(is.na(out))){ max_cov <- max(out, na.rm=TRUE) if(max_cov > 4){ warning(paste0("Covariates possibly not standardized (max value = ", max_cov, @@ -139,7 +139,7 @@ get_xlev <- function(data, model_frame){ get_reTrms <- function(formula, data, newdata=NULL){ fb <- lme4::findbars(formula) mf <- model.frame(lme4::subbars(formula), data, na.action=stats::na.pass) - if(is.null(newdata)) return(lme4::mkReTrms(fb, mf)) + if(is.null(newdata)) return(lme4::mkReTrms(fb, mf, drop.unused.levels=FALSE)) new_mf <- model.frame(stats::terms(mf), newdata, na.action=stats::na.pass, xlev=get_xlev(data, mf)) lme4::mkReTrms(fb, new_mf, drop.unused.levels=FALSE) diff --git a/R/ubmsFit-methods.R b/R/ubmsFit-methods.R index f54ff3c..2f0bac8 100644 --- a/R/ubmsFit-methods.R +++ b/R/ubmsFit-methods.R @@ -160,7 +160,7 @@ setMethod("getY", "ubmsFit", function(object){ #' @importFrom loo loo #' @export setMethod("loo", "ubmsFit", function(x, ..., cores=getOption("mc.cores", 1)){ - loglik <- loo::extract_log_lik(x@stanfit, merge_chains=FALSE) + loglik <- extract_log_lik(x, merge_chains=FALSE) r_eff <- loo::relative_eff(exp(loglik), cores=cores) loo::loo(loglik, r_eff=r_eff, cores=cores) }) @@ -178,7 +178,7 @@ setMethod("loo", "ubmsFit", function(x, ..., cores=getOption("mc.cores", 1)){ #' @importFrom loo waic #' @export setMethod("waic", "ubmsFit", function(x, ...){ - loglik <- loo::extract_log_lik(x@stanfit) + loglik <- extract_log_lik(x) loo::waic(loglik) }) @@ -285,3 +285,11 @@ was_parallel <- function(object){ setMethod("get_stancode", "ubmsFit", function(object, ...){ rstan::get_stancode(object@stanfit, ...) }) + +# Re-create inputs +rebuild_inputs <- function(object){ + inps <- build_stan_inputs(object@stanfit@model_name, object@response, + object@submodels, loglik_saved(object)) + inps$submodels <- object@submodels + inps +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 748f0ca..1625f2d 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -20,6 +20,7 @@ reference: - title: Model selection contents: - fitList + - kfold,ubmsFit-method - loo,ubmsFit-method - modSel,ubmsFitList-method - waic,ubmsFit-method diff --git a/inst/stan/colext.stan b/inst/stan/colext.stan index e1aecf5..8d599da 100644 --- a/inst/stan/colext.stan +++ b/inst/stan/colext.stan @@ -42,15 +42,15 @@ real lp_colext(int[] y, int[] Tsamp, int[] J, row_vector psi, matrix phi_raw, if(T > 1){ for (t in 1:(T-1)){ phi = get_phi(phi_raw, Tsamp[t], Tsamp[t+1]); - end = idx + J[t] - 1; - Dpt = get_pY(y[idx:end], logit_p[idx:end], nd[t]); + end = idx + J[Tsamp[t]] - 1; + Dpt = get_pY(y[idx:end], logit_p[idx:end], nd[Tsamp[t]]); phi_prod *= diag_pre_multiply(Dpt, phi); - idx += J[t]; + idx += J[Tsamp[t]]; } } - end = idx + J[T] - 1; - Dpt = get_pY(y[idx:end], logit_p[idx:end], nd[T]); + end = idx + J[Tsamp[T]] - 1; + Dpt = get_pY(y[idx:end], logit_p[idx:end], nd[Tsamp[T]]); return log(dot_product(psi * phi_prod, Dpt)); } diff --git a/man/extract_log_lik.Rd b/man/extract_log_lik.Rd new file mode 100644 index 0000000..6a2e861 --- /dev/null +++ b/man/extract_log_lik.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/loglik.R +\name{extract_log_lik} +\alias{extract_log_lik} +\alias{extract_log_lik,ubmsFit-method} +\alias{extract_log_lik,ubmsFitDistsamp-method} +\title{Extract Pointwise Log-likelihood From Model} +\usage{ +extract_log_lik(object, parameter_name = "log_lik", merge_chains = TRUE) +} +\arguments{ +\item{object}{A \code{ubmsFit} or \code{stanfit} object} + +\item{parameter_name}{The name of the log-likelihood parameter in the +Stan output; ignored when \code{object} is a \code{ubmsFit}} + +\item{merge_chains}{If \code{TRUE} (the default), all Markov chains are +merged together (i.e., stacked) and a matrix is returned. If ‘FALSE’ +they are kept separate and an array is returned.} +} +\value{ +A matrix (samples x sites) or array (samples x chains x sites) +} +\description{ +Extracts the pointwise log-likelihood matrix or array from a model. +This is useful as an input for functions in the \code{loo} package. +If called on a \code{ubmsFit} object, the log-likelihood matrix or array +is calculated using the posterior distributions of the parameters; the +log-likelihood values are not actually saved inside the model object. +If called on a \code{stanfit} object, \code{loo::extract_log_lik} is used. +In this case, the log-likelihood must be saved as one of the output +parameters in Stan. +} diff --git a/man/kfold-ubmsFit-method.Rd b/man/kfold-ubmsFit-method.Rd new file mode 100644 index 0000000..c5cd1e7 --- /dev/null +++ b/man/kfold-ubmsFit-method.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/kfold.R +\name{kfold,ubmsFit-method} +\alias{kfold,ubmsFit-method} +\title{K-fold Cross-validation of a ubmsFit Model} +\usage{ +\S4method{kfold}{ubmsFit}(x, K = 10, folds = NULL, quiet = FALSE, ...) +} +\arguments{ +\item{x}{A \code{ubmsFit} model} + +\item{K}{Number of folds into which the data will be partitioned} + +\item{folds}{An optional vector with length equal to the number of sites in the data and containing integers from 1 to K, to manually assign sites to folds. You should use this if you plan to compare multiple models, since the folds for each model should be identical. You can use \code{loo::kfold_split_random} to generate this vector} + +\item{quiet}{If \code{TRUE}, suppress progress bar} + +\item{...}{Currently ignored} +} +\value{ +An object of class \code{elpd_generic} that is compatible with \code{loo::loo_compare} +} +\description{ +Randomly partition data into K subsets of equal size (by site). Re-fit the model +K times, each time leaving out one of the subsets. Calculate the log-likelihood +for each of the sites that was left out. This function is an alternative +to \code{loo} (leave-one-out cross validation). +} diff --git a/man/stan_colext.Rd b/man/stan_colext.Rd index 8ce5688..8482a6c 100644 --- a/man/stan_colext.Rd +++ b/man/stan_colext.Rd @@ -19,6 +19,7 @@ stan_colext( prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ... ) } @@ -60,6 +61,10 @@ the detection model} \item{prior_sigma}{Prior distribution on random effect standard deviations} +\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values +in the output. This can greatly increase the size of the model. If +\code{FALSE}, the values are calculated post-hoc from the posteriors} + \item{...}{Arguments passed to the \code{\link{stan}} call, such as number of chains \code{chains} or iterations \code{iter}} } diff --git a/man/stan_multinomPois.Rd b/man/stan_multinomPois.Rd index ff34467..19d7281 100644 --- a/man/stan_multinomPois.Rd +++ b/man/stan_multinomPois.Rd @@ -12,6 +12,7 @@ stan_multinomPois( prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ... ) } @@ -35,6 +36,10 @@ the detection model} \item{prior_sigma}{Prior distribution on random effect standard deviations} +\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values +in the output. This can greatly increase the size of the model. If +\code{FALSE}, the values are calculated post-hoc from the posteriors} + \item{...}{Arguments passed to the \code{\link{stan}} call, such as number of chains \code{chains} or iterations \code{iter}} } diff --git a/man/stan_occu.Rd b/man/stan_occu.Rd index 512bc64..3092691 100644 --- a/man/stan_occu.Rd +++ b/man/stan_occu.Rd @@ -12,6 +12,7 @@ stan_occu( prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ... ) } @@ -35,6 +36,10 @@ the detection model} \item{prior_sigma}{Prior distribution on random effect standard deviations} +\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values +in the output. This can greatly increase the size of the model. If +\code{FALSE}, the values are calculated post-hoc from the posteriors} + \item{...}{Arguments passed to the \code{\link{stan}} call, such as number of chains \code{chains} or iterations \code{iter}} } diff --git a/man/stan_occuRN.Rd b/man/stan_occuRN.Rd index fcf0f02..8e820e6 100644 --- a/man/stan_occuRN.Rd +++ b/man/stan_occuRN.Rd @@ -13,6 +13,7 @@ stan_occuRN( prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ... ) } @@ -40,6 +41,10 @@ the detection model} \item{prior_sigma}{Prior distribution on random effect standard deviations} +\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values +in the output. This can greatly increase the size of the model. If +\code{FALSE}, the values are calculated post-hoc from the posteriors} + \item{...}{Arguments passed to the \code{\link{stan}} call, such as number of chains \code{chains} or iterations \code{iter}} } diff --git a/man/stan_occuTTD.Rd b/man/stan_occuTTD.Rd index 3583848..4589a90 100644 --- a/man/stan_occuTTD.Rd +++ b/man/stan_occuTTD.Rd @@ -18,6 +18,7 @@ stan_occuTTD( prior_coef_det = normal(0, 2.5), prior_intercept_shape = normal(0, 2.5), prior_sigma = gamma(1, 1), + log_lik = TRUE, ... ) } @@ -61,6 +62,10 @@ shape parameter (i.e., log(shape)) for Weibull TTD models} \item{prior_sigma}{Prior distribution on random effect standard deviations} +\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values +in the output. This can greatly increase the size of the model. If +\code{FALSE}, the values are calculated post-hoc from the posteriors} + \item{...}{Arguments passed to the \code{\link{stan}} call, such as number of chains \code{chains} or iterations \code{iter}} } diff --git a/man/stan_pcount.Rd b/man/stan_pcount.Rd index 2f5a635..4bd52db 100644 --- a/man/stan_pcount.Rd +++ b/man/stan_pcount.Rd @@ -14,6 +14,7 @@ stan_pcount( prior_intercept_det = logistic(0, 1), prior_coef_det = logistic(0, 1), prior_sigma = gamma(1, 1), + log_lik = TRUE, ... ) } @@ -43,6 +44,10 @@ the detection model} \item{prior_sigma}{Prior distribution on random effect standard deviations} +\item{log_lik}{If \code{TRUE}, Stan will save pointwise log-likelihood values +in the output. This can greatly increase the size of the model. If +\code{FALSE}, the values are calculated post-hoc from the posteriors} + \item{...}{Arguments passed to the \code{\link{stan}} call, such as number of chains \code{chains} or iterations \code{iter}} } diff --git a/src/Makevars b/src/Makevars index a9f7369..e7a95c7 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,10 +1,19 @@ STANHEADERS_SRC = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "message()" -e "cat(system.file('include', 'src', package = 'StanHeaders', mustWork = TRUE))" -e "message()" | grep "StanHeaders") STANC_FLAGS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "cat(ifelse(utils::packageVersion('rstan') >= 2.26, '-DUSE_STANC3',''))") + +# COMMENT THIS OUT WHEN TESTING CLANG PKG_CPPFLAGS = -I"../inst/include" -I"$(STANHEADERS_SRC)" -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -DBOOST_MATH_OVERFLOW_ERROR_POLICY=errno_on_error $(STANC_FLAGS) + +# UNCOMMENT THIS WHEN TESTING CLANG +#PKG_CPPFLAGS = -I"../inst/include" -isystem"$(STANHEADERS_SRC)" -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -DBOOST_MATH_OVERFLOW_ERROR_POLICY=errno_on_error $(STANC_FLAGS) + PKG_CXXFLAGS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "RcppParallel::CxxFlags()") $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "StanHeaders:::CxxFlags()") + PKG_LIBS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "RcppParallel::RcppParallelLibs()") $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "StanHeaders:::LdFlags()") +PKG_LIBS+= $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +# COMMENT THIS OUT WHEN TESTING CLANG COMPILER CXX_STD = CXX14 diff --git a/src/Makevars.win b/src/Makevars.win index eb698ce..a16dabe 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -3,7 +3,9 @@ STANHEADERS_SRC = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "message()" -e STANC_FLAGS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "cat(ifelse(utils::packageVersion('rstan') >= 2.26, '-DUSE_STANC3',''))") PKG_CPPFLAGS = -I"../inst/include" -I"$(STANHEADERS_SRC)" -DBOOST_DISABLE_ASSERTS -DEIGEN_NO_DEBUG -DRCPP_PARALLEL_USE_TBB=1 $(STANC_FLAGS) PKG_CXXFLAGS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "RcppParallel::CxxFlags()") $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "StanHeaders:::CxxFlags()") + PKG_LIBS = $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "RcppParallel::RcppParallelLibs()") $(shell "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "StanHeaders:::LdFlags()") +PKG_LIBS+= $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) CXX_STD = CXX14 diff --git a/src/exp_counts_occu.cpp b/src/exp_counts_occu.cpp index a79df67..fc37099 100644 --- a/src/exp_counts_occu.cpp +++ b/src/exp_counts_occu.cpp @@ -28,7 +28,6 @@ arma::vec exp_counts_occu(arma::mat obs, arma::ivec no_detects, int n_eh = obs.n_cols; vec counts_expect = zeros(n_eh); - int idx = 0; int pstart, pend; for (int i=0; i < n_eh; i++){ pstart = 0; @@ -73,7 +72,6 @@ arma::vec exp_counts_occuRN(arma::mat obs, arma::ivec Kmin, int n_eh = obs.n_cols; vec counts_expect = zeros(n_eh); - int idx = 0; int pstart, pend; for (int i=0; i < n_eh; i++){ pstart = 0; diff --git a/src/loglik.cpp b/src/loglik.cpp new file mode 100644 index 0000000..4c92f93 --- /dev/null +++ b/src/loglik.cpp @@ -0,0 +1,345 @@ +// [[Rcpp::depends(RcppArmadillo)]] +#include <RcppArmadillo.h> + +using namespace arma; + +// [[Rcpp::export]] +arma::mat get_loglik_occu(arma::vec y, int M, arma::imat si, arma::mat psimat, + arma::mat pmat, arma::ivec Kmin){ + + int S = psimat.n_cols; + mat out(S,M); + vec p_s(pmat.n_rows); + int J; + double cp; + vec psub, ysub; + + for (int s = 0; s < S; s++){ + + p_s = pmat.col(s); + + for (int m = 0; m < M; m++){ + + ysub = y.subvec(si(m,0), si(m,1)); + psub = p_s.subvec(si(m,0), si(m,1)); + J = psub.size(); + cp = 1.0; + + for (int j = 0; j < J; j++){ + cp *= pow(psub(j), ysub(j)) * pow(1-psub(j), 1-ysub(j)); + } + + if(Kmin(m) == 1){ + out(s,m) = log(cp * psimat(m,s) + DBL_MIN); + } else if(Kmin(m) == 0){ + out(s,m) = log(cp * psimat(m,s) + (1-psimat(m,s)) + DBL_MIN); + } + } + } + + return out; +} + + +// [[Rcpp::export]] +arma::mat get_loglik_pcount(arma::vec y, int M, arma::imat si, + arma::mat lammat, arma::mat pmat, int K, + arma::ivec Kmin){ + + int S = lammat.n_cols; + mat out(S,M); + vec p_s(pmat.n_rows); + vec psub, ysub; + int J; + double fac, ff, N, ky, numN; + + for (int s = 0; s < S; s++){ + + p_s = pmat.col(s); + + for (int m = 0; m < M; m++){ + ysub = y.subvec(si(m,0), si(m,1)); + J = ysub.size(); + psub = p_s.subvec(si(m,0), si(m,1)); + fac = 1.0; + ff = lammat(m,s) * prod(1 - psub); + numN = K - Kmin(m); + + for (int i = 1; i < (numN + 1); i++){ + N = K - i + 1; + ky = 1.0; + for (int j = 0; j < J; j++){ + ky *= N / (N - ysub(j)); + } + fac = 1 + fac * ff * ky / N; + } + + out(s,m) = log(fac) + Rf_dpois(Kmin(m), lammat(m,s), true); + for (int j = 0; j < J; j ++){ + out(s,m) += Rf_dbinom(ysub(j), Kmin(m), psub(j), true); + } + } + } + + return out; +} + +// [[Rcpp::export]] +arma::mat get_loglik_occuRN(arma::vec y, int M, arma::imat si, + arma::mat lammat, arma::mat rmat, int K, + arma::ivec Kmin){ + + int S = lammat.n_cols; + mat out(S,M); + vec r_s(rmat.n_rows); + vec ysub, qsub; + int J; + double f, g, p, lik; + + for (int s = 0; s < S; s++){ + + r_s = rmat.col(s); + + for (int m = 0; m < M; m++){ + lik = 0.0; + ysub = y.subvec(si(m,0), si(m,1)); + qsub = 1 - r_s.subvec(si(m,0), si(m,1)); + J = ysub.size(); + for (int k=Kmin(m); k < (K+1); k++){ + f = Rf_dpois(k, lammat(m, s), false); + g = 0.0; + for (int j = 0; j < J; j++){ + p = 1 - pow(qsub(j), k); + g += Rf_dbinom(ysub(j), 1, p, true); + } + lik += f * exp(g); + } + out(s,m) = log(lik + DBL_MIN); + } + } + + return out; + +} + +// pi functions +arma::vec pi_removal(arma::vec p){ + int J = p.size(); + vec pi_out(J); + pi_out(0) = p(0); + for (int j = 1; j < J; j++){ + pi_out(j) = pi_out(j-1) / p(j-1) * (1-p(j-1)) * p(j); + } + return pi_out; +} + +arma::vec pi_double(arma::vec p){ + vec pi_out(3); + pi_out(0) = p(0) * (1 - p(1)); + pi_out(1) = p(1) * (1 - p(0)); + pi_out(2) = p(0) * p(1); + return pi_out; +} + +arma::vec pi_fun(int pi_type, arma::vec p, int J){ + vec out(J); + if(pi_type == 0){ + out = pi_double(p); + } else if(pi_type == 1){ + out = pi_removal(p); + } else{ + Rcpp::stop("Invalid pi function"); + } + return out; +} + +// [[Rcpp::export]] +arma::mat get_loglik_multinomPois(arma::vec y, int M, arma::imat si, + arma::mat lammat, arma::mat pmat, int pi_type){ + + int S = lammat.n_cols; + mat out = zeros(S,M); + vec p_s(pmat.n_rows); + vec ysub, psub; + int R = pmat.n_rows / M; + int pstart, pend, J; + vec cp; + + for (int s = 0; s < S; s++){ + + p_s = pmat.col(s); + pstart = 0; + + for (int m = 0; m < M; m++){ + pend = pstart + R - 1; + ysub = y.subvec(si(m,0), si(m,1)); + psub = p_s.subvec(pstart, pend); + J = ysub.size(); + cp = pi_fun(pi_type, psub, J); + + for (int j = 0; j < J; j++){ + out(s, m) += Rf_dpois(ysub(j), lammat(m, s) * cp(j), true); + } + pstart += R; + } + } + return out; +} + + +arma::vec get_pY(arma::vec y, arma::vec p, int nd){ + + vec out(2); + int J = y.size(); + out(0) = nd; + out(1) = 1.0; + for (int j = 0; j < J; j++){ + out(1) *= Rf_dbinom(y(j), 1, p(j), false); + } + return out; +} + +arma::mat phi_matrix(arma::rowvec phi_raw){ + mat out(2,2); + out(0,0) = phi_raw(0); + out(0,1) = phi_raw(1); + out(1,0) = phi_raw(2); + out(1,1) = phi_raw(3); + return out; +} + +arma::mat get_phi(arma::mat phi_raw, int Tstart, int Tnext){ + mat phi = eye(2,2); + int delta = Tnext - Tstart; + if(delta == 1){ + return phi_matrix(phi_raw.row(Tstart)); + } + + for (int d = 1; d < (delta + 1); d++){ + phi = phi * phi_matrix(phi_raw.row(Tstart + d - 1)); + } + return phi; +} + + + +// [[Rcpp::export]] +arma::mat get_loglik_colext(arma::vec y, int M, arma::ivec Tsamp, arma::imat J, + arma::imat si, arma::cube psicube, arma::cube phicube, + arma::mat pmat, arma::imat nd){ + + int S = pmat.n_cols; + mat out(S,M); + vec p_s(pmat.n_rows); + mat psi_s(psicube.n_rows, psicube.n_cols); + mat phi_s(phicube.n_rows, phicube.n_cols); + + rowvec psi_m; + mat phi_m; + ivec Tsamp_m; + irowvec J_m, nd_m; + vec y_m, p_m; + int T; + mat phi_prod(2,2); + mat phi(2,2); + vec Dpt(2); + int idx, end; + + for (int s = 0; s < S; s++){ + + p_s = pmat.col(s); + psi_s = psicube.slice(s); + phi_s = phicube.slice(s); + + for (int m = 0; m < M; m++){ + idx = 0; + y_m = y.subvec(si(m,0), si(m,1)); + Tsamp_m = Tsamp.subvec(si(m,2), si(m,3)); + J_m = J.row(m); + + psi_m = psi_s.row(m); + phi_m = phi_s.rows(si(m,4),si(m,5)); + p_m = p_s.subvec(si(m,0), si(m,1)); + nd_m = nd.row(m); + + T = Tsamp_m.size(); + phi_prod = eye(2,2); + + if(T > 1){ + for (int t = 0; t < (T-1); t++){ + phi = get_phi(phi_m, Tsamp_m(t), Tsamp_m(t+1)); + end = idx + J_m(Tsamp_m(t)) - 1; + Dpt = get_pY(y_m.subvec(idx, end), p_m.subvec(idx, end), + nd_m(Tsamp_m(t))); + + phi_prod *= diagmat(Dpt) * phi; + idx += J_m(Tsamp_m(t)); + } + } + + end = idx + J_m(Tsamp_m(T-1)) - 1; + Dpt = get_pY(y_m.subvec(idx,end), p_m.subvec(idx,end), nd_m(Tsamp_m(T-1))); + out(s, m) = log(dot(psi_m * phi_prod, Dpt)); + } + } + + return out; +} + + +arma::vec ttd_prob_exp(arma::vec y, arma::vec lam, arma::ivec delta){ + int J = y.size(); + vec e_lamt(J); + + for (int j=0; j < J; j++){ + e_lamt(j) = pow(lam(j), delta(j)) * exp(-lam(j)*y(j)); + } + return e_lamt; +} + +arma::vec ttd_prob_weib(arma::vec y, arma::vec lam, arma::ivec delta, double k){ + int J = y.size(); + vec e_lamt(J); + for (int j = 0; j < J; j++){ + e_lamt(j) = pow(k*lam(j)*pow(lam(j)*y(j), (k-1)), delta(j)) * + exp(-1*pow(lam(j)*y(j), k)); + } + return e_lamt; +} + + +// [[Rcpp::export]] +arma::mat get_loglik_occuTTD(arma::vec y, int M, arma::imat si, arma::mat psimat, + arma::mat lammat, arma::vec k, + arma::ivec delta, int ydist){ + + int S = psimat.n_cols; + mat out(S,M); + vec lam_s(lammat.n_rows); + vec e_lamt; + vec y_m, lam_m; + ivec delta_m; + double lik; + + for (int s = 0; s < S; s++){ + + lam_s = lammat.col(s); + + for (int m = 0; m < M; m++){ + y_m = y.subvec(si(m,0), si(m,1)); + lam_m = lam_s.subvec(si(m,0), si(m,1)); + delta_m = delta.subvec(si(m,0), si(m,1)); + + if(ydist == 1){ + e_lamt = ttd_prob_exp(y_m, lam_m, delta_m); + } else if(ydist == 3){ + e_lamt = ttd_prob_weib(y_m, lam_m, delta_m, k(s)); + } + + lik = psimat(m,s) * prod(e_lamt) + (1-psimat(m,s)) * (1-max(delta_m)); + out(s,m) = log(lik); + } + } + + return out; +} diff --git a/tests/testthat/test_colext.R b/tests/testthat/test_colext.R index 618714c..4882a6e 100644 --- a/tests/testthat/test_colext.R +++ b/tests/testthat/test_colext.R @@ -82,6 +82,13 @@ test_that("stan_colext handles NA values",{ expect_RMSE(coef(fit_na), coef(fit), 1.5) }) +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(100/2 * 2, numSites(fit@data))) + expect_between(sum(ll), -5700, -5500) +}) + test_that("ubmsFitColext gof method works",{ set.seed(123) g <- gof(fit, draws=5, quiet=TRUE) diff --git a/tests/testthat/test_distsamp.R b/tests/testthat/test_distsamp.R index 8a8c22a..c26afa0 100644 --- a/tests/testthat/test_distsamp.R +++ b/tests/testthat/test_distsamp.R @@ -121,6 +121,13 @@ test_that("stan_distsamp handles NA values",{ expect_error(stan_distsamp(~1~1, ltUMF_na)) }) +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit_pt_hn) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(200/2 * 2, numSites(fit_pt_hn@data))) + expect_between(sum(ll), -39000, -38000) +}) + test_that("ubmsFitDistsamp gof method works",{ set.seed(123) g <- lapply(line_mods, function(x) gof(x, draws=5, quiet=TRUE)) @@ -311,3 +318,8 @@ test_that("distsamp spatial works", { ps <- plot_spatial(fit_spat) expect_is(ps, "gg") }) + +test_that("kfold errors when used on a stan_distsamp model", { + expect_error(kfold(fit_line_hn), + "kfold method not yet supported for stan_distsamp models") +}) diff --git a/tests/testthat/test_fit.R b/tests/testthat/test_fit.R index 636dd01..bfd2497 100644 --- a/tests/testthat/test_fit.R +++ b/tests/testthat/test_fit.R @@ -23,7 +23,7 @@ resp <- ubmsResponse(umf@y,"binomial","binomial",max_primary=1) #Build a stanfit object set.seed(123) -inp <- build_stan_inputs("occu", resp, sl) +inp <- build_stan_inputs("occu", resp, sl, log_lik=FALSE) good_fit <- TRUE tryCatch({ @@ -36,8 +36,9 @@ test <- process_stanfit(sf, sl) skip_if(!good_fit, "Test setup failed") test_that("ubmsFit object is constructed correctly",{ - - ufit <- suppressWarnings(ubmsFit("occu", call("dummy"), umf, resp, sl, + ufit <- suppressWarnings(ubmsFit("occu", + as.call(str2lang("stan_occu(formula = ~1 ~ 1, data = umf, chains = 2, iter = 40)")), + umf, resp, sl, chains=2, iter=40, refresh=0)) expect_true(inherits(ufit, "ubmsFit")) expect_true(inherits(ufit@stanfit, "stanfit")) @@ -58,14 +59,14 @@ test_that("remove_placeholders removes placeholder submodels from list",{ expect_equal(list_remove, sl) }) -test_that("get_loo generates loo object from stanfit",{ - loo_obj <- suppressWarnings(get_loo(sf)) - expect_true(inherits(loo_obj, "psis_loo")) -}) +#test_that("get_loo generates loo object from stanfit",{ +# loo_obj <- suppressWarnings(get_loo(sf)) +# expect_true(inherits(loo_obj, "psis_loo")) +#}) test_that("fit_model builds model correctly",{ ufit <- suppressWarnings( - fit_model("occu", resp, sl, chains=2, iter=20, refresh=0)) + fit_model("occu", resp, sl, log_lik=FALSE, chains=2, iter=20, refresh=0)) expect_true(inherits(ufit, "stanfit")) nms <- stanfit_names(sl) expect_equal(ufit@sim$fnames_oi[1:length(nms)], nms) @@ -80,7 +81,7 @@ test_that("fit_model builds model correctly",{ test_that("specific model name is shown in console output",{ # e.g. 'occu' instead of 'single_season' out <- capture.output(ufit <- suppressWarnings( - fit_model("occu", resp, sl, chains=2, iter=20))) + fit_model("occu", resp, sl, log_lik=FALSE, chains=2, iter=20))) expect_true(any(grepl("occu", out))) expect_false(any(grepl("single_season", out))) }) diff --git a/tests/testthat/test_inputs.R b/tests/testthat/test_inputs.R index 2a221a7..9bf5e9e 100644 --- a/tests/testthat/test_inputs.R +++ b/tests/testthat/test_inputs.R @@ -10,7 +10,7 @@ test_that("stan inputs are built correctly", { y <- matrix(c(1,0,0,1,1,1,0,0,1), nrow=3, byrow=T) resp <- ubmsResponse(y, "binomial", "P") - inp <- build_stan_inputs("occu", resp, sl) + inp <- build_stan_inputs("occu", resp, sl, log_lik=FALSE) expect_is(inp, "list") expect_equal(names(inp), c("stan_data", "pars")) @@ -33,6 +33,10 @@ test_that("parameter list for stan is generated correctly",{ sl <- ubmsSubmodelList(state, det) expect_equal(get_pars(det), "beta_det") expect_equal(get_pars(state), c("beta_state", "b_state", "sigma_state")) + expect_equal(get_pars(sl, "occu", log_lik=FALSE), + c("beta_state", "b_state", "sigma_state", "beta_det")) + expect_equal(get_pars(sl, "occu", log_lik=TRUE), + c("beta_state", "b_state", "sigma_state", "beta_det", "log_lik")) }) test_that("get_stan_data pulls necessary info from response object",{ diff --git a/tests/testthat/test_multinomPois.R b/tests/testthat/test_multinomPois.R index 6959786..c9bb896 100644 --- a/tests/testthat/test_multinomPois.R +++ b/tests/testthat/test_multinomPois.R @@ -96,6 +96,13 @@ test_that("stan_multinomPois handles NA values",{ expect_is(coef(fit_rem_na), "numeric") }) +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit_double) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(100/2 * 2, numSites(fit_double@data))) + expect_between(sum(ll), -5000, -4600) +}) + test_that("ubmsFitMultinomPois gof method works",{ ##here set.seed(123) g <- gof(fit_double, draws=5, quiet=TRUE) diff --git a/tests/testthat/test_occu.R b/tests/testthat/test_occu.R index bc4737a..558c2f7 100644 --- a/tests/testthat/test_occu.R +++ b/tests/testthat/test_occu.R @@ -68,6 +68,41 @@ test_that("stan_occu handles NA values",{ expect_RMSE(coef(fit), coef(fit_na), 2) }) +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(100/2 * 2, numSites(fit@data))) + expect_between(sum(ll), -3500, -3200) +}) + +test_that("extract_log_lik works when there are missing values and random effects",{ + skip_on_cran() + skip_on_ci() + umf3 <- umf2 + umf3@siteCovs$group <- sample(letters[1:5], nrow(umf2@siteCovs), replace=TRUE) + fit_na <- suppressWarnings(stan_occu(~x2+(1|group)~1, umf3[1:10,], chains=2, + iter=50, refresh=0)) + expect_is(fit_na, "ubmsFitOccu") + ll <- extract_log_lik(fit_na) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(50,9)) +}) + +test_that("log_lik argument controls saving log_lik parameter", { + skip_on_cran() + set.seed(123) + fit <- suppressWarnings(stan_occu(~x2~x1, umf[1:10,], chains=2, + iter=100, refresh=0)) + set.seed(123) + fit2 <- suppressWarnings(stan_occu(~x2~x1, umf[1:10,], chains=2, + iter=100, refresh=0, log_lik=FALSE)) + + expect_equal(fit@loo$estimates, fit2@loo$estimates) + expect_true("log_lik" %in% fit@stanfit@sim$pars_oi) + expect_false("log_lik" %in% fit2@stanfit@sim$pars_oi) + +}) + test_that("ubmsFitOccu gof method works",{ set.seed(123) g <- gof(fit, draws=5, quiet=TRUE) @@ -181,3 +216,15 @@ test_that("Fitted/residual methods work with ubmsFitOccu",{ expect_is(rp3, "gtable") expect_is(mp, "gg") }) + +test_that("stan_occu kfold method works",{ + +set.seed(123) +kf <- kfold(fit, K=2, quiet=TRUE) +expect_is(kf, "elpd_generic") +expect_equal(kf$estimates[1,1], -37.3802, tol=1e-4) + +expect_error(kfold(fit, K=2, folds=rep(3,10))) +expect_error(kfold(fit, K=2, folds=rep(1,5))) + +}) diff --git a/tests/testthat/test_occuRN.R b/tests/testthat/test_occuRN.R index 1c6ae39..5d01b90 100644 --- a/tests/testthat/test_occuRN.R +++ b/tests/testthat/test_occuRN.R @@ -68,6 +68,13 @@ test_that("stan_occuRN handles NA values",{ expect_is(coef(fit_na), "numeric") }) +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(100/2 * 2, numSites(fit@data))) + expect_between(sum(ll), -3000, -2700) +}) + test_that("ubmsFitOccuRN gof method works",{ set.seed(123) g <- gof(fit, draws=5, quiet=TRUE) diff --git a/tests/testthat/test_occuTTD.R b/tests/testthat/test_occuTTD.R index a36a77b..1eff257 100644 --- a/tests/testthat/test_occuTTD.R +++ b/tests/testthat/test_occuTTD.R @@ -85,6 +85,13 @@ test_that("stan_occuTTD handles NA values",{ expect_is(coef(fit_na), "numeric") }) +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(100/2 * 2, numSites(fit@data))) + expect_between(sum(ll), -700, -400) +}) + test_that("ubmsFitOccuTTD gof method gives error",{ expect_error(gof(fit, draws=5, quiet=TRUE)) }) diff --git a/tests/testthat/test_pcount.R b/tests/testthat/test_pcount.R index 32b0710..f0dc013 100644 --- a/tests/testthat/test_pcount.R +++ b/tests/testthat/test_pcount.R @@ -84,6 +84,13 @@ test_that("stan_pcount handles NA values",{ expect_true(is.numeric(coef(fit_na))) }) +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(100/2 * 2, numSites(fit@data))) + expect_between(sum(ll), -7000, -6000) +}) + test_that("ubmsFitPcount gof method works",{ set.seed(123) g <- gof(fit, draws=5, quiet=TRUE) diff --git a/tests/testthat/test_spatial.R b/tests/testthat/test_spatial.R index ed1307a..e4a4140 100644 --- a/tests/testthat/test_spatial.R +++ b/tests/testthat/test_spatial.R @@ -26,10 +26,11 @@ for (i in sample(1:500, 300, replace=FALSE)){ } umf <- unmarkedFrameOccu(y=y, siteCovs=dat_occ, obsCovs=dat_p) - +umf20 <- umf[1:20,] +umf10 <- umf[1:10,] fit <- suppressMessages(suppressWarnings(stan_occu(~1~cov1+RSR(x,y,1), - data=umf[1:20,], chains=2, iter=200, refresh=0))) -fit2 <- suppressWarnings(stan_occu(~1~1, data=umf[1:10,], chains=2, iter=200, refresh=0)) + data=umf20, chains=2, iter=200, refresh=0))) +fit2 <- suppressWarnings(stan_occu(~1~1, data=umf10, chains=2, iter=200, refresh=0)) test_that("spatial model output structure", { expect_is(fit, "ubmsFitOccu") @@ -125,6 +126,11 @@ test_that("has_spatial works on lists of formulas", { expect_error(has_spatial(list(det=~1,state=~RSR(x,y,1)),support=FALSE)) }) +test_that("has_spatial works on ubmsFit objects",{ + expect_true(has_spatial(fit)) + expect_false(has_spatial(fit2)) +}) + test_that("construction of ubmsSubmodelSpatial objects", { ex <- extract_missing_sites(umf) sm <- ubmsSubmodelSpatial("Test","test", ex$umf@siteCovs, ~1+RSR(x,y,1), "plogis", @@ -196,3 +202,14 @@ test_that("plot_spatial returns ggplot", { expect_error(plot_spatial(umf)) expect_error(plot_spatial(fit2)) }) + +test_that("extract_log_lik method works",{ + ll <- extract_log_lik(fit) + expect_is(ll, "matrix") + expect_equal(dim(ll), c(200/2 * 2, numSites(fit@data)-7)) + expect_between(sum(ll), -7000, -6500) +}) + +test_that("kfold errors when used on spatial model",{ + expect_error(kfold(fit)) +}) diff --git a/tests/testthat/test_ubmsFit_methods.R b/tests/testthat/test_ubmsFit_methods.R index 8acc733..3237ce2 100644 --- a/tests/testthat/test_ubmsFit_methods.R +++ b/tests/testthat/test_ubmsFit_methods.R @@ -158,3 +158,28 @@ test_that("get_stancode method works",{ out <- get_stancode(fit) expect_is(out, "character") }) + +test_that("rebuild_inputs can get Stan inputs from fitted model",{ + inps <- rebuild_inputs(fit) + cl <- fit@call + cl$return_inputs <- TRUE + refit <- eval(cl) + expect_equal(names(inps), c("stan_data", "pars", "submodels")) + expect_equal(inps, refit) +}) + +# Functions in loglik.R +test_that("calculate_par works",{ + inps <- rebuild_inputs(fit) + cp <- calculate_par(fit, inps, 'state') + expect_equal(dim(cp), c(3,40)) + cp_mean <- apply(cp, 1, mean) + pr <- predict(fit, 'state') + expect_equivalent(cp_mean, pr$Predicted) +}) + +test_that("extract_posterior works",{ + post <- extract_posterior(fit, 'beta_state') + expect_is(post, "matrix") + expect_equal(dim(post), c(40, 2)) +}) diff --git a/vignettes/spatial-models.Rmd b/vignettes/spatial-models.Rmd index 8e6a780..63cc2f7 100644 --- a/vignettes/spatial-models.Rmd +++ b/vignettes/spatial-models.Rmd @@ -250,21 +250,22 @@ fit_ubms ## ## Occupancy (logit-scale): ## Estimate SD 2.5% 97.5% n_eff Rhat -## (Intercept) -1.038 0.2286 -1.5060 -0.601 1269 1.00 -## habCov12 -0.212 0.3231 -0.8503 0.427 2404 1.00 -## habCov13 -0.258 0.3273 -0.9035 0.374 2058 1.00 -## habCov2 1.080 0.1708 0.7683 1.439 1429 1.00 -## RSR [tau] 0.155 0.0589 0.0704 0.293 253 1.01 +## (Intercept) -1.054 0.2324 -1.5121 -0.611 1177 1.01 +## habCov12 -0.216 0.3201 -0.8831 0.383 2499 1.00 +## habCov13 -0.261 0.3223 -0.8950 0.356 2262 1.00 +## habCov2 1.089 0.1801 0.7416 1.461 1295 1.01 +## RSR [tau] 0.144 0.0579 0.0675 0.283 264 1.02 ## ## Detection (logit-scale): ## Estimate SD 2.5% 97.5% n_eff Rhat -## (Intercept) 1.727 0.370 1.004 2.481 1610 1 -## detCov1 0.113 0.210 -0.310 0.533 4962 1 -## detCov21 1.185 0.555 0.122 2.320 2189 1 -## detCov22 1.372 0.557 0.356 2.496 2455 1 -## detCov23 -0.242 0.484 -1.197 0.719 2013 1 +## (Intercept) 1.734 0.371 1.035 2.488 1756 1 +## detCov1 0.108 0.203 -0.291 0.497 4281 1 +## detCov21 1.183 0.556 0.102 2.270 2319 1 +## detCov22 1.379 0.561 0.307 2.533 2471 1 +## detCov23 -0.245 0.486 -1.204 0.704 2218 1 ## -## LOOIC: 668.02 +## LOOIC: 668.475 +## Runtime: 88.788 sec ``` The precision term (tau, $\tau$) associated with the spatial random effect is the final row of the occupancy model estimates. @@ -291,12 +292,12 @@ head(psi) ``` ## Predicted SD 2.5% 97.5% -## 1 0.1441475 0.04413909 0.07266996 0.2439409 -## 2 0.2360008 0.08352206 0.10501601 0.4271224 -## 3 0.3578073 0.11260628 0.17170766 0.6029410 -## 4 0.4527547 0.12485601 0.23020833 0.7192525 -## 5 0.4522272 0.11820159 0.23762047 0.6955722 -## 6 0.3847294 0.11279317 0.18753530 0.6292943 +## 1 0.1441804 0.04624230 0.07020541 0.2517342 +## 2 0.2390063 0.08739192 0.10297331 0.4388321 +## 3 0.3640516 0.11781411 0.16512449 0.6239752 +## 4 0.4593940 0.13070577 0.22866630 0.7265760 +## 5 0.4570641 0.12368794 0.23700396 0.7078097 +## 6 0.3865731 0.11876489 0.18351350 0.6434938 ``` ## Plotting @@ -349,8 +350,8 @@ round(modSel(fl),2) ``` ## elpd nparam elpd_diff se_diff weight -## fit_ubms -334.01 44.80 0.00 0.0 1 -## fit_nonspatial -380.54 8.43 -46.53 8.1 0 +## fit_ubms -334.24 46.25 0.00 0.00 1 +## fit_nonspatial -380.46 8.36 -46.22 8.26 0 ``` A measure of predictive accuracy (`elpd`) is shown for each model, and the model with the highest predictive accuracy is listed first. @@ -376,7 +377,7 @@ names <- list( fit_stocc <- spatial.occupancy( detection.model = ~ detCov1 + detCov2, occupancy.model = ~ habCov1 + habCov2, - spatial.model = list(model="rsr", threshold=1, moran.cut=0.1*nrow(habData)), + spatial.model = list(model="rsr", threshold=1.01, moran.cut=0.1*nrow(habData)), so.data = make.so.data(visitData, habData, names), prior = list(a.tau=0.5, b.tau=0.00005, Q.b=0.1, Q.g=0.1), control = list(burnin=1000/5, iter=4000/5, thin=5) diff --git a/vignettes/spatial-models.Rmd.orig b/vignettes/spatial-models.Rmd.orig index d467b30..1382c7a 100644 --- a/vignettes/spatial-models.Rmd.orig +++ b/vignettes/spatial-models.Rmd.orig @@ -300,7 +300,7 @@ names <- list( fit_stocc <- spatial.occupancy( detection.model = ~ detCov1 + detCov2, occupancy.model = ~ habCov1 + habCov2, - spatial.model = list(model="rsr", threshold=1, moran.cut=0.1*nrow(habData)), + spatial.model = list(model="rsr", threshold=1.01, moran.cut=0.1*nrow(habData)), so.data = make.so.data(visitData, habData, names), prior = list(a.tau=0.5, b.tau=0.00005, Q.b=0.1, Q.g=0.1), control = list(burnin=1000/5, iter=4000/5, thin=5) @@ -320,7 +320,7 @@ if(file.exists("fit_stocc.Rds")){ nul <- capture.output(fit_stocc <- spatial.occupancy( detection.model = ~ detCov1 + detCov2, occupancy.model = ~ habCov1 + habCov2, - spatial.model = list(model="rsr", threshold=1, moran.cut=0.1*nrow(habData)), + spatial.model = list(model="rsr", threshold=1.01, moran.cut=0.1*nrow(habData)), so.data = make.so.data(visitData, habData, names), prior = list(a.tau=0.5, b.tau=0.00005, Q.b=0.1, Q.g=0.1), control = list(burnin=1000/5, iter=4000/5, thin=5) |