Update a spOccupancy or spAbundance model run with more MCMC iterations
updateMCMC.RdFunction for updating a previously run spOccupancy or spAbundance model with additional MCMC iterations. This function is useful for situations where a model is run for a long time but convergence/adequate mixing of the MCMC chains is not reached. Instead of re-running the entire model again, this function allows you to pick up where you left off. This function is currently in development, and only currently works with the following spOccupancy and spAbundance model objects: msAbund and sfJSDM. Note that cross-validation is not possible when updating the model.
Usage
updateMCMC(object, n.batch, n.samples, n.burn = 0, n.thin,
keep.orig = TRUE, verbose = TRUE, n.report = 100, ...)Arguments
- object
a
spOccupancyorspAbundancemodel object. Currently supports objects of classmsAbundandsfJSDM.- n.batch
the number of additional MCMC batches in each chain to run for the adaptive MCMC sampler. Only valid for model types fit with an adaptive MCMC sampler
- n.samples
the number of posterior samples to collect in each chain. Only valid for model types that are run with a fully Gibbs sampler and have
n.samplesas an argument in the original model fitting function.- n.burn
the number of samples out of the total
n.batch * batchlengthto discard as burn-in for each chain from the updated samples. Note this argument does not discard samples from the previous model run, and rather only applies to the samples in the updated run of the model. Defaults to 0- n.thin
the thinning interval for collection of MCMC samples in the updated model run. The thinning occurs after the
n.burnsamples are discarded. Default value is set to 1.- keep.orig
A logical value indicating whether or not the samples from the original run of the model should be kept or discarded.
- verbose
if
TRUE, messages about data preparation, model specification, and progress of the sampler are printed to the screen. Otherwise, no messages are printed.- n.report
the interval to report Metropolis sampler acceptance and MCMC progress.
- ...
currently no additional arguments
Author
Jeffrey W. Doser doserjef@msu.edu,
Value
An object of the same class as the original model fit provided in the argument object. See the manual page for the original model type for complete details.
Examples
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep<- sample(2:4, size = J, replace = TRUE)
N <- 6
# Community-level covariate effects
# Occurrence
beta.mean <- c(0.2)
p.occ <- length(beta.mean)
tau.sq.beta <- c(0.6)
# Detection
alpha.mean <- c(0)
tau.sq.alpha <- c(1)
p.det <- length(alpha.mean)
# Random effects
psi.RE <- list()
p.RE <- list()
# Draw species-level effects from community means.
beta <- matrix(NA, nrow = N, ncol = p.occ)
alpha <- matrix(NA, nrow = N, ncol = p.det)
for (i in 1:p.occ) {
beta[, i] <- rnorm(N, beta.mean[i], sqrt(tau.sq.beta[i]))
}
for (i in 1:p.det) {
alpha[, i] <- rnorm(N, alpha.mean[i], sqrt(tau.sq.alpha[i]))
}
alpha.true <- alpha
n.factors <- 3
phi <- rep(3 / .7, n.factors)
sigma.sq <- rep(2, n.factors)
nu <- rep(2, n.factors)
dat <- simMsOcc(J.x = J.x, J.y = J.y, n.rep = n.rep, N = N, beta = beta, alpha = alpha,
psi.RE = psi.RE, p.RE = p.RE, sp = TRUE, sigma.sq = sigma.sq,
phi = phi, nu = nu, cov.model = 'matern', factor.model = TRUE,
n.factors = n.factors)
#> sigma.sq is specified but will be set to 1 for spatial latent factor model
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y <- dat$y[, -pred.indx, , drop = FALSE]
# Occupancy covariates
X <- dat$X[-pred.indx, , drop = FALSE]
coords <- as.matrix(dat$coords[-pred.indx, , drop = FALSE])
# Prediction covariates
X.0 <- dat$X[pred.indx, , drop = FALSE]
coords.0 <- as.matrix(dat$coords[pred.indx, , drop = FALSE])
# Detection covariates
X.p <- dat$X.p[-pred.indx, , , drop = FALSE]
y <- apply(y, c(1, 2), max, na.rm = TRUE)
data.list <- list(y = y, coords = coords)
# Priors
prior.list <- list(beta.comm.normal = list(mean = 0, var = 2.72),
tau.sq.beta.ig = list(a = 0.1, b = 0.1),
nu.unif = list(0.5, 2.5))
# Starting values
inits.list <- list(beta.comm = 0,
beta = 0,
fix = TRUE,
tau.sq.beta = 1)
# Tuning
tuning.list <- list(phi = 1, nu = 0.25)
batch.length <- 25
n.batch <- 2
n.report <- 100
formula <- ~ 1
out <- sfJSDM(formula = formula,
data = data.list,
inits = inits.list,
n.batch = n.batch,
batch.length = batch.length,
accept.rate = 0.43,
priors = prior.list,
cov.model = "matern",
tuning = tuning.list,
n.factors = 3,
n.omp.threads = 1,
verbose = TRUE,
NNGP = TRUE,
n.neighbors = 5,
search.type = 'cb',
n.report = 10,
n.burn = 0,
n.thin = 1,
n.chains = 2)
#> ----------------------------------------
#> Preparing to run the model
#> ----------------------------------------
#> covariates (covs) not specified in data.
#> Assuming intercept only model.
#> No prior specified for phi.unif.
#> Setting uniform bounds based on the range of observed spatial coordinates.
#> phi is not specified in initial values.
#> Setting initial value to random values from the prior distribution
#> lambda is not specified in initial values.
#> Setting initial values of the lower triangle to 0
#> nu is not specified in initial values.
#> Setting initial values to random values from the prior distribution
#> w is not specified in initial values.
#> Setting initial value to 0
#> Fixing initial values across all chains
#> ----------------------------------------
#> Building the neighbor list
#> ----------------------------------------
#> ----------------------------------------
#> Building the neighbors of neighbors list
#> ----------------------------------------
#> ----------------------------------------
#> Model description
#> ----------------------------------------
#> Spatial Factor NNGP JSDM with Polya-Gamma latent
#> variable fit with 48 sites and 6 species.
#>
#> Samples per chain: 50 (2 batches of length 25)
#> Burn-in: 0
#> Thinning Rate: 1
#> Number of Chains: 2
#> Total Posterior Samples: 100
#>
#> Using the matern spatial correlation model.
#>
#> Using 3 latent spatial factors.
#> Using 5 nearest neighbors.
#>
#> Source compiled with OpenMP support and model fit using 1 thread(s).
#>
#> Adaptive Metropolis with target acceptance rate: 43.0
#> ----------------------------------------
#> Chain 1
#> ----------------------------------------
#> Sampling ...
#> Batch: 2 of 2, 100.00%
#> ----------------------------------------
#> Chain 2
#> ----------------------------------------
#> Sampling ...
#> Batch: 2 of 2, 100.00%
summary(out)
#>
#> Call:
#> sfJSDM(formula = formula, data = data.list, inits = inits.list,
#> priors = prior.list, tuning = tuning.list, cov.model = "matern",
#> NNGP = TRUE, n.neighbors = 5, search.type = "cb", n.factors = 3,
#> n.batch = n.batch, batch.length = batch.length, accept.rate = 0.43,
#> n.omp.threads = 1, verbose = TRUE, n.report = 10, n.burn = 0,
#> n.thin = 1, n.chains = 2)
#>
#> Samples per Chain: 50
#> Burn-in: 0
#> Thinning Rate: 1
#> Number of Chains: 2
#> Total Posterior Samples: 100
#> Run Time (min): 0.0044
#>
#> ----------------------------------------
#> Community Level
#> ----------------------------------------
#> Means (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) 0.2322 0.2556 -0.2771 0.2453 0.6704 1.1063 30
#>
#> Variances (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) 0.3482 0.3397 0.0589 0.2063 1.2993 1.3084 43
#>
#> ----------------------------------------
#> Species Level
#> ----------------------------------------
#> Estimates (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept)-sp1 0.3622 0.3254 -0.2586 0.3328 1.0234 1.5419 23
#> (Intercept)-sp2 0.0555 0.3471 -0.6236 0.0471 0.6463 1.0451 41
#> (Intercept)-sp3 0.3972 0.3998 -0.2835 0.3856 1.1926 1.3783 22
#> (Intercept)-sp4 -0.2545 0.3477 -0.9843 -0.2749 0.3344 0.9963 32
#> (Intercept)-sp5 0.3955 0.3726 -0.2382 0.3961 1.1989 1.0224 21
#> (Intercept)-sp6 0.4831 0.3029 -0.0982 0.4817 1.0232 1.0899 40
#>
#> ----------------------------------------
#> Spatial Covariance
#> ----------------------------------------
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> phi-1 12.5684 3.1266 7.6594 11.5839 19.9881 1.2173 10
#> phi-2 11.3220 4.5255 2.9456 10.6578 18.4969 1.7669 9
#> phi-3 16.2968 2.7835 11.4508 16.1413 20.4122 1.8401 8
#> nu-1 1.5832 0.2628 1.1442 1.5479 2.0639 6.4558 2
#> nu-2 0.8560 0.1528 0.6664 0.8105 1.2290 3.5835 6
#> nu-3 1.8787 0.1380 1.6133 1.8818 2.0921 3.8551 10
# Update the initial model fit
out.new <- updateMCMC(out, n.batch = 1, keep.orig = TRUE,
verbose = TRUE, n.report = 1)
#> ----------------------------------------
#> Preparing to run the model
#> ----------------------------------------
#> covariates (covs) not specified in data.
#> Assuming intercept only model.
#> No prior specified for phi.unif.
#> Setting uniform bounds based on the range of observed spatial coordinates.
#> ----------------------------------------
#> Building the neighbor list
#> ----------------------------------------
#> ----------------------------------------
#> Building the neighbors of neighbors list
#> ----------------------------------------
#> ----------------------------------------
#> Model description
#> ----------------------------------------
#> Spatial Factor NNGP JSDM with Polya-Gamma latent
#> variable fit with 48 sites and 6 species.
#>
#> Samples per chain: 25 (1 batches of length 25)
#> Burn-in: 0
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 25
#>
#> Using the matern spatial correlation model.
#>
#> Using 3 latent spatial factors.
#> Using 5 nearest neighbors.
#>
#> Source compiled with OpenMP support and model fit using 1 thread(s).
#>
#> Adaptive Metropolis with target acceptance rate: 43.0
#> ----------------------------------------
#> Chain 1
#> ----------------------------------------
#> Sampling ...
#> Batch: 1 of 1, 100.00%
#> ----------------------------------------
#> Preparing to run the model
#> ----------------------------------------
#> covariates (covs) not specified in data.
#> Assuming intercept only model.
#> No prior specified for phi.unif.
#> Setting uniform bounds based on the range of observed spatial coordinates.
#> ----------------------------------------
#> Building the neighbor list
#> ----------------------------------------
#> ----------------------------------------
#> Building the neighbors of neighbors list
#> ----------------------------------------
#> ----------------------------------------
#> Model description
#> ----------------------------------------
#> Spatial Factor NNGP JSDM with Polya-Gamma latent
#> variable fit with 48 sites and 6 species.
#>
#> Samples per chain: 25 (1 batches of length 25)
#> Burn-in: 0
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 25
#>
#> Using the matern spatial correlation model.
#>
#> Using 3 latent spatial factors.
#> Using 5 nearest neighbors.
#>
#> Source compiled with OpenMP support and model fit using 1 thread(s).
#>
#> Adaptive Metropolis with target acceptance rate: 43.0
#> ----------------------------------------
#> Chain 1
#> ----------------------------------------
#> Sampling ...
#> Batch: 1 of 1, 100.00%
summary(out.new)
#>
#> Call:
#> sfJSDM(formula = formula, data = data.list, inits = inits.list,
#> priors = prior.list, tuning = tuning.list, cov.model = "matern",
#> NNGP = TRUE, n.neighbors = 5, search.type = "cb", n.factors = 3,
#> n.batch = n.batch, batch.length = batch.length, accept.rate = 0.43,
#> n.omp.threads = 1, verbose = TRUE, n.report = 10, n.burn = 0,
#> n.thin = 1, n.chains = 2)
#>
#> Samples per Chain: 75
#> Burn-in: 0
#> Thinning Rate: 1
#> Number of Chains: 2
#> Total Posterior Samples: 150
#> Run Time (min): 0.0077
#>
#> ----------------------------------------
#> Community Level
#> ----------------------------------------
#> Means (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) 0.2401 0.262 -0.29 0.2476 0.6886 1.3649 54
#>
#> Variances (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept) 0.3506 0.3568 0.0594 0.2343 1.2924 1.6777 50
#>
#> ----------------------------------------
#> Species Level
#> ----------------------------------------
#> Estimates (logit scale):
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> (Intercept)-sp1 0.3221 0.3554 -0.3099 0.3159 1.0233 2.3092 32
#> (Intercept)-sp2 0.1151 0.3677 -0.5017 0.0942 0.7948 1.1260 49
#> (Intercept)-sp3 0.4239 0.3744 -0.2423 0.4364 1.1843 1.6156 38
#> (Intercept)-sp4 -0.2585 0.3602 -1.0239 -0.2551 0.3480 1.0451 72
#> (Intercept)-sp5 0.4280 0.3457 -0.1491 0.4208 1.1378 1.0361 38
#> (Intercept)-sp6 0.4532 0.3227 -0.1227 0.4602 1.0081 1.4523 56
#>
#> ----------------------------------------
#> Spatial Covariance
#> ----------------------------------------
#> Mean SD 2.5% 50% 97.5% Rhat ESS
#> phi-1 11.4096 3.2402 6.3683 11.0621 18.6425 1.3188 16
#> phi-2 10.3523 5.1303 2.6140 10.3988 18.4381 3.2532 8
#> phi-3 16.2921 2.5379 11.5285 16.2769 20.2804 1.4989 15
#> nu-1 1.6034 0.2401 1.2044 1.5479 2.0638 7.1522 2
#> nu-2 0.8075 0.1516 0.6171 0.7851 1.2288 2.1567 5
#> nu-3 1.9362 0.1859 1.6263 1.9308 2.2593 4.0330 8