From f3f651b2d0aac0608b83f858c523d7718da55cfc Mon Sep 17 00:00:00 2001 From: Jeff Date: Wed, 28 Aug 2024 15:51:55 -0400 Subject: [PATCH] updates --- .gitignore | 2 + NEWS.md | 5 +- R/PGOcc.R | 16 +- R/intMsPGOcc.R | 462 +++++++++++++++++----------------- R/intPGOcc.R | 32 +-- R/lfJSDM.R | 8 +- R/lfMsPGOcc.R | 16 +- R/msPGOcc.R | 16 +- R/postHocLM.R | 4 +- R/sfJSDM.R | 8 +- R/sfMsPGOcc.R | 16 +- R/spIntPGOcc.R | 28 +-- R/spMsPGOcc.R | 28 +-- R/spPGOcc.R | 28 +-- R/stIntPGOcc.R | 12 +- R/stMsPGOcc.R | 12 +- R/stPGOcc.R | 16 +- R/svcMsPGOcc.R | 12 +- R/svcPGBinom.R | 8 +- R/svcPGOcc.R | 16 +- R/svcTIntPGOcc.R | 12 +- R/svcTMsPGOcc.R | 12 +- R/svcTPGBinom.R | 8 +- R/svcTPGOcc.R | 16 +- R/tIntPGOcc.R | 12 +- R/tPGOcc.R | 16 +- man/PGOcc.Rd | 2 + man/intMsPGOcc.Rd | 21 +- man/intPGOcc.Rd | 2 + man/lfJSDM.Rd | 2 + man/lfMsPGOcc.Rd | 2 + man/msPGOcc.Rd | 2 + man/predict.PGOcc.Rd | 3 +- man/predict.intMsPGOcc.Rd | 60 ++--- man/predict.intPGOcc.Rd | 2 + man/predict.lfJSDM.Rd | 3 +- man/predict.lfMsPGOcc.Rd | 3 +- man/predict.msPGOcc.Rd | 3 +- man/predict.sfJSDM.Rd | 3 +- man/predict.sfMsPGOcc.Rd | 3 +- man/predict.spIntPGOcc.Rd | 3 +- man/predict.spMsPGOcc.Rd | 37 +-- man/predict.spPGOcc.Rd | 3 +- man/predict.stMsPGOcc.Rd | 17 +- man/predict.stPGOcc.Rd | 2 + man/predict.svcMsPGOcc.Rd | 12 +- man/predict.svcPGBinom.Rd | 2 + man/predict.svcPGOcc.Rd | 2 + man/predict.svcTMsPGOcc.Rd | 18 +- man/predict.svcTPGBinom.Rd | 2 + man/predict.svcTPGOcc.Rd | 2 + man/predict.tMsPGOcc.Rd | 2 + man/predict.tPGOcc.Rd | 2 + man/sfJSDM.Rd | 2 + man/sfMsPGOcc.Rd | 2 + man/spIntPGOcc.Rd | 2 + man/spMsPGOcc.Rd | 2 + man/spPGOcc.Rd | 2 + man/stMsPGOcc.Rd | 2 + man/stPGOcc.Rd | 2 + man/svcMsPGOcc.Rd | 2 + man/svcPGBinom.Rd | 2 + man/svcPGOcc.Rd | 2 + man/svcTMsPGOcc.Rd | 2 + man/svcTPGBinom.Rd | 2 + man/svcTPGOcc.Rd | 2 + man/tMsPGOcc.Rd | 2 + man/tPGOcc.Rd | 2 + src/intMsPGOcc.cpp | 2 - src/intPGOcc.cpp | 2 - vignettes/parallelization.Rmd | 216 ++++++++++++++++ 71 files changed, 775 insertions(+), 511 deletions(-) create mode 100644 vignettes/parallelization.Rmd diff --git a/.gitignore b/.gitignore index ce4f7a8..db7cf3e 100644 --- a/.gitignore +++ b/.gitignore @@ -32,6 +32,8 @@ vignettes/svcModels_cache vignettes/svcModels_files vignettes/identifiability_cache vignettes/identifiability_files +vignettes/parallelization_cache +vignettes/parallelization_files *.swp vignettes/modelConsiderations_cache vignettes/modelConsiderations_files diff --git a/NEWS.md b/NEWS.md index af9db82..f8cc791 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,9 +4,9 @@ + New functionality for fitting multi-season, single-species integrated occupancy models. The function `tIntPGOcc()` fits a non-spatial multi-season integrated occupancy model, `stIntPGOcc()` fits a spatial multi-season integrated occupancy model, and `svcTIntPGOcc()` fits a spatially-varying coefficient multi-season occupancy model. Random intercepts are supported in both the occurrence and detection formulas for both model types. + Added in functionality for both occupancy and detection random intercepts in single-species single-season integrated models (`intPGOcc()` and `spIntPGOcc()`) using `lme4` syntax (e.g., `(1 | observer)` for a random effect of observer). + `simTIntPGOcc()` is a new function that allows simulation of single-species multi-season detection-nondetection data from multiple data sources. -+ Update `simMsIntPGOcc()` to now include simulation of data sets with spatially-varying coefficients and unstructured random effects on both occurrence and detection. ++ Updated `simMsIntPGOcc()` to now include simulation of data sets with spatially-varying coefficients and unstructured random effects on both occurrence and detection. + Fixed a bug in the k-fold cross-validation for spatial integrated occupancy models (NNGP models only) that could lead to incorrect model deviance results under certain situations depending on how the spatial coordinates were ordered on the user-side relative to how they are re-ordered when fitting the model. If using `spIntPGOcc()` with `NNGP = TRUE` and using cross-validation, I suggest re-running the analysis. Apologies for the inconvenience. -+ Added in a `residuals()` function to extract occupancy and detection residuals following the approach of [Wilson et al. (2019)](https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1002/ecy.2703) for single-season, single-species occupancy models (functions `PGOcc()`, `spPGOcc()`, and `svcPGOcc()`). ++ Added in a `residuals()` function to extract occupancy and detection residuals following the approach of [Wilson et al. (2019)](https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1002/ecy.2703) for single-season, single-species occupancy models (functions `PGOcc()`, `spPGOcc()`, and `svcPGOcc()`). I'm hoping to implement this for all model functions and improve GoF functionality. If anyone has any interest in helping out with this, then please let me know! + `waicOcc()` for integrated single-species models is now substantially faster. + `updateMCMC()` now works with lfJSDM. + Fixed a bug in `updateMCMC()` that prevented it from working with `spAbundance::msAbund()` when there were random effects in the model. Also added the `save.fitted` argument to `updateMCMC()` to allow it to work with `msAbund()` and not save the replicate/fitted data values in cases where the amount of RAM is an important consideration. @@ -37,6 +37,7 @@ + Wrote a new "vignette" (really more of a blog post) on some recommendations to help improve interpretability of inferences in SVC models. + Fixed a few typos in the MCMC sampler vignettes for factor models and SVC models. + Fixed a bug that prevented cross-validation from working properly in multi-species models when setting `k.fold.only = TRUE`. Thanks to Zack Steel for pointing this out. ++ Fixed a typo in the generation of initial values for latent unstructured random effects in all model functions. The typo had no major ramifications, if anything it would have just led to slower convergence, as it resulted in very large (or very small) initial values for the latent random effects that are not really viable on the logit scale. # spOccupancy 0.7.2 diff --git a/R/PGOcc.R b/R/PGOcc.R index 27a5604..dde9a2f 100644 --- a/R/PGOcc.R +++ b/R/PGOcc.R @@ -551,7 +551,7 @@ PGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -581,7 +581,7 @@ PGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -677,14 +677,14 @@ PGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits.list[[i]] <- 1 beta.star.inits.list[[i]] <- 1 } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits.list[[i]] <- 1 alpha.star.inits.list[[i]] <- 1 @@ -721,11 +721,11 @@ PGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } storage.mode(chain.info) <- "integer" @@ -880,7 +880,7 @@ PGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) p.re.level.names.fit <- list() for (t in 1:p.det.re) { @@ -900,7 +900,7 @@ PGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/R/intMsPGOcc.R b/R/intMsPGOcc.R index f8a6483..9bf0ad7 100644 --- a/R/intMsPGOcc.R +++ b/R/intMsPGOcc.R @@ -1,8 +1,7 @@ intMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, n.samples, n.omp.threads = 1, verbose = TRUE, n.report = 100, n.burn = round(.10 * n.samples), n.thin = 1, n.chains = 1, - parallel.chains = FALSE, k.fold, k.fold.threads = 1, - k.fold.seed = 100, k.fold.only = FALSE, ...){ + parallel.chains = FALSE, ...){ ptm <- proc.time() @@ -102,10 +101,6 @@ intMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } } - # TODO: add in cross-validation eventually - if (!missing(k.fold)) { - message("k-fold cross-validation is not currently supported for integrated multi-species occupancy models. Please keep an eye out for future updates with this functionality.") - } # Check if all detection covariates are at site level for a given # data set. @@ -910,245 +905,244 @@ intMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, # Fit the model ------------------------------------------------------- out.tmp <- list() out <- list() - if (!k.fold.only) { - if (parallel.chains) { - if (verbose) { - cat("\n----------------------------------------\n"); - cat("\tRunning the model\n"); - cat("----------------------------------------\n"); - message("MCMC chains are running in parallel. Model progress output is suppressed.") - } - beta.comm.inits.list <- list() - alpha.comm.inits.list <- list() - tau.sq.beta.inits.list <- list() - tau.sq.alpha.inits.list <- list() - beta.inits.list <- list() - alpha.inits.list <- list() - sigma.sq.psi.inits.list <- list() - beta.star.inits.list <- list() - # sigma.sq.p.inits.list <- list() - # alpha.star.inits.list <- list() - for (i in 1:n.chains) { - beta.comm.inits.list[[i]] <- beta.comm.inits - alpha.comm.inits.list[[i]] <- alpha.comm.inits - tau.sq.beta.inits.list[[i]] <- tau.sq.beta.inits - tau.sq.alpha.inits.list[[i]] <- tau.sq.alpha.inits - beta.inits.list[[i]] <- beta.inits - alpha.inits.list[[i]] <- alpha.inits - sigma.sq.psi.inits.list[[i]] <- sigma.sq.psi.inits - beta.star.inits.list[[i]] <- beta.star.inits - # sigma.sq.p.inits.list[[i]] <- sigma.sq.p.inits - # alpha.star.inits.list[[i]] <- alpha.star.inits - } - for (i in 2:n.chains) { - if ((!fix.inits)) { - beta.comm.inits.list[[i]] <- rnorm(p.occ, mu.beta.comm, sqrt(sigma.beta.comm)) - alpha.comm.inits.list[[i]] <- rnorm(p.det, mu.alpha.comm, sqrt(sigma.alpha.comm)) - tau.sq.beta.inits.list[[i]] <- runif(p.occ, 0.5, 10) - tau.sq.alpha.inits.list[[i]] <- runif(p.det, 0.5, 10) - beta.inits.list[[i]] <- matrix(rnorm(N * p.occ, beta.comm.inits, - sqrt(tau.sq.beta.inits.list[[i]])), N, p.occ) - alpha.inits.list[[i]] <- matrix(rnorm(N * p.det, alpha.comm.inits, - sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) - if (p.occ.re > 0) { - sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, - sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) - beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) - } - # if (p.det.re > 0) { - # sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - # alpha.star.inits.list[[i]] <- rnorm(n.det.re, - # sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) - # alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) - # } + if (parallel.chains) { + if (verbose) { + cat("\n----------------------------------------\n"); + cat("\tRunning the model\n"); + cat("----------------------------------------\n"); + message("MCMC chains are running in parallel. Model progress output is suppressed.") + } + beta.comm.inits.list <- list() + alpha.comm.inits.list <- list() + tau.sq.beta.inits.list <- list() + tau.sq.alpha.inits.list <- list() + beta.inits.list <- list() + alpha.inits.list <- list() + sigma.sq.psi.inits.list <- list() + beta.star.inits.list <- list() + # sigma.sq.p.inits.list <- list() + # alpha.star.inits.list <- list() + for (i in 1:n.chains) { + beta.comm.inits.list[[i]] <- beta.comm.inits + alpha.comm.inits.list[[i]] <- alpha.comm.inits + tau.sq.beta.inits.list[[i]] <- tau.sq.beta.inits + tau.sq.alpha.inits.list[[i]] <- tau.sq.alpha.inits + beta.inits.list[[i]] <- beta.inits + alpha.inits.list[[i]] <- alpha.inits + sigma.sq.psi.inits.list[[i]] <- sigma.sq.psi.inits + beta.star.inits.list[[i]] <- beta.star.inits + # sigma.sq.p.inits.list[[i]] <- sigma.sq.p.inits + # alpha.star.inits.list[[i]] <- alpha.star.inits + } + for (i in 2:n.chains) { + if ((!fix.inits)) { + beta.comm.inits.list[[i]] <- rnorm(p.occ, mu.beta.comm, sqrt(sigma.beta.comm)) + alpha.comm.inits.list[[i]] <- rnorm(p.det, mu.alpha.comm, sqrt(sigma.alpha.comm)) + tau.sq.beta.inits.list[[i]] <- runif(p.occ, 0.5, 10) + tau.sq.alpha.inits.list[[i]] <- runif(p.det, 0.5, 10) + beta.inits.list[[i]] <- matrix(rnorm(N * p.occ, beta.comm.inits, + sqrt(tau.sq.beta.inits.list[[i]])), N, p.occ) + alpha.inits.list[[i]] <- matrix(rnorm(N * p.det, alpha.comm.inits, + sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) + if (p.occ.re > 0) { + sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) + beta.star.inits.list[[i]] <- rnorm(n.occ.re, + sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) + beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } + # if (p.det.re > 0) { + # sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) + # alpha.star.inits.list[[i]] <- rnorm(n.det.re, + # sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) + # alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) + # } } - par.cl <- parallel::makePSOCKcluster(n.chains) - registerDoParallel(par.cl) - out.tmp <- foreach(i = 1:n.chains) %dorng% { - .Call("intMsPGOcc", y, X, X.p.all, X.re, consts, - n.occ.re.long, p.det.long, n.obs.long, N.long, - beta.inits.list[[i]], alpha.inits.list[[i]], z.inits, beta.comm.inits.list[[i]], - alpha.comm.inits.list[[i]], tau.sq.beta.inits.list[[i]], - tau.sq.alpha.inits.list[[i]], sigma.sq.psi.inits.list[[i]], - beta.star.inits.list[[i]], z.long.indx.c, - alpha.comm.indx.c, data.indx.c, alpha.sp.indx.c, - alpha.dat.indx.c, alpha.p.det.indx.c, - sp.long.indx.c, obs.long.indx.c, sp.site.indx.c, - sp.dat.long.indx.c, beta.star.indx, beta.level.indx, - mu.beta.comm, mu.alpha.comm, - Sigma.beta.comm, sigma.alpha.comm, - tau.sq.beta.a, tau.sq.beta.b, tau.sq.alpha.a, - tau.sq.alpha.b, sigma.sq.psi.a, sigma.sq.psi.b, - n.samples, n.omp.threads, - verbose, n.report, samples.info, chain.info) - } - parallel::stopCluster(par.cl) - } else { - for (i in 1:n.chains) { - # Change initial values if i > 1 - if ((i > 1) & (!fix.inits)) { - beta.comm.inits <- rnorm(p.occ, mu.beta.comm, sqrt(sigma.beta.comm)) - alpha.comm.inits <- rnorm(p.det, mu.alpha.comm, sqrt(sigma.alpha.comm)) - tau.sq.beta.inits <- runif(p.occ, 0.5, 10) - tau.sq.alpha.inits <- runif(p.det, 0.5, 10) - beta.inits <- matrix(rnorm(N * p.occ, beta.comm.inits, - sqrt(tau.sq.beta.inits)), N, p.occ) - alpha.inits <- rnorm(N * p.det, alpha.comm.inits, sqrt(tau.sq.alpha.inits)) - alpha.inits <- alpha.inits[alpha.indx == 1] - if (p.occ.re > 0) { - sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) - beta.star.inits <- rep(beta.star.inits, N) - } + } + par.cl <- parallel::makePSOCKcluster(n.chains) + registerDoParallel(par.cl) + out.tmp <- foreach(i = 1:n.chains) %dorng% { + .Call("intMsPGOcc", y, X, X.p.all, X.re, consts, + n.occ.re.long, p.det.long, n.obs.long, N.long, + beta.inits.list[[i]], alpha.inits.list[[i]], z.inits, beta.comm.inits.list[[i]], + alpha.comm.inits.list[[i]], tau.sq.beta.inits.list[[i]], + tau.sq.alpha.inits.list[[i]], sigma.sq.psi.inits.list[[i]], + beta.star.inits.list[[i]], z.long.indx.c, + alpha.comm.indx.c, data.indx.c, alpha.sp.indx.c, + alpha.dat.indx.c, alpha.p.det.indx.c, + sp.long.indx.c, obs.long.indx.c, sp.site.indx.c, + sp.dat.long.indx.c, beta.star.indx, beta.level.indx, + mu.beta.comm, mu.alpha.comm, + Sigma.beta.comm, sigma.alpha.comm, + tau.sq.beta.a, tau.sq.beta.b, tau.sq.alpha.a, + tau.sq.alpha.b, sigma.sq.psi.a, sigma.sq.psi.b, + n.samples, n.omp.threads, + verbose, n.report, samples.info, chain.info) + } + parallel::stopCluster(par.cl) + } else { + for (i in 1:n.chains) { + # Change initial values if i > 1 + if ((i > 1) & (!fix.inits)) { + beta.comm.inits <- rnorm(p.occ, mu.beta.comm, sqrt(sigma.beta.comm)) + alpha.comm.inits <- rnorm(p.det, mu.alpha.comm, sqrt(sigma.alpha.comm)) + tau.sq.beta.inits <- runif(p.occ, 0.5, 10) + tau.sq.alpha.inits <- runif(p.det, 0.5, 10) + beta.inits <- matrix(rnorm(N * p.occ, beta.comm.inits, + sqrt(tau.sq.beta.inits)), N, p.occ) + alpha.inits <- rnorm(N * p.det, alpha.comm.inits, sqrt(tau.sq.alpha.inits)) + alpha.inits <- alpha.inits[alpha.indx == 1] + if (p.occ.re > 0) { + sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) + beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rep(beta.star.inits, N) } - - storage.mode(chain.info) <- "integer" - out.tmp[[i]] <- .Call("intMsPGOcc", y, X, X.p.all, X.re, consts, - n.occ.re.long, p.det.long, n.obs.long, N.long, - beta.inits, alpha.inits, z.inits, beta.comm.inits, - alpha.comm.inits, tau.sq.beta.inits, tau.sq.alpha.inits, - sigma.sq.psi.inits, - beta.star.inits, z.long.indx.c, - alpha.comm.indx.c, data.indx.c, alpha.sp.indx.c, - alpha.dat.indx.c, alpha.p.det.indx.c, - sp.long.indx.c, obs.long.indx.c, sp.site.indx.c, - sp.dat.long.indx.c, beta.star.indx, beta.level.indx, - mu.beta.comm, mu.alpha.comm, - Sigma.beta.comm, sigma.alpha.comm, - tau.sq.beta.a, tau.sq.beta.b, tau.sq.alpha.a, - tau.sq.alpha.b, sigma.sq.psi.a, sigma.sq.psi.b, - n.samples, n.omp.threads, - verbose, n.report, samples.info, chain.info) - chain.info[1] <- chain.info[1] + 1 } + + storage.mode(chain.info) <- "integer" + out.tmp[[i]] <- .Call("intMsPGOcc", y, X, X.p.all, X.re, consts, + n.occ.re.long, p.det.long, n.obs.long, N.long, + beta.inits, alpha.inits, z.inits, beta.comm.inits, + alpha.comm.inits, tau.sq.beta.inits, tau.sq.alpha.inits, + sigma.sq.psi.inits, + beta.star.inits, z.long.indx.c, + alpha.comm.indx.c, data.indx.c, alpha.sp.indx.c, + alpha.dat.indx.c, alpha.p.det.indx.c, + sp.long.indx.c, obs.long.indx.c, sp.site.indx.c, + sp.dat.long.indx.c, beta.star.indx, beta.level.indx, + mu.beta.comm, mu.alpha.comm, + Sigma.beta.comm, sigma.alpha.comm, + tau.sq.beta.a, tau.sq.beta.b, tau.sq.alpha.a, + tau.sq.alpha.b, sigma.sq.psi.a, sigma.sq.psi.b, + n.samples, n.omp.threads, + verbose, n.report, samples.info, chain.info) + chain.info[1] <- chain.info[1] + 1 } - # Calculate R-Hat --------------- - out$rhat <- list() - if (n.chains > 1) { - # as.vector removes the "Upper CI" when there is only 1 variable. - out$rhat$beta.comm <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) - mcmc(t(a$beta.comm.samples)))), - autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) - out$rhat$alpha.comm <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) - mcmc(t(a$alpha.comm.samples)))), - autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) - out$rhat$tau.sq.beta <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) - mcmc(t(a$tau.sq.beta.samples)))), - autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) - out$rhat$tau.sq.alpha <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) - mcmc(t(a$tau.sq.alpha.samples)))), - autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) - out$rhat$beta <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) - mcmc(t(a$beta.samples)))), + } + # Calculate R-Hat --------------- + out$rhat <- list() + if (n.chains > 1) { + # as.vector removes the "Upper CI" when there is only 1 variable. + out$rhat$beta.comm <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) + mcmc(t(a$beta.comm.samples)))), + autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) + out$rhat$alpha.comm <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) + mcmc(t(a$alpha.comm.samples)))), + autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) + out$rhat$tau.sq.beta <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) + mcmc(t(a$tau.sq.beta.samples)))), + autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) + out$rhat$tau.sq.alpha <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) + mcmc(t(a$tau.sq.alpha.samples)))), + autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) + out$rhat$beta <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) + mcmc(t(a$beta.samples)))), + autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) + out$rhat$alpha <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) + mcmc(t(a$alpha.samples)))), + autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) + if (p.occ.re > 0) { + out$rhat$sigma.sq.psi <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) + mcmc(t(a$sigma.sq.psi.samples)))), autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) - out$rhat$alpha <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) - mcmc(t(a$alpha.samples)))), - autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) - if (p.occ.re > 0) { - out$rhat$sigma.sq.psi <- as.vector(gelman.diag(mcmc.list(lapply(out.tmp, function(a) - mcmc(t(a$sigma.sq.psi.samples)))), - autoburnin = FALSE, multivariate = FALSE)$psrf[, 2]) - } - } else { - out$rhat$beta.comm <- rep(NA, p.occ) - out$rhat$alpha.comm <- rep(NA, p.det) - out$rhat$tau.sq.beta <- rep(NA, p.occ) - out$rhat$tau.sq.alpha <- rep(NA, p.det) - out$rhat$beta <- rep(NA, p.occ * N) - out$rhat$alpha <- rep(NA, sum(p.det.long * N.long)) - if (p.occ.re > 0) { - out$rhat$sigma.sq.psi <- rep(NA, p.occ.re) - } } - # Put everything into MCMC objects - out$beta.comm.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$beta.comm.samples)))) - colnames(out$beta.comm.samples) <- x.names - out$alpha.comm.samples <- mcmc(do.call(rbind, - lapply(out.tmp, function(a) t(a$alpha.comm.samples)))) - colnames(out$alpha.comm.samples) <- x.p.names - out$tau.sq.beta.samples <- mcmc(do.call(rbind, - lapply(out.tmp, function(a) t(a$tau.sq.beta.samples)))) - colnames(out$tau.sq.beta.samples) <- x.names - out$tau.sq.alpha.samples <- mcmc(do.call(rbind, - lapply(out.tmp, function(a) t(a$tau.sq.alpha.samples)))) - colnames(out$tau.sq.alpha.samples) <- x.p.names - - if (is.null(sp.names)) { - sp.names <- paste('sp', 1:N, sep = '') - } - coef.names <- paste(rep(x.names, each = N), sp.names, sep = '-') - out$beta.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$beta.samples)))) - colnames(out$beta.samples) <- coef.names - out$alpha.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$alpha.samples)))) - coef.names.det <- paste(rep(x.p.names, each = N), sp.names, sep = '-') - coef.names.det <- coef.names.det[alpha.indx == 1] - colnames(out$alpha.samples) <- coef.names.det - out$z.samples <- do.call(abind, lapply(out.tmp, function(a) array(a$z.samples, - dim = c(N, J, n.post.samples)))) - out$z.samples <- aperm(out$z.samples, c(3, 1, 2)) - out$psi.samples <- do.call(abind, lapply(out.tmp, function(a) array(a$psi.samples, - dim = c(N, J, n.post.samples)))) - out$psi.samples <- aperm(out$psi.samples, c(3, 1, 2)) - out$like.samples <- do.call(abind, lapply(out.tmp, function(a) array(a$like.samples, - dim = c(N, J, n.post.samples)))) - out$like.samples <- aperm(out$like.samples, c(3, 1, 2)) - if (p.occ.re > 0) { - out$sigma.sq.psi.samples <- mcmc( - do.call(rbind, lapply(out.tmp, function(a) t(a$sigma.sq.psi.samples)))) - colnames(out$sigma.sq.psi.samples) <- x.re.names - out$beta.star.samples <- mcmc( - do.call(rbind, lapply(out.tmp, function(a) t(a$beta.star.samples)))) - tmp.names <- unlist(re.level.names) - beta.star.names <- paste(rep(x.re.names, n.occ.re.long), tmp.names, sep = '-') - beta.star.names <- paste(beta.star.names, rep(sp.names, each = n.occ.re), sep = '-') - colnames(out$beta.star.samples) <- beta.star.names - out$re.level.names <- re.level.names - } - # Calculate effective sample sizes - out$ESS <- list() - out$ESS$beta.comm <- effectiveSize(out$beta.comm.samples) - out$ESS$alpha.comm <- effectiveSize(out$alpha.comm.samples) - out$ESS$tau.sq.beta <- effectiveSize(out$tau.sq.beta.samples) - out$ESS$tau.sq.alpha <- effectiveSize(out$tau.sq.alpha.samples) - out$ESS$beta <- effectiveSize(out$beta.samples) - out$ESS$alpha <- effectiveSize(out$alpha.samples) - if (p.occ.re > 0) { - out$ESS$sigma.sq.psi <- effectiveSize(out$sigma.sq.psi.samples) - } - out$X <- X - out$X.p <- X.p - out$X.p.re <- X.p.re - out$X.re <- X.re - out$y <- y.big - out$call <- cl - out$n.samples <- n.samples - out$x.names <- x.names - out$sp.names <- sp.names - out$x.p.names <- x.p.names - out$n.post <- n.post.samples - out$n.thin <- n.thin - out$n.burn <- n.burn - out$n.chains <- n.chains - out$alpha.indx <- alpha.dat.indx.c + 1 - out$alpha.comm.indx <- alpha.comm.indx.r - out$N <- N - out$sites <- sites - out$species <- data$species - # TODO: add in detection random effects eventually - # if (p.det.re > 0) { - # out$pRE <- TRUE - # } else { - out$pRE <- FALSE - # } + } else { + out$rhat$beta.comm <- rep(NA, p.occ) + out$rhat$alpha.comm <- rep(NA, p.det) + out$rhat$tau.sq.beta <- rep(NA, p.occ) + out$rhat$tau.sq.alpha <- rep(NA, p.det) + out$rhat$beta <- rep(NA, p.occ * N) + out$rhat$alpha <- rep(NA, sum(p.det.long * N.long)) if (p.occ.re > 0) { - out$psiRE <- TRUE - } else { - out$psiRE <- FALSE + out$rhat$sigma.sq.psi <- rep(NA, p.occ.re) } } + # Put everything into MCMC objects + out$beta.comm.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$beta.comm.samples)))) + colnames(out$beta.comm.samples) <- x.names + out$alpha.comm.samples <- mcmc(do.call(rbind, + lapply(out.tmp, function(a) t(a$alpha.comm.samples)))) + colnames(out$alpha.comm.samples) <- x.p.names + out$tau.sq.beta.samples <- mcmc(do.call(rbind, + lapply(out.tmp, function(a) t(a$tau.sq.beta.samples)))) + colnames(out$tau.sq.beta.samples) <- x.names + out$tau.sq.alpha.samples <- mcmc(do.call(rbind, + lapply(out.tmp, function(a) t(a$tau.sq.alpha.samples)))) + colnames(out$tau.sq.alpha.samples) <- x.p.names + + if (is.null(sp.names)) { + sp.names <- paste('sp', 1:N, sep = '') + } + coef.names <- paste(rep(x.names, each = N), sp.names, sep = '-') + out$beta.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$beta.samples)))) + colnames(out$beta.samples) <- coef.names + out$alpha.samples <- mcmc(do.call(rbind, lapply(out.tmp, function(a) t(a$alpha.samples)))) + coef.names.det <- paste(rep(x.p.names, each = N), sp.names, sep = '-') + coef.names.det <- coef.names.det[alpha.indx == 1] + colnames(out$alpha.samples) <- coef.names.det + out$z.samples <- do.call(abind, lapply(out.tmp, function(a) array(a$z.samples, + dim = c(N, J, n.post.samples)))) + out$z.samples <- aperm(out$z.samples, c(3, 1, 2)) + out$psi.samples <- do.call(abind, lapply(out.tmp, function(a) array(a$psi.samples, + dim = c(N, J, n.post.samples)))) + out$psi.samples <- aperm(out$psi.samples, c(3, 1, 2)) + out$like.samples <- do.call(abind, lapply(out.tmp, function(a) array(a$like.samples, + dim = c(N, J, n.post.samples)))) + out$like.samples <- aperm(out$like.samples, c(3, 1, 2)) + if (p.occ.re > 0) { + out$sigma.sq.psi.samples <- mcmc( + do.call(rbind, lapply(out.tmp, function(a) t(a$sigma.sq.psi.samples)))) + colnames(out$sigma.sq.psi.samples) <- x.re.names + out$beta.star.samples <- mcmc( + do.call(rbind, lapply(out.tmp, function(a) t(a$beta.star.samples)))) + tmp.names <- unlist(re.level.names) + beta.star.names <- paste(rep(x.re.names, n.occ.re.long), tmp.names, sep = '-') + beta.star.names <- paste(beta.star.names, rep(sp.names, each = n.occ.re), sep = '-') + colnames(out$beta.star.samples) <- beta.star.names + out$re.level.names <- re.level.names + } + # Calculate effective sample sizes + out$ESS <- list() + out$ESS$beta.comm <- effectiveSize(out$beta.comm.samples) + out$ESS$alpha.comm <- effectiveSize(out$alpha.comm.samples) + out$ESS$tau.sq.beta <- effectiveSize(out$tau.sq.beta.samples) + out$ESS$tau.sq.alpha <- effectiveSize(out$tau.sq.alpha.samples) + out$ESS$beta <- effectiveSize(out$beta.samples) + out$ESS$alpha <- effectiveSize(out$alpha.samples) + if (p.occ.re > 0) { + out$ESS$sigma.sq.psi <- effectiveSize(out$sigma.sq.psi.samples) + } + out$X <- X + out$X.p <- X.p + out$X.p.re <- X.p.re + out$X.re <- X.re + out$y <- y.big + out$call <- cl + out$n.samples <- n.samples + out$x.names <- x.names + out$sp.names <- sp.names + out$x.p.names <- x.p.names + out$n.post <- n.post.samples + out$n.thin <- n.thin + out$n.burn <- n.burn + out$n.chains <- n.chains + out$alpha.indx <- alpha.dat.indx.c + 1 + out$alpha.comm.indx <- alpha.comm.indx.r + out$N <- N + out$sites <- sites + out$species <- data$species + # TODO: add in detection random effects eventually + # if (p.det.re > 0) { + # out$pRE <- TRUE + # } else { + out$pRE <- FALSE + # } + if (p.occ.re > 0) { + out$psiRE <- TRUE + } else { + out$psiRE <- FALSE + } class(out) <- "intMsPGOcc" out$run.time <- proc.time() - ptm out } + diff --git a/R/intPGOcc.R b/R/intPGOcc.R index 5ad7f07..411bebb 100644 --- a/R/intPGOcc.R +++ b/R/intPGOcc.R @@ -233,16 +233,16 @@ intPGOcc <- function(occ.formula, det.formula, data, inits, priors, tmp <- parseFormula(det.formula[[i]], data$det.covs[[i]]) X.p[[i]] <- as.matrix(tmp[[1]]) x.p.names[[i]] <- tmp[[2]] - if (ncol(tmp[[4]]) > 0) { + if (ncol(tmp[[4]]) > 0) { X.p.re[[i]] <- as.matrix(tmp[[4]]) x.p.re.names[[i]] <- colnames(X.p.re[[i]]) - p.re.level.names[[i]] <- lapply(data$det.covs[[i]][, x.p.re.names[[i]], drop = FALSE], - function (a) sort(unique(a))) - } else { - X.p.re[[i]] <- matrix(NA, 0, 0) - x.p.re.names[[i]] <- NULL - p.re.level.names[[i]] <- NULL - } + p.re.level.names[[i]] <- lapply(data$det.covs[[i]][, x.p.re.names[[i]], drop = FALSE], + function (a) sort(unique(a))) + } else { + X.p.re[[i]] <- matrix(NA, 0, 0) + x.p.re.names[[i]] <- NULL + p.re.level.names[[i]] <- NULL + } } else { stop(paste("error: det.formula for data source ", i, " is misspecified", sep = '')) } @@ -664,7 +664,7 @@ intPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -712,7 +712,7 @@ intPGOcc <- function(occ.formula, det.formula, data, inits, priors, } alpha.col.indx <- unlist(alpha.col.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -823,12 +823,12 @@ intPGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } @@ -857,11 +857,11 @@ intPGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } storage.mode(chain.info) <- "integer" @@ -1078,7 +1078,7 @@ intPGOcc <- function(occ.formula, det.formula, data, inits, priors, } alpha.col.indx.fit <- unlist(alpha.col.fit.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits.fit <- rnorm(n.det.re.fit, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.n.re.indx.fit <- apply(X.p.re.all.fit, 1, function(a) sum(!is.na(a))) } else { alpha.star.indx.fit <- alpha.star.indx @@ -1095,7 +1095,7 @@ intPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/R/lfJSDM.R b/R/lfJSDM.R index 1b30031..6e857a4 100644 --- a/R/lfJSDM.R +++ b/R/lfJSDM.R @@ -424,7 +424,7 @@ lfJSDM <- function(formula, data, inits, priors, n.factors, n.samples, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } else { sigma.sq.psi.inits <- 0 @@ -522,7 +522,7 @@ lfJSDM <- function(formula, data, inits, priors, n.factors, n.samples, sqrt(tau.sq.beta.inits.list[[i]])), N, p.occ) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } @@ -566,7 +566,7 @@ lfJSDM <- function(formula, data, inits, priors, n.factors, n.samples, } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } } @@ -750,7 +750,7 @@ lfJSDM <- function(formula, data, inits, priors, n.factors, n.samples, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) beta.star.inits.fit <- rep(beta.star.inits.fit, N) re.level.names.fit <- list() diff --git a/R/lfMsPGOcc.R b/R/lfMsPGOcc.R index d25b205..ae8ac2a 100644 --- a/R/lfMsPGOcc.R +++ b/R/lfMsPGOcc.R @@ -777,7 +777,7 @@ lfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } else { sigma.sq.psi.inits <- 0 @@ -809,7 +809,7 @@ lfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } else { sigma.sq.p.inits <- 0 @@ -939,13 +939,13 @@ lfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -996,12 +996,12 @@ lfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, lambda.inits <- c(lambda.inits) if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } } @@ -1234,7 +1234,7 @@ lfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.star.inits.fit <- rep(alpha.star.inits.fit, N) p.re.level.names.fit <- list() @@ -1255,7 +1255,7 @@ lfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) beta.star.inits.fit <- rep(beta.star.inits.fit, N) re.level.names.fit <- list() diff --git a/R/msPGOcc.R b/R/msPGOcc.R index 08a467f..4ee3307 100644 --- a/R/msPGOcc.R +++ b/R/msPGOcc.R @@ -717,7 +717,7 @@ msPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) # Starting values for all species beta.star.inits <- rep(beta.star.inits, N) } else { @@ -748,7 +748,7 @@ msPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } else { sigma.sq.p.inits <- 0 @@ -869,13 +869,13 @@ msPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -914,12 +914,12 @@ msPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits)), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } } @@ -1134,7 +1134,7 @@ msPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.star.inits.fit <- rep(alpha.star.inits.fit, N) p.re.level.names.fit <- list() @@ -1156,7 +1156,7 @@ msPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) beta.star.inits.fit <- rep(beta.star.inits.fit, N) re.level.names.fit <- list() diff --git a/R/postHocLM.R b/R/postHocLM.R index c6e7d1a..acddb69 100644 --- a/R/postHocLM.R +++ b/R/postHocLM.R @@ -284,7 +284,7 @@ postHocLM <- function(formula, data, inits, priors, verbose = FALSE, } } beta.star.indx <- rep(0:(p.re - 1), n.re.long) - beta.star.inits <- rnorm(n.re, sqrt(sigma.sq.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.re, 0, sqrt(sigma.sq.inits[beta.star.indx + 1])) } else { sigma.sq.inits <- 0 beta.star.indx <- 0 @@ -341,7 +341,7 @@ postHocLM <- function(formula, data, inits, priors, verbose = FALSE, tau.sq.inits <- runif(1, 0.5, 10) if (p.re > 0) { sigma.sq.inits <- runif(p.re, 0.5, 10) - beta.star.inits <- rnorm(n.re, sqrt(sigma.sq.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.re, 0, sqrt(sigma.sq.inits[beta.star.indx + 1])) } } storage.mode(chain.info) <- "integer" diff --git a/R/sfJSDM.R b/R/sfJSDM.R index 7cc8a7c..10640e1 100644 --- a/R/sfJSDM.R +++ b/R/sfJSDM.R @@ -655,7 +655,7 @@ sfJSDM <- function(formula, data, inits, priors, } beta.star.inits <- t(beta.star.inits) } else { - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) if (verbose) { message('beta.star is not specified in initial values.\nSetting initial values to random values from the random effects variance\n') @@ -960,7 +960,7 @@ sfJSDM <- function(formula, data, inits, priors, sqrt(tau.sq.beta.inits.list[[i]])), N, p.occ) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } @@ -1020,7 +1020,7 @@ sfJSDM <- function(formula, data, inits, priors, } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } } @@ -1360,7 +1360,7 @@ sfJSDM <- function(formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) beta.star.inits.fit <- rep(beta.star.inits.fit, N) re.level.names.fit <- list() diff --git a/R/sfMsPGOcc.R b/R/sfMsPGOcc.R index 735ef89..879a38e 100644 --- a/R/sfMsPGOcc.R +++ b/R/sfMsPGOcc.R @@ -1020,7 +1020,7 @@ sfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } else { sigma.sq.psi.inits <- 0 @@ -1052,7 +1052,7 @@ sfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } else { sigma.sq.p.inits <- 0 @@ -1298,13 +1298,13 @@ sfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -1367,12 +1367,12 @@ sfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } } @@ -1667,7 +1667,7 @@ sfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.star.inits.fit <- rep(alpha.star.inits.fit, N) p.re.level.names.fit <- list() @@ -1688,7 +1688,7 @@ sfMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) beta.star.inits.fit <- rep(beta.star.inits.fit, N) re.level.names.fit <- list() diff --git a/R/spIntPGOcc.R b/R/spIntPGOcc.R index 0cab015..a0a58bc 100644 --- a/R/spIntPGOcc.R +++ b/R/spIntPGOcc.R @@ -843,7 +843,7 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -891,7 +891,7 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } alpha.col.indx <- unlist(alpha.col.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -1093,12 +1093,12 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } @@ -1142,11 +1142,11 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } storage.mode(chain.info) <- "integer" @@ -1388,7 +1388,7 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } alpha.col.indx.fit <- unlist(alpha.col.fit.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits.fit <- rnorm(n.det.re.fit, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.n.re.indx.fit <- apply(X.p.re.all.fit, 1, function(a) sum(!is.na(a))) } else { alpha.star.indx.fit <- alpha.star.indx @@ -1405,7 +1405,7 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { @@ -1836,12 +1836,12 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } @@ -1886,11 +1886,11 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } storage.mode(chain.info) <- "integer" @@ -2135,7 +2135,7 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } alpha.col.indx.fit <- unlist(alpha.col.fit.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits.fit <- rnorm(n.det.re.fit, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.n.re.indx.fit <- apply(X.p.re.all.fit, 1, function(a) sum(!is.na(a))) } else { alpha.star.indx.fit <- alpha.star.indx @@ -2152,7 +2152,7 @@ spIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/R/spMsPGOcc.R b/R/spMsPGOcc.R index 7b18e90..66e31ad 100644 --- a/R/spMsPGOcc.R +++ b/R/spMsPGOcc.R @@ -973,7 +973,7 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } else { sigma.sq.psi.inits <- 0 @@ -1005,7 +1005,7 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } else { sigma.sq.p.inits <- 0 @@ -1217,13 +1217,13 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -1288,12 +1288,12 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } } @@ -1541,7 +1541,7 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.star.inits.fit <- rep(alpha.star.inits.fit, N) p.re.level.names.fit <- list() @@ -1562,7 +1562,7 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) beta.star.inits.fit <- rep(beta.star.inits.fit, N) re.level.names.fit <- list() @@ -1922,13 +1922,13 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -1994,12 +1994,12 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } } @@ -2270,7 +2270,7 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) alpha.star.inits.fit <- rep(alpha.star.inits.fit, N) p.re.level.names.fit <- list() @@ -2291,7 +2291,7 @@ spMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) beta.star.inits.fit <- rep(beta.star.inits.fit, N) re.level.names.fit <- list() diff --git a/R/spPGOcc.R b/R/spPGOcc.R index 0ff778c..9695d99 100644 --- a/R/spPGOcc.R +++ b/R/spPGOcc.R @@ -796,7 +796,7 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -826,7 +826,7 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -1031,14 +1031,14 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.p')]) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } if (p.occ.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.psi')]) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } @@ -1093,13 +1093,13 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.p')]) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } if (p.occ.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.psi')]) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } } @@ -1280,7 +1280,7 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) p.re.level.names.fit <- list() for (t in 1:p.det.re) { @@ -1300,7 +1300,7 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { @@ -1641,14 +1641,14 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.p')]) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } if (p.occ.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.psi')]) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } @@ -1704,13 +1704,13 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.p')]) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } } if (p.occ.re > 0) { if (!fixed.params[which(all.params == 'sigma.sq.psi')]) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } } @@ -1929,7 +1929,7 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) p.re.level.names.fit <- list() for (t in 1:p.det.re) { @@ -1949,7 +1949,7 @@ spPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/R/stIntPGOcc.R b/R/stIntPGOcc.R index b9cc930..67d5b2b 100644 --- a/R/stIntPGOcc.R +++ b/R/stIntPGOcc.R @@ -926,7 +926,7 @@ stIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -974,7 +974,7 @@ stIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } alpha.col.indx <- unlist(alpha.col.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -1269,12 +1269,12 @@ stIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (ar1) { @@ -1330,11 +1330,11 @@ stIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (ar1) { ar1.vals[5] <- runif(1, rho.a, rho.b) diff --git a/R/stMsPGOcc.R b/R/stMsPGOcc.R index d650a42..53fcbae 100644 --- a/R/stMsPGOcc.R +++ b/R/stMsPGOcc.R @@ -1026,7 +1026,7 @@ stMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } else { sigma.sq.psi.inits <- 0 @@ -1058,7 +1058,7 @@ stMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } else { sigma.sq.p.inits <- 0 @@ -1377,13 +1377,13 @@ stMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -1453,12 +1453,12 @@ stMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } if (ar1) { diff --git a/R/stPGOcc.R b/R/stPGOcc.R index ac14712..fabe8d6 100644 --- a/R/stPGOcc.R +++ b/R/stPGOcc.R @@ -855,7 +855,7 @@ stPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -885,7 +885,7 @@ stPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -1179,12 +1179,12 @@ stPGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (ar1) { @@ -1242,11 +1242,11 @@ stPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (ar1) { ar1.vals[5] <- runif(1, rho.a, rho.b) @@ -1493,7 +1493,7 @@ stPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) p.re.level.names.fit <- list() for (t in 1:p.det.re) { @@ -1513,7 +1513,7 @@ stPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/R/svcMsPGOcc.R b/R/svcMsPGOcc.R index d2068b0..82cd75a 100644 --- a/R/svcMsPGOcc.R +++ b/R/svcMsPGOcc.R @@ -974,7 +974,7 @@ svcMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } else { sigma.sq.psi.inits <- 0 @@ -1006,7 +1006,7 @@ svcMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } else { sigma.sq.p.inits <- 0 @@ -1299,13 +1299,13 @@ svcMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -1374,12 +1374,12 @@ svcMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } } diff --git a/R/svcPGBinom.R b/R/svcPGBinom.R index 70ecd5e..0f5c79c 100644 --- a/R/svcPGBinom.R +++ b/R/svcPGBinom.R @@ -555,7 +555,7 @@ svcPGBinom <- function(formula, data, inits, priors, tuning, } } beta.star.indx <- rep(0:(p.re - 1), n.re.long) - beta.star.inits <- rnorm(n.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -780,7 +780,7 @@ svcPGBinom <- function(formula, data, inits, priors, tuning, beta.inits.list[[i]] <- rnorm(p, mu.beta, sqrt(sigma.beta)) if (p.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.re, + beta.star.inits.list[[i]] <- rnorm(n.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (!fixed.params[which(all.params == 'sigma.sq')]) { @@ -834,7 +834,7 @@ svcPGBinom <- function(formula, data, inits, priors, tuning, } if (p.re > 0) { sigma.sq.psi.inits <- runif(p.re, 0.5, 10) - beta.star.inits <- rnorm(n.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } storage.mode(chain.info) <- "integer" @@ -996,7 +996,7 @@ svcPGBinom <- function(formula, data, inits, priors, tuning, if (p.re > 0) { beta.star.indx.fit <- rep(0:(p.re - 1), n.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.re.fit, + beta.star.inits.fit <- rnorm(n.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.re) { diff --git a/R/svcPGOcc.R b/R/svcPGOcc.R index 7f3659b..20b6ed5 100644 --- a/R/svcPGOcc.R +++ b/R/svcPGOcc.R @@ -847,7 +847,7 @@ svcPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -877,7 +877,7 @@ svcPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -1125,12 +1125,12 @@ svcPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (!fixed.params[which(all.params == 'sigma.sq')]) { @@ -1186,11 +1186,11 @@ svcPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } } storage.mode(chain.info) <- "integer" @@ -1423,7 +1423,7 @@ svcPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) p.re.level.names.fit <- list() for (t in 1:p.det.re) { @@ -1443,7 +1443,7 @@ svcPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/R/svcTIntPGOcc.R b/R/svcTIntPGOcc.R index b4754fe..e7a1dce 100644 --- a/R/svcTIntPGOcc.R +++ b/R/svcTIntPGOcc.R @@ -1006,7 +1006,7 @@ svcTIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -1054,7 +1054,7 @@ svcTIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } alpha.col.indx <- unlist(alpha.col.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -1358,12 +1358,12 @@ svcTIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (sigma.sq.ig) { @@ -1419,11 +1419,11 @@ svcTIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (ar1) { ar1.vals[5] <- runif(1, rho.a, rho.b) diff --git a/R/svcTMsPGOcc.R b/R/svcTMsPGOcc.R index 353aee9..3a0f548 100644 --- a/R/svcTMsPGOcc.R +++ b/R/svcTMsPGOcc.R @@ -1162,7 +1162,7 @@ svcTMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } else { sigma.sq.psi.inits <- 0 @@ -1194,7 +1194,7 @@ svcTMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } else { sigma.sq.p.inits <- 0 @@ -1533,13 +1533,13 @@ svcTMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, sqrt(tau.sq.alpha.inits.list[[i]])), N, p.det) if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) beta.star.inits.list[[i]] <- rep(beta.star.inits.list[[i]], N) } if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) alpha.star.inits.list[[i]] <- rep(alpha.star.inits.list[[i]], N) } @@ -1606,12 +1606,12 @@ svcTMsPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) alpha.star.inits <- rep(alpha.star.inits, N) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) beta.star.inits <- rep(beta.star.inits, N) } if (ar1) { diff --git a/R/svcTPGBinom.R b/R/svcTPGBinom.R index 9936630..e17cc4f 100644 --- a/R/svcTPGBinom.R +++ b/R/svcTPGBinom.R @@ -587,7 +587,7 @@ svcTPGBinom <- function(formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.re - 1), n.re.long) - beta.star.inits <- rnorm(n.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -870,7 +870,7 @@ svcTPGBinom <- function(formula, data, inits, priors, beta.inits.list[[i]] <- rnorm(p, mu.beta, sqrt(sigma.beta)) if (p.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.re, + beta.star.inits.list[[i]] <- rnorm(n.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (sigma.sq.ig) { @@ -922,7 +922,7 @@ svcTPGBinom <- function(formula, data, inits, priors, } if (p.re > 0) { sigma.sq.psi.inits <- runif(p.re, 0.5, 10) - beta.star.inits <- rnorm(n.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (ar1) { ar1.vals[5] <- runif(1, rho.a, rho.b) @@ -1103,7 +1103,7 @@ svcTPGBinom <- function(formula, data, inits, priors, if (p.re > 0) { beta.star.indx.fit <- rep(0:(p.re - 1), n.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.re.fit, + beta.star.inits.fit <- rnorm(n.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.re) { diff --git a/R/svcTPGOcc.R b/R/svcTPGOcc.R index 51a39e7..c594ab1 100644 --- a/R/svcTPGOcc.R +++ b/R/svcTPGOcc.R @@ -926,7 +926,7 @@ svcTPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -956,7 +956,7 @@ svcTPGOcc <- function(occ.formula, det.formula, data, inits, priors, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -1261,12 +1261,12 @@ svcTPGOcc <- function(occ.formula, det.formula, data, inits, priors, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (sigma.sq.ig) { @@ -1323,11 +1323,11 @@ svcTPGOcc <- function(occ.formula, det.formula, data, inits, priors, } if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (ar1) { ar1.vals[5] <- runif(1, rho.a, rho.b) @@ -1584,7 +1584,7 @@ svcTPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) p.re.level.names.fit <- list() for (t in 1:p.det.re) { @@ -1604,7 +1604,7 @@ svcTPGOcc <- function(occ.formula, det.formula, data, inits, priors, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/R/tIntPGOcc.R b/R/tIntPGOcc.R index ddcc1af..a5a0344 100644 --- a/R/tIntPGOcc.R +++ b/R/tIntPGOcc.R @@ -767,7 +767,7 @@ tIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -815,7 +815,7 @@ tIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, } alpha.col.indx <- unlist(alpha.col.list) # Index that indicates the data source the random effect corresponds to. - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -990,12 +990,12 @@ tIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (ar1) { @@ -1031,11 +1031,11 @@ tIntPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, alpha.inits <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (ar1) { ar1.vals[5] <- runif(1, rho.a, rho.b) diff --git a/R/tPGOcc.R b/R/tPGOcc.R index f147321..fcfb945 100644 --- a/R/tPGOcc.R +++ b/R/tPGOcc.R @@ -643,7 +643,7 @@ tPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, } } beta.star.indx <- rep(0:(p.occ.re - 1), n.occ.re.long) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } else { sigma.sq.psi.inits <- 0 beta.star.indx <- 0 @@ -673,7 +673,7 @@ tPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, } } alpha.star.indx <- rep(0:(p.det.re - 1), n.det.re.long) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } else { sigma.sq.p.inits <- 0 alpha.star.indx <- 0 @@ -840,12 +840,12 @@ tPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, alpha.inits.list[[i]] <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits.list[[i]] <- runif(p.det.re, 0.5, 10) - alpha.star.inits.list[[i]] <- rnorm(n.det.re, + alpha.star.inits.list[[i]] <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits.list[[i]][alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits.list[[i]] <- runif(p.occ.re, 0.5, 10) - beta.star.inits.list[[i]] <- rnorm(n.occ.re, + beta.star.inits.list[[i]] <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits.list[[i]][beta.star.indx + 1])) } if (ar1) { @@ -881,11 +881,11 @@ tPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, alpha.inits <- rnorm(p.det, mu.alpha, sqrt(sigma.alpha)) if (p.det.re > 0) { sigma.sq.p.inits <- runif(p.det.re, 0.5, 10) - alpha.star.inits <- rnorm(n.det.re, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) + alpha.star.inits <- rnorm(n.det.re, 0, sqrt(sigma.sq.p.inits[alpha.star.indx + 1])) } if (p.occ.re > 0) { sigma.sq.psi.inits <- runif(p.occ.re, 0.5, 10) - beta.star.inits <- rnorm(n.occ.re, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) + beta.star.inits <- rnorm(n.occ.re, 0, sqrt(sigma.sq.psi.inits[beta.star.indx + 1])) } if (ar1) { ar1.vals[5] <- runif(1, rho.a, rho.b) @@ -1075,7 +1075,7 @@ tPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, if (p.det.re > 0) { alpha.star.indx.fit <- rep(0:(p.det.re - 1), n.det.re.long.fit) alpha.level.indx.fit <- sort(unique(c(X.p.re.fit))) - alpha.star.inits.fit <- rnorm(n.det.re.fit, + alpha.star.inits.fit <- rnorm(n.det.re.fit, 0, sqrt(sigma.sq.p.inits[alpha.star.indx.fit + 1])) p.re.level.names.fit <- list() for (t in 1:p.det.re) { @@ -1095,7 +1095,7 @@ tPGOcc <- function(occ.formula, det.formula, data, inits, priors, tuning, if (p.occ.re > 0) { beta.star.indx.fit <- rep(0:(p.occ.re - 1), n.occ.re.long.fit) beta.level.indx.fit <- sort(unique(c(X.re.fit))) - beta.star.inits.fit <- rnorm(n.occ.re.fit, + beta.star.inits.fit <- rnorm(n.occ.re.fit, 0, sqrt(sigma.sq.psi.inits[beta.star.indx.fit + 1])) re.level.names.fit <- list() for (t in 1:p.occ.re) { diff --git a/man/PGOcc.Rd b/man/PGOcc.Rd index c166a94..b325493 100644 --- a/man/PGOcc.Rd +++ b/man/PGOcc.Rd @@ -238,6 +238,8 @@ inits.list <- list(alpha = 0, beta = 0, n.samples <- 5000 n.report <- 1000 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- PGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov, data = data.list, diff --git a/man/intMsPGOcc.Rd b/man/intMsPGOcc.Rd index ac2ce9d..b02145f 100644 --- a/man/intMsPGOcc.Rd +++ b/man/intMsPGOcc.Rd @@ -11,8 +11,7 @@ Function for fitting integrated multi-species occupancy models using Polya-Gamma intMsPGOcc(occ.formula, det.formula, data, inits, priors, n.samples, n.omp.threads = 1, verbose = TRUE, n.report = 100, n.burn = round(.10 * n.samples), n.thin = 1, n.chains = 1, - parallel.chains = FALSE, k.fold, k.fold.threads = 1, - k.fold.seed, k.fold.only = FALSE, ...) + parallel.chains = FALSE, ...) } \arguments{ @@ -131,26 +130,10 @@ intMsPGOcc(occ.formula, det.formula, data, inits, priors, n.samples, arugment, which implements \emph{within-chain} parallelization for spatially-explicit models. Note if set to \code{TRUE}, all model progress output will be suppressed.} -\item{k.fold}{cross-validation is not currently supported for integrated - multi-species occupancy models.} - -\item{k.fold.threads}{cross-validation is not currently supported for - integrated multi-species occupancy models.} - -\item{k.fold.seed}{cross-validation is not currently supported for integrated - multi-species occupancy models.} - -\item{k.fold.only}{cross-validation is not currently supported for integrated - multi-species occupancy models.} - \item{...}{currently no additional arguments} } -\note{ - Basic functionality of this function is stable, but some components are still in development and not currently available. Please create a GitHub issue on the package GitHub page if you use this function and encounter an error. -} - \references{ Polson, N.G., J.G. Scott, and J. Windle. (2013) Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables. @@ -324,6 +307,8 @@ inits.list <- list(alpha.comm = list(0, 0), beta = 0) # Fit the model. +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- intMsPGOcc(occ.formula = ~ occ.cov.1, det.formula = list(f.1 = ~ det.cov.1.1 + det.cov.1.2 + det.cov.1.3, f.2 = ~ det.cov.2.1 + det.cov.2.2), diff --git a/man/intPGOcc.Rd b/man/intPGOcc.Rd index 8600a7e..4997ee6 100644 --- a/man/intPGOcc.Rd +++ b/man/intPGOcc.Rd @@ -257,6 +257,8 @@ prior.list <- list(beta.normal = list(mean = 0, var = 2.72), alpha.normal = list(mean = list(0, 0, 0, 0), var = list(2.72, 2.72, 2.72, 2.72))) n.samples <- 5000 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- intPGOcc(occ.formula = ~ occ.cov, det.formula = list(f.1 = ~ det.cov.1.1, f.2 = ~ det.cov.2.1, diff --git a/man/lfJSDM.Rd b/man/lfJSDM.Rd index 39418a1..32d8a76 100644 --- a/man/lfJSDM.Rd +++ b/man/lfJSDM.Rd @@ -266,6 +266,8 @@ data.list <- list(y = y[, , 1], prior.list <- list(beta.comm.normal = list(mean = 0, var = 2.72), tau.sq.beta.ig = list(a = 0.1, b = 0.1)) inits.list <- list(beta.comm = 0, beta = 0, tau.sq.beta = 1) +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- lfJSDM(formula = ~ occ.cov.1 + occ.cov.2 + (1 | occ.re.1), data = data.list, inits = inits.list, diff --git a/man/lfMsPGOcc.Rd b/man/lfMsPGOcc.Rd index aeb1e13..6cc500d 100644 --- a/man/lfMsPGOcc.Rd +++ b/man/lfMsPGOcc.Rd @@ -323,6 +323,8 @@ n.samples <- 300 n.burn <- 200 n.thin <- 1 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- lfMsPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1 + det.cov.2 + (1 | det.re), data = data.list, diff --git a/man/msPGOcc.Rd b/man/msPGOcc.Rd index 5c88278..cd9628e 100644 --- a/man/msPGOcc.Rd +++ b/man/msPGOcc.Rd @@ -292,6 +292,8 @@ n.samples <- 3000 n.burn <- 2000 n.thin <- 1 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- msPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/predict.PGOcc.Rd b/man/predict.PGOcc.Rd index 68317be..75456e1 100644 --- a/man/predict.PGOcc.Rd +++ b/man/predict.PGOcc.Rd @@ -92,7 +92,8 @@ inits.list <- list(alpha = rep(0, p.det), n.samples <- 5000 n.report <- 1000 - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- PGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov, data = data.list, diff --git a/man/predict.intMsPGOcc.Rd b/man/predict.intMsPGOcc.Rd index 16f42fe..469b455 100644 --- a/man/predict.intMsPGOcc.Rd +++ b/man/predict.intMsPGOcc.Rd @@ -94,8 +94,8 @@ sp <- FALSE factor.model <- FALSE # Simulate occupancy data dat <- simIntMsOcc(n.data = n.data, J.x = J.x, J.y = J.y, - J.obs = J.obs, n.rep = n.rep, N = N, beta = beta, alpha = alpha, - psi.RE = psi.RE, p.RE = p.RE, sp = sp, factor.model = factor.model, + J.obs = J.obs, n.rep = n.rep, N = N, beta = beta, alpha = alpha, + psi.RE = psi.RE, p.RE = p.RE, sp = sp, factor.model = factor.model, n.factors = n.factors) J <- nrow(dat$coords.obs) y <- dat$y @@ -113,49 +113,49 @@ colnames(occ.covs) <- c('int', 'occ.cov.1') det.covs <- list() # Add covariates one by one det.covs[[1]] <- list(det.cov.1.1 = X.p[[1]][, , 2], - det.cov.1.2 = X.p[[1]][, , 3], - det.cov.1.3 = X.p[[1]][, , 4]) + det.cov.1.2 = X.p[[1]][, , 3], + det.cov.1.3 = X.p[[1]][, , 4]) det.covs[[2]] <- list(det.cov.2.1 = X.p[[2]][, , 2], - det.cov.2.2 = X.p[[2]][, , 3]) + det.cov.2.2 = X.p[[2]][, , 3]) data.list <- list(y = y, - occ.covs = occ.covs, - det.covs = det.covs, + occ.covs = occ.covs, + det.covs = det.covs, sites = sites, species = species) # Take a look at the data.list structure for integrated multi-species # occupancy models. # Priors -prior.list <- list(beta.comm.normal = list(mean = 0, - var = 2.73), - - alpha.comm.normal = list(mean = list(0, 0), - var = list(2.72, 2.72)), +prior.list <- list(beta.comm.normal = list(mean = 0, var = 2.73), + alpha.comm.normal = list(mean = list(0, 0), + var = list(2.72, 2.72)), tau.sq.beta.ig = list(a = 0.1, b = 0.1), tau.sq.alpha.ig = list(a = list(0.1, 0.1), - b = list(0.1, 0.1))) + b = list(0.1, 0.1))) inits.list <- list(alpha.comm = list(0, 0), - beta.comm = 0, - tau.sq.beta = 1, - tau.sq.alpha = list(1, 1), + beta.comm = 0, + tau.sq.beta = 1, + tau.sq.alpha = list(1, 1), alpha = list(a = matrix(rnorm(p.det.long[1] * N[1]), N[1], p.det.long[1]), - b = matrix(rnorm(p.det.long[2] * N[2]), N[2], p.det.long[2])), - beta = 0) + b = matrix(rnorm(p.det.long[2] * N[2]), N[2], p.det.long[2])), + beta = 0) -# Fit the model. +# Fit the model. +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- intMsPGOcc(occ.formula = ~ occ.cov.1, det.formula = list(f.1 = ~ det.cov.1.1 + det.cov.1.2 + det.cov.1.3, f.2 = ~ det.cov.2.1 + det.cov.2.2), - inits = inits.list, - priors = prior.list, - data = data.list, - n.samples = 100, - n.omp.threads = 1, - verbose = TRUE, - n.report = 10, - n.burn = 50, - n.thin = 1, - n.chains = 1) + inits = inits.list, + priors = prior.list, + data = data.list, + n.samples = 100, + n.omp.threads = 1, + verbose = TRUE, + n.report = 10, + n.burn = 50, + n.thin = 1, + n.chains = 1) #Predict at new locations. X.0 <- dat$X.pred psi.0 <- dat$psi.pred @@ -164,7 +164,7 @@ out.pred <- predict(out, X.0, ignore.RE = TRUE) # Create prediction for one species. curr.sp <- 2 psi.hat.quants <- apply(out.pred$psi.0.samples[,curr.sp, ], - 2, quantile, c(0.025, 0.5, 0.975)) + 2, quantile, c(0.025, 0.5, 0.975)) plot(psi.0[curr.sp, ], psi.hat.quants[2, ], pch = 19, xlab = 'True', ylab = 'Predicted', ylim = c(min(psi.hat.quants), max(psi.hat.quants)), main = paste("Species ", curr.sp, sep = '')) diff --git a/man/predict.intPGOcc.Rd b/man/predict.intPGOcc.Rd index df359a0..33872d1 100644 --- a/man/predict.intPGOcc.Rd +++ b/man/predict.intPGOcc.Rd @@ -100,6 +100,8 @@ inits.list <- list(alpha = list(0, 0, 0, 0), prior.list <- list(beta.normal = list(mean = 0, var = 2.72), alpha.normal = list(mean = list(0, 0, 0, 0), var = list(2.72, 2.72, 2.72, 2.72))) +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. n.samples <- 5000 out <- intPGOcc(occ.formula = ~ occ.cov, det.formula = list(f.1 = ~ det.cov.1.1, diff --git a/man/predict.lfJSDM.Rd b/man/predict.lfJSDM.Rd index 0edb30b..a3fef5a 100644 --- a/man/predict.lfJSDM.Rd +++ b/man/predict.lfJSDM.Rd @@ -112,7 +112,8 @@ inits.list <- list(alpha.comm = 0, beta = 0, tau.sq.beta = 1, lambda = lambda.inits) - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- lfJSDM(formula = ~ occ.cov, data = data.list, inits = inits.list, diff --git a/man/predict.lfMsPGOcc.Rd b/man/predict.lfMsPGOcc.Rd index d22397b..ccb87a5 100644 --- a/man/predict.lfMsPGOcc.Rd +++ b/man/predict.lfMsPGOcc.Rd @@ -130,7 +130,8 @@ inits.list <- list(alpha.comm = 0, tau.sq.alpha = 1, lambda = lambda.inits, z = apply(y, c(1, 2), max, na.rm = TRUE)) - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- lfMsPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/predict.msPGOcc.Rd b/man/predict.msPGOcc.Rd index 044823a..6394d79 100644 --- a/man/predict.msPGOcc.Rd +++ b/man/predict.msPGOcc.Rd @@ -112,7 +112,8 @@ inits.list <- list(alpha.comm = 0, tau.sq.beta = 1, tau.sq.alpha = 1, z = apply(y, c(1, 2), max, na.rm = TRUE)) - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- msPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/predict.sfJSDM.Rd b/man/predict.sfJSDM.Rd index f499487..6dc4803 100644 --- a/man/predict.sfJSDM.Rd +++ b/man/predict.sfJSDM.Rd @@ -145,7 +145,8 @@ inits.list <- list(beta.comm = 0, lambda = lambda.inits) # Tuning tuning.list <- list(phi = 1) - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- sfJSDM(formula = ~ occ.cov, data = data.list, inits = inits.list, diff --git a/man/predict.sfMsPGOcc.Rd b/man/predict.sfMsPGOcc.Rd index 19422b1..becfe96 100644 --- a/man/predict.sfMsPGOcc.Rd +++ b/man/predict.sfMsPGOcc.Rd @@ -167,7 +167,8 @@ inits.list <- list(alpha.comm = 0, z = apply(y, c(1, 2), max, na.rm = TRUE)) # Tuning tuning.list <- list(phi = 1) - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- sfMsPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/predict.spIntPGOcc.Rd b/man/predict.spIntPGOcc.Rd index 635dc91..fad49a3 100644 --- a/man/predict.spIntPGOcc.Rd +++ b/man/predict.spIntPGOcc.Rd @@ -149,7 +149,8 @@ tuning.list <- list(phi = 1) n.batch <- 40 # Batch length batch.length <- 25 - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- spIntPGOcc(occ.formula = ~ occ.cov, det.formula = list(f.1 = ~ det.cov.1.1, f.2 = ~ det.cov.2.1 + det.cov.2.2, diff --git a/man/predict.spMsPGOcc.Rd b/man/predict.spMsPGOcc.Rd index ef9f23b..60bec78 100644 --- a/man/predict.spMsPGOcc.Rd +++ b/man/predict.spMsPGOcc.Rd @@ -162,25 +162,26 @@ inits.list <- list(alpha.comm = 0, z = apply(y, c(1, 2), max, na.rm = TRUE)) # Tuning tuning.list <- list(phi = 1) - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- spMsPGOcc(occ.formula = ~ occ.cov, - det.formula = ~ det.cov.1 + det.cov.2, - data = data.list, - inits = inits.list, - n.batch = n.batch, - batch.length = batch.length, - accept.rate = 0.43, - priors = prior.list, - cov.model = "exponential", - tuning = tuning.list, - n.omp.threads = 1, - verbose = TRUE, - NNGP = TRUE, - n.neighbors = 5, - search.type = 'cb', - n.report = 10, - n.burn = 500, - n.thin = 1) + det.formula = ~ det.cov.1 + det.cov.2, + data = data.list, + inits = inits.list, + n.batch = n.batch, + batch.length = batch.length, + accept.rate = 0.43, + priors = prior.list, + cov.model = "exponential", + tuning = tuning.list, + n.omp.threads = 1, + verbose = TRUE, + NNGP = TRUE, + n.neighbors = 5, + search.type = 'cb', + n.report = 10, + n.burn = 500, + n.thin = 1) summary(out, level = 'both') diff --git a/man/predict.spPGOcc.Rd b/man/predict.spPGOcc.Rd index 8063c40..01fe2b0 100644 --- a/man/predict.spPGOcc.Rd +++ b/man/predict.spPGOcc.Rd @@ -141,7 +141,8 @@ inits.list <- list(alpha = 0, beta = 0, z = apply(y, 1, max, na.rm = TRUE)) # Tuning tuning.list <- list(phi = 1) - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- spPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1, data = data.list, diff --git a/man/predict.stMsPGOcc.Rd b/man/predict.stMsPGOcc.Rd index 3a0998d..7385d3a 100644 --- a/man/predict.stMsPGOcc.Rd +++ b/man/predict.stMsPGOcc.Rd @@ -181,7 +181,8 @@ batch.length <- 25 n.burn <- 25 n.thin <- 1 n.samples <- n.batch * batch.length - +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- stMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, @@ -189,19 +190,19 @@ out <- stMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, n.batch = n.batch, batch.length = batch.length, accept.rate = 0.43, - ar1 = TRUE, - NNGP = TRUE, - n.neighbors = 5, - n.factors = n.factors, - cov.model = 'exponential', + ar1 = TRUE, + NNGP = TRUE, + n.neighbors = 5, + n.factors = n.factors, + cov.model = 'exponential', priors = prior.list, tuning = tuning.list, n.omp.threads = 1, verbose = TRUE, n.report = 1, n.burn = n.burn, - n.thin = n.thin, - n.chains = 1) + n.thin = n.thin, + n.chains = 1) summary(out) diff --git a/man/predict.stPGOcc.Rd b/man/predict.stPGOcc.Rd index 6010ced..4cd804e 100644 --- a/man/predict.stPGOcc.Rd +++ b/man/predict.stPGOcc.Rd @@ -162,6 +162,8 @@ batch.length <- 25 n.iter <- n.batch * batch.length # Run the model +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- stPGOcc(occ.formula = ~ trend + occ.cov.1, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/predict.svcMsPGOcc.Rd b/man/predict.svcMsPGOcc.Rd index f29d415..1204925 100644 --- a/man/predict.svcMsPGOcc.Rd +++ b/man/predict.svcMsPGOcc.Rd @@ -176,18 +176,20 @@ n.burn <- 0 n.thin <- 1 n.samples <- n.batch * batch.length +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2 + occ.cov.3 + occ.cov.4, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, inits = inits.list, n.batch = n.batch, - n.factors = n.factors, + n.factors = n.factors, batch.length = batch.length, - std.by.sp = TRUE, + std.by.sp = TRUE, accept.rate = 0.43, priors = prior.list, - svc.cols = svc.cols, + svc.cols = svc.cols, cov.model = "spherical", tuning = tuning.list, n.omp.threads = 1, @@ -197,8 +199,8 @@ out <- svcMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2 + occ.cov.3 + search.type = 'cb', n.report = 10, n.burn = n.burn, - n.thin = n.thin, - n.chains = 1) + n.thin = n.thin, + n.chains = 1) summary(out) # Predict at new locations ------------------------------------------------ diff --git a/man/predict.svcPGBinom.Rd b/man/predict.svcPGBinom.Rd index 4c8c0dc..48791d3 100644 --- a/man/predict.svcPGBinom.Rd +++ b/man/predict.svcPGBinom.Rd @@ -140,6 +140,8 @@ batch.length <- 25 n.burn <- 100 n.thin <- 1 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcPGBinom(formula = ~ cov.1 + cov.2 + cov.3, svc.cols = c(1, 2), data = data.list, diff --git a/man/predict.svcPGOcc.Rd b/man/predict.svcPGOcc.Rd index 3da980a..75452a0 100644 --- a/man/predict.svcPGOcc.Rd +++ b/man/predict.svcPGOcc.Rd @@ -148,6 +148,8 @@ inits.list <- list(alpha = 0, beta = 0, # Tuning tuning.list <- list(phi = 1) +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1, data = data.list, diff --git a/man/predict.svcTMsPGOcc.Rd b/man/predict.svcTMsPGOcc.Rd index 5e714b8..f475c82 100644 --- a/man/predict.svcTMsPGOcc.Rd +++ b/man/predict.svcTMsPGOcc.Rd @@ -182,6 +182,8 @@ n.burn <- 25 n.thin <- 1 n.samples <- n.batch * batch.length +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcTMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, @@ -189,20 +191,20 @@ out <- svcTMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, n.batch = n.batch, batch.length = batch.length, accept.rate = 0.43, - ar1 = TRUE, - svc.cols = svc.cols, - NNGP = TRUE, - n.neighbors = 5, - n.factors = n.factors, - cov.model = 'exponential', + ar1 = TRUE, + svc.cols = svc.cols, + NNGP = TRUE, + n.neighbors = 5, + n.factors = n.factors, + cov.model = 'exponential', priors = prior.list, tuning = tuning.list, n.omp.threads = 1, verbose = TRUE, n.report = 1, n.burn = n.burn, - n.thin = n.thin, - n.chains = 1) + n.thin = n.thin, + n.chains = 1) summary(out) diff --git a/man/predict.svcTPGBinom.Rd b/man/predict.svcTPGBinom.Rd index 4ae9847..1248065 100644 --- a/man/predict.svcTPGBinom.Rd +++ b/man/predict.svcTPGBinom.Rd @@ -154,6 +154,8 @@ n.burn <- 0 n.thin <- 1 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcTPGBinom(formula = ~ trend + cov.1 + cov.2, svc.cols = svc.cols, data = data.list, diff --git a/man/predict.svcTPGOcc.Rd b/man/predict.svcTPGOcc.Rd index 36da862..d7b39a6 100644 --- a/man/predict.svcTPGOcc.Rd +++ b/man/predict.svcTPGOcc.Rd @@ -169,6 +169,8 @@ batch.length <- 25 n.iter <- n.batch * batch.length # Run the model +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcTPGOcc(occ.formula = ~ trend + occ.cov.1, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/predict.tMsPGOcc.Rd b/man/predict.tMsPGOcc.Rd index ff5c1fc..23d587e 100644 --- a/man/predict.tMsPGOcc.Rd +++ b/man/predict.tMsPGOcc.Rd @@ -139,6 +139,8 @@ n.burn <- 25 n.thin <- 1 n.samples <- n.batch * batch.length +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- tMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/predict.tPGOcc.Rd b/man/predict.tPGOcc.Rd index 6bf792c..7da4e59 100644 --- a/man/predict.tPGOcc.Rd +++ b/man/predict.tPGOcc.Rd @@ -122,6 +122,8 @@ n.burn <- 2000 n.thin <- 1 # Run the model +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- tPGOcc(occ.formula = ~ trend + occ.cov.1, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/sfJSDM.Rd b/man/sfJSDM.Rd index d6f8d2c..1881e65 100644 --- a/man/sfJSDM.Rd +++ b/man/sfJSDM.Rd @@ -419,6 +419,8 @@ n.batch <- 5 n.report <- 100 formula <- ~ 1 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- sfJSDM(formula = formula, data = data.list, inits = inits.list, diff --git a/man/sfMsPGOcc.Rd b/man/sfMsPGOcc.Rd index fcfd71a..6d045b7 100644 --- a/man/sfMsPGOcc.Rd +++ b/man/sfMsPGOcc.Rd @@ -438,6 +438,8 @@ inits.list <- list(alpha.comm = 0, # Tuning tuning.list <- list(phi = 1) +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- sfMsPGOcc(occ.formula = ~ occ.cov + (1 | occ.re), det.formula = ~ det.cov.1 + det.cov.2 + (1 | det.re), data = data.list, diff --git a/man/spIntPGOcc.Rd b/man/spIntPGOcc.Rd index 22b670d..f63a9c2 100644 --- a/man/spIntPGOcc.Rd +++ b/man/spIntPGOcc.Rd @@ -370,6 +370,8 @@ n.batch <- 2 # Batch length batch.length <- 25 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- spIntPGOcc(occ.formula = ~ occ.cov, det.formula = list(f.1 = ~ det.cov.1.1, f.2 = ~ det.cov.2.1 + det.cov.2.2, diff --git a/man/spMsPGOcc.Rd b/man/spMsPGOcc.Rd index a70b519..c2ff235 100644 --- a/man/spMsPGOcc.Rd +++ b/man/spMsPGOcc.Rd @@ -391,6 +391,8 @@ inits.list <- list(alpha.comm = 0, # Tuning tuning.list <- list(phi = 1) +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- spMsPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/spPGOcc.Rd b/man/spPGOcc.Rd index 1ce37f8..6a0bdf2 100644 --- a/man/spPGOcc.Rd +++ b/man/spPGOcc.Rd @@ -348,6 +348,8 @@ inits.list <- list(alpha = 0, beta = 0, # Tuning tuning.list <- list(phi = 1) +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- spPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/stMsPGOcc.Rd b/man/stMsPGOcc.Rd index 45f480d..74c3f98 100644 --- a/man/stMsPGOcc.Rd +++ b/man/stMsPGOcc.Rd @@ -437,6 +437,8 @@ n.burn <- 25 n.thin <- 1 n.samples <- n.batch * batch.length +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- stMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/stPGOcc.Rd b/man/stPGOcc.Rd index b8f6841..3572e9c 100644 --- a/man/stPGOcc.Rd +++ b/man/stPGOcc.Rd @@ -388,6 +388,8 @@ batch.length <- 25 n.iter <- n.batch * batch.length # Run the model +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- stPGOcc(occ.formula = ~ trend + occ.cov.1, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/svcMsPGOcc.Rd b/man/svcMsPGOcc.Rd index 0cdc543..55c7e93 100644 --- a/man/svcMsPGOcc.Rd +++ b/man/svcMsPGOcc.Rd @@ -435,6 +435,8 @@ n.burn <- 0 n.thin <- 1 n.samples <- n.batch * batch.length +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2 + occ.cov.3 + occ.cov.4 + (1 | occ.factor.1), det.formula = ~ det.cov.1 + det.cov.2 + (1 | det.factor.1), diff --git a/man/svcPGBinom.Rd b/man/svcPGBinom.Rd index 870cb48..84f2915 100644 --- a/man/svcPGBinom.Rd +++ b/man/svcPGBinom.Rd @@ -334,6 +334,8 @@ batch.length <- 25 n.burn <- 100 n.thin <- 1 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcPGBinom(formula = ~ cov.1 + cov.2 + cov.3, svc.cols = c(1, 2), data = data.list, diff --git a/man/svcPGOcc.Rd b/man/svcPGOcc.Rd index 0f39720..0bf0fa4 100644 --- a/man/svcPGOcc.Rd +++ b/man/svcPGOcc.Rd @@ -367,6 +367,8 @@ inits.list <- list(alpha = 0, beta = 0, # Tuning tuning.list <- list(phi = 1) +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcPGOcc(occ.formula = ~ occ.cov, det.formula = ~ det.cov.1, data = data.list, diff --git a/man/svcTMsPGOcc.Rd b/man/svcTMsPGOcc.Rd index 0ccefc1..20cf608 100644 --- a/man/svcTMsPGOcc.Rd +++ b/man/svcTMsPGOcc.Rd @@ -478,6 +478,8 @@ n.burn <- 25 n.thin <- 1 n.samples <- n.batch * batch.length +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcTMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/svcTPGBinom.Rd b/man/svcTPGBinom.Rd index c384ddf..b20066d 100644 --- a/man/svcTPGBinom.Rd +++ b/man/svcTPGBinom.Rd @@ -382,6 +382,8 @@ n.batch <- 2 n.burn <- 0 n.thin <- 1 +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcTPGBinom(formula = ~ trend + cov.1 + cov.2, svc.cols = svc.cols, data = data.list, diff --git a/man/svcTPGOcc.Rd b/man/svcTPGOcc.Rd index 568e3b2..9ecf506 100644 --- a/man/svcTPGOcc.Rd +++ b/man/svcTPGOcc.Rd @@ -414,6 +414,8 @@ n.burn <- 0 n.thin <- 1 # Run the model +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- svcTPGOcc(occ.formula = ~ trend + occ.cov.1 + occ.cov.2, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/tMsPGOcc.Rd b/man/tMsPGOcc.Rd index 6e8a2ad..da001e8 100644 --- a/man/tMsPGOcc.Rd +++ b/man/tMsPGOcc.Rd @@ -333,6 +333,8 @@ n.burn <- 25 n.thin <- 1 n.samples <- n.batch * batch.length +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- tMsPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/man/tPGOcc.Rd b/man/tPGOcc.Rd index 15db5f6..c2e2529 100644 --- a/man/tPGOcc.Rd +++ b/man/tPGOcc.Rd @@ -328,6 +328,8 @@ n.burn <- 100 n.thin <- 1 # Run the model +# Note that this is just a test case and more iterations/chains may need to +# be run to ensure convergence. out <- tPGOcc(occ.formula = ~ trend + occ.cov.1, det.formula = ~ det.cov.1 + det.cov.2, data = data.list, diff --git a/src/intMsPGOcc.cpp b/src/intMsPGOcc.cpp index 4b46377..82c098e 100644 --- a/src/intMsPGOcc.cpp +++ b/src/intMsPGOcc.cpp @@ -787,5 +787,3 @@ extern "C" { return(result_r); } } - - diff --git a/src/intPGOcc.cpp b/src/intPGOcc.cpp index 494b85e..a62795e 100755 --- a/src/intPGOcc.cpp +++ b/src/intPGOcc.cpp @@ -350,8 +350,6 @@ extern "C" { GetRNGstate(); for (s = 0; s < nSamples; s++) { - // for (s = 0; s < 1; s++) { - // for (s = 0; s < 1; s++) { /******************************************************************** *Update Occupancy Auxiliary Variables *******************************************************************/ diff --git a/vignettes/parallelization.Rmd b/vignettes/parallelization.Rmd new file mode 100644 index 0000000..301d8cc --- /dev/null +++ b/vignettes/parallelization.Rmd @@ -0,0 +1,216 @@ +--- +title: "Parallelization in spOccupancy and spAbundance" +author: "Jeffrey W. Doser" +date: "2024 (last update August 5, 2024)" +description: Learn about how to use parallelization to speed up your models +output: + rmarkdown::html_vignette: + toc: true + toc_depth: 3 +vignette: > + %\VignetteIndexEntry{parallelization} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +options(rmarkdown.html_vignette.check_title = FALSE) +knitr::opts_chunk$set( + comment = "", cache = TRUE +) +``` + +\newcommand{\bm}{\boldsymbol} + +# Introduction + +This vignette describes options for fitting `spOccupancy` and `spAbundance` models in parallel. As of `spOccupancy` v0.8.0 and `spAbundance` v0.2.0, all model-fitting functions have two arguments related to parallelization: `n.omp.threads` and `parallel.chains`. This document will describe what the two arguments do and when to use one vs. the other. I wrote this short article just to give a bit more detail on the two forms of parallelization and in what situations one approach might be better than the other. + +# `parallel.chains` vs. `n.omp.threads` + +The arguments `parallel.chains` and `n.omp.threads` both focus on fitting models in parallel to speed up model run times, but they do so in completely different ways. The newer `parallel.chains` argument is likely the form of parallelization that most users of JAGS, Stan, and/or NIMBLE are most familiar with. When setting `parallel.chains = TRUE`, the `n.chains` MCMC chains specified in the function arguments will be run simultaneously in parallel across different cores. Thus, if you set `n.chains = 3` and `parallel.chains = TRUE`, your computer will attempt to initiate three cores to run each of the individual MCMC chains. This results in substantial improvements in run time as opposed to running the MCMC chains sequentially. The downside to using `parallel.chains = TRUE` is that the model progress that is normally printed to the screen when `verbose = TRUE` will no longer show up. To get a sense of how long a model run in parallel will take, you can pretty simply just run a single MCMC chain of the model for a few iterations, see how long that takes, and then multiply that number accordingly based on the total number of MCMC iterations you plan on running for the full model run. This will not be a perfect representation of model run time, but it should give a general sense of how long the parallel model will take to run. + +The `n.omp.threads` argument implements *within-chain* parallelization using OpenMP. The situations in which `n.omp.threads` will lead to improvements in model performance are a bit more nuanced than when using `parallel.chains`. **First, `n.omp.threads` only results in changes in model run times for spatial models (i.e., models that include spatial random effects).** The `n.omp.threads` argument is still included in all non-spatial models in both `spOccupancy` and `spAbundance`, but it currently does not do anything, so setting this to any value above 1 will not actually result in any differences in run time. For spatial models, setting `n.omp.threads` to a number greater than 1 can result in fairly substantial decreases in model run time for models that use a Nearest Neighbor Gaussian Process (`NNGP = TRUE`). However, the benefits that this provides will depend on three things: (1) the number of spatial locations; (2) how many neighbors are used in the NNGP approximation; and (3) the number of species (for multi-species models). Generally, increasing `n.omp.threads` to a value greater than 1 (usually a value between 5-10 is the sweet spot) will improve performance more when working with a larger number of spatial locations. This becomes clear when considering how the parallelization is implemented. The within-chain parallellization is done when updating each of the $J$ spatial random effects. Essentially, if running a model with $J$ sites and setting `n.omp.threads = 5`, the model will update $\frac{J}{5}$ of the spatial random effects simultaneously across 5 different cores. The larger $J$ is, the more benefits this provides. Similarly, I have found that the within-chain parallelization will only result in improvements when the number of neighbors is 10 or higher. So, if you're fitting a model with less than 10 neighbors, setting `n.omp.threads` to a value greater than 1 may not result in much (if any) improvements in run time. Lastly, for multi-species models, the more species (or the more factors in spatial factor models) you have, generally the more benefits you will see from using `n.omp.threads` with a value greater than 1. + +# Should I set `parallel.chains = TRUE` and `n.omp.threads` to a value greater than 1? + +While perhaps a bit counterintuitive, setting `parallel.chains = TRUE` and `n.omp.threads` > 1 will actually not give improved computational performance, and the resulting model fit may actually take longer than if you did not use any form of parallelization! This is a result of how the parallelization is implemented in the two approaches. I generally do not recommend using both of these arguments together. However, as I will point out later in the document, one can fairly easily run multiple chains in parallel by establishing different R sessions (either via RStudio or via the command line) and running a separate chian in each session, and then within each of those chains `n.omp.threads` can be set. Such an approach does allow you to get the best of both worlds (i.e., maximize both within and across chain parallelization). + +# Which type of parallelization should I use? + +For the large majority of users, `parallel.chains` will provide more decreases in model run time than `n.omp.threads`. For fitting any type of non-spatial model, users should use `parallel.chains`. When fitting spatially-explicit models with tens of thousands of locations, using `n.omp.threads` may also provide substantial improvements in run time, in which case one might want to compare some shorter runs using `parallel.chains` vs. `n.omp.threads` to determine which provides the best speed up. For most users, I recommend first fitting some initial small model runs without using any form of parallelization to just make sure things are running properly. Then, setting `parallel.chains = TRUE` for the full model run is a very easy way to substantially improve model run times for nearly all situations when fitting models in `spOccupancy` or `spAbundance`. + +# A simple example with a single-species spatial occupancy model + +Below, I provide some code (without a whole lot of explanation) for fitting and comparing run times for a spatial occupancy model using a simulated data set: (1) no parallelizaition; (2) setting `n.omp.threads = 5`; (3) setting `parallel.chains = TRUE`; and (4) setting `n.omp.threads = 5` and `parallel.chains = TRUE`. Note that I do not recommend doing the final option, and my reason for doing this here is just to show how it can actually substantially slow things down. Note that run times may be different on your machine, but that the relative performance of the four models should be similar across devices. + +```{r} +library(spOccupancy) +set.seed(11111) +# Simulate the data (625 spatial locations) +J.x <- 25 +J.y <- 25 +J <- J.x * J.y +n.rep <- sample(1, J, replace = TRUE) +beta <- c(-1, 0.2, 0.3, -0.3) +p.occ <- length(beta) +alpha <- c(0.3, 0.5) +p.det <- length(alpha) +phi <- 3 / .5 +sigma.sq <- 1.5 +p.RE <- list() +psi.RE <- list() +sp <- TRUE +cov.model = 'exponential' +x.positive <- FALSE +svc.cols <- 1 +n.rep.max <- max(n.rep) + +dat <- simOcc(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha, + p.RE = p.RE, psi.RE = psi.RE, sigma.sq = sigma.sq, phi = phi, sp = TRUE, + cov.model = cov.model) + +# Package everything up into list for spOccupancy +y <- dat$y +X.p <- dat$X.p +X <- dat$X +coords <- dat$coords + +# Package all data into a list +occ.covs <- cbind(X) +colnames(occ.covs) <- c('int', 'occ.cov.1', 'occ.cov.2', 'occ.cov.3') +det.covs <- list(int = X.p[, , 1], + det.cov.1 = X.p[, , 2]) +data.list <- list(y = y, + occ.covs = occ.covs, + det.covs = det.covs, + coords = coords) + +# MCMC settings +n.batch <- 400 +batch.length <- 25 +n.burn <- 5000 +n.thin <- 5 +n.chains <- 3 + +# Priors +prior.list <- list(sigma.sq.ig = c(2, 1), + phi.unif = c(3 / 1, 3 / .1)) + +# Starting values +inits.list <- list(alpha = 0, + beta = 0, + phi = 3 / .5, + sigma.sq = 1, + nu = 1, + z = apply(y, 1, max, na.rm = TRUE)) +# Tuning +tuning.list <- list(phi = 0.5) + + +# Model 1: no parallelization --------------------------------------------- +out.1 <- spPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2 + occ.cov.3, + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + cov.model = "exponential", + tuning = tuning.list, + verbose = FALSE, + NNGP = TRUE, + n.neighbors = 15, + n.report = 5, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains, + parallel.chains = FALSE, + n.omp.threads = 1) + +# Model 2: within-chain parallelization ----------------------------------- +out.2 <- spPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2 + occ.cov.3, + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + cov.model = "exponential", + tuning = tuning.list, + verbose = FALSE, + NNGP = TRUE, + n.neighbors = 15, + n.report = 5, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains, + parallel.chains = FALSE, + n.omp.threads = 5) + +# Model 3: across chain parallelization ----------------------------------- +out.3 <- spPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2 + occ.cov.3, + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + cov.model = "exponential", + tuning = tuning.list, + verbose = FALSE, + NNGP = TRUE, + n.neighbors = 15, + n.report = 5, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains, + parallel.chains = TRUE, + n.omp.threads = 1) +# Model 4: within and across chain parallelization ------------------------ +out.4 <- spPGOcc(occ.formula = ~ occ.cov.1 + occ.cov.2 + occ.cov.3, + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + cov.model = "exponential", + tuning = tuning.list, + verbose = FALSE, + NNGP = TRUE, + n.neighbors = 15, + n.report = 5, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains, + parallel.chains = TRUE, + n.omp.threads = 5) + +# Compare model run times ------------------------------------------------- +# Model 1: no parallelization +out.1$run.time +# Model 2: within-chain parallelization (n.omp.threads) +out.2$run.time +# Model 3: across-chain parallelization (parallel.chains = TRUE) +out.3$run.time +# Model 4: within and across chain parallelization +# NOTE: this is far slower than even the model with no parallelization!!! +out.4$run.time +``` + +We can see in this simple example, using `parallel.chains = TRUE` improved model run time a bit more than using `n.omp.threads = 5`, while using both of them simultaneously did not result in any improvements in run time. + +# Using `n.omp.threads` while running chains in parallel + +While using `n.omp.threads` and `parallel.chains` generally does not improve model run time (like we saw in the previous example), that does not mean we cannot use another approach to leverage both forms of parallelization. Doing this will just require us to do a bit more as opposed to using a simple argument in the function. + +Say for example we want to run three chains of our model. Instead of running chains in parallel using `parallel.chains`, we can instead simply run multiple instances of our R script in different R sessions, where in each instance of the script we only run a single chain (i.e., we set `n.chains = 1`). We can then run three different instances of the script simultaneously, which will effectively run the chains in parallel. If we use this approach, we can also set `n.omp.threads` to a value greater than 1, and we will get the benefits of both *within-chain* parallelization and across chain parallelization, which we did not see when trying to use `parallel.chains` and `n.omp.threads` simultaneously. The three instances of the R script should save out the resulting object from the model-fitting function, saving it in a form that allows you to distinguish that chain from the other chains. You will then have to subsequently combine the results from the different chains after running the chains, which admittedly can be a bit tricky, depending on what your objectives our with the analysis. The chains can be run simultaneously in a few different ways. Perhaps the most broadly applicable approach is to run the chains directly from a terminal window using `Rscript` or `R CMD BATCH`, which will be very suitable for people trying to run `spOccupancy` or `spAbundance` models on High Performance Computers. Altneratively, if using `RStudio`, one could simply open up three different RStudio windows, or could leverage the "Background Jobs" tab in a single Rstudio script to easily run multiple instances of the script. + +[This GitHub repository](https://github.com/doserjef/Switzerland24-Spatial-Workshop) contains materials from a workshop in which there is an example of how to do this. In particular, the files that start with `12` give an example of how to do this with a multi-species spatial occupancy model. + +