From c873061de3cf3613be848e1cb33fc6904a8a11be Mon Sep 17 00:00:00 2001 From: doserjef Date: Sat, 14 Oct 2023 09:47:42 -0400 Subject: [PATCH] lots of monotonous edits to documentation --- man/DS.Rd | 18 +- man/abund.Rd | 67 +- man/dataNMixSim.rda.Rd | 30 +- man/fitted.DS.Rd | 2 +- man/fitted.abund.Rd | 2 +- man/fitted.lfMsAbund.Rd | 2 +- man/fitted.lfMsDS.Rd | 2 +- man/fitted.msAbund.Rd | 2 +- man/fitted.msDS.Rd | 2 +- man/fitted.sfMsAbund.Rd | 2 +- man/fitted.sfMsDS.Rd | 2 +- man/fitted.spAbund.Rd | 2 +- man/fitted.spDS.Rd | 2 +- man/fitted.svcAbund.Rd | 2 +- man/fitted.svcMsAbund.Rd | 2 +- man/hbefCount2015.rda.Rd | 2 +- man/lfMsAbund.Rd | 85 +- man/lfMsDS.Rd | 62 +- man/msAbund.Rd | 68 +- man/msDS.Rd | 69 +- man/neonDWP.rda.Rd | 15 +- man/ppcAbund.Rd | 50 +- man/predict.DS.Rd | 2 +- man/predict.NMix.Rd | 28 +- man/predict.abund.Rd | 39 +- man/predict.lfMsAbund.Rd | 44 +- man/predict.lfMsDS.Rd | 54 +- man/predict.lfMsNMix.Rd | 34 +- man/predict.msAbund.Rd | 42 +- man/predict.msDS.Rd | 45 +- man/predict.msNMix.Rd | 29 +- man/predict.sfMsAbund.Rd | 55 +- man/predict.sfMsDS.Rd | 66 +- man/predict.sfMsNMix.Rd | 44 +- man/predict.spAbund.Rd | 45 +- man/predict.spDS.Rd | 14 +- man/predict.spNMix.Rd | 61 +- man/predict.svcAbund.Rd | 24 +- man/predict.svcMsAbund.Rd | 85 +- man/sfMsAbund.Rd | 87 +- man/sfMsDS.Rd | 75 +- man/simAbund.Rd | 18 +- man/simDS.Rd | 16 +- man/simMsAbund.Rd | 17 +- man/simMsDS.Rd | 12 +- man/simMsNMix.Rd | 4 +- man/simNMix.Rd | 16 +- man/spAbund.Rd | 92 +- man/spDS.Rd | 12 +- man/summary.abund.Rd | 2 +- man/summary.lfMsAbund.Rd | 2 +- man/summary.msAbund.Rd | 2 +- man/summary.sfMsAbund.Rd | 2 +- man/summary.spAbund.Rd | 2 +- man/summary.spDS.Rd | 2 +- man/summary.svcAbund.Rd | 2 +- man/summary.svcMsAbund.Rd | 2 +- man/svcAbund.Rd | 76 +- man/svcMsAbund.Rd | 217 ++-- man/waicAbund.Rd | 56 +- vignettes/distanceSampling.html | 11 +- vignettes/nMixtureModels.Rmd | 2 +- vignettes/nMixtureModels.html | 1799 ++++++++++++++++--------------- 63 files changed, 1852 insertions(+), 1875 deletions(-) diff --git a/man/DS.Rd b/man/DS.Rd index e70f560..b65bcc6 100644 --- a/man/DS.Rd +++ b/man/DS.Rd @@ -4,11 +4,11 @@ \usage{ DS(abund.formula, det.formula, data, inits, priors, tuning, - n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', - transect = 'line', det.func = 'halfnormal', - n.omp.threads = 1, verbose = TRUE, - n.report = 100, n.burn = round(.10 * n.batch * batch.length), n.thin = 1, - n.chains = 1, ...) + n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', + transect = 'line', det.func = 'halfnormal', + n.omp.threads = 1, verbose = TRUE, + n.report = 100, n.burn = round(.10 * n.batch * batch.length), n.thin = 1, + n.chains = 1, ...) } \description{ @@ -73,7 +73,7 @@ DS(abund.formula, det.formula, data, inits, priors, tuning, Gamma distribution. The hyperparameters of the inverse-Gamma distribution are passed as a list of length two with first and second elements corresponding to the shape and scale parameters, respectively, which are each specified as - vectors of length equal to the number of random intercepts or of length one + vectors of length equal to the number of random intercepts/slopes or of length one if priors are the same for all random effect variances.} \item{tuning}{a single numeric value representing the initial variance of the @@ -81,10 +81,10 @@ DS(abund.formula, det.formula, data, inits, priors, tuning, random effect values), \code{alpha.star} (the detection random effect values), and \code{kappa}. See Roberts and Rosenthal (2009) for details.} -\item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC +\item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} -\item{batch.length}{the length of each MCMC batch in each chain to run for the Adaptive +\item{batch.length}{the number of MCMC samples in each batch in each chain to run for the Adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} \item{accept.rate}{target acceptance rate for Adaptive MCMC. Default is @@ -238,7 +238,7 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, covs = covs, diff --git a/man/abund.Rd b/man/abund.Rd index 73632ef..008f25a 100644 --- a/man/abund.Rd +++ b/man/abund.Rd @@ -1,17 +1,17 @@ \name{abund} \alias{abund} -\title{Function for Fitting Single-Species Abundance GLMs} +\title{Function for Fitting Univariate Abundance GLMMs} \usage{ abund(formula, data, inits, priors, tuning, - n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', - n.omp.threads = 1, verbose = TRUE, - n.report = 100, n.burn = round(.10 * n.batch * batch.length), n.thin = 1, - n.chains = 1, save.fitted = TRUE, ...) + n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', + n.omp.threads = 1, verbose = TRUE, + n.report = 100, n.burn = round(.10 * n.batch * batch.length), n.thin = 1, + n.chains = 1, save.fitted = TRUE, ...) } \description{ - Function for fitting single-species abundance generalized linear models + Function for fitting univariate abundance generalized linear (mixed) models } \arguments{ @@ -21,7 +21,7 @@ abund(formula, data, inits, priors, tuning, and slopes are allowed using lme4 syntax (Bates et al. 2015).} \item{data}{a list containing data necessary for model fitting. - Valid tags are \code{y}, \code{covs}, and \code{offset}. \code{y} + Valid tags are \code{y}, \code{covs}, \code{z}, and \code{offset}. \code{y} is a vector, matrix, or data frame of the observed count values. If a vector, the values represent the observed counts at each site. If multiple replicate observations are obtained at the sites (e.g., sub-samples, repeated sampling over @@ -38,7 +38,7 @@ abund(formula, data, inits, priors, tuning, given site. For zero-inflated Gaussian models, the tag \code{z} is used to specify the binary component of the zero-inflated model and should have the same length as \code{y}. \code{offset} is an offset to use in the abundance model (e.g., an area offset). - This can be either a single value, a vector with an offset for each site (e.g., if survey area differed in size), or a six x replicate matrix if more than one counts are available at a given site.} + This can be either a single value, a vector with an offset for each site (e.g., if survey area differed in size), or a site x replicate matrix if more than one count is available at a given site.} \item{inits}{a list with each tag corresponding to a parameter name. Valid tags are \code{beta}, \code{kappa}, \code{sigma.sq.mu}, and \code{tau.sq}. @@ -71,20 +71,20 @@ abund(formula, data, inits, priors, tuning, Gamma distribution. The hyperparameters of the inverse-Gamma distribution are passed as a list of length two with first and second elements corresponding to the shape and scale parameters, respectively, which are each specified as - vectors of length equal to the number of random intercepts or of length one + vectors of length equal to the number of random effects or of length one if priors are the same for all random effect variances. \code{tau.sq} is the residual variance for Gaussian (or zero-inflated Gaussian) models, and it is assigned an inverse-Gamma prior. The hyperparameters of the inverse-Gamma are passed as a vector of length two, with the first and second element corresponding to the shape and scale parameters, respectively.} -\item{tuning}{a list with each tag corresponding to a parameter name, whose +\item{tuning}{a list with each tag corresponding to a parameter name, whose value defines the initial variance of the adaptive sampler. Valid tags are \code{beta}, \code{beta.star} (the abundance random effect values), and \code{kappa}. See Roberts and Rosenthal (2009) for details. Note that no tuning is necessary for Gaussian or zero-inflated Gaussian models.} -\item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC +\item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} \item{batch.length}{the length of each MCMC batch in each chain to run for the Adaptive @@ -204,27 +204,25 @@ kappa <- 0.5 sp <- FALSE family <- 'NB' dat <- simAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, - kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') + kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') y <- dat$y X <- dat$X X.re <- dat$X.re covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.cov.2 = X[, , 3], - abund.cov.3 = X[, , 4], - abund.factor.1 = X.re[, , 1]) + abund.cov.1 = X[, , 2], + abund.cov.2 = X[, , 3], + abund.cov.3 = X[, , 4], + abund.factor.1 = X.re[, , 1]) -data.list <- list(y = y, - covs = covs) +data.list <- list(y = y, covs = covs) # Priors prior.list <- list(beta.normal = list(mean = 0, var = 100), kappa.unif = c(0.001, 10)) # Starting values -inits.list <- list(beta = 0, - kappa = kappa) +inits.list <- list(beta = 0, kappa = kappa) tuning <- list(kappa = 0.2, beta = 0.1, beta.star = 0.2) n.batch <- 5 @@ -233,19 +231,20 @@ n.burn <- 0 n.thin <- 1 n.chains <- 1 -out <- abund(formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3 + (1 | abund.factor.1), - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - tuning = tuning, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) +out <- abund(formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3 + + (1 | abund.factor.1), + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + tuning = tuning, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) } diff --git a/man/dataNMixSim.rda.Rd b/man/dataNMixSim.rda.Rd index a7fe4d0..6ddad28 100644 --- a/man/dataNMixSim.rda.Rd +++ b/man/dataNMixSim.rda.Rd @@ -9,7 +9,7 @@ \description{ A simulated data set of repeated count data for 6 species across 225 sites with a maximum -of 5 replicate surveys performed at a given site. +of 3 replicate surveys performed at a given site. } \usage{ @@ -20,19 +20,19 @@ data(dataNMixSim) \code{dataNMixSim} is a list with four elements: \code{y}: a three-dimensional array of count data with - dimensions of species (6), sites (225) and replicates (5). + dimensions of species (6), sites (225) and replicates (3). - \code{abund.covs}: a numeric matrix with 373 rows and two columns consisting + \code{abund.covs}: a numeric matrix with 225 rows and two columns consisting of a continuous covariate and a categorical variable which may both influence abundance of the different species. - \code{det.covs}: a list of two numeric matrices with 225 rows and 5 columns. + \code{det.covs}: a list of two numeric matrices with 225 rows and 3 columns. Both matrices contain a continuous covariate that may affect detection probability of the species \code{coords}: a numeric matrix with 225 rows and two columns containing the - site coordinates (Easting and Northing). Note the data are generated across - a a unit square (i.e., the x and y coordinates are both between 0 and 1). + site coordinates (X and Y). Note the data are generated across + a unit square (i.e., the x and y coordinates are both between 0 and 1). } \keyword{datasets} @@ -62,8 +62,8 @@ data(dataNMixSim) # Random effects mu.RE <- list() mu.RE <- list(levels = c(10), - sigma.sq.mu = c(0.5), - beta.indx = list(1)) + sigma.sq.mu = c(0.5), + beta.indx = list(1)) p.RE <- list() # Draw species-level effects from community means. beta <- matrix(NA, nrow = n.sp, ncol = p.abund) @@ -83,11 +83,9 @@ data(dataNMixSim) family <- 'Poisson' dat <- simMsNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, alpha = alpha, - mu.RE = mu.RE, p.RE = p.RE, sp = sp, kappa = kappa, family = family, - factor.model = factor.model, phi = phi, - cov.model = 'exponential', n.factors = n.factors) - table(dat$N) - apply(dat$N, 1, sum) + mu.RE = mu.RE, p.RE = p.RE, sp = sp, kappa = kappa, family = family, + factor.model = factor.model, phi = phi, + cov.model = 'exponential', n.factors = n.factors) y <- dat$y X <- dat$X @@ -102,10 +100,10 @@ data(dataNMixSim) colnames(abund.covs) <- c('int', 'abund.cov.1', 'abund.factor.1') abund.covs <- abund.covs[, -1] det.covs <- list(det.cov.1 = X.p[, , 2], - det.cov.2 = X.p[, , 3]) + det.cov.2 = X.p[, , 3]) dataNMixSim <- list(y = y, - abund.covs = abund.covs, - det.covs = det.covs, + abund.covs = abund.covs, + det.covs = det.covs, coords = coords) } } diff --git a/man/fitted.DS.Rd b/man/fitted.DS.Rd index d6bc3eb..02bb41f 100644 --- a/man/fitted.DS.Rd +++ b/man/fitted.DS.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for DS Object} \description{ - Method for extracting model fitted values and cell-specific detection probabilities from a distance sampling (\code{DS}) model. + Method for extracting model fitted values and cell-specific detection probabilities from a hierarchical distance sampling (\code{DS}) model. } \usage{ diff --git a/man/fitted.abund.Rd b/man/fitted.abund.Rd index a9e54c4..97aaf05 100644 --- a/man/fitted.abund.Rd +++ b/man/fitted.abund.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for abund Object} \description{ - Method for extracting model fitted values from a fitted generalized linear abundance (\code{abund}) model. + Method for extracting model fitted values from a fitted GLMM (\code{abund}). } \usage{ diff --git a/man/fitted.lfMsAbund.Rd b/man/fitted.lfMsAbund.Rd index fe41cac..62e97b5 100644 --- a/man/fitted.lfMsAbund.Rd +++ b/man/fitted.lfMsAbund.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for lfMsAbund Object} \description{ - Method for extracting model fitted values from a fitted latent factor multi-species abundance GLM (\code{lfMsAbund}). + Method for extracting model fitted values from a fitted latent factor multivariate GLMM (\code{lfMsAbund}). } \usage{ diff --git a/man/fitted.lfMsDS.Rd b/man/fitted.lfMsDS.Rd index 2098dba..df3542e 100644 --- a/man/fitted.lfMsDS.Rd +++ b/man/fitted.lfMsDS.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for lfMsDS Object} \description{ - Method for extracting model fitted values and cell-specific detection probabilities from a latent factor multi-species distance sampling (\code{lfMsDS}) model. + Method for extracting model fitted values and cell-specific detection probabilities from a latent factor multi-species hierarchical distance sampling (\code{lfMsDS}) model. } \usage{ diff --git a/man/fitted.msAbund.Rd b/man/fitted.msAbund.Rd index 59dc7a6..88e707f 100644 --- a/man/fitted.msAbund.Rd +++ b/man/fitted.msAbund.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for msAbund Object} \description{ - Method for extracting model fitted values from a fitted multi-species abundance GLM (\code{msAbund}). + Method for extracting model fitted values from a fitted multivariate GLMM (\code{msAbund}). } \usage{ diff --git a/man/fitted.msDS.Rd b/man/fitted.msDS.Rd index ddac640..4e28981 100644 --- a/man/fitted.msDS.Rd +++ b/man/fitted.msDS.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for msDS Object} \description{ - Method for extracting model fitted values and cell-specific detection probabilities from a multi-species distance sampling (\code{msDS}) model. + Method for extracting model fitted values and cell-specific detection probabilities from a multi-species hierarchical distance sampling (\code{msDS}) model. } \usage{ diff --git a/man/fitted.sfMsAbund.Rd b/man/fitted.sfMsAbund.Rd index bb010a6..4a88a9a 100644 --- a/man/fitted.sfMsAbund.Rd +++ b/man/fitted.sfMsAbund.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for sfMsAbund Object} \description{ - Method for extracting model fitted values from a fitted spatial factor multi-species abundance GLM (\code{sfMsAbund}). + Method for extracting model fitted values from a fitted spatial factor multivariate GLMM (\code{sfMsAbund}). } \usage{ diff --git a/man/fitted.sfMsDS.Rd b/man/fitted.sfMsDS.Rd index 4649750..dc31e09 100644 --- a/man/fitted.sfMsDS.Rd +++ b/man/fitted.sfMsDS.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for sfMsDS Object} \description{ - Method for extracting model fitted values and cell-specific detection probabilities from a spatial factor multi-species distance sampling (\code{sfMsDS}) model. + Method for extracting model fitted values and cell-specific detection probabilities from a spatial factor multi-species hierarchical distance sampling (\code{sfMsDS}) model. } \usage{ diff --git a/man/fitted.spAbund.Rd b/man/fitted.spAbund.Rd index e8632d2..6e7343e 100644 --- a/man/fitted.spAbund.Rd +++ b/man/fitted.spAbund.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for spAbund Object} \description{ - Method for extracting model fitted values from a fitted spatial generalized linear abundance (\code{spAbund}) model. + Method for extracting model fitted values from a fitted spatial GLMM (\code{spAbund}). } \usage{ diff --git a/man/fitted.spDS.Rd b/man/fitted.spDS.Rd index e28ad26..7498ce9 100644 --- a/man/fitted.spDS.Rd +++ b/man/fitted.spDS.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for spDS Object} \description{ - Method for extracting model fitted values and cell-specific detection probabilities from a spatial distance sampling (\code{spDS}) model. + Method for extracting model fitted values and cell-specific detection probabilities from a spatial hierarchical distance sampling (\code{spDS}) model. } \usage{ diff --git a/man/fitted.svcAbund.Rd b/man/fitted.svcAbund.Rd index 323a483..e4ab531 100644 --- a/man/fitted.svcAbund.Rd +++ b/man/fitted.svcAbund.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for svcAbund Object} \description{ - Method for extracting model fitted values from a fitted spatially-varying coefficient generalized linear abundance (\code{svcAbund}) model. + Method for extracting model fitted values from a fitted spatially-varying coefficient GLMM (\code{svcAbund}). } \usage{ diff --git a/man/fitted.svcMsAbund.Rd b/man/fitted.svcMsAbund.Rd index ff1fa57..b064328 100644 --- a/man/fitted.svcMsAbund.Rd +++ b/man/fitted.svcMsAbund.Rd @@ -5,7 +5,7 @@ \title{Extract Model Fitted Values for svcMsAbund Object} \description{ - Method for extracting model fitted values from a fitted spatial factor multi-species abundance GLM (\code{svcMsAbund}). + Method for extracting model fitted values from a fitted multivatiate spatially-varying coefficient GLMM (\code{svcMsAbund}). } \usage{ diff --git a/man/hbefCount2015.rda.Rd b/man/hbefCount2015.rda.Rd index 37ba391..baa548d 100644 --- a/man/hbefCount2015.rda.Rd +++ b/man/hbefCount2015.rda.Rd @@ -40,7 +40,7 @@ data(hbefCount2015) \code{y}: a three-dimensional array of count data with dimensions of species (12), sites (373) and replicates (3). - \code{abund.covs}: a numeric matrix with 373 rows and one column consisting of the + \code{abund.covs}: a data frame with 373 rows and one column consisting of the elevation at each site. \code{det.covs}: a list of two numeric matrices with 373 rows and 3 columns. diff --git a/man/lfMsAbund.Rd b/man/lfMsAbund.Rd index 988800b..7e9ed11 100644 --- a/man/lfMsAbund.Rd +++ b/man/lfMsAbund.Rd @@ -1,9 +1,9 @@ \name{lfMsAbund} \alias{lfMsAbund} -\title{Function for Fitting Latent Factor Multi-Species Abundance GLMs} +\title{Function for Fitting Latent Factor Multivariate Abundance GLMMs} \description{ - Function for fitting multi-species generalized linear (mixed) models with species + Function for fitting multivariate generalized linear (mixed) models with species correlations (i.e., an abundance-based joint species distribution model). We use a factor modeling approach for dimension reduction. } @@ -24,14 +24,16 @@ lfMsAbund(formula, data, inits, priors, tuning, n.factors, and slopes are allowed using lme4 syntax (Bates et al. 2015).} \item{data}{a list containing data necessary for model fitting. - Valid tags are \code{y}, \code{covs}, \code{coords}, and \code{offset}. + Valid tags are \code{y}, \code{covs}, \code{coords}, \code{z}, and \code{offset}. \code{y} is a two or three-dimensional array of observed count data. The first dimension of the array is equal to the number of species and the second dimension is equal to the number of sites. If specified as a three-dimensional array, the third dimension corresponds to replicate observations at each site (e.g., sub-samples, repeated sampling over multiple seasons). \code{covs} is a list or data frame - containing the variables used in the model. Each list element is a different + containing the variables used in the model. If a data frame, each row + of \code{covs} is a site and each column is a variable. + If specified as a list, each list element is a different covariate, which can be site-level or observation-level. Site-level covariates are specified as a vector of length \eqn{J}{J}, while observation-level covariates are specified as a matrix or data frame with the number of rows equal to \eqn{J}{J} @@ -39,8 +41,8 @@ lfMsAbund(formula, data, inits, priors, tuning, n.factors, given site. \code{coords} is a matrix or data frame with two columns that contain the spatial coordinates of each site. Note that \code{spAbundance} assumes coordinates are specified - in a projected coordinate system. For Gaussian hurdle models, the tag \code{z} is used to specify the - binary component of the hurdle model and should have the same dimensions as \code{y}. \code{offset} + in a projected coordinate system. For zi-Gaussian models, the tag \code{z} is used to specify the + binary component of the zi-Gaussian model and should have the same dimensions as \code{y}. \code{offset} is an offset to use in the abundance model (e.g., an area offset). This can be either a single value, a vector with an offset for each site (e.g., if survey area differed in size), or a site x replicate matrix if more than one count is available @@ -50,7 +52,7 @@ lfMsAbund(formula, data, inits, priors, tuning, n.factors, Valid tags are \code{beta.comm}, \code{beta}, \code{tau.sq.beta}, \code{sigma.sq.mu}, \code{kappa}, \code{lambda}, \code{w}, \code{tau.sq}. \code{kappa} is only specified if \code{family = 'NB'}, \code{tau.sq} - is only speified for Gaussian and Gaussian-hurdle models, + is only specified for Gaussian and zi-Gaussian models, and \code{sigma.sq.mu} is only specified if random effects are included in \code{formula}. The value portion of each tag is the parameter's initial value. See \code{priors} description for definition @@ -89,7 +91,7 @@ lfMsAbund(formula, data, inits, priors, tuning, n.factors, lower and upper bounds of the uniform distribution, respectively, which are each specified as vectors of length equal to the number of species or of length one if priors are the same for all species-specific dispersion parameters. \code{tau.sq} is the - species-specific residual variance for Gaussian (or Gaussian hurdle) models, and it is assigned + species-specific residual variance for Gaussian (or zi-Gaussian) models, and it is assigned an inverse-Gamma prior. The hyperparameters of the inverse-Gamma are passed as a list of length two, with the first and second element corresponding to the shape and scale parameters, respectively, which are each specified as vectors of length @@ -100,7 +102,7 @@ lfMsAbund(formula, data, inits, priors, tuning, n.factors, Valid tags are \code{beta}, \code{beta.star} (the abundance random effect values), \code{kappa}, \code{lambda} (the latent factor loadings), and \code{w} (the latent factors). See Roberts and Rosenthal (2009) for details. - Note that no tuning is necessary for Gaussian or Gaussian hurdle models.} + Note that no tuning is necessary for Gaussian or zi-Gaussian models.} \item{n.factors}{the number of factors to use in the latent factor model approach. Typically, the number of factors is set to be small (e.g., 4-5) relative to the @@ -108,18 +110,18 @@ lfMsAbund(formula, data, inits, priors, tuning, n.factors, decreases in computation time. However, the value can be anywhere between 1 and N (the number of species in the community).} - \item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC + \item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{batch.length}{the length of each MCMC batch to run for the Adaptive + \item{batch.length}{the length of each MCMC batch to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{accept.rate}{target acceptance rate for Adaptive MCMC. Defaul is + \item{accept.rate}{target acceptance rate for adaptive MCMC. Defaul is 0.43. See Roberts and Rosenthal (2009) for details.} \item{family}{the distribution to use for the latent abundance process. Currently supports \code{'NB'} (negative binomial), \code{'Poisson'}, \code{'Gaussian'}, - and \code{'Gaussian-hurdle'}.} + and \code{'zi-Gaussian'}.} \item{n.omp.threads}{a positive integer indicating the number of threads to use for SMP parallel processing. The package must @@ -189,13 +191,14 @@ lfMsAbund(formula, data, inits, priors, tuning, n.factors, \item{tau.sq.samples}{a \code{coda} object of posterior samples for the Gaussian residual variance parameter. Only included when - \code{family = 'Gaussian'} or \code{family = 'Gaussian-hurdle'}.} + \code{family = 'Gaussian'} or \code{family = 'zi-Gaussian'}.} \item{lambda.samples}{a \code{coda} object of posterior samples for the latent factor loadings.} \item{w.samples}{a three-dimensional array of posterior samples for - the latent effects for each latent factor.} + the latent effects for each latent factor. Array dimensions correspond + to MCMC sample, latent factor, then site.} \item{y.rep.samples}{a three or four-dimensional array of posterior samples for the fitted (replicate) values for each species with dimensions corresponding @@ -242,7 +245,7 @@ p.abund <- length(beta.mean) tau.sq.beta <- c(0.2, 1.2) # Random effects (two random intercepts) mu.RE <- list(levels = c(10, 15), - sigma.sq.mu = c(0.43, 0.5)) + sigma.sq.mu = c(0.43, 0.5)) # Draw species-level effects from community means. beta <- matrix(NA, nrow = n.sp, ncol = p.abund) for (i in 1:p.abund) { @@ -253,8 +256,9 @@ kappa <- runif(n.sp, 0.1, 1) factor.model <- TRUE n.factors <- 3 -dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB') +dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, + n.sp = n.sp, beta = beta, mu.RE = mu.RE, + sp = sp, kappa = kappa, family = 'NB') y <- dat$y X <- dat$X @@ -263,19 +267,15 @@ coords <- dat$coords # Package all data into a list covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.factor.1 = X.re[, , 1], - abund.factor.2 = X.re[, , 2]) -data.list <- list(y = y, - covs = covs, - coords = coords) + abund.cov.1 = X[, , 2], + abund.factor.1 = X.re[, , 1], + abund.factor.2 = X.re[, , 2]) +data.list <- list(y = y, covs = covs, coords = coords) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), kappa.unif = list(a = 0, b = 10), tau.sq.beta.ig = list(a = .1, b = .1)) -inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - tau.sq.beta = 1) +inits.list <- list(beta.comm = 0, beta = 0, kappa = 0.5, + tau.sq.beta = 1) tuning.list <- list(kappa = 0.3, beta = 0.1, beta.star = 0.1, lambda = 0.5, w = 0.5) @@ -286,19 +286,20 @@ n.burn <- 20 n.thin <- 1 n.chains <- 1 -out <- lfMsAbund(formula = ~ abund.cov.1 + (1 | abund.factor.1) + (1 | abund.factor.2), - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - tuning = tuning.list, - batch.length = batch.length, - n.factors = n.factors, - n.omp.threads = 3, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) +out <- lfMsAbund(formula = ~ abund.cov.1 + (1 | abund.factor.1) + + (1 | abund.factor.2), + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + tuning = tuning.list, + batch.length = batch.length, + n.factors = n.factors, + n.omp.threads = 3, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) } diff --git a/man/lfMsDS.Rd b/man/lfMsDS.Rd index 320a645..bbe5ad8 100644 --- a/man/lfMsDS.Rd +++ b/man/lfMsDS.Rd @@ -34,10 +34,10 @@ lfMsDS(abund.formula, det.formula, data, inits, priors, of species, second dimension equal to the number of sites, and third dimension equal to the maximum number of replicates at a given site. \code{covs} is a matrix or data frame containing the variables - used in the the abundance and/or the detection portion of the model, with + used in the abundance and/or the detection portion of the model, with \eqn{J}{J} rows for each column (variable). \code{dist.breaks} is a vector of distances that denote the breakpoints of the distance bands. \code{dist.breaks} should - have length equal to the number of columns in \code{y} plus one. \code{offset} is an + have length equal to the third dimension of \code{y} plus one. \code{offset} is an offset that can be used to scale estimates from abundance per transect to density per some desired unit of measure. This can be either a single value or a vector with an offset value for each site (e.g., if transects differ in length). \code{coords} is a matrix @@ -93,7 +93,7 @@ lfMsDS(abund.formula, det.formula, data, inits, priors, if priors are the same for all species-specific dispersion parameters.} \item{tuning}{a list with each tag corresponding to a parameter name, whose - whose value defines the initial variance of the adaptive sampler. + value defines the initial variance of the adaptive sampler. Valid tags are \code{beta}, \code{alpha}, \code{lambda} (the latent factor loadings), \code{w} (the latent factors), \code{beta.star} (the abundance random effect values), \code{alpha.star} (the detection random effect values), and @@ -204,12 +204,13 @@ lfMsDS(abund.formula, det.formula, data, inits, priors, \item{N.samples}{a three-dimensional array of posterior samples for the latent abundance values for each species. Note that these values always represent transect-level abundance, even when an offset is - supplied.} + supplied. Array dimensions correspond to MCMC sample, species, and site.} \item{mu.samples}{a three-dimensional array of posterior samples for the latent expected abundance values for each species. When an offset is supplied in the \code{data} object, these correspond to expected - abundance per unit area (i.e., density).} + abundance per unit area (i.e., density). Array dimensions correspond to + MCMC sample, species, and site.} \item{sigma.sq.mu.samples}{a \code{coda} object of posterior samples for variances of random effects included in the abundance portion @@ -288,7 +289,7 @@ factor.model <- TRUE n.factors <- 3 dat <- simMsDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect, factor.model = factor.model, @@ -302,24 +303,23 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, - covs = covs, - dist.breaks = dist.breaks, - coords = coords, + covs = covs, + dist.breaks = dist.breaks, + coords = coords, offset = offset) # Priors prior.list <- list(beta.comm.normal = list(mean = 0, var = 10), - alpha.comm.normal = list(mean = 0, - var = 10), + alpha.comm.normal = list(mean = 0, var = 10), kappa.unif = list(0, 100), tau.sq.beta.ig = list(a = 0.1, b = 0.1), tau.sq.alpha.ig = list(a = 0.1, b = 0.1)) # Starting values inits.list <- list(alpha.comm = 0, beta.comm = 0, beta = 0, - alpha = 0, kappa = 1) + alpha = 0, kappa = 1) tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1, kappa = 0.8, lambda = 1, w = 1) @@ -331,23 +331,23 @@ n.thin <- 1 n.chains <- 1 out <- lfMsDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - det.formula = ~ det.cov.1, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - family = 'Poisson', - det.func = 'halfnormal', - transect = transect, - tuning = tuning, - n.factors = n.factors, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 10, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + family = 'Poisson', + det.func = 'halfnormal', + transect = transect, + tuning = tuning, + n.factors = n.factors, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 10, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out, level = 'community') } diff --git a/man/msAbund.Rd b/man/msAbund.Rd index f4fa0f7..c3f8395 100644 --- a/man/msAbund.Rd +++ b/man/msAbund.Rd @@ -1,15 +1,15 @@ \name{msAbund} \alias{msAbund} -\title{Function for Fitting Multi-Species Abundance GLMs} +\title{Function for Fitting Multivariate Abundance GLMMs} \description{ - The function \code{msAbund} fits multi-species abundance GLMs. + The function \code{msAbund} fits multivariate abundance GLMMs. } \usage{ msAbund(formula, data, inits, priors, tuning, n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', - n.omp.threads = 1, verbose = TRUE, n.report = 100, + n.omp.threads = 1, verbose = TRUE, n.report = 100, n.burn = round(.10 * n.batch * batch.length), n.thin = 1, n.chains = 1, save.fitted = TRUE, ...) } @@ -22,15 +22,17 @@ msAbund(formula, data, inits, priors, tuning, and slopes are allowed using lme4 syntax (Bates et al. 2015).} \item{data}{a list containing data necessary for model fitting. - Valid tags are \code{y}, \code{covs}, and \code{offset}. + Valid tags are \code{y}, \code{covs}, \code{z}, and \code{offset}. \code{y} is a two or three-dimensional array of observed count data. The first dimension of the array is equal to the number of species and the second dimension is equal to the number of sites. If specified as a three-dimensional array, the third dimension corresponds to replicate observations at each site (e.g., sub-samples, repeated sampling over multiple seasons). \code{covs} is a list or data frame - containing the variables used in the model. Each list element is a different - covariate, which can be site-level or observation-level. Site-level covariates + containing the variables used in the model. If a data frame, each row + of \code{covs} is a site and each column is a variable. If a list, + each list element is a different covariate, which + can be site-level or observation-level. Site-level covariates are specified as a vector of length \eqn{J}{J}, while observation-level covariates are specified as a matrix or data frame with the number of rows equal to \eqn{J}{J} and number of columns equal to the maximum number of replicate observations at a @@ -89,19 +91,19 @@ msAbund(formula, data, inits, priors, tuning, scale parameters, respectively, which are each specified as vectors of length equal to the number of species or a single value if priors are the same for all species.} -\item{tuning}{a list with each tag corresponding to a parameter name, whose +\item{tuning}{a list with each tag corresponding to a parameter name, whose value defines the initial variance of the adaptive sampler. Valid tags are \code{beta}, \code{beta.star} (the abundance random effect values), and \code{kappa}. See Roberts and Rosenthal (2009) for details. Note that no tuning is necessary for Gaussian or zero-inflated Gaussian models.} - \item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC + \item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{batch.length}{the length of each MCMC batch to run for the Adaptive + \item{batch.length}{the length of each MCMC batch to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{accept.rate}{target acceptance rate for Adaptive MCMC. Defaul is + \item{accept.rate}{target acceptance rate for adaptive MCMC. Defaul is 0.43. See Roberts and Rosenthal (2009) for details.} \item{family}{the distribution to use for the latent abundance process. Currently @@ -223,7 +225,7 @@ p.abund <- length(beta.mean) tau.sq.beta <- c(0.2, 1.2) # Random effects (two random intercepts) mu.RE <- list(levels = c(10, 15), - sigma.sq.mu = c(0.43, 0.5)) + sigma.sq.mu = c(0.43, 0.5)) # Draw species-level effects from community means. beta <- matrix(NA, nrow = n.sp, ncol = p.abund) for (i in 1:p.abund) { @@ -233,7 +235,7 @@ sp <- FALSE kappa <- runif(n.sp, 0.1, 1) dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB') + mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB') y <- dat$y X <- dat$X @@ -241,18 +243,17 @@ X.re <- dat$X.re # Package all data into a list covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.factor.1 = X.re[, , 1], - abund.factor.2 = X.re[, , 2]) -data.list <- list(y = y, - covs = covs) + abund.cov.1 = X[, , 2], + abund.factor.1 = X.re[, , 1], + abund.factor.2 = X.re[, , 2]) +data.list <- list(y = y, covs = covs) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), kappa.unif = list(a = 0, b = 10), tau.sq.beta.ig = list(a = .1, b = .1)) inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - tau.sq.beta = 1) + beta = 0, + kappa = 0.5, + tau.sq.beta = 1) tuning.list <- list(kappa = 0.3, beta = 0.1, beta.star = 0.1) # Small @@ -262,18 +263,19 @@ n.burn <- 20 n.thin <- 1 n.chains <- 1 -out <- msAbund(formula = ~ abund.cov.1 + (1 | abund.factor.1) + (1 | abund.factor.2), - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - tuning = tuning.list, - batch.length = batch.length, - n.omp.threads = 3, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) +out <- msAbund(formula = ~ abund.cov.1 + (1 | abund.factor.1) + + (1 | abund.factor.2), + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + tuning = tuning.list, + batch.length = batch.length, + n.omp.threads = 3, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) } diff --git a/man/msDS.Rd b/man/msDS.Rd index 812c3cc..855f1a5 100644 --- a/man/msDS.Rd +++ b/man/msDS.Rd @@ -8,11 +8,11 @@ Function for fitting multi-species hierarchical distance sampling models. \usage{ msDS(abund.formula, det.formula, data, inits, priors, - tuning, n.batch, batch.length, accept.rate = 0.43, - family = 'Poisson', transect = 'line', det.func = 'halfnormal', - n.omp.threads = 1, verbose = TRUE, n.report = 100, - n.burn = round(.10 * n.batch * batch.length), n.thin = 1, - n.chains = 1, ...) + tuning, n.batch, batch.length, accept.rate = 0.43, + family = 'Poisson', transect = 'line', det.func = 'halfnormal', + n.omp.threads = 1, verbose = TRUE, n.report = 100, + n.burn = round(.10 * n.batch * batch.length), n.thin = 1, + n.chains = 1, ...) } \arguments{ @@ -33,10 +33,10 @@ msDS(abund.formula, det.formula, data, inits, priors, of species, second dimension equal to the number of sites, and third dimension equal to the maximum number of replicates at a given site. \code{covs} is a matrix or data frame containing the variables - used in the the abundance and/or the detection portion of the model, with + used in the abundance and/or the detection portion of the model, with \eqn{J}{J} rows for each column (variable). \code{dist.breaks} is a vector of distances that denote the breakpoints of the distance bands. \code{dist.breaks} should - have length equal to the number of columns in \code{y} plus one. \code{offset} is an + have length equal to the third dimension of \code{y} plus one. \code{offset} is an offset that can be used to scale estimates from abundance per transect to density per some desired unit of measure. This can be either a single value or a vector with an offset value for each site (e.g., if transects differ in length)} @@ -79,7 +79,7 @@ msDS(abund.formula, det.formula, data, inits, priors, Gamma distribution. The hyperparameters of the inverse-Gamma distribution are passed as a list of length two with first and second elements corresponding to the shape and scale parameters, respectively, which are each specified as - vectors of length equal to the number of random intercepts or of length one + vectors of length equal to the number of random effects or of length one if priors are the same for all random effect variances. \code{kappa} is the negative binomial dispersion parameter for each species and is assumed to follow a uniform distribution. The hyperparameters of the uniform distribution @@ -187,12 +187,13 @@ msDS(abund.formula, det.formula, data, inits, priors, \item{N.samples}{a three-dimensional array of posterior samples for the latent abundance values for each species. Note that these values always represent transect-level abundance, even when an offset is - supplied.} + supplied. Array dimensions correspond to MCMC sample, species, and site.} \item{mu.samples}{a three-dimensional array of posterior samples for the latent expected abundance values for each species. When an offset is supplied in the \code{data} object, these correspond to expected - abundance per unit area (i.e., density).} + abundance per unit area (i.e., density). Array dimensions correspond to + MCMC samples, species, and site.} \item{sigma.sq.mu.samples}{a \code{coda} object of posterior samples for variances of random effects included in the abundance portion @@ -270,7 +271,7 @@ transect <- 'line' factor.model <- FALSE dat <- simMsDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, mu.RE = mu.RE, p.RE = p.RE, sp = sp, cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect, factor.model = factor.model) @@ -284,23 +285,23 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, - covs = covs, - dist.breaks = dist.breaks, + covs = covs, + dist.breaks = dist.breaks, offset = offset) # Priors prior.list <- list(beta.comm.normal = list(mean = 0, var = 10), - alpha.comm.normal = list(mean = 0, - var = 10), + alpha.comm.normal = list(mean = 0, + var = 10), kappa.unif = list(0, 100), tau.sq.beta.ig = list(a = 0.1, b = 0.1), tau.sq.alpha.ig = list(a = 0.1, b = 0.1)) # Starting values inits.list <- list(alpha.comm = 0, beta.comm = 0, beta = 0, - alpha = 0, kappa = 1) + alpha = 0, kappa = 1) tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1, kappa = 0.8) @@ -312,22 +313,22 @@ n.thin <- 1 n.chains <- 1 out <- msDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - det.formula = ~ det.cov.1, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - family = 'Poisson', - det.func = 'halfnormal', - transect = transect, - tuning = tuning, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 10, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + family = 'Poisson', + det.func = 'halfnormal', + transect = transect, + tuning = tuning, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 10, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out, level = 'community') } diff --git a/man/neonDWP.rda.Rd b/man/neonDWP.rda.Rd index a575264..8eceb94 100644 --- a/man/neonDWP.rda.Rd +++ b/man/neonDWP.rda.Rd @@ -15,11 +15,12 @@ the number of all bird species observed during a six minute, unlimited radius po during the breeding season. Distance of each individual bird to the observer was recorded using a laser rangefinder. For modeling, we binned the distance measurements into 4 distance bins, and only used observations within 250m. The 16 species included in the data set are as follows: -(1) BACS (Baird's Sparrow); (2) CARW (Carolina Wren); (3) COGD (Common Ground-Dove); -(4) COYE (Common Yellowthroat); (5) EABL (Eastern Bluebird); (6) EAME (Eastern Meadowlark); -(7) EATO (Eastern Towhee); (8) GCFL (Great-crested Flycatcher); (9) MODO (Mourning Dove); -(10) NOCA (Northern Cardinal); (11) NOMO (Northern Mockingbird); -(12) RBWO (Red-bellied Woodpecker); (13) WEVI (White-eyed Vireo). +(1) EATO (Eastern Towhee); (2) EAME (Eastern Meadowlark); (3) AMCR (American Crow); +(4) BACS (Bachman's Sparrow); (5) CARW (Carolina Wren); (6) COGD (Common Ground Dove); +(7) CONI (Common Nighthawk); (8) COYE (Common Yellowthroat); (9) EABL (Eastern Bluebird); +(10) GCFL (Great-crested Flycatcher); (11) MODO (Mourning Dover); (12) NOCA (Northern Cardinal); +(13) NOMO (Northern Mockingbird); (14) RBWO (Red-bellied Woodpecker); (15) RHWO (Red-headed Woodpecker); +(16) WEVI (White-eyed Vireo). } \usage{ @@ -41,7 +42,7 @@ data(neonDWP) \code{neonDWP} is a list with five elements: \code{y}: a three-dimensional array of distance sampling data with - dimensions of species (16), sites (90) and distance bin (5). + dimensions of species (16), sites (90) and distance bin (4). \code{covs}: a data frame with 90 rows and four columns consisting of covariates for use in modeling abundance and/or detection. @@ -49,7 +50,7 @@ data(neonDWP) \code{dist.breaks}: a vector of five values indicating the break points of the four distance bands in the data. - \code{offset}: an offset used to scale the 125m radius point count surveys + \code{offset}: an offset used to scale the 250m radius point count surveys to ha, such that resulting estimates will be the number of individuals per ha. diff --git a/man/ppcAbund.Rd b/man/ppcAbund.Rd index 1a15b54..4ae6542 100644 --- a/man/ppcAbund.Rd +++ b/man/ppcAbund.Rd @@ -27,7 +27,7 @@ abundance data for the posterior predictive check. Value 1 will group values by row (site), and value 2 will group values by column (replicate).} -\item{type}{a character string indicating whether fitted values should be generated conditional on the estimated latent abundance values (\code{type = 'conditional'}) estimated during the model or based on the marginal expected abundance values (\code{type = 'marginal'}).} +\item{type}{a character string indicating whether fitted values should be generated conditional on the estimated latent abundance values (\code{type = 'conditional'}) estimated during the model or based on the marginal expected abundance values (\code{type = 'marginal'}). This is only relevant for N-mixture models.} \item{...}{currently no additional arguments} } @@ -107,7 +107,7 @@ sp <- FALSE cov.model <- 'exponential' dist <- 'NB' dat <- simNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha, - kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, + kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, phi = phi, sigma.sq = sigma.sq, cov.model = cov.model, family = 'NB') @@ -118,24 +118,20 @@ X.p <- dat$X.p abund.covs <- X colnames(abund.covs) <- c('int', 'abund.cov.1') -det.covs <- list(det.cov.1 = X.p[, , 2], - det.cov.2 = X.p[, , 3]) +det.covs <- list(det.cov.1 = X.p[, , 2], det.cov.2 = X.p[, , 3]) -data.list <- list(y = y, - abund.covs = abund.covs, - det.covs = det.covs) +data.list <- list(y = y, abund.covs = abund.covs, + det.covs = det.covs) # Priors prior.list <- list(beta.normal = list(mean = rep(0, p.abund), - var = rep(100, p.abund)), - alpha.normal = list(mean = rep(0, p.det), - var = rep(2.72, p.det)), + var = rep(100, p.abund)), + alpha.normal = list(mean = rep(0, p.det), + var = rep(2.72, p.det)), kappa.unif = c(0, 10)) # Starting values -inits.list <- list(alpha = 0, - beta = 0, - kappa = kappa, - N = apply(y, 1, max, na.rm = TRUE)) +inits.list <- list(alpha = 0, beta = 0, kappa = kappa, + N = apply(y, 1, max, na.rm = TRUE)) tuning <- 0.5 @@ -146,19 +142,19 @@ n.thin <- 1 n.chains <- 1 out <- NMix(abund.formula = ~ abund.cov.1, - det.formula = ~ det.cov.1 + det.cov.2, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1 + det.cov.2, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Posterior predictive check ppc.out <- ppcAbund(out, fit.stat = 'chi-squared', group = 0) diff --git a/man/predict.DS.Rd b/man/predict.DS.Rd index 1ae8fe2..41d2fac 100644 --- a/man/predict.DS.Rd +++ b/man/predict.DS.Rd @@ -93,7 +93,7 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, covs = covs, diff --git a/man/predict.NMix.Rd b/man/predict.NMix.Rd index adb5843..c1104df 100644 --- a/man/predict.NMix.Rd +++ b/man/predict.NMix.Rd @@ -89,7 +89,7 @@ prior.list <- list(beta.normal = list(mean = rep(0, p.abund), # Initial values inits.list <- list(alpha = rep(0, p.det), beta = rep(0, p.abund), - kappa = 0.5, + kappa = 0.5, N = apply(y, 1, max, na.rm = TRUE)) n.batch <- 10 @@ -99,19 +99,19 @@ n.thin <- 1 n.chains <- 1 out <- NMix(abund.formula = ~ abund.cov, - det.formula = ~ det.cov, - data = data.list, - inits = inits.list, - n.batch = n.batch, - batch.length = batch.length, - family = 'Poisson', - priors = prior.list, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov, + data = data.list, + inits = inits.list, + n.batch = n.batch, + batch.length = batch.length, + family = 'Poisson', + priors = prior.list, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) diff --git a/man/predict.abund.Rd b/man/predict.abund.Rd index 8f9478a..0671a89 100644 --- a/man/predict.abund.Rd +++ b/man/predict.abund.Rd @@ -1,6 +1,6 @@ \name{predict.abund} \alias{predict.abund} -\title{Function for prediction at new locations for single-species abundance GLMs} +\title{Function for prediction at new locations for univariate GLMMs} \description{ The function \code{predict} collects posterior predictive samples for a set of new locations given an object of class `abund`. @@ -37,7 +37,7 @@ \item{mu.0.samples}{a three-dimensional object of posterior predictive samples for the expected abundance values with dimensions corresponding to posterior predictive sample, site, and replicate. When there is no replication, this will be a two-dimensional - matrix. Note if an offset was used when fitting the model with \code{abund}, thes abundance + matrix. Note if an offset was used when fitting the model with \code{abund}, the abundance values are reported per unit of the offset.} \item{y.0.samples}{a three-dimensional object of posterior predictive samples for the @@ -62,7 +62,7 @@ kappa <- 0.5 sp <- FALSE family <- 'NB' dat <- simAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, - kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') + kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') # Split into fitting and prediction data set pred.indx <- sample(1:J, round(J * .25), replace = FALSE) @@ -75,9 +75,9 @@ coords <- as.matrix(dat$coords[-pred.indx, ]) coords.0 <- as.matrix(dat$coords[pred.indx, ]) abund.covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.cov.2 = X[, , 3], - abund.cov.3 = X[, , 4]) + abund.cov.1 = X[, , 2], + abund.cov.2 = X[, , 3], + abund.cov.3 = X[, , 4]) data.list <- list(y = y, covs = abund.covs) @@ -85,8 +85,7 @@ data.list <- list(y = y, covs = abund.covs) prior.list <- list(beta.normal = list(mean = 0, var = 100), kappa.unif = c(0.001, 10)) # Starting values -inits.list <- list(beta = 0, - kappa = kappa) +inits.list <- list(beta = 0, kappa = kappa) n.batch <- 5 batch.length <- 25 @@ -95,18 +94,18 @@ n.thin <- 1 n.chains <- 1 out <- abund(formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Predict at new locations ------------------------------------------------ colnames(X.0) <- c('intercept', 'abund.cov', 'abund.cov.2', 'abund.cov.3') diff --git a/man/predict.lfMsAbund.Rd b/man/predict.lfMsAbund.Rd index 0d8ce6e..0c6b7c5 100644 --- a/man/predict.lfMsAbund.Rd +++ b/man/predict.lfMsAbund.Rd @@ -1,6 +1,6 @@ \name{predict.lfMsAbund} \alias{predict.lfMsAbund} -\title{Function for prediction at new locations for latent factor multi-species abundance GLMs} +\title{Function for prediction at new locations for latent factor multivariate GLMMs} \description{ The function \code{predict} collects posterior predictive samples for a set of new locations given an object of class `lfMsAbund`. @@ -22,7 +22,7 @@ \item{ignore.RE}{logical value that specifies whether or not to remove unstructured random effects from the subsequent predictions. If \code{TRUE}, unstructured random effects will be included. If \code{FALSE}, unstructured random effects will be set to 0 and predictions will only be generated from the fixed effects.} - \item{z.0.samples}{a three-dimensional array with dimensions corresponding to MCMC samples, species, and prediction locations. The array contains the full posterior samples of the predicted binary portion of a Gaussian hurdle model. In the context of abundance models, this typically corresponds to estimates of the presence or absence of each species at the location. When using \code{spOccupancy} to generate the first stage samples of the Gaussian-hurdle model, this is the object contained in the \code{z.0.samples} object of the predition function for the spOccupancy object. Ignored for all model types other than Gaussian-hurdle.} + \item{z.0.samples}{a three-dimensional array with dimensions corresponding to MCMC samples, species, and prediction locations. The array contains the full posterior samples of the predicted binary portion of a Gaussian zero-inflated model. In the context of abundance models, this typically corresponds to estimates of the presence or absence of each species at the location. When using \code{spOccupancy} to generate the first stage samples of the zi-Gaussian model, this is the object contained in the \code{z.0.samples} object of the predition function for the spOccupancy object. Ignored for all model types other than zi-Gaussian.} \item{include.w}{a logical value used to indicate whether the latent random effects should be included in the predictions. By default, this is set to \code{TRUE}. If set to \code{FALSE}, predictions are given using the covariates and any unstructured random effects in the model. If \code{FALSE}, the \code{coords.0} argument is not required.} @@ -79,7 +79,7 @@ factor.model <- TRUE n.factors <- 3 dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', + mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', factor.model = factor.model, n.factors = n.factors) # Split into fitting and prediction data set @@ -96,17 +96,15 @@ coords.0 <- dat$coords[pred.indx, ] # Package all data into a list covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2]) -data.list <- list(y = y, - covs = covs, - coords = coords) + abund.cov.1 = X[, , 2]) +data.list <- list(y = y, covs = covs, coords = coords) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), kappa.unif = list(a = 0, b = 10), tau.sq.beta.ig = list(a = .1, b = .1)) inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - tau.sq.beta = 1) + beta = 0, + kappa = 0.5, + tau.sq.beta = 1) tuning.list <- list(kappa = 0.3, beta = 0.1, lambda = 0.5, w = 0.5) # Small @@ -117,19 +115,19 @@ n.thin <- 1 n.chains <- 1 out <- lfMsAbund(formula = ~ abund.cov.1, - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - tuning = tuning.list, - batch.length = batch.length, - n.omp.threads = 1, - n.factors = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + tuning = tuning.list, + batch.length = batch.length, + n.omp.threads = 1, + n.factors = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Predict at new locations out.pred <- predict(out, X.0, coords.0) diff --git a/man/predict.lfMsDS.Rd b/man/predict.lfMsDS.Rd index 65be58c..f29abf4 100644 --- a/man/predict.lfMsDS.Rd +++ b/man/predict.lfMsDS.Rd @@ -8,7 +8,7 @@ \usage{ \method{predict}{lfMsDS}(object, X.0, coords.0, ignore.RE = FALSE, - type = 'abundance', include.w = TRUE, ...) + type = 'abundance', include.w = TRUE, ...) } \arguments{ @@ -98,7 +98,7 @@ factor.model <- TRUE n.factors <- 3 dat <- simMsDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect, factor.model = factor.model, @@ -120,24 +120,24 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, - covs = covs, - dist.breaks = dist.breaks, - coords = coords, + covs = covs, + dist.breaks = dist.breaks, + coords = coords, offset = offset) # Priors prior.list <- list(beta.comm.normal = list(mean = 0, var = 10), - alpha.comm.normal = list(mean = 0, - var = 10), + alpha.comm.normal = list(mean = 0, + var = 10), kappa.unif = list(0, 100), tau.sq.beta.ig = list(a = 0.1, b = 0.1), tau.sq.alpha.ig = list(a = 0.1, b = 0.1)) # Starting values inits.list <- list(alpha.comm = 0, beta.comm = 0, beta = 0, - alpha = 0, kappa = 1) + alpha = 0, kappa = 1) tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1, kappa = 0.8, lambda = 1, w = 1) @@ -149,24 +149,24 @@ n.thin <- 1 n.chains <- 1 out <- lfMsDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - det.formula = ~ det.cov.1, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - family = 'Poisson', - det.func = 'halfnormal', - transect = transect, - tuning = tuning, - n.factors = n.factors, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 10, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + family = 'Poisson', + det.func = 'halfnormal', + transect = transect, + tuning = tuning, + n.factors = n.factors, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 10, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out, level = 'community') # Predict at new locations ------------------------------------------------ diff --git a/man/predict.lfMsNMix.Rd b/man/predict.lfMsNMix.Rd index b58e31f..8f5e416 100644 --- a/man/predict.lfMsNMix.Rd +++ b/man/predict.lfMsNMix.Rd @@ -8,7 +8,7 @@ \usage{ \method{predict}{lfMsNMix}(object, X.0, coords.0, ignore.RE = FALSE, - type = 'abundance', include.w = TRUE, ...) + type = 'abundance', include.w = TRUE, ...) } \arguments{ @@ -85,8 +85,8 @@ family <- 'Poisson' n.factors <- 3 dat <- simMsNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, alpha = alpha, - sp = FALSE, family = 'Poisson', factor.model = TRUE, - n.factors = n.factors) + sp = FALSE, family = 'Poisson', factor.model = TRUE, + n.factors = n.factors) # Split into fitting and prediction data set pred.indx <- sample(1:J, round(J * .25), replace = FALSE) y <- dat$y[, -pred.indx, ] @@ -108,7 +108,7 @@ det.covs <- list(det.cov.1 = X.p[, , 2], data.list <- list(y = y, abund.covs = abund.covs, det.covs = det.covs, - coords = coords) + coords = coords) # Occupancy initial values prior.list <- list(beta.comm.normal = list(mean = 0, var = 2.72), @@ -130,19 +130,19 @@ batch.length <- 25 accept.rate <- 0.43 out <- lfMsNMix(abund.formula = ~ abund.cov, - det.formula = ~ det.cov.1 + det.cov.2, - data = data.list, - inits = inits.list, - family = 'Poisson', - n.factors = n.factors, - n.batch = n.batch, - batch.length = batch.length, - accept.rate = 0.43, - tuning = tuning, - priors = prior.list, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1) + det.formula = ~ det.cov.1 + det.cov.2, + data = data.list, + inits = inits.list, + family = 'Poisson', + n.factors = n.factors, + n.batch = n.batch, + batch.length = batch.length, + accept.rate = 0.43, + tuning = tuning, + priors = prior.list, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1) summary(out, level = 'community') diff --git a/man/predict.msAbund.Rd b/man/predict.msAbund.Rd index 8d98678..a2aed5f 100644 --- a/man/predict.msAbund.Rd +++ b/man/predict.msAbund.Rd @@ -1,6 +1,6 @@ \name{predict.msAbund} \alias{predict.msAbund} -\title{Function for prediction at new locations for multi-species abundance GLMs} +\title{Function for prediction at new locations for multivariate GLMMs} \description{ The function \code{predict} collects posterior predictive samples for a set of new locations given an object of class `msAbund`. @@ -18,7 +18,7 @@ \item{ignore.RE}{logical value that specifies whether or not to remove unstructured random effects from the subsequent predictions. If \code{TRUE}, unstructured random effects will be included. If \code{FALSE}, unstructured random effects will be set to 0 and predictions will only be generated from the fixed effects.} - \item{z.0.samples}{a three-dimensional array with dimensions corresponding to MCMC samples, species, and prediction locations. The array contains the full posterior samples of the predicted binary portion of a Gaussian hurdle model. In the context of abundance models, this typically corresponds to estimates of the presence or absence of each species at the location. When using \code{spOccupancy} to generate the first stage samples of the Gaussian-hurdle model, this is the object contained in the \code{z.0.samples} object of the predition function for the spOccupancy object. Ignored for all model types other than Gaussian-hurdle.} + \item{z.0.samples}{a three-dimensional array with dimensions corresponding to MCMC samples, species, and prediction locations. The array contains the full posterior samples of the predicted binary portion of a Gaussian zero-inflated model. In the context of abundance models, this typically corresponds to estimates of the presence or absence of each species at the location. When using \code{spOccupancy} to generate the first stage samples of the zi-Gaussian model, this is the object contained in the \code{z.0.samples} object of the predition function for the spOccupancy object. Ignored for all model types other than zi-Gaussian.} \item{...}{currently no additional arguments} } @@ -68,7 +68,7 @@ sp <- FALSE kappa <- runif(n.sp, 0.1, 1) dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB') + mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB') # Split into fitting and prediction data set pred.indx <- sample(1:J, round(J * .25), replace = FALSE) @@ -80,17 +80,15 @@ y.0 <- dat$y[, pred.indx, , drop = FALSE] X.0 <- dat$X[pred.indx, , , drop = FALSE] # Package all data into a list -covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2]) -data.list <- list(y = y, - covs = covs) +covs <- list(int = X[, , 1], abund.cov.1 = X[, , 2]) +data.list <- list(y = y, covs = covs) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), kappa.unif = list(a = 0, b = 10), tau.sq.beta.ig = list(a = .1, b = .1)) inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - tau.sq.beta = 1) + beta = 0, + kappa = 0.5, + tau.sq.beta = 1) tuning.list <- list(kappa = 0.3, beta = 0.1) # Small @@ -101,18 +99,18 @@ n.thin <- 1 n.chains <- 1 out <- msAbund(formula = ~ abund.cov.1, - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - tuning = tuning.list, - batch.length = batch.length, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + tuning = tuning.list, + batch.length = batch.length, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Predict at new locations out.pred <- predict(out, X.0) diff --git a/man/predict.msDS.Rd b/man/predict.msDS.Rd index 51819c5..c72d795 100644 --- a/man/predict.msDS.Rd +++ b/man/predict.msDS.Rd @@ -7,8 +7,7 @@ } \usage{ -\method{predict}{msDS}(object, X.0, ignore.RE = FALSE, type = 'abundance', - ...) +\method{predict}{msDS}(object, X.0, ignore.RE = FALSE, type = 'abundance', ...) } \arguments{ @@ -92,7 +91,7 @@ transect <- 'line' factor.model <- FALSE dat <- simMsDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect, factor.model = factor.model) @@ -113,11 +112,11 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, - covs = covs, - dist.breaks = dist.breaks, + covs = covs, + dist.breaks = dist.breaks, offset = offset) # Priors @@ -141,23 +140,23 @@ n.thin <- 1 n.chains <- 1 out <- msDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - det.formula = ~ det.cov.1, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - family = 'Poisson', - det.func = 'halfnormal', - transect = transect, - tuning = tuning, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 10, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + family = 'Poisson', + det.func = 'halfnormal', + transect = transect, + tuning = tuning, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 10, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out, level = 'community') # Predict at new locations ------------------------------------------------ diff --git a/man/predict.msNMix.Rd b/man/predict.msNMix.Rd index d9b2b1f..129a849 100644 --- a/man/predict.msNMix.Rd +++ b/man/predict.msNMix.Rd @@ -7,8 +7,7 @@ } \usage{ -\method{predict}{msNMix}(object, X.0, ignore.RE = FALSE, type = 'abundance', - ...) +\method{predict}{msNMix}(object, X.0, ignore.RE = FALSE, type = 'abundance', ...) } \arguments{ @@ -80,7 +79,7 @@ for (i in 1:p.det) { family <- 'Poisson' dat <- simMsNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, alpha = alpha, - sp = FALSE, family = 'Poisson') + sp = FALSE, family = 'Poisson') # Split into fitting and prediction data set pred.indx <- sample(1:J, round(J * .25), replace = FALSE) y <- dat$y[, -pred.indx, ] @@ -120,18 +119,18 @@ batch.length <- 25 accept.rate <- 0.43 out <- msNMix(abund.formula = ~ abund.cov, - det.formula = ~ det.cov.1 + det.cov.2, - data = data.list, - inits = inits.list, - family = 'Poisson', - n.batch = n.batch, - batch.length = batch.length, - accept.rate = 0.43, - tuning = tuning, - priors = prior.list, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1) + det.formula = ~ det.cov.1 + det.cov.2, + data = data.list, + inits = inits.list, + family = 'Poisson', + n.batch = n.batch, + batch.length = batch.length, + accept.rate = 0.43, + tuning = tuning, + priors = prior.list, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1) summary(out, level = 'community') diff --git a/man/predict.sfMsAbund.Rd b/man/predict.sfMsAbund.Rd index cf4e080..cf81198 100644 --- a/man/predict.sfMsAbund.Rd +++ b/man/predict.sfMsAbund.Rd @@ -1,6 +1,6 @@ \name{predict.sfMsAbund} \alias{predict.sfMsAbund} -\title{Function for prediction at new locations for spatial factor multi-species abundance GLMs} +\title{Function for prediction at new locations for spatial factor multivariate GLMMs} \description{ The function \code{predict} collects posterior predictive samples for a set of new locations given an object of class `sfMsAbund`. @@ -8,8 +8,8 @@ \usage{ \method{predict}{sfMsAbund}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, ignore.RE = FALSE, - z.0.samples, include.sp = TRUE, ...) + verbose = TRUE, n.report = 100, ignore.RE = FALSE, + z.0.samples, include.sp = TRUE, ...) } \arguments{ @@ -65,7 +65,7 @@ and replicate. These will be in the same units as \code{mu.0.samples}.} \item{w.0.samples}{a three-dimensional array of posterior predictive samples for the - latent factors.} + latent factors. Array dimensions correspond to MCMC sample, latent factor, and site.} The return object will include additional objects used for standard extractor functions. @@ -96,7 +96,7 @@ cov.model <- 'spherical' phi <- runif(n.factors, 3 / 1, 3 / .1) dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', + mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', factor.model = factor.model, n.factors = n.factors, phi = phi, cov.model = cov.model) @@ -113,20 +113,17 @@ X.0 <- dat$X[pred.indx, , , drop = FALSE] coords.0 <- dat$coords[pred.indx, ] # Package all data into a list -covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2]) -data.list <- list(y = y, - covs = covs, - coords = coords) +covs <- list(int = X[, , 1], abund.cov.1 = X[, , 2]) +data.list <- list(y = y, covs = covs, coords = coords) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), kappa.unif = list(a = 0, b = 10), - phi.unif = list(a = 3 / 1, b = 3 / .1), + phi.unif = list(a = 3 / 1, b = 3 / .1), tau.sq.beta.ig = list(a = .1, b = .1)) inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - phi = 3 / .5, - tau.sq.beta = 1) + beta = 0, + kappa = 0.5, + phi = 3 / .5, + tau.sq.beta = 1) tuning.list <- list(kappa = 0.3, beta = 0.1, lambda = 0.5, w = 0.5, phi = 1) @@ -138,20 +135,20 @@ n.thin <- 1 n.chains <- 1 out <- sfMsAbund(formula = ~ abund.cov.1, - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - tuning = tuning.list, - batch.length = batch.length, - n.omp.threads = 1, - n.factors = 1, - verbose = TRUE, - n.neighbors = 5, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + tuning = tuning.list, + batch.length = batch.length, + n.omp.threads = 1, + n.factors = 1, + verbose = TRUE, + n.neighbors = 5, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Predict at new locations out.pred <- predict(out, X.0, coords.0) diff --git a/man/predict.sfMsDS.Rd b/man/predict.sfMsDS.Rd index 3689d22..c539da4 100644 --- a/man/predict.sfMsDS.Rd +++ b/man/predict.sfMsDS.Rd @@ -8,8 +8,8 @@ \usage{ \method{predict}{sfMsDS}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, ignore.RE = FALSE, - type = 'abundance', include.sp = TRUE, ...) + verbose = TRUE, n.report = 100, ignore.RE = FALSE, + type = 'abundance', include.sp = TRUE, ...) } \arguments{ @@ -48,8 +48,7 @@ } \author{ - Jeffrey W. Doser \email{doserjef@msu.edu}, \cr - Andrew O. Finley \email{finleya@msu.edu} + Jeffrey W. Doser \email{doserjef@msu.edu} \cr } \value{ @@ -117,7 +116,7 @@ phi <- runif(n.factors, 3 / 1, 3 / .2) cov.model <- 'exponential' dat <- simMsDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect, factor.model = factor.model, @@ -140,25 +139,24 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, - covs = covs, - dist.breaks = dist.breaks, - coords = coords, + covs = covs, + dist.breaks = dist.breaks, + coords = coords, offset = offset) # Priors prior.list <- list(beta.comm.normal = list(mean = 0, var = 10), - alpha.comm.normal = list(mean = 0, - var = 10), + alpha.comm.normal = list(mean = 0, var = 10), kappa.unif = list(0, 100), - phi.unif = list(3 / 1, 3 / .1), + phi.unif = list(3 / 1, 3 / .1), tau.sq.beta.ig = list(a = 0.1, b = 0.1), tau.sq.alpha.ig = list(a = 0.1, b = 0.1)) # Starting values inits.list <- list(alpha.comm = 0, beta.comm = 0, beta = 0, - alpha = 0, kappa = 1, phi = 3 / .5) + alpha = 0, kappa = 1, phi = 3 / .5) tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1, kappa = 0.8, lambda = 1, w = 1, phi = 0.8) @@ -170,27 +168,27 @@ n.thin <- 1 n.chains <- 1 out <- sfMsDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - det.formula = ~ det.cov.1, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - family = 'Poisson', - det.func = 'halfnormal', - transect = transect, - tuning = tuning, - cov.model = 'exponential', - NNGP = TRUE, - n.neighbors = 5, - n.factors = n.factors, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 10, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + family = 'Poisson', + det.func = 'halfnormal', + transect = transect, + tuning = tuning, + cov.model = 'exponential', + NNGP = TRUE, + n.neighbors = 5, + n.factors = n.factors, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 10, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out, level = 'community') # Predict at new locations ------------------------------------------------ diff --git a/man/predict.sfMsNMix.Rd b/man/predict.sfMsNMix.Rd index 9be9295..0fc91fb 100644 --- a/man/predict.sfMsNMix.Rd +++ b/man/predict.sfMsNMix.Rd @@ -8,9 +8,9 @@ \usage{ \method{predict}{sfMsNMix}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, - ignore.RE = FALSE, type = 'abundance', - include.sp = TRUE, ...) + verbose = TRUE, n.report = 100, + ignore.RE = FALSE, type = 'abundance', + include.sp = TRUE, ...) } \arguments{ @@ -105,8 +105,8 @@ n.factors <- 3 phi <- runif(n.factors, 3 / 1, 3 / .1) dat <- simMsNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, - beta = beta, alpha = alpha, sp = TRUE, - family = 'Poisson', factor.model = TRUE, + beta = beta, alpha = alpha, sp = TRUE, + family = 'Poisson', factor.model = TRUE, n.factors = n.factors, phi = phi, cov.model = 'exponential') # Split into fitting and prediction data set pred.indx <- sample(1:J, round(J * .25), replace = FALSE) @@ -129,7 +129,7 @@ det.covs <- list(det.cov.1 = X.p[, , 2], data.list <- list(y = y, abund.covs = abund.covs, det.covs = det.covs, - coords = coords) + coords = coords) # Initial values prior.list <- list(beta.comm.normal = list(mean = 0, var = 2.72), @@ -142,7 +142,7 @@ inits.list <- list(alpha.comm = 0, beta.comm = 0, beta = 0, alpha = 0, - phi = 3 / .5, + phi = 3 / .5, tau.sq.beta = 1, tau.sq.alpha = 1, N = apply(y, c(1, 2), max, na.rm = TRUE)) @@ -153,21 +153,21 @@ batch.length <- 25 accept.rate <- 0.43 out <- sfMsNMix(abund.formula = ~ abund.cov, - det.formula = ~ det.cov.1 + det.cov.2, - data = data.list, - inits = inits.list, - family = 'Poisson', - n.factors = n.factors, - n.batch = n.batch, - batch.length = batch.length, - accept.rate = 0.43, - cov.model = 'exponential', - n.neighbors = 5, - tuning = tuning, - priors = prior.list, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1) + det.formula = ~ det.cov.1 + det.cov.2, + data = data.list, + inits = inits.list, + family = 'Poisson', + n.factors = n.factors, + n.batch = n.batch, + batch.length = batch.length, + accept.rate = 0.43, + cov.model = 'exponential', + n.neighbors = 5, + tuning = tuning, + priors = prior.list, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1) summary(out, level = 'community') diff --git a/man/predict.spAbund.Rd b/man/predict.spAbund.Rd index 70ffa87..6482f2f 100644 --- a/man/predict.spAbund.Rd +++ b/man/predict.spAbund.Rd @@ -1,6 +1,6 @@ \name{predict.spAbund} \alias{predict.spAbund} -\title{Function for prediction at new locations for single-species spatial abundance GLMs} +\title{Function for prediction at new locations for univariate spatial GLMMs} \description{ The function \code{predict} collects posterior predictive samples for a set of new locations given an object of class `spAbund`. @@ -8,8 +8,8 @@ \usage{ \method{predict}{spAbund}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, - ignore.RE = FALSE, z.0.samples, include.sp = TRUE, ...) + verbose = TRUE, n.report = 100, + ignore.RE = FALSE, z.0.samples, include.sp = TRUE, ...) } \arguments{ @@ -86,7 +86,7 @@ phi <- 3 / .5 family <- 'NB' cov.model = 'exponential' dat <- simAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, - kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB', + kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB', sigma.sq = sigma.sq, phi = phi, cov.model = cov.model) # Split into fitting and prediction data set @@ -100,9 +100,9 @@ coords <- as.matrix(dat$coords[-pred.indx, ]) coords.0 <- as.matrix(dat$coords[pred.indx, ]) abund.covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.cov.2 = X[, , 3], - abund.cov.3 = X[, , 4]) + abund.cov.1 = X[, , 2], + abund.cov.2 = X[, , 3], + abund.cov.3 = X[, , 4]) data.list <- list(y = y, covs = abund.covs, coords = coords) @@ -110,8 +110,7 @@ data.list <- list(y = y, covs = abund.covs, coords = coords) prior.list <- list(beta.normal = list(mean = 0, var = 100), kappa.unif = c(0.001, 10)) # Starting values -inits.list <- list(beta = 0, - kappa = kappa) +inits.list <- list(beta = 0, kappa = kappa) n.batch <- 5 batch.length <- 25 @@ -120,20 +119,20 @@ n.thin <- 1 n.chains <- 1 out <- spAbund(formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - priors = prior.list, - accept.rate = 0.43, - n.neighbors = 5, - cov.model = cov.model, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + n.neighbors = 5, + cov.model = cov.model, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Predict at new locations ------------------------------------------------ colnames(X.0) <- c('intercept', 'abund.cov', 'abund.cov.2', 'abund.cov.3') diff --git a/man/predict.spDS.Rd b/man/predict.spDS.Rd index 53cdca8..13b381b 100644 --- a/man/predict.spDS.Rd +++ b/man/predict.spDS.Rd @@ -8,8 +8,8 @@ \usage{ \method{predict}{spDS}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, ignore.RE = FALSE, - type = 'abundance', include.sp = TRUE, ...) + verbose = TRUE, n.report = 100, ignore.RE = FALSE, + type = 'abundance', include.sp = TRUE, ...) } \arguments{ @@ -121,12 +121,12 @@ coords.0 <- dat$coords[pred.indx, ] covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, covs = covs, dist.breaks = dist.breaks, - coords = coords, + coords = coords, offset = offset) # Priors @@ -155,9 +155,9 @@ out <- spDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, family = 'NB', det.func = 'halfnormal', transect = 'point', - cov.model = 'exponential', - NNGP = TRUE, - n.neighbors = 5, + cov.model = 'exponential', + NNGP = TRUE, + n.neighbors = 5, tuning = tuning, priors = prior.list, accept.rate = 0.43, diff --git a/man/predict.spNMix.Rd b/man/predict.spNMix.Rd index d308db4..33a3054 100644 --- a/man/predict.spNMix.Rd +++ b/man/predict.spNMix.Rd @@ -8,8 +8,8 @@ \usage{ \method{predict}{spNMix}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, ignore.RE = FALSE, - type = 'abundance', include.sp = TRUE, ...) + verbose = TRUE, n.report = 100, ignore.RE = FALSE, + type = 'abundance', include.sp = TRUE, ...) } \arguments{ @@ -92,7 +92,7 @@ kappa <- 0.5 sp <- TRUE cov.model <- 'exponential' dat <- simNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha, - kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, + kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, phi = phi, sigma.sq = sigma.sq, cov.model = cov.model, family = 'NB') @@ -113,27 +113,26 @@ w.0 <- dat$w[pred.indx] abund.covs <- X colnames(abund.covs) <- c('int', 'abund.cov.1') -det.covs <- list(det.cov.1 = X.p[, , 2], - det.cov.2 = X.p[, , 3]) +det.covs <- list(det.cov.1 = X.p[, , 2], det.cov.2 = X.p[, , 3]) data.list <- list(y = y, - abund.covs = abund.covs, - det.covs = det.covs, + abund.covs = abund.covs, + det.covs = det.covs, coords = coords) # Priors prior.list <- list(beta.normal = list(mean = rep(0, p.abund), - var = rep(100, p.abund)), - alpha.normal = list(mean = rep(0, p.det), - var = rep(2.72, p.det)), + var = rep(100, p.abund)), + alpha.normal = list(mean = rep(0, p.det), + var = rep(2.72, p.det)), kappa.unif = c(0, 10)) # Starting values inits.list <- list(alpha = alpha, - beta = beta, - kappa = kappa, - phi = 3 / 0.5, - sigma.sq = 1, - N = apply(y, 1, max, na.rm = TRUE)) + beta = beta, + kappa = kappa, + phi = 3 / 0.5, + sigma.sq = 1, + N = apply(y, 1, max, na.rm = TRUE)) # Tuning values tuning.list <- list(phi = 0.5, kappa = 0.5, beta = 0.1, alpha = 0.1, w = 0.1) @@ -145,22 +144,22 @@ n.thin <- 1 n.chains <- 1 out <- spNMix(abund.formula = ~ abund.cov.1, - det.formula = ~ det.cov.1 + det.cov.2, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - priors = prior.list, - NNGP = TRUE, - cov.model = 'spherical', - n.neighbors = 10, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1 + det.cov.2, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + NNGP = TRUE, + cov.model = 'spherical', + n.neighbors = 10, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) # Predict at new locations ------------------------------------------------ diff --git a/man/predict.svcAbund.Rd b/man/predict.svcAbund.Rd index fa3e701..3798c81 100644 --- a/man/predict.svcAbund.Rd +++ b/man/predict.svcAbund.Rd @@ -8,8 +8,8 @@ \usage{ \method{predict}{svcAbund}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, ignore.RE = FALSE, - z.0.samples, ...) + verbose = TRUE, n.report = 100, ignore.RE = FALSE, + z.0.samples, ...) } \arguments{ @@ -35,7 +35,7 @@ \item{ignore.RE}{logical value that specifies whether or not to remove unstructured random effects from the subsequent predictions. If \code{TRUE}, random effects will be included. If \code{FALSE}, random effects will be set to 0 and predictions will only be generated from the fixed effects.} - \item{z.0.samples}{a matrix with rows corresponding to MCMC samples and columns corresponding to prediction locations containing the full posterior samples of the predicted binary portion of a Gaussian hurdle model. In the context of abundance models, this typically corresponds to estimates of the presence or absence of the species at the location. When using \code{spOccupancy} to generate the first stage samples of the zero-inflated Gaussian model, this is the object contained in the \code{z.0.samples} object of the predition function for th spOccupancy object. Ignored for all model types other than zero-inflated Gaussian.} + \item{z.0.samples}{a matrix with rows corresponding to MCMC samples and columns corresponding to prediction locations containing the full posterior samples of the predicted binary portion of a zero-inflated Gaussian model. In the context of abundance models, this typically corresponds to estimates of the presence or absence of the species at the location. When using \code{spOccupancy} to generate the first stage samples of the zero-inflated Gaussian model, this is the object contained in the \code{z.0.samples} object of the predition function for th spOccupancy object. Ignored for all model types other than zero-inflated Gaussian.} \item{...}{currently no additional arguments} } @@ -59,7 +59,7 @@ abundance values.} \item{w.0.samples}{a three-dimensional array of posterior predictive samples - for the spatial random effects, with dimensions corresponding to MCMC iteration, + for the spatially-varying coefficients, with dimensions corresponding to MCMC iteration, coefficient, and site.} The return object will include additional objects used for standard @@ -87,8 +87,8 @@ tau.sq <- 2 # Get all the data dat <- simAbund(J.x = J.x, J.y = J.y, beta = beta, tau.sq = tau.sq, - mu.RE = mu.RE, sp = sp, svc.cols = svc.cols, family = 'Gaussian', - cov.model = cov.model, sigma.sq = sigma.sq, phi = phi) + mu.RE = mu.RE, sp = sp, svc.cols = svc.cols, family = 'Gaussian', + cov.model = cov.model, sigma.sq = sigma.sq, phi = phi) # Prep the data for spAbundance ------------------------------------------- y <- dat$y @@ -109,18 +109,16 @@ covs <- cbind(X) colnames(covs) <- c('int', 'cov.1', 'cov.2', 'cov.3') # Data list bundle -data.list <- list(y = y, - covs = covs, - coords = coords) +data.list <- list(y = y, covs = covs, coords = coords) # Priors prior.list <- list(beta.normal = list(mean = 0, var = 1000), - sigma.sq.ig = list(a = 2, b = 1), tau.sq = c(2, 1), - sigma.sq.mu.ig = list(a = 2, b = 1), + sigma.sq.ig = list(a = 2, b = 1), tau.sq = c(2, 1), + sigma.sq.mu.ig = list(a = 2, b = 1), phi.unif = list(a = 3 / 1, b = 3 / 0.1)) # Starting values inits.list <- list(beta = 0, alpha = 0, - sigma.sq = 1, phi = phi, tau.sq = 2, sigma.sq.mu = 0.5) + sigma.sq = 1, phi = phi, tau.sq = 2, sigma.sq.mu = 0.5) # Tuning tuning.list <- list(phi = 1) @@ -138,7 +136,7 @@ out <- svcAbund(formula = ~ cov.1 + cov.2 + cov.3, inits = inits.list, priors = prior.list, accept.rate = 0.43, - family = 'Gaussian', + family = 'Gaussian', cov.model = "exponential", tuning = tuning.list, n.omp.threads = 1, diff --git a/man/predict.svcMsAbund.Rd b/man/predict.svcMsAbund.Rd index ea65b73..aa9d20f 100644 --- a/man/predict.svcMsAbund.Rd +++ b/man/predict.svcMsAbund.Rd @@ -1,6 +1,6 @@ \name{predict.svcMsAbund} \alias{predict.svcMsAbund} -\title{Function for prediction at new locations for spatial factor multi-species abundance GLMs} +\title{Function for prediction at new locations for multivariate spatially-varying coefficient GLMMs} \description{ The function \code{predict} collects posterior predictive samples for a set of new locations given an object of class `svcMsAbund`. @@ -8,8 +8,8 @@ \usage{ \method{predict}{svcMsAbund}(object, X.0, coords.0, n.omp.threads = 1, - verbose = TRUE, n.report = 100, ignore.RE = FALSE, - z.0.samples, ...) + verbose = TRUE, n.report = 100, ignore.RE = FALSE, + z.0.samples, ...) } \arguments{ @@ -61,8 +61,9 @@ abundance values with dimensions corresponding to posterior predictive sample, species, site, and replicate.} - \item{w.0.samples}{a three-dimensional array of posterior predictive samples for the - latent factors.} + \item{w.0.samples}{a four-dimensional array of posterior predictive samples for the + spatial factors for each spatially-varying coefficient. Dimensions correspond to MCMC sample, + spatial factor, site, and spatially varying coefficient.} The return object will include additional objects used for standard extractor functions. @@ -73,7 +74,7 @@ set.seed(408) J.x <- 8 J.y <- 8 J <- J.x * J.y -n.rep <- sample(3, size = J, replace = TRUE) +n.rep <- rep(1, J) n.sp <- 6 # Community-level covariate effects beta.mean <- c(-2, 0.5) @@ -86,44 +87,43 @@ for (i in 1:p.abund) { beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i])) } sp <- TRUE -kappa <- runif(n.sp, 0.1, 1) factor.model <- TRUE -n.factors <- 3 +n.factors <- 2 +svc.cols <- c(1, 2) cov.model <- 'spherical' -phi <- runif(n.factors, 3 / 1, 3 / .1) +tau.sq <- runif(n.sp, 0.1, 2) +phi <- runif(n.factors * length(svc.cols), 3 / 1, 3 / .1) dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', + mu.RE = mu.RE, sp = sp, family = 'Gaussian', tau.sq = tau.sq, factor.model = factor.model, n.factors = n.factors, - phi = phi, cov.model = cov.model) + phi = phi, cov.model = cov.model, svc.cols = svc.cols) # Split into fitting and prediction data set pred.indx <- sample(1:J, round(J * .25), replace = FALSE) -y <- dat$y[, -pred.indx, , drop = FALSE] +y <- dat$y[, -pred.indx, drop = FALSE] # Occupancy covariates -X <- dat$X[-pred.indx, , , drop = FALSE] +X <- dat$X[-pred.indx, , drop = FALSE] # Coordinates coords <- dat$coords[-pred.indx, ] # Prediction values -y.0 <- dat$y[, pred.indx, , drop = FALSE] -X.0 <- dat$X[pred.indx, , , drop = FALSE] +y.0 <- dat$y[, pred.indx, drop = FALSE] +X.0 <- dat$X[pred.indx, , drop = FALSE] coords.0 <- dat$coords[pred.indx, ] # Package all data into a list -covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2]) -data.list <- list(y = y, - covs = covs, - coords = coords) +covs <- data.frame(abund.cov.1 = X[, 2]) +data.list <- list(y = y, covs = covs, coords = coords) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), - kappa.unif = list(a = 0, b = 10), - phi.unif = list(a = 3 / 1, b = 3 / .1), + tau.sq.ig = list(a = .01, b = .01), + phi.unif = list(a = 3 / 1, b = 3 / .1), tau.sq.beta.ig = list(a = .1, b = .1)) inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - phi = 3 / .5, - tau.sq.beta = 1) + beta = 0, + kappa = 0.5, + tau.sq = 1, + phi = 3 / .5, + tau.sq.beta = 1) tuning.list <- list(kappa = 0.3, beta = 0.1, lambda = 0.5, w = 0.5, phi = 1) @@ -134,21 +134,24 @@ n.burn <- 20 n.thin <- 1 n.chains <- 1 -out <- sfMsAbund(formula = ~ abund.cov.1, - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - tuning = tuning.list, - batch.length = batch.length, - n.omp.threads = 1, - n.factors = 1, - verbose = TRUE, - n.neighbors = 5, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) +out <- svcMsAbund(formula = ~ abund.cov.1, + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + tuning = tuning.list, + batch.length = batch.length, + n.omp.threads = 1, + svc.cols = c(1, 2), + n.factors = n.factors, + cov.model = 'exponential', + family = 'Gaussian', + verbose = TRUE, + n.neighbors = 5, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Predict at new locations out.pred <- predict(out, X.0, coords.0) diff --git a/man/sfMsAbund.Rd b/man/sfMsAbund.Rd index d980ff6..08db9f4 100644 --- a/man/sfMsAbund.Rd +++ b/man/sfMsAbund.Rd @@ -1,19 +1,19 @@ \name{sfMsAbund} \alias{sfMsAbund} -\title{Function for Fitting Spatial Factor Multi-Species Abundance GLMs} +\title{Function for Fitting Spatial Factor Multivariate Abundance GLMMs} \description{ - The function \code{sfMsAbund} fits multi-species spatial abundance GLMs with species correlations (i.e., a spatially-explicit abundace-based joint species distribution model). We use a spatial factor modeling approach. Currently, models are implemented using a Nearest Neighbor Gaussian Process. Future development may allow for running the models using full Gaussian Processes. + The function \code{sfMsAbund} fits multivariate spatial abundance GLMMs with species correlations (i.e., a spatially-explicit abundace-based joint species distribution model). We use a spatial factor modeling approach. Currently, models are implemented using a Nearest Neighbor Gaussian Process. Future development may allow for running the models using full Gaussian Processes. } \usage{ sfMsAbund(formula, data, inits, priors, - tuning, cov.model = 'exponential', NNGP = TRUE, - n.neighbors = 15, search.type = 'cb', n.factors, - n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', - n.omp.threads = 1, verbose = TRUE, n.report = 100, - n.burn = round(.10 * n.batch * batch.length), n.thin = 1, n.chains = 1, - save.fitted = TRUE, ...) + tuning, cov.model = 'exponential', NNGP = TRUE, + n.neighbors = 15, search.type = 'cb', n.factors, + n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', + n.omp.threads = 1, verbose = TRUE, n.report = 100, + n.burn = round(.10 * n.batch * batch.length), n.thin = 1, n.chains = 1, + save.fitted = TRUE, ...) } \arguments{ @@ -24,14 +24,16 @@ sfMsAbund(formula, data, inits, priors, and slopes are allowed using lme4 syntax (Bates et al. 2015).} \item{data}{a list containing data necessary for model fitting. - Valid tags are \code{y}, \code{covs}, \code{coords}, and \code{offset}. + Valid tags are \code{y}, \code{covs}, \code{z}, \code{coords}, and \code{offset}. \code{y} is a two or three-dimensional array of observed count data. The first dimension of the array is equal to the number of species and the second dimension is equal to the number of sites. If specified as a three-dimensional array, the third dimension corresponds to replicate observations at each site (e.g., sub-samples, repeated sampling over multiple seasons). \code{covs} is a list - containing the variables used in the model. Each list element is a different + containing the variables used in the model. If a data frame, each row + of \code{covs} is a site and each column is a variable. + If a list, each list element is a different covariate, which can be site-level or observation-level. Site-level covariates are specified as a vector of length \eqn{J}{J}, while observation-level covariates are specified as a matrix or data frame with the number of rows equal to \eqn{J}{J} @@ -150,13 +152,13 @@ sfMsAbund(formula, data, inits, priors, decreases in computation time. However, the value can be anywhere between 1 and the number of species in the modeled community.} - \item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC + \item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{batch.length}{the length of each MCMC batch to run for the Adaptive + \item{batch.length}{the length of each MCMC batch to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{accept.rate}{target acceptance rate for Adaptive MCMC. Defaul is + \item{accept.rate}{target acceptance rate for adaptive MCMC. Defaul is 0.43. See Roberts and Rosenthal (2009) for details.} \item{family}{the distribution to use for the abundance. Currently @@ -263,7 +265,8 @@ sfMsAbund(formula, data, inits, priors, corresponding to MCMC samples, species, site, and replicate.} \item{w.samples}{a three-dimensional array of posterior samples for - the latent spatial random effects for each latent factor.} + the latent effects for each latent factor. Array dimensions correspond + to MCMC sample, latent factor, then site.} \item{sigma.sq.mu.samples}{a \code{coda} object of posterior samples for variances of random effects included in the abundance portion @@ -302,7 +305,7 @@ p.abund <- length(beta.mean) tau.sq.beta <- c(0.2, 1.2) # Random effects (two random intercepts) mu.RE <- list(levels = c(10, 15), - sigma.sq.mu = c(0.43, 0.5)) + sigma.sq.mu = c(0.43, 0.5)) # Draw species-level effects from community means. beta <- matrix(NA, nrow = n.sp, ncol = p.abund) for (i in 1:p.abund) { @@ -315,7 +318,7 @@ phi <- runif(n.factors, 3/1, 3 / .1) kappa <- runif(n.sp, 0.1, 1) dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', + mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', factor.model = factor.model, phi = phi, cov.model = 'exponential', n.factors = n.factors) @@ -326,21 +329,16 @@ coords <- dat$coords # Package all data into a list covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.factor.1 = X.re[, , 1], - abund.factor.2 = X.re[, , 2]) -data.list <- list(y = y, - covs = covs, - coords = coords) + abund.cov.1 = X[, , 2], + abund.factor.1 = X.re[, , 1], + abund.factor.2 = X.re[, , 2]) +data.list <- list(y = y, covs = covs, coords = coords) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), kappa.unif = list(a = 0, b = 10), - phi.unif = list(a = 3 / 1, b = 3 / .1), + phi.unif = list(a = 3 / 1, b = 3 / .1), tau.sq.beta.ig = list(a = .1, b = .1)) -inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - tau.sq.beta = 1, - phi = 3 / 0.5) +inits.list <- list(beta.comm = 0, beta = 0, kappa = 0.5, + tau.sq.beta = 1, phi = 3 / 0.5) # Small n.batch <- 2 @@ -349,22 +347,23 @@ n.burn <- 20 n.thin <- 1 n.chains <- 1 -out <- sfMsAbund(formula = ~ abund.cov.1 + (1 | abund.factor.1) + (1 | abund.factor.2), - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - NNGP = TRUE, - cov.model = 'exponential', - n.neighbors = 5, - n.factors = n.factors, - batch.length = batch.length, - n.omp.threads = 3, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) +out <- sfMsAbund(formula = ~ abund.cov.1 + (1 | abund.factor.1) + + (1 | abund.factor.2), + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + NNGP = TRUE, + cov.model = 'exponential', + n.neighbors = 5, + n.factors = n.factors, + batch.length = batch.length, + n.omp.threads = 3, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) } diff --git a/man/sfMsDS.Rd b/man/sfMsDS.Rd index cceb9cd..6f19a6a 100644 --- a/man/sfMsDS.Rd +++ b/man/sfMsDS.Rd @@ -36,10 +36,10 @@ sfMsDS(abund.formula, det.formula, data, inits, priors, of species, second dimension equal to the number of sites, and third dimension equal to the maximum number of replicates at a given site. \code{covs} is a matrix or data frame containing the variables - used in the the abundance and/or the detection portion of the model, with + used in the abundance and/or the detection portion of the model, with \eqn{J}{J} rows for each column (variable). \code{dist.breaks} is a vector of distances that denote the breakpoints of the distance bands. \code{dist.breaks} should - have length equal to the number of columns in \code{y} plus one. \code{offset} is an + have length equal to the third dimension of \code{y} plus one. \code{offset} is an offset that can be used to scale estimates from abundance per transect to density per some desired unit of measure. This can be either a single value or a vector with an offset value for each site (e.g., if transects differ in length). \code{coords} is a matrix @@ -88,7 +88,7 @@ sfMsDS(abund.formula, det.formula, data, inits, priors, Gamma distribution. The hyperparameters of the inverse-Gamma distribution are passed as a list of length two with first and second elements corresponding to the shape and scale parameters, respectively, which are each specified as - vectors of length equal to the number of random intercepts or of length one + vectors of length equal to the number of random effects or of length one if priors are the same for all random effect variances. \code{kappa} is the negative binomial dispersion parameter for each species and is assumed to follow a uniform distribution. The hyperparameters of the uniform distribution @@ -109,7 +109,7 @@ sfMsDS(abund.formula, det.formula, data, inits, priors, diagonal elements are fixed at 1. The lower triangular elements are assigned a standard normal prior (i.e., mean 0 and variance 1).} -\item{tuning}{a list with each tag corresponding to a parameter name, whose +\item{tuning}{a list with each tag corresponding to a parameter name, whose value defines the initial variance of the adaptive sampler. Valid tags are \code{beta}, \code{alpha}, \code{lambda} (the latent factor loadings), \code{w} (the latent factors), \code{beta.star} (the abundance @@ -259,17 +259,19 @@ of neighbors using WAIC.} for the spatial factor loadings.} \item{w.samples}{a three-dimensional array of posterior samples for - the latent effects for each spatial factor.} + the latent effects for each spatial factor. Array dimensions correspond + to MCMC sample, spatial factor, then site.} \item{N.samples}{a three-dimensional array of posterior samples for the latent abundance values for each species. Note that these values always represent transect-level abundance, even when an offset is - supplied.} + supplied. Array dimensions correspond to MCMC sample, species, then site.} \item{mu.samples}{a three-dimensional array of posterior samples for the latent expected abundance values for each species. When an offset is supplied in the \code{data} object, these correspond to expected - abundance per unit area (i.e., density).} + abundance per unit area (i.e., density). Array dimensions correspond to + MCMC sample, species, then site.} \item{sigma.sq.mu.samples}{a \code{coda} object of posterior samples for variances of random effects included in the abundance portion @@ -350,7 +352,7 @@ phi <- runif(n.factors, 3 / 1, 3 / .2) cov.model <- 'exponential' dat <- simMsDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, mu.RE = mu.RE, p.RE = p.RE, sp = sp, cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect, factor.model = factor.model, @@ -364,25 +366,24 @@ dist.breaks <- dat$dist.breaks covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, - covs = covs, - dist.breaks = dist.breaks, - coords = coords, + covs = covs, + dist.breaks = dist.breaks, + coords = coords, offset = offset) # Priors prior.list <- list(beta.comm.normal = list(mean = 0, var = 10), - alpha.comm.normal = list(mean = 0, - var = 10), + alpha.comm.normal = list(mean = 0, var = 10), kappa.unif = list(0, 100), - phi.unif = list(3 / 1, 3 / .1), + phi.unif = list(3 / 1, 3 / .1), tau.sq.beta.ig = list(a = 0.1, b = 0.1), tau.sq.alpha.ig = list(a = 0.1, b = 0.1)) # Starting values inits.list <- list(alpha.comm = 0, beta.comm = 0, beta = 0, - alpha = 0, kappa = 1, phi = 3 / .5) + alpha = 0, kappa = 1, phi = 3 / .5) tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1, kappa = 0.8, lambda = 1, w = 1, phi = 0.8) @@ -394,26 +395,26 @@ n.thin <- 1 n.chains <- 1 out <- sfMsDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, - det.formula = ~ det.cov.1, - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - family = 'Poisson', - det.func = 'halfnormal', - transect = transect, - tuning = tuning, - cov.model = 'exponential', - NNGP = TRUE, - n.neighbors = 5, - n.factors = n.factors, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 10, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + det.formula = ~ det.cov.1, + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + family = 'Poisson', + det.func = 'halfnormal', + transect = transect, + tuning = tuning, + cov.model = 'exponential', + NNGP = TRUE, + n.neighbors = 5, + n.factors = n.factors, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 10, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out, level = 'community') } diff --git a/man/simAbund.Rd b/man/simAbund.Rd index 90c4dfc..10d12e9 100644 --- a/man/simAbund.Rd +++ b/man/simAbund.Rd @@ -1,15 +1,15 @@ \name{simAbund} \alias{simAbund} -\title{Simulate Single-Species Count Data without Imperfect Detection} +\title{Simulate Univariate Data for Testing GLMMs} \description{ - The function \code{simAbund} simulates single-species count data without imperfect detection for simulation studies, power assessments, or function testing. Data can be optionally simulated with a spatial Gaussian Process in the model. Non-spatial random intercepts can also be included in the model. + The function \code{simAbund} simulates univariate data without imperfect detection for simulation studies, power assessments, or function testing related to GLMMs. Data can be optionally simulated with a spatial Gaussian Process in the model. Non-spatial random effects can also be included in the model. } \usage{ simAbund(J.x, J.y, n.rep, n.rep.max, beta, kappa, tau.sq, mu.RE = list(), - offset = 1, sp = FALSE, svc.cols = 1, cov.model, sigma.sq, phi, nu, - family = 'Poisson', z, x.positive = FALSE, ...) + offset = 1, sp = FALSE, svc.cols = 1, cov.model, sigma.sq, phi, nu, + family = 'Poisson', z, x.positive = FALSE, ...) } \arguments{ @@ -58,8 +58,7 @@ simAbund(J.x, J.y, n.rep, n.rep.max, beta, kappa, tau.sq, mu.RE = list(), } \author{ - Jeffrey W. Doser \email{doserjef@msu.edu}, \cr - Andrew O. Finley \email{finleya@msu.edu} + Jeffrey W. Doser \email{doserjef@msu.edu} \cr } \value{ @@ -67,7 +66,7 @@ simAbund(J.x, J.y, n.rep, n.rep.max, beta, kappa, tau.sq, mu.RE = list(), \item{X}{a three-dimensional numeric design array of covariates with dimensions corresponding to sites, replicates, and number of covariates (including an intercept) for the model.} \item{coords}{a \eqn{J \times 2}{J x 2} numeric matrix of coordinates of each site. Required for spatial models.} - \item{w}{a \eqn{J \times 1}{J x 1} matrix of the spatial random effects. Only used to simulate data when \code{sp = TRUE}.} + \item{w}{a matrix of the spatial random effects. Only used to simulate data when \code{sp = TRUE}. If simulating data with spatially-varying coefficients, the number of columns equals the number of spatially-varying coefficients and each row corresponds to a site.} \item{mu}{a \code{J x max(n.rep)} matrix of the expected abundance values for each site and replicate survey.} \item{y}{a \code{J x max(n.rep)} matrix of the raw count data for each site and replicate combination.} \item{X.re}{a numeric three-dimensional array containing the levels of any abundance random effect included in the model. Only relevant when abundance random effects are specified in \code{mu.RE}.} @@ -82,11 +81,10 @@ J <- J.x * J.y n.rep <- sample(3, J, replace = TRUE) beta <- c(0, -1.5, 0.3, -0.8) p.abund <- length(beta) -mu.RE <- list(levels = c(30), - sigma.sq.mu = c(1.3)) +mu.RE <- list(levels = c(30), sigma.sq.mu = c(1.3)) kappa <- 0.5 sp <- FALSE family <- 'NB' dat <- simAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, - kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') + kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') } diff --git a/man/simDS.Rd b/man/simDS.Rd index 38d8f00..8794620 100644 --- a/man/simDS.Rd +++ b/man/simDS.Rd @@ -3,7 +3,7 @@ \title{Simulate Single-Species Distance Sampling Data} \description{ - The function \code{simDS} simulates single-species distance sampling data for simulation studies, power assessments, or function testing. Data can be optionally simulated with a spatial Gaussian Process in the abundance portion of the model. Non-spatial random effects can also be included in the detection or abundance portions of the distance sampling. + The function \code{simDS} simulates single-species distance sampling data for simulation studies, power assessments, or function testing. Data can be optionally simulated with a spatial Gaussian Process in the abundance portion of the model. Non-spatial random effects can also be included in the detection or abundance portions of the distance sampling model. } \usage{ @@ -33,7 +33,7 @@ simDS(J.x, J.y, n.bins, bin.width, beta, alpha, det.func, transect = 'line', \item{transect}{the type of transect. Currently supports line transects (\code{'line'}) or circular transects (i.e., point counts; \code{'point'}).} -\item{kappa}{a single numeric value containing the dispersion parameter for the abundance portion of the N-mixture model. Only relevant when \code{family = 'NB'}.} +\item{kappa}{a single numeric value containing the dispersion parameter for the abundance portion of the hierarchical distance sampling model. Only relevant when \code{family = 'NB'}.} \item{mu.RE}{a list used to specify the non-spatial random intercepts included in the abundance portion of the model. The list must have two tags: \code{levels} and \code{sigma.sq.mu}. \code{levels} is a vector of length equal to the number of distinct random intercepts to include in the model and contains the number of levels there are in each intercept. \code{sigma.sq.mu} is a vector of length equal to the number of distinct random intercepts to include in the model and contains the variances for each random effect. If not specified, no random effects are included in the abundance portion of the model.} @@ -41,7 +41,7 @@ simDS(J.x, J.y, n.bins, bin.width, beta, alpha, det.func, transect = 'line', \item{offset}{either a single numeric value or a vector of length \code{J} that contains the offset for each location in the data set.} -\item{sp}{a logical value indicating whether to simulate a spatially-explicit N-mixture model with a Gaussian process. By default set to \code{FALSE}.} +\item{sp}{a logical value indicating whether to simulate a spatially-explicit HDS model with a Gaussian process. By default set to \code{FALSE}.} \item{cov.model}{a quoted keyword that specifies the covariance function used to model the spatial dependence structure among the latent abundance values. Supported covariance model key words are: \code{"exponential"}, \code{"matern"}, \code{"spherical"}, and \code{"gaussian"}.} @@ -72,12 +72,12 @@ simDS(J.x, J.y, n.bins, bin.width, beta, alpha, det.func, transect = 'line', \item{mu}{a \eqn{J \times 1}{J x 1} matrix of the expected abundance values for each site.} \item{N}{a length \eqn{J}{J} vector of the latent abundances at each site.} \item{p}{ a length J vector of the detection probabilities at each site.} - \item{pi.full}{a \code{J x n.bins + 1} vector of the cell-specific detection probabilities for each site, where the last column indicates the probability of not detecting the individual at that site.} + \item{pi.full}{a \code{J x n.bins + 1} vector of the cell-specific detection probabilities for each site, where the last column indicates the probability of not detecting an individual at that site.} \item{y}{a \code{J x max(n.bins)} matrix of the raw count data for each site and distance bin.} \item{X.p.re}{a numeric matrix containing the levels of any detection random effect included in the model. Only relevant when detection random effects are specified in \code{p.RE}.} \item{X.re}{a numeric matrix containing the levels of any abundance random effect included in the model. Only relevant when abundance random effects are specified in \code{mu.RE}.} \item{alpha.star}{a numeric vector that contains the simulated detection random effects for each given level of the random effects included in the detection model. Only relevant when detection random effects are included in the model.} - \item{beta.star}{a numeric vector that contains the simulated abundance random effects for each given level of the random effects included in the N-mixture model. Only relevant when abundance random effects are included in the model.} + \item{beta.star}{a numeric vector that contains the simulated abundance random effects for each given level of the random effects included in the HDS model. Only relevant when abundance random effects are included in the model.} } \examples{ @@ -106,8 +106,8 @@ offset <- 1.8 transect <- 'point' dat <- simDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, - mu.RE = mu.RE, p.RE = p.RE, sp = sp, - sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, + beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + mu.RE = mu.RE, p.RE = p.RE, sp = sp, + sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect) } diff --git a/man/simMsAbund.Rd b/man/simMsAbund.Rd index c948fe3..571d22f 100644 --- a/man/simMsAbund.Rd +++ b/man/simMsAbund.Rd @@ -1,9 +1,9 @@ \name{simMsAbund} \alias{simMsAbund} -\title{Simulate Multi-Species Count Data without Imperfect Detection} +\title{Simulate Multivariate Data for Testing GLMMs} \description{ - The function \code{simMsAbund} simulates multi-species count data without imperfect detection for simulation studies, power assessments, or function testing. Data can be optionally simulated with a spatial Gaussian Process in the model, as well as an option to allow for species correlations using a factor modeling approach. Non-spatial random intercepts can also be included in the abundance portions of the model. + The function \code{simMsAbund} simulates multivariate data without imperfect detection for simulation studies, power assessments, or function testing related to GLMMs. Data can be optionally simulated with a spatial Gaussian Process in the model, as well as an option to allow for species correlations using a factor modeling approach. Non-spatial random effects can also be included in the abundance portions of the model. } \usage{ @@ -29,7 +29,7 @@ simMsAbund(J.x, J.y, n.rep, n.rep.max, n.sp, beta, kappa, tau.sq, mu.RE = list() \item{kappa}{a numeric vector of length \code{n.sp} containing the dispersion parameter for the model for each species. Only relevant when \code{family = 'NB'}.} -\item{tau.sq}{a numeric vector of length \code{n.sp} containing the residual variance parameters for the model for each species. Only relevant for Gaussian or Gaussian-hurdle models.} +\item{tau.sq}{a numeric vector of length \code{n.sp} containing the residual variance parameters for the model for each species. Only relevant for Gaussian or zero-inflated Gaussian models.} \item{mu.RE}{a list used to specify the non-spatial random intercepts included in the model. The list must have two tags: \code{levels} and \code{sigma.sq.mu}. \code{levels} is a vector of length equal to the number of distinct random intercepts to include in the model and contains the number of levels there are in each intercept. \code{sigma.sq.mu} is a vector of length equal to the number of distinct random intercepts to include in the model and contains the variances for each random effect. If not specified, no random effects are included in the model.} @@ -57,14 +57,13 @@ simMsAbund(J.x, J.y, n.rep, n.rep.max, n.sp, beta, kappa, tau.sq, mu.RE = list() \item{family}{the distribution to use for the latent abundance process. Currently supports \code{'NB'} (negative binomial) and \code{'Poisson'}.} -\item{z}{a matrix with \code{n.sp} rows and \code{J} columns containing the binary presence/absence portion of a Gaussian-hurdle model for each species. Only relevant when \code{family = 'Gaussian-hurdle'}.} +\item{z}{a matrix with \code{n.sp} rows and \code{J} columns containing the binary presence/absence portion of a zero-inflated Gaussian model for each species. Only relevant when \code{family = 'zi-Gaussian'}.} \item{...}{currently no additional arguments} } \author{ - Jeffrey W. Doser \email{doserjef@msu.edu}, \cr - Andrew O. Finley \email{finleya@msu.edu} + Jeffrey W. Doser \email{doserjef@msu.edu} \cr } \value{ @@ -72,7 +71,7 @@ simMsAbund(J.x, J.y, n.rep, n.rep.max, n.sp, beta, kappa, tau.sq, mu.RE = list() \item{X}{a three-dimensional numeric design array of covariates with dimensions corresponding to sites, replicates, and number of covariates (including an intercept) for the model.} \item{coords}{a \eqn{J \times 2}{J x 2} numeric matrix of coordinates of each site. Required for spatial models.} - \item{w}{a \eqn{N \times J}{N x J} matrix of the spatial random effects for each species. Only used to simulate data when \code{sp = TRUE}. If \code{factor.model = TRUE}, the first dimension is \code{n.factors}.} + \item{w}{a list of \eqn{N \times J}{N x J} matrices of the spatially-varying coefficients for each species. Each element of the list corresponds to a different spatially-varying coefficient. Only used to simulate data when \code{sp = TRUE}. If \code{factor.model = TRUE}, the first dimension of each matrix is \code{n.factors}.} \item{mu}{a \code{n.sp x J} matrix of the mean abundances for each species at each site.} \item{y}{a \code{n.sp x J x max(n.rep)} array of the raw count data for each species at each site and replicate combination. Sites with fewer than \code{max(n.rep)} replicates will contain \code{NA} values.} \item{X.re}{a numeric matrix containing the levels of any abundance random effect included in the model. Only relevant when abundance random effects are specified in \code{mu.RE}.} @@ -92,7 +91,7 @@ p.abund <- length(beta.mean) tau.sq.beta <- c(0.2, 1.2) # Random effects (two random intercepts) mu.RE <- list(levels = c(10, 15), - sigma.sq.mu = c(0.43, 0.5)) + sigma.sq.mu = c(0.43, 0.5)) # Draw species-level effects from community means. beta <- matrix(NA, nrow = n.sp, ncol = p.abund) for (i in 1:p.abund) { @@ -105,7 +104,7 @@ phi <- runif(n.factors, 3/1, 3 / .1) kappa <- runif(n.sp, 0.1, 1) dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', + mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', factor.model = factor.model, phi = phi, cov.model = 'spherical', n.factors = n.factors) } diff --git a/man/simMsDS.Rd b/man/simMsDS.Rd index 94c6b74..8aef368 100644 --- a/man/simMsDS.Rd +++ b/man/simMsDS.Rd @@ -8,7 +8,7 @@ \usage{ simMsDS(J.x, J.y, n.bins, bin.width, n.sp, beta, alpha, - det.func, transect = 'line', kappa, mu.RE = list(), + det.func, transect = 'line', kappa, mu.RE = list(), p.RE = list(), offset = 1, sp = FALSE, cov.model, sigma.sq, phi, nu, family = 'Poisson', factor.model = FALSE, n.factors, ...) @@ -26,9 +26,9 @@ simMsDS(J.x, J.y, n.bins, bin.width, n.sp, beta, alpha, \item{n.sp}{a single numeric value indicating the number of species to simulate count data.} -\item{beta}{a numeric matrix with \code{n.sp} rows containing the intercept and regression coefficient parameters for the abundance portion of the multi-species N-mixture model. Each row corresponds to the regression coefficients for a given species.} +\item{beta}{a numeric matrix with \code{n.sp} rows containing the intercept and regression coefficient parameters for the abundance portion of the multi-species hierarchical distance sampling (HDS) model. Each row corresponds to the regression coefficients for a given species.} -\item{alpha}{a numeric matrix with \code{n.sp} rows containing the intercept and regression coefficient parameters for the detection portion of the multi-species N-mixture model. Each row corresponds to the regression coefficients for a given species.} +\item{alpha}{a numeric matrix with \code{n.sp} rows containing the intercept and regression coefficient parameters for the detection portion of the multi-species HDS model. Each row corresponds to the regression coefficients for a given species.} \item{det.func}{the detection model used to describe how detection probability varies with distance. In other software, this is often referred to as the key function. Currently @@ -38,7 +38,7 @@ simMsDS(J.x, J.y, n.bins, bin.width, n.sp, beta, alpha, \item{transect}{the type of transect. Currently supports line transects (\code{'line'}) or circular transects (i.e., point counts; \code{'point'}).} -\item{kappa}{a numeric vector of length \code{n.sp} containing the dispersion parameter for the abundance portion of the N-mixture model for each species. Only relevant when \code{family = 'NB'}.} +\item{kappa}{a numeric vector of length \code{n.sp} containing the dispersion parameter for the abundance portion of the HDS model for each species. Only relevant when \code{family = 'NB'}.} \item{mu.RE}{a list used to specify the non-spatial random effects included in the abundance portion of the model. The list must have two tags: \code{levels} and \code{sigma.sq.mu}. \code{levels} is a vector of length equal to the number of distinct random effects to include in the model and contains the number of levels there are in each effect. \code{sigma.sq.mu} is a vector of length equal to the number of distinct random effects to include in the model and contains the variances for each random effect. If not specified, no random effects are included in the abundance portion of the model. An optional third tag, \code{beta.indx}, is a list that contains integers denoting the corresponding value of \code{beta} that each random effect corresponds to. This allows specification of random intercepts as well as slopes. By default, all effects are assumed to be random intercepts.} @@ -67,7 +67,7 @@ simMsDS(J.x, J.y, n.bins, bin.width, n.sp, beta, alpha, } \author{ - Jeffrey W. Doser \email{doserjef@msu.edu}, \cr + Jeffrey W. Doser \email{doserjef@msu.edu} \cr } \value{ @@ -126,7 +126,7 @@ transect <- 'line' factor.model <- FALSE dat <- simMsDS(J.x = J.x, J.y = J.y, n.bins = n.bins, bin.width = bin.width, - n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, + n.sp = n.sp, beta = beta, alpha = alpha, det.func = det.func, kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = sp, cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, nu = nu, family = family, offset = offset, transect = transect, factor.model = factor.model) diff --git a/man/simMsNMix.Rd b/man/simMsNMix.Rd index 756d351..9707830 100644 --- a/man/simMsNMix.Rd +++ b/man/simMsNMix.Rd @@ -94,8 +94,8 @@ alpha.mean <- c(0.5, 0.2) tau.sq.alpha <- c(0.2, 0.3) p.det <- length(alpha.mean) mu.RE <- list(levels = c(10, 12), - sigma.sq.mu = c(1.5, 0.3), - beta.indx = list(1, 2)) + sigma.sq.mu = c(1.5, 0.3), + beta.indx = list(1, 2)) p.RE <- list(levels = c(15, 10), sigma.sq.p = c(0.8, 0.5), alpha.indx = list(1, 2)) diff --git a/man/simNMix.Rd b/man/simNMix.Rd index 6668a6d..0e8b54e 100644 --- a/man/simNMix.Rd +++ b/man/simNMix.Rd @@ -8,8 +8,8 @@ \usage{ simNMix(J.x, J.y, n.rep, n.rep.max, beta, alpha, kappa, mu.RE = list(), - p.RE = list(), offset = 1, sp = FALSE, cov.model, sigma.sq, phi, nu, - family = 'Poisson', ...) + p.RE = list(), offset = 1, sp = FALSE, cov.model, sigma.sq, phi, nu, + family = 'Poisson', ...) } \arguments{ @@ -50,7 +50,7 @@ simNMix(J.x, J.y, n.rep, n.rep.max, beta, alpha, kappa, mu.RE = list(), } \author{ - Jeffrey W. Doser \email{doserjef@msu.edu}, \cr + Jeffrey W. Doser \email{doserjef@msu.edu} \cr } \value{ @@ -80,12 +80,10 @@ alpha <- c(0.7, 0.4) kappa <- 0.5 phi <- 3 / .6 sigma.sq <- 2 -mu.RE <- list(levels = 10, - sigma.sq.mu = 1.2) -p.RE <- list(levels = 15, - sigma.sq.p = 0.8) +mu.RE <- list(levels = 10, sigma.sq.mu = 1.2) +p.RE <- list(levels = 15, sigma.sq.p = 0.8) dat <- simNMix(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, alpha = alpha, - kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = TRUE, - cov.model = 'spherical', sigma.sq = sigma.sq, phi = phi, + kappa = kappa, mu.RE = mu.RE, p.RE = p.RE, sp = TRUE, + cov.model = 'spherical', sigma.sq = sigma.sq, phi = phi, family = 'NB') } diff --git a/man/spAbund.Rd b/man/spAbund.Rd index 5801ab0..c50b199 100644 --- a/man/spAbund.Rd +++ b/man/spAbund.Rd @@ -1,19 +1,19 @@ \name{spAbund} \alias{spAbund} -\title{Function for Fitting Single-Species Spatial Abundance GLMs} +\title{Function for Fitting Univariate Spatial Abundance GLMs} \description{ - The function \code{spAbund} fits single-species spatial abundance GLMs. + The function \code{spAbund} fits univariate spatial abundance GLMs. } \usage{ spAbund(formula, data, inits, priors, tuning, - cov.model = 'exponential', NNGP = TRUE, - n.neighbors = 15, search.type = 'cb', - n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', - n.omp.threads = 1, verbose = TRUE, n.report = 100, + cov.model = 'exponential', NNGP = TRUE, + n.neighbors = 15, search.type = 'cb', + n.batch, batch.length, accept.rate = 0.43, family = 'Poisson', + n.omp.threads = 1, verbose = TRUE, n.report = 100, n.burn = round(.10 * n.batch * batch.length), n.thin = 1, - n.chains = 1, save.fitted = TRUE, ...) + n.chains = 1, save.fitted = TRUE, ...) } \arguments{ @@ -24,15 +24,18 @@ spAbund(formula, data, inits, priors, tuning, and slopes are allowed using lme4 syntax (Bates et al. 2015).} \item{data}{a list containing data necessary for model fitting. - Valid tags are \code{y} and \code{covs}. \code{y} + Valid tags are \code{y}, \code{covs}, \code{z}, \code{coords}, and \code{offset}. \code{y} is a vector, matrix, or data frame of the observed count values. If a vector, the values represent the observed counts at each site. If multiple replicate observations are obtained at the sites (e.g., sub-samples, repeated sampling over multiple seasons), \code{y} can be specified as a matrix or data frame with first dimension equal to the number of sites (\eqn{J}{J}) and second dimension equal to the maximum number of - replicates at a given site. \code{covs} is a list - containing the variables used in the model. Each list element is a different + replicates at a given site. \code{covs} is either a data frame or list + containing the variables used in the model. When only fitting a model with site-level + data, \code{covs} can be specified as a data frame, with each row corresponding to + site and each column corresponding to a variable. When multiple abundance values + are available at a site, \code{covs} is specified as a list, where each list element is a different covariate, which can be site-level or observation-level. Site-level covariates are specified as a vector of length \eqn{J}{J}, while observation-level covariates are specified as a matrix or data frame with the number of rows equal to \eqn{J}{J} @@ -126,13 +129,13 @@ spAbund(formula, data, inits, priors, tuning, might produce different, but equally valid, neighbor sets, e.g., if data are on a grid. } - \item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC + \item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{batch.length}{the length of each MCMC batch in each chain to run for the Adaptive + \item{batch.length}{the length of each MCMC batch in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{accept.rate}{target acceptance rate for Adaptive MCMC. Default is + \item{accept.rate}{target acceptance rate for adaptive MCMC. Default is 0.43. See Roberts and Rosenthal (2009) for details.} \item{family}{the distribution to use for the latent abundance process. Currently @@ -270,8 +273,8 @@ sp <- TRUE cov.model <- 'exponential' family <- 'NB' dat <- simAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, - kappa = kappa, mu.RE = mu.RE, sp = sp, phi = phi, - sigma.sq = sigma.sq, cov.model = cov.model, family = 'NB') + kappa = kappa, mu.RE = mu.RE, sp = sp, phi = phi, + sigma.sq = sigma.sq, cov.model = cov.model, family = 'NB') y <- dat$y X <- dat$X @@ -279,26 +282,21 @@ X.re <- dat$X.re coords <- dat$coords covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.cov.2 = X[, , 3], - abund.cov.3 = X[, , 4], - abund.factor.1 = X.re[, , 1], - abund.factor.2 = X.re[, , 2]) + abund.cov.1 = X[, , 2], + abund.cov.2 = X[, , 3], + abund.cov.3 = X[, , 4], + abund.factor.1 = X.re[, , 1], + abund.factor.2 = X.re[, , 2]) -data.list <- list(y = y, - covs = covs, - coords = coords) +data.list <- list(y = y, covs = covs, coords = coords) # Priors prior.list <- list(beta.normal = list(mean = 0, var = 100), - phi.unif = c(3 / 1, 3 / .1), - sigma.sq.ig = c(2, 1), + phi.unif = c(3 / 1, 3 / .1), + sigma.sq.ig = c(2, 1), kappa.unif = c(0.001, 10)) # Starting values -inits.list <- list(beta = beta, - kappa = kappa, - sigma.sq = sigma.sq, - phi = phi) +inits.list <- list(beta = beta, kappa = kappa, sigma.sq = sigma.sq, phi = phi) tuning <- list(phi = 0.3, kappa = 0.05, beta = 0.1, beta.star = 0.1, w = 0.1) n.batch <- 4 @@ -308,23 +306,23 @@ n.thin <- 1 n.chains <- 1 out <- spAbund(formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3 + - (1 | abund.factor.1) + (abund.cov.1 | abund.factor.2), - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - tuning = tuning, - priors = prior.list, - NNGP = TRUE, - cov.model = 'exponential', - search.type = 'cb', - n.neighbors = 5, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) + (1 | abund.factor.1) + (abund.cov.1 | abund.factor.2), + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + tuning = tuning, + priors = prior.list, + NNGP = TRUE, + cov.model = 'exponential', + search.type = 'cb', + n.neighbors = 5, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) } diff --git a/man/spDS.Rd b/man/spDS.Rd index cfd4d44..9f34bb9 100644 --- a/man/spDS.Rd +++ b/man/spDS.Rd @@ -34,7 +34,7 @@ spDS(abund.formula, det.formula, data, inits, priors, tuning, with first dimension equal to the number of sites (\eqn{J}{J}) and second dimension equal to the number of distance bins. \code{covs} is a matrix or data frame - containing the variables used in the the abundance and/or the detection + containing the variables used in the abundance and/or the detection portion of the model, with \eqn{J}{J} rows for each column (variable). \code{dist.breaks} is a vector of distances that denote the breakpoints of the distance bands. \code{dist.breaks} should @@ -87,7 +87,7 @@ spDS(abund.formula, det.formula, data, inits, priors, tuning, Gamma distribution. The hyperparameters of the inverse-Gamma distribution are passed as a list of length two with first and second elements corresponding to the shape and scale parameters, respectively, which are each specified as - vectors of length equal to the number of random intercepts or of length one + vectors of length equal to the number of random intercepts/slopes or of length one if priors are the same for all random effect variances.} \item{cov.model}{a quoted keyword that specifies the covariance @@ -302,12 +302,12 @@ coords <- dat$coords covs <- cbind(X, X.p) colnames(covs) <- c('int.abund', 'abund.cov.1', 'abund.cov.2', 'abund.cov.3', - 'int.det', 'det.cov.1') + 'int.det', 'det.cov.1') data.list <- list(y = y, covs = covs, dist.breaks = dist.breaks, - coords = coords, + coords = coords, offset = offset) # Priors @@ -336,8 +336,8 @@ out <- spDS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3, family = 'NB', det.func = 'halfnormal', transect = 'point', - cov.model = 'exponential', - NNGP = TRUE, + cov.model = 'exponential', + NNGP = TRUE, n.neighbors = 5, tuning = tuning, priors = prior.list, diff --git a/man/summary.abund.Rd b/man/summary.abund.Rd index 53eb744..d90bf6d 100644 --- a/man/summary.abund.Rd +++ b/man/summary.abund.Rd @@ -7,7 +7,7 @@ \title{Methods for abund Object} \description{ - Methods for extracting information from fitted abundance GLMs (\code{abund}). + Methods for extracting information from fitted univariate GLMMs (\code{abund}). } \usage{ diff --git a/man/summary.lfMsAbund.Rd b/man/summary.lfMsAbund.Rd index 4acd256..e6fc8f8 100644 --- a/man/summary.lfMsAbund.Rd +++ b/man/summary.lfMsAbund.Rd @@ -7,7 +7,7 @@ \title{Methods for lfMsAbund Object} \description{ - Methods for extracting information from fitted latent factor multi-species abundance GLMs (\code{lfMsAbund}). + Methods for extracting information from fitted latent factor multivariate abundance GLMMs (\code{lfMsAbund}). } \usage{ diff --git a/man/summary.msAbund.Rd b/man/summary.msAbund.Rd index cde23af..7e05e41 100644 --- a/man/summary.msAbund.Rd +++ b/man/summary.msAbund.Rd @@ -7,7 +7,7 @@ \title{Methods for msAbund Object} \description{ - Methods for extracting information from fitted multi-species abundance GLMs (\code{msAbund}). + Methods for extracting information from fitted multivariate abundance GLMMs (\code{msAbund}). } \usage{ diff --git a/man/summary.sfMsAbund.Rd b/man/summary.sfMsAbund.Rd index 4dbd5f4..228d220 100644 --- a/man/summary.sfMsAbund.Rd +++ b/man/summary.sfMsAbund.Rd @@ -7,7 +7,7 @@ \title{Methods for sfMsAbund Object} \description{ - Methods for extracting information from fitted spatial factor multi-species abundance GLMs (\code{sfMsAbund}). + Methods for extracting information from fitted spatial factor multivariate abundance GLMMs (\code{sfMsAbund}). } \usage{ diff --git a/man/summary.spAbund.Rd b/man/summary.spAbund.Rd index 50bdbe7..449d3f8 100644 --- a/man/summary.spAbund.Rd +++ b/man/summary.spAbund.Rd @@ -7,7 +7,7 @@ \title{Methods for spAbund Object} \description{ - Methods for extracting information from fitted spatially-explicit abundance GLMs (\code{spAbund}). + Methods for extracting information from fitted univariate spatially-explicit GLMMs (\code{spAbund}). } \usage{ diff --git a/man/summary.spDS.Rd b/man/summary.spDS.Rd index 67b1646..6833c2e 100644 --- a/man/summary.spDS.Rd +++ b/man/summary.spDS.Rd @@ -7,7 +7,7 @@ \title{Methods for spDS Object} \description{ - Methods for extracting information from fitted single-species hierarchiacl distance sampling (\code{spDS}) models. + Methods for extracting information from fitted single-species spatial hierarchiacl distance sampling (\code{spDS}) models. } \usage{ diff --git a/man/summary.svcAbund.Rd b/man/summary.svcAbund.Rd index 94f5991..0308c3e 100644 --- a/man/summary.svcAbund.Rd +++ b/man/summary.svcAbund.Rd @@ -7,7 +7,7 @@ \title{Methods for svcAbund Object} \description{ - Methods for extracting information from fitted spatially-varying coefficient abundance GLMs (\code{svcAbund}). + Methods for extracting information from fitted univariate spatially-varying coefficient GLMMs (\code{svcAbund}). } \usage{ diff --git a/man/summary.svcMsAbund.Rd b/man/summary.svcMsAbund.Rd index 1ea410e..7eb443c 100644 --- a/man/summary.svcMsAbund.Rd +++ b/man/summary.svcMsAbund.Rd @@ -7,7 +7,7 @@ \title{Methods for svcMsAbund Object} \description{ - Methods for extracting information from fitted spatial factor multi-species abundance GLMs (\code{svcMsAbund}). + Methods for extracting information from fitted multivariate spatially-varying coefficient abundance GLMMs (\code{svcMsAbund}). } \usage{ diff --git a/man/svcAbund.Rd b/man/svcAbund.Rd index 07a1fbb..85509fb 100644 --- a/man/svcAbund.Rd +++ b/man/svcAbund.Rd @@ -1,19 +1,19 @@ \name{svcAbund} \alias{svcAbund} -\title{Function for Fitting Single-Species Spatialy-Varying Coefficient GLMs} +\title{Function for Fitting Univariate Spatialy-Varying Coefficient GLMMs} \description{ - The function \code{svcAbund} fits univariate spatially-varying coefficient GLMs. + The function \code{svcAbund} fits univariate spatially-varying coefficient GLMMs. } \usage{ svcAbund(formula, data, inits, priors, tuning, - svc.cols = 1, cov.model = 'exponential', NNGP = TRUE, - n.neighbors = 15, search.type = 'cb', n.batch, - batch.length, accept.rate = 0.43, family = 'Gaussian', - n.omp.threads = 1, verbose = TRUE, n.report = 100, - n.burn = round(.10 * n.batch * batch.length), n.thin = 1, - n.chains = 1, ...) + svc.cols = 1, cov.model = 'exponential', NNGP = TRUE, + n.neighbors = 15, search.type = 'cb', n.batch, + batch.length, accept.rate = 0.43, family = 'Gaussian', + n.omp.threads = 1, verbose = TRUE, n.report = 100, + n.burn = round(.10 * n.batch * batch.length), n.thin = 1, + n.chains = 1, ...) } \arguments{ @@ -31,9 +31,9 @@ svcAbund(formula, data, inits, priors, tuning, frame of covariates used in the model, where each column (or list element) represents a different covariate. \code{coords} is a \eqn{J \times 2}{J x 2} matrix of the observation coordinates. Note that \code{svcAbundance} assumes coordinates are specified - in a projected coordinate system. \code{z} is used for fitting a two-stage - Gaussian hurdle model. It is a vector where each value indicates the binary - component of the hurdle model. In the context of abundance models, this can be + in a projected coordinate system. \code{z} is used for fitting a zero-inflated + Gaussian model. It is a vector where each value indicates the binary + component of the model. In the context of abundance models, this can be thought of as the component of the model that indicates whether the species is present at each location, and then the supplied values in \code{y} are the observed abundance values at those locations where \code{z = 1}.} @@ -65,11 +65,16 @@ svcAbund(formula, data, inits, priors, tuning, inverse-Gamma distribution. The spatial decay \code{phi} and spatial smoothness \code{nu}, parameters are assumed to follow Uniform distributions. The hyperparameters of the inverse-Gamma for \code{sigma.sq} - and \code{tau.sq} are passed as a vector of length two, with the first and second + is passed as a list of length two with the first and second elements corresponding + to the shape and scale parameters of the inverse-Gamma distribution either for + each spatially-varying coefficient, or a single value if assumign the same values + for all spatially-varying coefficients. The hyperparameters of the inverse-Gamma for + \code{tau.sq} is passed as a vector of length two, with the first and second elements corresponding to the \emph{shape} and \emph{scale}, respectively. - The hyperparameters of the Uniform are also passed as a vector of + The hyperparameters of the Uniform are also passed as a list of length two with the first and second elements corresponding to - the lower and upper support, respectively. \code{sigma.sq.mu} + the lower and upper support, respectively, for each SVC or a single value + if giving the same prior for each SVC. \code{sigma.sq.mu} are the random effect variances for any random effects, and are assumed to follow an inverse-Gamma distribution. The hyperparameters of the inverse-Gamma distribution are passed as a list of length two with the @@ -84,7 +89,7 @@ svcAbund(formula, data, inits, priors, tuning, in the model formula (with 1 being the intercept if specified), or it can be specified as a character vector with names corresponding to variable names in \code{occ.covs} (for the intercept, use '(Intercept)'). \code{svc.cols} - default argument of 1 results in a spatial occupancy model analogous to + default argument of 1 results in a univariate spatial GLMM analogous to \code{spAbund} (assuming an intercept is included in the model).} \item{cov.model}{a quoted keyword that specifies the covariance @@ -119,13 +124,13 @@ svcAbund(formula, data, inits, priors, tuning, might produce different, but equally valid, neighbor sets, e.g., if data are on a grid. } - \item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC + \item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{batch.length}{the length of each MCMC batch in each chain to run for the Adaptive + \item{batch.length}{the length of each MCMC batch in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{accept.rate}{target acceptance rate for Adaptive MCMC. Default is + \item{accept.rate}{target acceptance rate for adaptive MCMC. Default is 0.43. See Roberts and Rosenthal (2009) for details.} \item{family}{the distribution to use for abundance. Currently, spatially-varying @@ -191,8 +196,7 @@ svcAbund(formula, data, inits, priors, tuning, for the abundance regression coefficients.} \item{tau.sq.samples}{a \code{coda} object of posterior samples - for the residual variance parameter. Only included when - \code{family = 'Gaussian'} or \code{family = 'zi-Gaussian'}.} + for the residual variance parameter.} \item{y.rep.samples}{a \code{coda} object of posterior samples for the abundance replicate (fitted) values with dimensions @@ -205,8 +209,9 @@ svcAbund(formula, data, inits, priors, tuning, \item{theta.samples}{a \code{coda} object of posterior samples for spatial covariance parameters.} - \item{w.samples}{a \code{coda} object of posterior samples - for latent spatial random effects.} + \item{w.samples}{a three-dimensional array of posterior samples + for the spatially-varying coefficients with dimensions corresponding + to MCMC sample, SVC, and site.} \item{sigma.sq.mu.samples}{a \code{coda} object of posterior samples for variances of random effects included in the model. @@ -238,12 +243,12 @@ set.seed(1000) J.x <- 10 J.y <- 10 J <- J.x * J.y -# Occurrence -------------------------- +# Abundance --------------------------- beta <- c(5, 0.5, -0.2, 0.75) p <- length(beta) mu.RE <- list() mu.RE <- list(levels = c(35, 40), - sigma.sq.mu = c(0.7, 1.5), + sigma.sq.mu = c(0.7, 1.5), beta.indx = list(1, 1)) # Spatial parameters ------------------ sp <- TRUE @@ -257,9 +262,9 @@ z <- rbinom(J, 1, 0.5) # Get all the data dat <- simAbund(J.x = J.x, J.y = J.y, beta = beta, tau.sq = tau.sq, - mu.RE = mu.RE, sp = sp, svc.cols = svc.cols, family = 'zi-Gaussian', - cov.model = cov.model, sigma.sq = sigma.sq, phi = phi, - z = z) + mu.RE = mu.RE, sp = sp, svc.cols = svc.cols, + family = 'zi-Gaussian', cov.model = cov.model, + sigma.sq = sigma.sq, phi = phi, z = z) # Get data in format for spAbundance -------------------------------------- y <- dat$y X <- dat$X @@ -271,19 +276,17 @@ covs <- cbind(X, X.re) colnames(covs) <- c('int', 'cov.1', 'cov.2', 'cov.3', 'factor.1', 'factor.2') # Data list bundle -data.list <- list(y = y, - covs = covs, - coords = coords, - z = z) +data.list <- list(y = y, covs = covs, coords = coords, z = z) # Priors prior.list <- list(beta.normal = list(mean = 0, var = 1000), - sigma.sq.ig = list(a = 2, b = 1), tau.sq = c(2, 1), - sigma.sq.mu.ig = list(a = 2, b = 1), + sigma.sq.ig = list(a = 2, b = 1), tau.sq = c(2, 1), + sigma.sq.mu.ig = list(a = 2, b = 1), phi.unif = list(a = 3 / 1, b = 3 / 0.1)) # Starting values inits.list <- list(beta = 0, alpha = 0, - sigma.sq = 1, phi = 3 / 0.5, tau.sq = 2, sigma.sq.mu = 0.5) + sigma.sq = 1, phi = 3 / 0.5, + tau.sq = 2, sigma.sq.mu = 0.5) # Tuning tuning.list <- list(phi = 1) @@ -292,7 +295,8 @@ batch.length <- 25 n.burn <- 100 n.thin <- 1 -out <- svcAbund(formula = ~ cov.1 + cov.2 + cov.3 + (1 | factor.1) + (1 | factor.2), +out <- svcAbund(formula = ~ cov.1 + cov.2 + cov.3 + + (1 | factor.1) + (1 | factor.2), svc.cols = c(1, 2), data = data.list, n.batch = n.batch, @@ -300,7 +304,7 @@ out <- svcAbund(formula = ~ cov.1 + cov.2 + cov.3 + (1 | factor.1) + (1 | factor inits = inits.list, priors = prior.list, accept.rate = 0.43, - family = 'zi-Gaussian', + family = 'zi-Gaussian', cov.model = "exponential", tuning = tuning.list, n.omp.threads = 1, diff --git a/man/svcMsAbund.Rd b/man/svcMsAbund.Rd index a88c666..8e92ae8 100644 --- a/man/svcMsAbund.Rd +++ b/man/svcMsAbund.Rd @@ -1,19 +1,19 @@ \name{svcMsAbund} \alias{svcMsAbund} -\title{Function for Fitting Spatially-Varying Coefficient Multi-Species Abundance GLMs} +\title{Function for Fitting Spatially-Varying Coefficient Multivariate Abundance GLMMs} \description{ - The function \code{svcMsAbund} fits multi-species spatially-varying coefficient GLMs with species correlations (i.e., a spatially-explicit abundace-based joint species distribution model). We use a spatial factor modeling approach. Currently, models are implemented using a Nearest Neighbor Gaussian Process. Future development may allow for running the models using full Gaussian Processes. + The function \code{svcMsAbund} fits multivariate spatially-varying coefficient GLMs with species correlations (i.e., a spatially-explicit abundace-based joint species distribution model). We use a spatial factor modeling approach. Models are implemented using a Nearest Neighbor Gaussian Process. } \usage{ svcMsAbund(formula, data, inits, priors, tuning, - svc.cols = 1, cov.model = 'exponential', NNGP = TRUE, - n.neighbors = 15, search.type = 'cb', n.factors, - n.batch, batch.length, accept.rate = 0.43, family = 'Gaussian', - n.omp.threads = 1, verbose = TRUE, n.report = 100, - n.burn = round(.10 * n.batch * batch.length), n.thin = 1, n.chains = 1, - ...) + svc.cols = 1, cov.model = 'exponential', NNGP = TRUE, + n.neighbors = 15, search.type = 'cb', n.factors, + n.batch, batch.length, accept.rate = 0.43, family = 'Gaussian', + n.omp.threads = 1, verbose = TRUE, n.report = 100, + n.burn = round(.10 * n.batch * batch.length), n.thin = 1, n.chains = 1, + ...) } \arguments{ @@ -24,19 +24,11 @@ svcMsAbund(formula, data, inits, priors, tuning, and slopes are allowed using lme4 syntax (Bates et al. 2015).} \item{data}{a list containing data necessary for model fitting. - Valid tags are \code{y}, \code{covs}, \code{coords}. - \code{y} is a two or three-dimensional array of observed count data. The - first dimension of the array is equal to the - number of species and the second dimension is equal to the number of sites. If - specified as a three-dimensional array, the third dimension corresponds to - replicate observations at each site (e.g., sub-samples, repeated sampling - over multiple seasons). \code{covs} is a list - containing the variables used in the model. Each list element is a different - covariate, which can be site-level or observation-level. Site-level covariates - are specified as a vector of length \eqn{J}{J}, while observation-level covariates - are specified as a matrix or data frame with the number of rows equal to \eqn{J}{J} - and number of columns equal to the maximum number of replicate observations at a - given site. \code{coords} is a + Valid tags are \code{y}, \code{covs}, \code{coords}, and \code{z}. + \code{y} is a matrix with sites corresponding to species and columns + corresponding to sites. \code{covs} is a list, matrix, or data + frame of covariates used in the model, where each column (or list element) + represents a different covariate. \code{coords} is a \eqn{J \times 2}{J x 2} matrix of the observation coordinates. Note that \code{spAbundance} assumes coordinates are specified in a projected coordinate system. For zero-inflated Gaussian models, the tag \code{z} is used to specify the @@ -44,10 +36,10 @@ svcMsAbund(formula, data, inits, priors, tuning, \item{inits}{a list with each tag corresponding to a parameter name. Valid tags are \code{beta.comm}, \code{beta}, - \code{tau.sq.beta}, \code{sigma.sq.mu}, \code{kappa}, + \code{tau.sq.beta}, \code{sigma.sq.mu}, \code{phi}, \code{lambda}, \code{nu}, and \code{tau.sq}. \code{nu} is only specified if - \code{cov.model = "matern"}, \code{kappa} is only specified if \code{family = 'NB'}, - \code{tau.sq} is only specified for Gaussian and zero-inflated Gaussian models, + \code{cov.model = "matern"}, \code{tau.sq} is only + specified for Gaussian and zero-inflated Gaussian models, and \code{sigma.sq.mu} is only specified if random effects are included in \code{formula}. The value portion of each tag is the parameter's initial value. See \code{priors} description for definition @@ -57,7 +49,7 @@ svcMsAbund(formula, data, inits, priors, tuning, \item{priors}{a list with each tag corresponding to a parameter name. Valid tags are \code{beta.comm.normal}, \code{tau.sq.beta.ig}, \code{sigma.sq.mu}, - \code{kappa.unif}, \code{phi.unif}, \code{nu.unif}, and \code{tau.sq.ig}. + \code{phi.unif}, \code{nu.unif}, and \code{tau.sq.ig}. Community-level (\code{beta.comm}) regression coefficients are assumed to follow a normal distribution. The hyperparameters of the normal distribution are passed as a list of length two with the first and second elements @@ -73,16 +65,27 @@ svcMsAbund(formula, data, inits, priors, tuning, which are each specified as vectors of length equal to the number of coefficients to be estimated or a single value if priors are the same for all parameters. If not specified, prior shape and scale - parameters are set to 0.1. The spatial factor model fits \code{n.factors} independent + parameters are set to 0.1. If desired, the species-specific regression coefficients + (\code{beta}) can also be estimated indepdendently by specifying the + tag \code{independent.betas = TRUE}. If specified, this will not estimate species-specific + coefficients as random effects from a common-community-level distribution, and rather + the values of \code{beta.comm} and \code{tau.sq.beta} will be fixed at the + specified initial values. This is equivalent to specifying a Gaussian, independent + prior for each of the species-specific effects. + The spatial factor model fits \code{n.factors} independent spatial processes. The spatial decay \code{phi} and smoothness \code{nu} parameters - for each latent factor are assumed to follow Uniform distributions. + for each latent factor and spatially-varying coefficient + are assumed to follow Uniform distributions. The hyperparameters of the Uniform are passed as a list with two elements, - with both elements being vectors of length \code{n.factors} corresponding to the lower and + with both elements being vectors of length equal to the number of spatial factors + times the number of spatially-varying coefficients corresponding to the lower and upper support, respectively, or as a single value if the same value is assigned - for all factors. The priors for the factor loadings matrix \code{lambda} are fixed + for all factors and spatially-varying coefficients. + The priors for the factor loadings matrix \code{lambda} are fixed following the standard spatial factor model to ensure parameter identifiability (Christensen and Amemlya 2002). The - upper triangular elements of the \code{n.sp x n.factors} matrix are fixed at 0 and the + upper triangular elements of the \code{n.sp x n.factors} matrix + for each spatially-varying coefficient are fixed at 0 and the diagonal elements are fixed at 1. The lower triangular elements are assigned a standard normal prior (i.e., mean 0 and variance 1). \code{sigma.sq.mu} are the random @@ -91,13 +94,7 @@ svcMsAbund(formula, data, inits, priors, tuning, are passed as a list of length two with first and second elements corresponding to the shape and scale parameters, respectively, which are each specified as vectors of length equal to the number of random intercepts or of length one - if priors are the same for all random effect variances. \code{kappa} is the - negative binomial dispersion parameter for each species and is assumed to - follow a uniform distribution. The hyperparameters of the uniform distribution - are passed as a list of length two with first and second elements corresponding to the - lower and upper bounds of the uniform distribution, respectively, which are each - specified as vectors of length equal to the number of species or of length one - if priors are the same for all species-specific dispersion parameters. \code{tau.sq} is the + if priors are the same for all random effect variances. \code{tau.sq} is the species-specific residual variance for Gaussian (or zero-inflated Gaussian) models, and it is assigned an inverse-Gamma prior. The hyperparameters of the inverse-Gamma are passed as a list of length two, with the first and second element corresponding to the shape and @@ -105,10 +102,8 @@ svcMsAbund(formula, data, inits, priors, tuning, equal to the number of species or a single value if priors are the same for all species.} \item{tuning}{a single numeric value representing the initial variance of the - adaptive sampler for \code{beta}, \code{alpha}, \code{beta.star} (the abundance - random effect values), \code{kappa}, \code{phi}, \code{lambda}. - See Roberts and Rosenthal (2009) for details. Note that only \code{phi} and \code{nu} - are tuned for Gaussian or zero-inflated Gaussian models.} + adaptive sampler for \code{phi} and \code{nu}. + See Roberts and Rosenthal (2009) for details.} \item{svc.cols}{a vector indicating the variables whose effects will be estimated as spatially-varying coefficients. \code{svc.cols} can be an @@ -149,24 +144,26 @@ svcMsAbund(formula, data, inits, priors, tuning, might produce different, but equally valid, neighbor sets, e.g., if data are on a grid. } - \item{n.factors}{the number of factors to use in the spatial factor model approach. + \item{n.factors}{the number of factors to use in the spatial factor + model approach for each spatially-varying coefficient. Typically, the number of factors is set to be small (e.g., 4-5) relative to the total number of species in the community, which will lead to substantial decreases in computation time. However, the value can be anywhere between 1 and the number of species in the modeled community.} - \item{n.batch}{the number of MCMC batches in each chain to run for the Adaptive MCMC + \item{n.batch}{the number of MCMC batches in each chain to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{batch.length}{the length of each MCMC batch to run for the Adaptive + \item{batch.length}{the length of each MCMC batch to run for the adaptive MCMC sampler. See Roberts and Rosenthal (2009) for details.} - \item{accept.rate}{target acceptance rate for Adaptive MCMC. Defaul is + \item{accept.rate}{target acceptance rate for adaptive MCMC. Default is 0.43. See Roberts and Rosenthal (2009) for details.} - \item{family}{the distribution to use for the abundance. Currently - supports \code{'NB'} (negative binomial) and \code{'Poisson'}.} - + \item{family}{the distribution to use for abundance. Currently, spatially-varying + coefficient models are available for \code{family = 'Gaussian'} and + \code{family = 'zi-Gaussian'}.} + \item{n.omp.threads}{a positive integer indicating the number of threads to use for SMP parallel processing. The package must be compiled for OpenMP support. For most Intel-based machines, we @@ -235,19 +232,15 @@ svcMsAbund(formula, data, inits, priors, tuning, \item{beta.samples}{a \code{coda} object of posterior samples for the species level abundance regression coefficients.} - \item{kappa.samples}{a \code{coda} object of posterior samples - for the species level abundance dispersion parameters. Only included - when \code{family = 'NB'}.} - \item{tau.sq.samples}{a \code{coda} object of posterior samples - for the Gaussian residual variance parameter. Only included when - \code{family = 'Gaussian'} or \code{family = 'zi-Gaussian'}.} + for the Gaussian residual variance parameter.} \item{theta.samples}{a \code{coda} object of posterior samples for the spatial correlation parameters.} \item{lambda.samples}{a \code{coda} object of posterior samples - for the latent spatial factor loadings.} + for the latent spatial factor loadings for each spatially-varying + coefficient.} \item{y.rep.samples}{a three or four-dimensional array of posterior samples for the fitted (replicate) values for each species with dimensions corresponding @@ -257,8 +250,10 @@ svcMsAbund(formula, data, inits, priors, tuning, the expected abundance values for each species with dimensions corresponding to MCMC samples, species, site, and replicate.} - \item{w.samples}{a three-dimensional array of posterior samples for - the latent spatial random effects for each latent factor.} + \item{w.samples}{a four-dimensional array of posterior samples for + the latent spatial random effects for each spatial factor within each + spatially-varying coefficient. Dimensions correspond to MCMC sample, + factor, site, and spatially-varying coefficient.} \item{sigma.sq.mu.samples}{a \code{coda} object of posterior samples for variances of random effects included in the abundance portion @@ -285,81 +280,81 @@ svcMsAbund(formula, data, inits, priors, tuning, } \examples{ -set.seed(408) -J.x <- 8 -J.y <- 8 +set.seed(332) +J.x <- 10 +J.y <- 10 J <- J.x * J.y -n.rep <- sample(3, size = J, replace = TRUE) +n.rep <- rep(1, J) n.sp <- 6 # Community-level covariate effects -beta.mean <- c(-2, 0.5) +beta.mean <- c(0, 0.25, 0.6) p.abund <- length(beta.mean) -tau.sq.beta <- c(0.2, 1.2) -# Random effects (two random intercepts) -mu.RE <- list(levels = c(10, 15), - sigma.sq.mu = c(0.43, 0.5)) +tau.sq.beta <- c(0.2, 1.2, 0.4) +# Random effects +mu.RE <- list() # Draw species-level effects from community means. beta <- matrix(NA, nrow = n.sp, ncol = p.abund) for (i in 1:p.abund) { beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i])) } -sp <- TRUE +sp <- TRUE +svc.cols <- c(1, 2) n.factors <- 2 +q.p.svc <- length(svc.cols) * n.factors factor.model <- TRUE -phi <- runif(n.factors, 3/1, 3 / .1) -kappa <- runif(n.sp, 0.1, 1) +phi <- runif(q.p.svc, 3/1, 3 / .4) +tau.sq <- runif(n.sp, 0.1, 5) +cov.model <- 'exponential' +family <- 'Gaussian' -dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, - mu.RE = mu.RE, sp = sp, kappa = kappa, family = 'NB', - factor.model = factor.model, phi = phi, - cov.model = 'exponential', n.factors = n.factors) +dat <- simMsAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, n.sp = n.sp, beta = beta, + mu.RE = mu.RE, sp = sp, tau.sq = tau.sq, family = family, + factor.model = factor.model, phi = phi, + cov.model = cov.model, n.factors = n.factors, + svc.cols = svc.cols) y <- dat$y X <- dat$X -X.re <- dat$X.re coords <- dat$coords -# Package all data into a list -covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.factor.1 = X.re[, , 1], - abund.factor.2 = X.re[, , 2]) -data.list <- list(y = y, - covs = covs, - coords = coords) +covs <- data.frame(abund.cov.1 = X[, 2], + abund.cov.2 = X[, 3]) +data.list <- list(y = y, covs = covs, coords = coords) prior.list <- list(beta.comm.normal = list(mean = 0, var = 100), - kappa.unif = list(a = 0, b = 10), - phi.unif = list(a = 3 / 1, b = 3 / .1), - tau.sq.beta.ig = list(a = .1, b = .1)) -inits.list <- list(beta.comm = 0, - beta = 0, - kappa = 0.5, - tau.sq.beta = 1, - phi = 3 / 0.5) - -# Small -n.batch <- 2 + tau.sq.ig = list(a = 2, b = 2), + phi.unif = list(a = 3 / 1, b = 3 / .1), + tau.sq.beta.ig = list(a = .1, b = .1)) +inits.list <- list(beta.comm = 0, + beta = 0, + tau.sq = 1, + tau.sq.beta = 1, + phi = 3 / 0.5) +tuning.list <- list(phi = 0.5) + +n.batch <- 5 batch.length <- 25 -n.burn <- 20 +n.burn <- 0 n.thin <- 1 n.chains <- 1 -out <- sfMsAbund(formula = ~ abund.cov.1 + (1 | abund.factor.1) + (1 | abund.factor.2), - data = data.list, - n.batch = n.batch, - inits = inits.list, - priors = prior.list, - NNGP = TRUE, - cov.model = 'exponential', - n.neighbors = 5, - n.factors = n.factors, - batch.length = batch.length, - n.omp.threads = 3, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) +out <- svcMsAbund(formula = ~ abund.cov.1 + abund.cov.2, + data = data.list, + n.batch = n.batch, + inits = inits.list, + priors = prior.list, + tuning = tuning.list, + NNGP = TRUE, + svc.cols = c(1, 2), + family = 'Gaussian', + cov.model = 'exponential', + n.neighbors = 5, + n.factors = n.factors, + batch.length = batch.length, + n.omp.threads = 1, + verbose = TRUE, + n.report = 20, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) summary(out) - } diff --git a/man/waicAbund.Rd b/man/waicAbund.Rd index 1130972..590e34a 100644 --- a/man/waicAbund.Rd +++ b/man/waicAbund.Rd @@ -18,7 +18,7 @@ waicAbund(object, N.max, by.species = FALSE, ...) \code{msDS}, \code{lfMsDS}, \code{sfMsDS}.} \item{N.max}{values indicating the upper limit on the latent abundance - values when calculating WAIC for N-mixture models or distance sampling models. + values when calculating WAIC for N-mixture models or hierarchical distance sampling models. For single-species models, this can be a single value or a vector of different values for each site. For multi-species models, this can be a single value, a vector of values for each species, or a species by site matrix @@ -36,11 +36,11 @@ waicAbund(object, N.max, by.species = FALSE, ...) } \note{ - When fitting Gaussian-hurdle models, the WAIC is only calculated for the - values that are past the hurdle. If fitting a first stage model with + When fitting zero-inflated Gaussian models, the WAIC is only calculated for the + non-zero values. If fitting a first stage model with \code{spOccupancy} to the binary portion of the - hurdle model, you can use the \code{spOccupancy::waicOcc} function to - calculate WAIC for the binary portion of the hurdle model. + zero-inflated model, you can use the \code{spOccupancy::waicOcc} function to + calculate WAIC for the binary component. } \references{ @@ -83,33 +83,30 @@ J <- J.x * J.y n.rep <- sample(3, J, replace = TRUE) beta <- c(0, -1.5, 0.3, -0.8) p.abund <- length(beta) -mu.RE <- list(levels = c(30), - sigma.sq.mu = c(1.3)) +mu.RE <- list(levels = c(30), sigma.sq.mu = c(1.3)) kappa <- 0.5 sp <- FALSE family <- 'NB' dat <- simAbund(J.x = J.x, J.y = J.y, n.rep = n.rep, beta = beta, - kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') + kappa = kappa, mu.RE = mu.RE, sp = sp, family = 'NB') y <- dat$y X <- dat$X X.re <- dat$X.re abund.covs <- list(int = X[, , 1], - abund.cov.1 = X[, , 2], - abund.cov.2 = X[, , 3], - abund.cov.3 = X[, , 4], - abund.factor.1 = X.re[, , 1]) + abund.cov.1 = X[, , 2], + abund.cov.2 = X[, , 3], + abund.cov.3 = X[, , 4], + abund.factor.1 = X.re[, , 1]) -data.list <- list(y = y, - covs = abund.covs) +data.list <- list(y = y, covs = abund.covs) # Priors prior.list <- list(beta.normal = list(mean = 0, var = 100), kappa.unif = c(0.001, 10)) # Starting values -inits.list <- list(beta = 0, - kappa = kappa) +inits.list <- list(beta = 0, kappa = kappa) n.batch <- 5 batch.length <- 25 @@ -117,19 +114,20 @@ n.burn <- 0 n.thin <- 1 n.chains <- 1 -out <- abund(formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3 + (1 | abund.factor.1), - data = data.list, - n.batch = n.batch, - batch.length = batch.length, - inits = inits.list, - priors = prior.list, - accept.rate = 0.43, - n.omp.threads = 1, - verbose = TRUE, - n.report = 1, - n.burn = n.burn, - n.thin = n.thin, - n.chains = n.chains) +out <- abund(formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3 + + (1 | abund.factor.1), + data = data.list, + n.batch = n.batch, + batch.length = batch.length, + inits = inits.list, + priors = prior.list, + accept.rate = 0.43, + n.omp.threads = 1, + verbose = TRUE, + n.report = 1, + n.burn = n.burn, + n.thin = n.thin, + n.chains = n.chains) # Calculate WAIC waicAbund(out) diff --git a/vignettes/distanceSampling.html b/vignettes/distanceSampling.html index 4055e2c..543d146 100644 --- a/vignettes/distanceSampling.html +++ b/vignettes/distanceSampling.html @@ -687,7 +687,7 @@

Fitting single-species HDS models with DS()

Thinning Rate: 30 Number of Chains: 3 Total Posterior Samples: 3000 -Run Time (min): 0.602 +Run Time (min): 0.6826 Abundance (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS @@ -1106,7 +1106,7 @@

Fitting single-species spatial HDS models with spDS()

Thinning Rate: 30 Number of Chains: 3 Total Posterior Samples: 3000 -Run Time (min): 3.6429 +Run Time (min): 4.6714 Abundance (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS @@ -1688,7 +1688,7 @@

Fitting multi-species HDS models with msDS()

Thinning Rate: 30 Number of Chains: 3 Total Posterior Samples: 3000 -Run Time (min): 11.764 +Run Time (min): 14.5699 ---------------------------------------- Community Level @@ -1709,6 +1709,7 @@

Fitting multi-species HDS models with msDS()

Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) -2.5335 0.1017 -2.7243 -2.5372 -2.3301 1.0072 1616 scale(wind) -0.0360 0.0563 -0.1461 -0.0364 0.0778 1.0002 3000 + Detection Variances (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS (Intercept) 0.1397 0.0750 0.0533 0.1214 0.3411 1.0049 1124 @@ -2449,7 +2450,7 @@

Fitting latent factor multi-species HDS models with lfMsDS()Fitting spatial factor multi-species HDS models with sfMsDS()Fitting single-species N-mixture models with NMix()

n.omp.threads = 1, verbose = TRUE, n.report = 100, n.burn = round(.10 * n.batch * batch.length), n.thin = 1, n.chains = 1, ...) -

The first two arguments, abund.formula and det.formula, use standard R model syntax to denote the covariates to be included in the abundance and detection portions of the model, respectively. Only the right hand side of the formulas are included. Random intercepts and slopes can be included in both the abundance and detection portions of the single-species N-mixture model using lme4 syntax (Bates et al. 2015). For example, to include a random intercept for different observers in the detection portion of the model, we would include (1 | observer) in det.formula, where observer indicates the specific observer for each data point. The names of variables given in the formulas should correspond to those found in data, which is a list consisting of the following tags: y (count data), abund.covs (abundance covariates), det.covs (detection covariates). y should be stored as a sites x replicate matrix, abund.covs is a matrix or data frame with site-specific covariate values, and det.covs as a list with each list element corresponding to a covariate to include in the detection portion of the model. Covariates on detection can vary by site and/or survey, and so these covariates may be specified as a site by survey matrix for “observation”-level covariates (i.e., covariates that vary by site and survey) or as a one-dimensional vector for site-level covariates.

+

The first two arguments, abund.formula and det.formula, use standard R model syntax to denote the covariates to be included in the abundance and detection portions of the model, respectively. Only the right hand side of the formulas are included. Random intercepts and slopes can be included in both the abundance and detection portions of the single-species N-mixture model using lme4 syntax (Bates et al. 2015). For example, to include a random intercept for different observers in the detection portion of the model, we would include (1 | observer) in det.formula, where observer indicates the specific observer for each data point. The names of variables given in the formulas should correspond to those found in data, which is a list consisting of the following tags: y (count data), abund.covs (abundance covariates), det.covs (detection covariates). y should be stored as a sites x replicate matrix, abund.covs is a matrix or data frame with site-specific covariate values, and det.covs as a list with each list element corresponding to a covariate to include in the detection portion of the model. Covariates on detection can vary by site and/or survey, and so these covariates may be specified as a site by survey matrix for “observation”-level covariates (i.e., covariates that vary by site and survey) or as a one-dimensional vector for site-level covariates. Note the tag offset can also be specified to include an offset in the model.

spAbundance can handle “imbalanced” data sets (i.e., each site may not have the same number of repeat visits). This is the case with our simulated example, where some sites have three repeat visits, while others are only visited once. To accommodate such an imbalanced design, the detection-nondetection data matrix y should have a value of NA in any site/visit combination where there is not a survey. For example, let’s take a look at the first 10 rows of our detection-nondetection matrix in data.one.sp

head(data.one.sp$y, 10)
      [,1] [,2] [,3]
@@ -697,7 +697,7 @@ 
Fitting single-species N-mixture models with NMix()

 Thinning Rate: 10
 Number of Chains: 3
 Total Posterior Samples: 3000
-Run Time (min): 0.3006
+Run Time (min): 0.5099
 
 Abundance (log scale): 
                       Mean     SD    2.5%     50%   97.5%   Rhat  ESS
@@ -729,7 +729,7 @@ 
Fitting single-species N-mixture models with NMix()

 Thinning Rate: 10
 Number of Chains: 3
 Total Posterior Samples: 3000
-Run Time (min): 0.3006
+Run Time (min): 0.5099
 
 Abundance (log scale): 
                       Mean     SD     25%     50%     75%   Rhat  ESS
@@ -778,7 +778,7 @@ 
Posterior predictive checks

 Number of Chains: 3
 Total Posterior Samples: 3000
 
-Bayesian p-value:  0.0523 
+Bayesian p-value:  0.3103 
 Fit statistic:  freeman-tukey

The Bayesian p-value here is very close to 0, indicating inadequate model fit. In particular, a value close to 0 indicates there is more variability in the observed data points than the replicate data points generated by our model. This could be a result of (1) missing sources of variability in true abundance; (2) missing sources of variability in detection probability; or (3) missing sources of variability in both abundance and detection. We will see later in the vignette that our low Bayesian p-value in this case is a result of additional spatial variability in abundance. See the introductory spOccupancy vignette for ways to further explore resulting objects from posterior predictive checks.

@@ -801,19 +801,19 @@

A brief note on generation of fitted values for posterior predictive checks< load(url("https://www.jeffdoser.com/files/misc/nmix-sim/nmix-sim-results.rda")) # Conditional Bayesian p-values round(conditional.bps, 2)

-
 [1] 0.01 0.35 0.00 0.00 0.48 0.17 0.06 0.00 0.00 0.22 0.12 0.32 0.00 0.11 0.09
-[16] 0.18 0.30 0.00 0.27 0.00
+
 [1] 0.09 0.03 0.00 0.00 0.04 0.01 0.03 0.00 0.01 0.07 0.09 0.03 0.00 0.12 0.00
+[16] 0.00 0.02 0.01 0.02 0.08
# Marginal Bayesian p-values
 round(marginal.bps, 2)
-
 [1] 0.00 0.06 0.00 0.00 0.13 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
-[16] 0.00 0.01 0.00 0.15 0.00
+
 [1] 0.00 0.01 0.00 0.00 0.02 0.00 0.00 0.00 0.00 0.00 0.08 0.06 0.00 0.15 0.00
+[16] 0.00 0.03 0.00 0.00 0.06
# How many conditional values are not rejected (e.g., > 0.1 and < 0.9)
 sum(conditional.bps > 0.1 & conditional.bps < 0.9)
-
[1] 10
+
[1] 1
# How many marginal values are not reject (e.g., > 0.1 and < 0.9)
 sum(marginal.bps > 0.1 & marginal.bps < 0.9)
-
[1] 2
-

Here we see the marginal approach fails to reject the model (using cutoff values of 0.1 and 0.9) 2 times out of 20, while the conditional approach fails to reject the model 10 times out of 20, suggesting that in this very small and simplified simulation study the marginal approach is better at detecting unmodeled spatial heterogeneity relative to the conditional approach.

+
[1] 1
+

Here we see the marginal approach fails to reject the model (using cutoff values of 0.1 and 0.9) 1 times out of 20, while the conditional approach fails to reject the model 1 times out of 20, suggesting that in this very small and simplified simulation study the marginal approach is better at detecting unmodeled spatial heterogeneity relative to the conditional approach.

This concept of what we called “marginal” and “conditional” fitted values is not new to the statistical ecology literature. In particular, Conn et al. (2018) discuss this concept (and many more aspects of Goodness-of-Fit checking in hierarchical Bayesian models), and give a nice visual summary of this idea in their Figure 4 in the context of hierarchical spatial models. We believe additional exploration of posterior predictive tests and other associated Goodness-of-Fit tests in hierarchical models like N-mixture models, hierarchical distance sampling models, and occupancy models is an important avenue of future research. Here we provide options to perform posterior predictive checks using both types (with the default being marginal), and encourage users to explore multiple forms of performing posterior predictive checks to increase one’s comfort that the fitted model is an adequate representation of the data-generating process. Our very small and simple simulation study provided us with a little intuition on the performance of the two approaches in detecting unmodeled spatial heterogeneity. In particular, we found that the marginal approach had a higher ability to detect unmodeled spatial heterogeneity in abundance relative to the the conditional approach, although a much more in-depth simulation study is needed to assign any large amount of weight to this finding. This discussion applies to all N-mixture models discussed in this vignette.

@@ -853,7 +853,7 @@

Model selection using WAIC

N.max not specified. Setting upper index of integration of N to 10 plus
 the largest estimated abundance value at each site in object$N.samples
      elpd         pD       WAIC 
--160.36696   18.22889  357.19171
+-160.52411 17.93976 356.92773

From this simple model comparison exercise, we see the negative binomial model outperforms the Poisson model, indicating there is substantial overdispersion in latent abundance that is not captured by the covariate/random effect included in the model.

@@ -879,7 +879,7 @@

Prediction

out.pred <- predict(out.2, X.0, ignore.RE = TRUE) str(out.pred)
List of 3
- $ mu.0.samples: 'mcmc' num [1:3000, 1:100] 0.01403 0.0052 0.0053 0.00168 0.00545 ...
+ $ mu.0.samples: 'mcmc' num [1:3000, 1:100] 0.00768 0.01046 0.00448 0.00424 0.00111 ...
   ..- attr(*, "mcpar")= num [1:3] 1 3000 1
  $ N.0.samples : 'mcmc' num [1:3000, 1:100] 0 0 0 0 0 0 0 0 0 0 ...
   ..- attr(*, "mcpar")= num [1:3] 1 3000 1
@@ -898,7 +898,7 @@ 
Prediction

   geom_line() + 
   theme_bw() + 
   labs(x = 'Covariate', y = 'Expected abundance')
-

+

Here we see a strong positive relationship between the covariate and expected abundance.

@@ -1040,66 +1040,66 @@

Fitting single-species spatial N-mixture models with spNMix()Fitting single-species spatial N-mixture models with spNMix()Fitting single-species spatial N-mixture models with spNMix() 
summary(out.sp)
@@ -1253,31 +1253,31 @@

Fitting single-species spatial N-mixture models with spNMix() +sigma.sq 0.7613 0.3243 0.2926 0.7067 1.5438 1.0162 212 +phi 27.0302 9.4049 8.7160 27.5246 41.1774 1.0289 287

Looking at the model summary we see adequate convergence of all model parameters. The summary() output looks the same as what we saw previously, with the additional section titled “Spatial Covariance”. There we see the estimate of the spatial variance (sigma.sq) and spatial decay (phi) parameters. The spatial variance is near 1, which indicates substantial spatial variability in abundance across space. Interpretation of the spatial variance parameter can follow interpretation of variance parameters in “regular” (i.e., unstructured) random effect variances: when the variance is close to 0, that indicates little support for inclusion of the spatial random effect. When the variance is large, it indicates substantial support for the spatial random effect. What indicates “large” vs. “small” isn’t necessarily straightforward, but here the estimate is clearly above 1, which is a substantially large variance parameter on the log scale. We can also look at the magnitude of the spatial random effect estimates themselves to give an indication as to how much residual spatial autocorrelation there is in the abundance estimates. The posterior samples for the spatial random effects are stored in the w.samples tag of the resulting model fit list. Here we calculate the means of the spatial random effects and plot a histogram of their values.

w.means <- apply(out.sp$w.samples, 2, mean)
 hist(w.means)
-

+

We see a fair amount of large magnitude random effect values, both above and below zero. If there was limited support for spatial autocorrelation, these values would all be very close to zero.

@@ -1295,7 +1295,7 @@

Posterior predictive checks

Number of Chains: 3 Total Posterior Samples: 3000 -Bayesian p-value: 0.46 +Bayesian p-value: 0.4727 Fit statistic: freeman-tukey

Here we see a striking contrast to the Bayesian p-value from the non-spatial Poisson N-mixture model, which was essentially 0. Here, our estimate is almost exactly 0.5, indicating our model is adequately representing the variability in the data with the addition of the spatial random effect in abundance. If you take a look at the manual page for dataNMixSim, you can see how we simulated this data set and why the above result makes sense.

@@ -1307,13 +1307,13 @@

Model selection using WAIC

N.max not specified. Setting upper index of integration of N to 10 plus
 the largest estimated abundance value at each site in object$N.samples
      elpd         pD       WAIC 
--160.36696   18.22889  357.19171
+-160.52411 17.93976 356.92773 
# Poisson spatial model
 waicAbund(out.sp)
N.max not specified. Setting upper index of integration of N to 10 plus
 the largest estimated abundance value at each site in object$N.samples
      elpd         pD       WAIC 
--137.65357   29.05015  333.40744
+-137.44483 29.61766 334.12498

Here we see the spatial Poisson model outperforms the negative binomial non-spatial model. Both of these models can be viewed as a “standard” N-mixture model that account for overdispersion in two different ways: the negative binomial model has a parameter kappa that controls the amount of overdispersion, which is not assumed to have any spatial structure, while the spatial Poisson model can be thought of as a Poisson-lognormal model, where the log-normal random effects are assigned a spatial structure. Here, the latter model does a more adequate job of accounting for overdispersion in the latent abundance values. We could also fit a spatial negative binomial model, which essentially includes two mechanisms to account for overdispersion. We leave this to the interested reader to fit this model and see how it compares for this data set. We will note that spatial negative binomial models can often be quite difficult to successfully fit. This is particularly the case when there is no medium to large range spatial autocorrelation in the abundance values. When this is the case, the spatial random effects will be very difficult to identify from the negative binomial overdispersion parameter, and we may seem some unidentifiability occurring between the negative binomial overdispersion parameter and the parameters controlling the spatial dependence. In such a case, a spatial Poisson model or a negative binomial non-spatial model will usually suffice. Alternatively, we could set an informative prior on the spatial range to only allow the spatial random effects to explain long-range (i.e., broad-scale) spatial autocorrelation (if any exists) and let the negative binomial overdispersion parameter account for short-range (i.e., fine-scale) overdispersion. We could then compare that model to simpler versions of the model using WAIC.

@@ -1328,7 +1328,7 @@

Prediction

out.pred <- predict(out.sp, X.0, ignore.RE = TRUE, include.sp = FALSE)
 str(out.pred)
List of 3
- $ mu.0.samples: 'mcmc' num [1:3000, 1:100] 0.00573 0.0022 0.00471 0.00945 0.00319 ...
+ $ mu.0.samples: 'mcmc' num [1:3000, 1:100] 0.00488 0.00856 0.01105 0.00356 0.0024 ...
   ..- attr(*, "mcpar")= num [1:3] 1 3000 1
  $ N.0.samples : 'mcmc' int [1:3000, 1:100] 0 0 0 0 0 0 0 0 0 0 ...
   ..- attr(*, "mcpar")= num [1:3] 1 3000 1
@@ -1346,7 +1346,7 @@ 
Prediction

   geom_line() +
   theme_bw() +
   labs(x = 'Covariate', y = 'Expected abundance')
-

+

@@ -1452,48 +1452,48 @@

Fitting multi-species N-mixture models with msNMix()

Sampling ... Batch: 200 of 800, 25.00% Species Parameter Acceptance Tuning - 1 beta[1] 44.0 0.27972 - 1 alpha[1] 40.0 0.41729 - 2 beta[1] 44.0 0.16301 - 2 alpha[1] 32.0 0.27418 - 3 beta[1] 44.0 0.09127 - 3 alpha[1] 36.0 0.15047 - 4 beta[1] 32.0 0.33488 - 4 alpha[1] 56.0 0.63510 - 5 beta[1] 40.0 0.12320 - 5 alpha[1] 52.0 0.22003 - 6 beta[1] 28.0 0.16966 - 6 alpha[1] 52.0 0.29701 + 1 beta[1] 36.0 0.25821 + 1 alpha[1] 24.0 0.42572 + 2 beta[1] 52.0 0.16301 + 2 alpha[1] 48.0 0.26875 + 3 beta[1] 40.0 0.09127 + 3 alpha[1] 44.0 0.15351 + 4 beta[1] 44.0 0.34855 + 4 alpha[1] 52.0 0.63510 + 5 beta[1] 48.0 0.13615 + 5 alpha[1] 56.0 0.21141 + 6 beta[1] 60.0 0.18379 + 6 alpha[1] 52.0 0.30914 ------------------------------------------------- Batch: 400 of 800, 50.00% Species Parameter Acceptance Tuning - 1 beta[1] 28.0 0.25821 - 1 alpha[1] 44.0 0.44309 - 2 beta[1] 36.0 0.17308 - 2 alpha[1] 56.0 0.27972 - 3 beta[1] 48.0 0.08946 - 3 alpha[1] 36.0 0.14457 - 4 beta[1] 32.0 0.30302 - 4 alpha[1] 40.0 0.62252 - 5 beta[1] 44.0 0.12320 - 5 alpha[1] 60.0 0.21141 - 6 beta[1] 40.0 0.17308 - 6 alpha[1] 32.0 0.32175 + 1 beta[1] 52.0 0.27972 + 1 alpha[1] 40.0 0.40903 + 2 beta[1] 48.0 0.16966 + 2 alpha[1] 56.0 0.26343 + 3 beta[1] 44.0 0.08095 + 3 alpha[1] 36.0 0.15047 + 4 beta[1] 32.0 0.30914 + 4 alpha[1] 44.0 0.62252 + 5 beta[1] 52.0 0.14749 + 5 alpha[1] 40.0 0.22901 + 6 beta[1] 64.0 0.16966 + 6 alpha[1] 40.0 0.29701 ------------------------------------------------- Batch: 600 of 800, 75.00% Species Parameter Acceptance Tuning - 1 beta[1] 56.0 0.29113 - 1 alpha[1] 52.0 0.38521 - 2 beta[1] 32.0 0.16301 - 2 alpha[1] 52.0 0.25310 - 3 beta[1] 28.0 0.08595 - 3 alpha[1] 52.0 0.16630 - 4 beta[1] 56.0 0.32175 - 4 alpha[1] 52.0 0.58627 - 5 beta[1] 32.0 0.14749 - 5 alpha[1] 36.0 0.21141 - 6 beta[1] 40.0 0.16630 - 6 alpha[1] 40.0 0.30302 + 1 beta[1] 36.0 0.25310 + 1 alpha[1] 40.0 0.40093 + 2 beta[1] 32.0 0.16630 + 2 alpha[1] 24.0 0.30302 + 3 beta[1] 44.0 0.08946 + 3 alpha[1] 32.0 0.15351 + 4 beta[1] 20.0 0.32825 + 4 alpha[1] 44.0 0.61020 + 5 beta[1] 28.0 0.13890 + 5 alpha[1] 16.0 0.22003 + 6 beta[1] 56.0 0.18015 + 6 alpha[1] 48.0 0.31538 ------------------------------------------------- Batch: 800 of 800, 100.00% ---------------------------------------- @@ -1502,48 +1502,48 @@

Fitting multi-species N-mixture models with msNMix()

Sampling ... Batch: 200 of 800, 25.00% Species Parameter Acceptance Tuning - 1 beta[1] 68.0 0.25821 - 1 alpha[1] 40.0 0.36277 - 2 beta[1] 40.0 0.17658 - 2 alpha[1] 56.0 0.28537 - 3 beta[1] 36.0 0.08425 - 3 alpha[1] 44.0 0.14457 - 4 beta[1] 44.0 0.32825 - 4 alpha[1] 36.0 0.64793 - 5 beta[1] 32.0 0.13615 - 5 alpha[1] 40.0 0.21568 - 6 beta[1] 48.0 0.15978 - 6 alpha[1] 36.0 0.32175 + 1 beta[1] 48.0 0.23836 + 1 alpha[1] 64.0 0.40093 + 2 beta[1] 36.0 0.17658 + 2 alpha[1] 44.0 0.27418 + 3 beta[1] 40.0 0.08595 + 3 alpha[1] 40.0 0.14457 + 4 beta[1] 36.0 0.33488 + 4 alpha[1] 64.0 0.58627 + 5 beta[1] 48.0 0.13081 + 5 alpha[1] 52.0 0.21568 + 6 beta[1] 36.0 0.15661 + 6 alpha[1] 40.0 0.27972 ------------------------------------------------- Batch: 400 of 800, 50.00% Species Parameter Acceptance Tuning - 1 beta[1] 48.0 0.28537 - 1 alpha[1] 60.0 0.37758 - 2 beta[1] 44.0 0.15661 - 2 alpha[1] 44.0 0.27972 - 3 beta[1] 52.0 0.08258 - 3 alpha[1] 32.0 0.14171 - 4 beta[1] 52.0 0.32825 - 4 alpha[1] 40.0 0.64793 - 5 beta[1] 44.0 0.14749 - 5 alpha[1] 32.0 0.22448 - 6 beta[1] 28.0 0.16301 - 6 alpha[1] 60.0 0.29701 + 1 beta[1] 40.0 0.25821 + 1 alpha[1] 52.0 0.40903 + 2 beta[1] 28.0 0.16966 + 2 alpha[1] 56.0 0.25821 + 3 beta[1] 40.0 0.08946 + 3 alpha[1] 52.0 0.15351 + 4 beta[1] 20.0 0.32825 + 4 alpha[1] 36.0 0.67437 + 5 beta[1] 32.0 0.14171 + 5 alpha[1] 20.0 0.23364 + 6 beta[1] 48.0 0.16966 + 6 alpha[1] 36.0 0.28537 ------------------------------------------------- Batch: 600 of 800, 75.00% Species Parameter Acceptance Tuning - 1 beta[1] 40.0 0.26875 - 1 alpha[1] 40.0 0.40093 - 2 beta[1] 44.0 0.16630 - 2 alpha[1] 48.0 0.25310 - 3 beta[1] 40.0 0.08425 - 3 alpha[1] 40.0 0.14749 - 4 beta[1] 36.0 0.32825 - 4 alpha[1] 40.0 0.64793 - 5 beta[1] 48.0 0.14171 - 5 alpha[1] 48.0 0.23364 - 6 beta[1] 44.0 0.17308 - 6 alpha[1] 36.0 0.29701 + 1 beta[1] 32.0 0.25821 + 1 alpha[1] 48.0 0.39299 + 2 beta[1] 52.0 0.16630 + 2 alpha[1] 36.0 0.26875 + 3 beta[1] 32.0 0.08095 + 3 alpha[1] 40.0 0.13346 + 4 beta[1] 32.0 0.32825 + 4 alpha[1] 28.0 0.67437 + 5 beta[1] 40.0 0.12320 + 5 alpha[1] 52.0 0.21568 + 6 beta[1] 48.0 0.17308 + 6 alpha[1] 48.0 0.30302 ------------------------------------------------- Batch: 800 of 800, 100.00% ---------------------------------------- @@ -1552,48 +1552,48 @@

Fitting multi-species N-mixture models with msNMix()

Sampling ... Batch: 200 of 800, 25.00% Species Parameter Acceptance Tuning - 1 beta[1] 56.0 0.27418 - 1 alpha[1] 40.0 0.41729 - 2 beta[1] 40.0 0.18015 - 2 alpha[1] 44.0 0.26343 - 3 beta[1] 40.0 0.08258 - 3 alpha[1] 24.0 0.15978 - 4 beta[1] 56.0 0.30302 + 1 beta[1] 48.0 0.23836 + 1 alpha[1] 44.0 0.40093 + 2 beta[1] 48.0 0.16630 + 2 alpha[1] 44.0 0.25310 + 3 beta[1] 52.0 0.08425 + 3 alpha[1] 20.0 0.15661 + 4 beta[1] 52.0 0.34165 4 alpha[1] 36.0 0.67437 - 5 beta[1] 44.0 0.12822 - 5 alpha[1] 40.0 0.19910 - 6 beta[1] 52.0 0.17308 - 6 alpha[1] 56.0 0.29701 + 5 beta[1] 28.0 0.13081 + 5 alpha[1] 48.0 0.22901 + 6 beta[1] 28.0 0.16966 + 6 alpha[1] 40.0 0.29701 ------------------------------------------------- Batch: 400 of 800, 50.00% Species Parameter Acceptance Tuning - 1 beta[1] 40.0 0.27972 - 1 alpha[1] 32.0 0.40903 - 2 beta[1] 36.0 0.16966 - 2 alpha[1] 32.0 0.25821 - 3 beta[1] 40.0 0.08595 - 3 alpha[1] 56.0 0.14749 - 4 beta[1] 52.0 0.36277 - 4 alpha[1] 36.0 0.61020 - 5 beta[1] 52.0 0.12822 - 5 alpha[1] 36.0 0.21141 - 6 beta[1] 44.0 0.17658 - 6 alpha[1] 28.0 0.29701 + 1 beta[1] 52.0 0.25821 + 1 alpha[1] 56.0 0.45205 + 2 beta[1] 24.0 0.16630 + 2 alpha[1] 32.0 0.27972 + 3 beta[1] 28.0 0.08769 + 3 alpha[1] 28.0 0.15351 + 4 beta[1] 28.0 0.31538 + 4 alpha[1] 28.0 0.66102 + 5 beta[1] 44.0 0.13346 + 5 alpha[1] 44.0 0.23836 + 6 beta[1] 28.0 0.16301 + 6 alpha[1] 36.0 0.33488 ------------------------------------------------- Batch: 600 of 800, 75.00% Species Parameter Acceptance Tuning - 1 beta[1] 36.0 0.29113 - 1 alpha[1] 36.0 0.40903 - 2 beta[1] 48.0 0.15978 - 2 alpha[1] 52.0 0.25310 - 3 beta[1] 48.0 0.08595 - 3 alpha[1] 20.0 0.14749 - 4 beta[1] 52.0 0.32825 - 4 alpha[1] 48.0 0.66102 - 5 beta[1] 44.0 0.13081 - 5 alpha[1] 28.0 0.21568 - 6 beta[1] 48.0 0.18379 - 6 alpha[1] 56.0 0.30302 + 1 beta[1] 48.0 0.26343 + 1 alpha[1] 44.0 0.39299 + 2 beta[1] 40.0 0.16630 + 2 alpha[1] 32.0 0.28537 + 3 beta[1] 44.0 0.08258 + 3 alpha[1] 40.0 0.14749 + 4 beta[1] 40.0 0.34165 + 4 alpha[1] 32.0 0.67437 + 5 beta[1] 32.0 0.13081 + 5 alpha[1] 36.0 0.22448 + 6 beta[1] 52.0 0.17308 + 6 alpha[1] 48.0 0.31538 ------------------------------------------------- Batch: 800 of 800, 100.00%

The resulting object out.ms is a list of class msNMix consisting primarily of posterior samples of all community and species-level parameters, as well as some additional objects that are used for summaries, prediction, and model fit evaluation. We can display a nice summary of these results using the summary() function as before. For multi-species objects, when using summary we need to specify the level of parameters we want to summarize. We do this using the argument level, which takes values community, species, or both to print results for community-level parameters, species-level parameters, or all parameters. By default, level is set to both to display both species and community-level parameters.

@@ -1611,74 +1611,75 @@

Fitting multi-species N-mixture models with msNMix()

Thinning Rate: 10 Number of Chains: 3 Total Posterior Samples: 3000 -Run Time (min): 2.4082 +Run Time (min): 3.4674 ---------------------------------------- Community Level ---------------------------------------- Abundance Means (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) -0.6793 0.5408 -1.7682 -0.6712 0.4235 1.0109 3000 -scale(abund.cov.1) 0.3577 0.3767 -0.3879 0.3620 1.1255 1.0002 3000 +(Intercept) -0.6681 0.5421 -1.7512 -0.6696 0.4090 1.0027 1888 +scale(abund.cov.1) 0.3447 0.3866 -0.4206 0.3415 1.1174 1.0024 3000 Abundance Variances (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) 1.7606 2.3011 0.3832 1.2013 6.2662 1.0694 2830 -scale(abund.cov.1) 0.9262 1.1253 0.2094 0.6408 3.6262 1.0084 2403 +(Intercept) 1.7226 2.2140 0.3449 1.1721 6.5096 1.0083 2083 +scale(abund.cov.1) 0.9200 1.1334 0.2100 0.6348 3.1639 1.0400 2714 Abundance Random Effect Variances (log scale): - Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept)-abund.factor.1 0.744 0.1813 0.4495 0.7229 1.1648 1.0019 1492 + Mean SD 2.5% 50% 97.5% Rhat ESS +(Intercept)-abund.factor.1 0.7446 0.1875 0.4549 0.7161 1.1891 1.0121 1841 Detection Means (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) 0.2780 0.3898 -0.5561 0.2866 1.0505 0.9997 2376 -scale(det.cov.1) 0.5824 0.2515 0.0666 0.5924 1.0476 1.0016 2742 -scale(det.cov.2) 0.7042 0.5973 -0.5598 0.6980 1.8996 1.0010 2728 +(Intercept) 0.2743 0.3880 -0.5387 0.2791 1.0514 1.0015 3000 +scale(det.cov.1) 0.5869 0.2453 0.0997 0.5972 1.0389 1.0082 2729 +scale(det.cov.2) 0.6879 0.6038 -0.5740 0.7027 1.8776 1.0023 3000 + Detection Variances (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) 1.0167 1.2804 0.2030 0.6974 3.5672 1.0388 1159 -scale(det.cov.1) 0.3467 0.4305 0.0662 0.2335 1.3697 1.0150 2724 -scale(det.cov.2) 2.7573 3.3470 0.5307 1.8997 10.0553 1.0046 2761 +(Intercept) 0.9177 0.8501 0.2096 0.6714 3.2584 1.0159 2383 +scale(det.cov.1) 0.3304 0.4554 0.0597 0.2251 1.2220 1.0335 2422 +scale(det.cov.2) 2.7012 3.0959 0.5506 1.8827 10.3661 1.0204 2841 ---------------------------------------- Species Level ---------------------------------------- Abundance (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept)-sp1 -1.6735 0.3287 -2.3222 -1.6650 -1.0488 1.0160 463 -(Intercept)-sp2 -0.5628 0.2745 -1.1066 -0.5636 -0.0053 1.0462 214 -(Intercept)-sp3 1.0694 0.2676 0.5669 1.0730 1.5822 1.5827 60 -(Intercept)-sp4 -1.6473 0.3297 -2.3043 -1.6538 -1.0071 1.0085 577 -(Intercept)-sp5 -0.5719 0.3172 -1.1902 -0.5802 0.0474 1.0653 109 -(Intercept)-sp6 -0.7157 0.2936 -1.2813 -0.7077 -0.1272 1.1062 200 -scale(abund.cov.1)-sp1 1.0300 0.1363 0.7634 1.0267 1.2961 1.0070 1204 -scale(abund.cov.1)-sp2 -0.7193 0.0770 -0.8670 -0.7193 -0.5667 1.0222 1249 -scale(abund.cov.1)-sp3 0.2557 0.0405 0.1765 0.2554 0.3358 1.0040 2070 -scale(abund.cov.1)-sp4 0.2662 0.1625 -0.0523 0.2671 0.5875 1.0005 1892 -scale(abund.cov.1)-sp5 1.2825 0.0863 1.1175 1.2815 1.4521 1.0106 627 -scale(abund.cov.1)-sp6 -0.0100 0.0766 -0.1585 -0.0109 0.1407 0.9997 2544 +(Intercept)-sp1 -1.6540 0.3338 -2.3148 -1.6529 -1.0181 1.0161 388 +(Intercept)-sp2 -0.5394 0.2877 -1.1184 -0.5276 0.0086 1.0292 221 +(Intercept)-sp3 1.0239 0.2731 0.5083 1.0411 1.5912 1.1371 51 +(Intercept)-sp4 -1.6233 0.3413 -2.2931 -1.6208 -0.9511 1.0155 573 +(Intercept)-sp5 -0.5496 0.3066 -1.1587 -0.5529 0.0732 1.0820 107 +(Intercept)-sp6 -0.6855 0.2726 -1.1980 -0.6842 -0.1578 1.0104 195 +scale(abund.cov.1)-sp1 1.0248 0.1317 0.7740 1.0249 1.2840 0.9998 1227 +scale(abund.cov.1)-sp2 -0.7174 0.0782 -0.8710 -0.7195 -0.5648 0.9999 1290 +scale(abund.cov.1)-sp3 0.2579 0.0391 0.1844 0.2582 0.3328 1.0000 2052 +scale(abund.cov.1)-sp4 0.2668 0.1639 -0.0514 0.2684 0.5918 1.0018 1713 +scale(abund.cov.1)-sp5 1.2799 0.0866 1.1138 1.2801 1.4482 1.0064 539 +scale(abund.cov.1)-sp6 -0.0130 0.0754 -0.1596 -0.0136 0.1322 1.0053 2411 Detection (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept)-sp1 0.2238 0.2750 -0.3003 0.2267 0.7434 1.0023 659 -(Intercept)-sp2 0.1996 0.1745 -0.1521 0.1968 0.5398 1.0099 790 -(Intercept)-sp3 0.6227 0.1076 0.4172 0.6221 0.8346 1.0014 695 -(Intercept)-sp4 -0.1993 0.3424 -0.8663 -0.1997 0.4452 1.0176 967 -(Intercept)-sp5 -0.6907 0.2745 -1.2514 -0.6843 -0.1733 1.0356 134 -(Intercept)-sp6 1.5000 0.1995 1.1254 1.4917 1.9105 1.0085 1029 -scale(det.cov.1)-sp1 0.1518 0.2171 -0.2867 0.1584 0.5782 1.0055 1524 -scale(det.cov.1)-sp2 0.7398 0.1428 0.4651 0.7399 1.0244 1.0045 1946 -scale(det.cov.1)-sp3 1.1425 0.0881 0.9731 1.1417 1.3195 1.0017 1247 -scale(det.cov.1)-sp4 0.1772 0.3222 -0.4953 0.1913 0.7815 1.0123 1883 -scale(det.cov.1)-sp5 0.6868 0.1304 0.4382 0.6829 0.9498 1.0085 790 -scale(det.cov.1)-sp6 0.6634 0.1345 0.4045 0.6619 0.9346 1.0035 1818 -scale(det.cov.2)-sp1 0.6336 0.1961 0.2658 0.6262 1.0403 1.0008 1972 -scale(det.cov.2)-sp2 1.2872 0.1834 0.9443 1.2844 1.6596 1.0030 1873 -scale(det.cov.2)-sp3 0.9603 0.0811 0.8040 0.9588 1.1246 1.0012 1907 -scale(det.cov.2)-sp4 2.9437 0.6506 1.8404 2.8813 4.3336 1.0036 1660 -scale(det.cov.2)-sp5 -0.1191 0.1099 -0.3398 -0.1162 0.0912 1.0125 1151 -scale(det.cov.2)-sp6 -0.8564 0.1537 -1.1620 -0.8551 -0.5559 1.0007 1746 +(Intercept)-sp1 0.2073 0.2769 -0.3514 0.2118 0.7388 1.0135 667 +(Intercept)-sp2 0.1889 0.1731 -0.1505 0.1909 0.5290 1.0374 827 +(Intercept)-sp3 0.6277 0.1069 0.4218 0.6298 0.8267 1.0095 729 +(Intercept)-sp4 -0.2009 0.3332 -0.8395 -0.1961 0.4650 1.0064 1173 +(Intercept)-sp5 -0.6691 0.2699 -1.2348 -0.6482 -0.1826 1.1159 158 +(Intercept)-sp6 1.4959 0.2065 1.0980 1.4937 1.9043 1.0002 711 +scale(det.cov.1)-sp1 0.1599 0.2156 -0.2708 0.1584 0.5630 1.0025 1911 +scale(det.cov.1)-sp2 0.7384 0.1425 0.4667 0.7341 1.0237 1.0163 2014 +scale(det.cov.1)-sp3 1.1437 0.0872 0.9793 1.1437 1.3170 1.0051 1279 +scale(det.cov.1)-sp4 0.1785 0.3136 -0.4596 0.1906 0.7692 1.0066 2147 +scale(det.cov.1)-sp5 0.6884 0.1340 0.4383 0.6871 0.9617 1.0124 749 +scale(det.cov.1)-sp6 0.6536 0.1372 0.3851 0.6513 0.9237 1.0014 1892 +scale(det.cov.2)-sp1 0.6313 0.2038 0.2501 0.6227 1.0418 1.0000 1747 +scale(det.cov.2)-sp2 1.2839 0.1806 0.9382 1.2803 1.6497 1.0082 1853 +scale(det.cov.2)-sp3 0.9588 0.0832 0.8031 0.9578 1.1341 1.0006 1942 +scale(det.cov.2)-sp4 2.9508 0.6496 1.8662 2.8821 4.3887 1.0082 1510 +scale(det.cov.2)-sp5 -0.1246 0.1070 -0.3333 -0.1253 0.0846 1.0211 1305 +scale(det.cov.2)-sp6 -0.8546 0.1525 -1.1598 -0.8506 -0.5611 1.0009 1414 
# Or more explicitly
 # summary(out.ms, level = 'both')

We see adequate convergence and mixing of all model parameters, with ESS values fairly large and Rhat values less (or very close to) 1.1. Looking at the community-level variance parameters, we can see there is fairly minimal variation among average abundance of each species (i.e., the intercept), with larger variable in the effect of the covariate on abundance, which is further seen by looking at the species-specific effects of the covariate, which range from highly positive to highly negative.

@@ -1707,17 +1708,17 @@

Posterior predictive checks

---------------------------------------- Community Level ---------------------------------------- -Bayesian p-value: 0.5848 +Bayesian p-value: 0.5508 ---------------------------------------- Species Level ---------------------------------------- -sp1 Bayesian p-value: 0.5687 -sp2 Bayesian p-value: 0.4857 -sp3 Bayesian p-value: 0.4977 -sp4 Bayesian p-value: 0.674 -sp5 Bayesian p-value: 0.7103 -sp6 Bayesian p-value: 0.5727 +sp1 Bayesian p-value: 0.5547 +sp2 Bayesian p-value: 0.4673 +sp3 Bayesian p-value: 0.506 +sp4 Bayesian p-value: 0.59 +sp5 Bayesian p-value: 0.6497 +sp6 Bayesian p-value: 0.537 Fit statistic: chi-squared

Here we see mixed results of the posterior predictive check, with the first and fourth species having a fairly small Bayesian p-value, suggesting the model does not fit very well for those species, which as with the single-species case, is likely related to additional spatial variation in abundance that we are not accounting for.

@@ -1734,13 +1735,13 @@

Model selection using WAIC

Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
-
       elpd         pD      WAIC
-1 -163.6622  14.618585  356.5615
-2 -324.0871  23.530929  695.2361
-3 -951.7242 101.630446 2106.7093
-4 -113.5199   8.560084  244.1600
-5 -269.0655  12.665373  563.4617
-6 -322.3626  13.934981  672.5952
+
       elpd        pD      WAIC
+1 -163.6351 14.762143  356.7944
+2 -324.0913 23.360959  694.9045
+3 -953.3861 99.038745 2104.8498
+4 -113.5324  8.466194  243.9972
+5 -269.1548 12.425717  563.1610
+6 -322.4475 13.752504  672.4000
# Single-species Poisson N-mixture model for species 1
 waicAbund(out)
N.max not specified. Setting upper index of integration of N to 10 plus
@@ -1787,66 +1788,66 @@ 
Model selection using WAIC

 Sampling ... 
 Batch: 200 of 800, 25.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     44.0        0.32825
-    1   alpha[1]    44.0        0.41729
-    2   beta[1]     52.0        0.26875
-    2   alpha[1]    40.0        0.25310
-    3   beta[1]     48.0        0.25821
-    3   alpha[1]    24.0        0.15047
-    4   beta[1]     32.0        0.31538
-    4   alpha[1]    48.0        0.64793
-    5   beta[1]     32.0        0.16301
-    5   alpha[1]    40.0        0.21568
-    6   beta[1]     52.0        0.20312
-    6   alpha[1]    52.0        0.32825
-    1   kappa       24.0        2.00743
-    2   kappa       28.0        0.54709
-    3   kappa       48.0        0.32525
-    4   kappa       40.0        3.05522
-    5   kappa       36.0        2.82033
-    6   kappa       28.0        2.60349
+    1   beta[1]     44.0        0.31538
+    1   alpha[1]    56.0        0.38521
+    2   beta[1]     60.0        0.29701
+    2   alpha[1]    36.0        0.23836
+    3   beta[1]     28.0        0.26875
+    3   alpha[1]    44.0        0.13890
+    4   beta[1]     40.0        0.32825
+    4   alpha[1]    36.0        0.67437
+    5   beta[1]     36.0        0.15978
+    5   alpha[1]    32.0        0.22003
+    6   beta[1]     48.0        0.17658
+    6   alpha[1]    60.0        0.27418
+    1   kappa       60.0        1.74517
+    2   kappa       24.0        0.55814
+    3   kappa       48.0        0.31250
+    4   kappa       52.0        3.37654
+    5   kappa       40.0        2.99473
+    6   kappa       64.0        2.60349
 -------------------------------------------------
 Batch: 400 of 800, 50.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     48.0        0.30302
-    1   alpha[1]    44.0        0.41729
-    2   beta[1]     40.0        0.26343
-    2   alpha[1]    52.0        0.27418
-    3   beta[1]     52.0        0.27418
-    3   alpha[1]    48.0        0.14457
-    4   beta[1]     44.0        0.32175
-    4   alpha[1]    36.0        0.62252
-    5   beta[1]     44.0        0.15351
-    5   alpha[1]    36.0        0.22448
-    6   beta[1]     32.0        0.20312
-    6   alpha[1]    40.0        0.28537
-    1   kappa       28.0        2.26337
-    2   kappa       28.0        0.56941
-    3   kappa       48.0        0.32525
-    4   kappa       32.0        4.20743
-    5   kappa       44.0        3.37654
-    6   kappa       40.0        3.17991
+    1   beta[1]     20.0        0.29701
+    1   alpha[1]    40.0        0.40903
+    2   beta[1]     36.0        0.26343
+    2   alpha[1]    60.0        0.24809
+    3   beta[1]     36.0        0.24809
+    3   alpha[1]    40.0        0.14171
+    4   beta[1]     24.0        0.32825
+    4   alpha[1]    20.0        0.63510
+    5   beta[1]     68.0        0.15661
+    5   alpha[1]    48.0        0.20312
+    6   beta[1]     36.0        0.18015
+    6   alpha[1]    36.0        0.29113
+    1   kappa       72.0        2.04798
+    2   kappa       24.0        0.56941
+    3   kappa       60.0        0.31250
+    4   kappa       20.0        4.04246
+    5   kappa       28.0        3.24415
+    6   kappa       24.0        3.17991
 -------------------------------------------------
 Batch: 600 of 800, 75.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     56.0        0.29113
-    1   alpha[1]    48.0        0.41729
-    2   beta[1]     44.0        0.28537
-    2   alpha[1]    36.0        0.25310
-    3   beta[1]     40.0        0.23836
-    3   alpha[1]    56.0        0.14171
-    4   beta[1]     44.0        0.32825
-    4   alpha[1]    48.0        0.64793
-    5   beta[1]     44.0        0.15978
-    5   alpha[1]    36.0        0.22448
-    6   beta[1]     24.0        0.18750
-    6   alpha[1]    40.0        0.29701
-    1   kappa       40.0        2.50141
-    2   kappa       16.0        0.56941
-    3   kappa       64.0        0.29430
-    4   kappa       28.0        3.88395
-    5   kappa       44.0        3.17991
-    6   kappa       76.0        3.24415
+    1   beta[1]     44.0        0.31538
+    1   alpha[1]    36.0        0.38521
+    2   beta[1]     40.0        0.27418
+    2   alpha[1]    28.0        0.24809
+    3   beta[1]     44.0        0.23364
+    3   alpha[1]    48.0        0.13615
+    4   beta[1]     48.0        0.35559
+    4   alpha[1]    52.0        0.67437
+    5   beta[1]     28.0        0.15047
+    5   alpha[1]    52.0        0.22003
+    6   beta[1]     56.0        0.18750
+    6   alpha[1]    56.0        0.30302
+    1   kappa       28.0        2.26337
+    2   kappa       24.0        0.56941
+    3   kappa       60.0        0.29430
+    4   kappa       60.0        3.88395
+    5   kappa       40.0        2.99473
+    6   kappa       28.0        3.17991
 -------------------------------------------------
 Batch: 800 of 800, 100.00%
 ----------------------------------------
@@ -1855,66 +1856,66 @@ 
Model selection using WAIC

 Sampling ... 
 Batch: 200 of 800, 25.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     24.0        0.30302
-    1   alpha[1]    32.0        0.37758
-    2   beta[1]     40.0        0.26875
-    2   alpha[1]    48.0        0.23364
-    3   beta[1]     40.0        0.24318
-    3   alpha[1]    52.0        0.15047
-    4   beta[1]     52.0        0.32825
-    4   alpha[1]    44.0        0.70189
-    5   beta[1]     40.0        0.15661
-    5   alpha[1]    56.0        0.22901
-    6   beta[1]     56.0        0.19515
-    6   alpha[1]    32.0        0.29701
-    1   kappa       32.0        2.13156
-    2   kappa       56.0        0.56941
-    3   kappa       52.0        0.30631
-    4   kappa       52.0        4.04246
-    5   kappa       52.0        3.37654
-    6   kappa       48.0        3.44476
+    1   beta[1]     40.0        0.29701
+    1   alpha[1]    40.0        0.39299
+    2   beta[1]     24.0        0.29701
+    2   alpha[1]    20.0        0.25310
+    3   beta[1]     60.0        0.24809
+    3   alpha[1]    36.0        0.13346
+    4   beta[1]     32.0        0.32175
+    4   alpha[1]    24.0        0.64793
+    5   beta[1]     40.0        0.15978
+    5   alpha[1]    36.0        0.22448
+    6   beta[1]     36.0        0.17658
+    6   alpha[1]    52.0        0.30914
+    1   kappa       60.0        2.40332
+    2   kappa       32.0        0.60462
+    3   kappa       28.0        0.30631
+    4   kappa       52.0        3.80704
+    5   kappa       44.0        3.37654
+    6   kappa       80.0        3.44476
 -------------------------------------------------
 Batch: 400 of 800, 50.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     40.0        0.31538
-    1   alpha[1]    52.0        0.40903
-    2   beta[1]     72.0        0.26343
-    2   alpha[1]    48.0        0.22901
-    3   beta[1]     40.0        0.24318
-    3   alpha[1]    44.0        0.14457
-    4   beta[1]     56.0        0.37010
-    4   alpha[1]    52.0        0.67437
-    5   beta[1]     28.0        0.13081
-    5   alpha[1]    48.0        0.22448
-    6   beta[1]     44.0        0.18750
-    6   alpha[1]    40.0        0.30302
-    1   kappa       56.0        2.40332
-    2   kappa       28.0        0.60462
-    3   kappa       56.0        0.27716
-    4   kappa       40.0        4.20743
-    5   kappa       40.0        3.30968
-    6   kappa       36.0        3.37654
+    1   beta[1]     44.0        0.31538
+    1   alpha[1]    56.0        0.40903
+    2   beta[1]     48.0        0.26343
+    2   alpha[1]    36.0        0.25310
+    3   beta[1]     52.0        0.24809
+    3   alpha[1]    52.0        0.13890
+    4   beta[1]     32.0        0.33488
+    4   alpha[1]    24.0        0.61020
+    5   beta[1]     48.0        0.15661
+    5   alpha[1]    44.0        0.22448
+    6   beta[1]     40.0        0.17308
+    6   alpha[1]    40.0        0.29113
+    1   kappa       52.0        2.08935
+    2   kappa       44.0        0.58092
+    3   kappa       60.0        0.29430
+    4   kappa       40.0        3.65777
+    5   kappa       36.0        3.24415
+    6   kappa       24.0        3.24415
 -------------------------------------------------
 Batch: 600 of 800, 75.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     28.0        0.31538
-    1   alpha[1]    60.0        0.40093
-    2   beta[1]     32.0        0.29701
-    2   alpha[1]    36.0        0.26343
-    3   beta[1]     48.0        0.25310
-    3   alpha[1]    36.0        0.14171
-    4   beta[1]     40.0        0.32175
-    4   alpha[1]    40.0        0.66102
-    5   beta[1]     32.0        0.14457
-    5   alpha[1]    28.0        0.22003
-    6   beta[1]     40.0        0.17658
-    6   alpha[1]    44.0        0.29113
-    1   kappa       56.0        2.17462
+    1   beta[1]     44.0        0.33488
+    1   alpha[1]    44.0        0.40093
+    2   beta[1]     32.0        0.27418
+    2   alpha[1]    32.0        0.24809
+    3   beta[1]     68.0        0.24809
+    3   alpha[1]    40.0        0.13081
+    4   beta[1]     48.0        0.32175
+    4   alpha[1]    32.0        0.61020
+    5   beta[1]     32.0        0.14749
+    5   alpha[1]    36.0        0.23364
+    6   beta[1]     32.0        0.19129
+    6   alpha[1]    36.0        0.28537
+    1   kappa       56.0        2.40332
     2   kappa       40.0        0.58092
-    3   kappa       36.0        0.28847
-    4   kappa       52.0        3.96241
-    5   kappa       60.0        2.99473
-    6   kappa       52.0        3.30968
+    3   kappa       44.0        0.32525
+    4   kappa       44.0        4.04246
+    5   kappa       40.0        3.51434
+    6   kappa       56.0        3.24415
 -------------------------------------------------
 Batch: 800 of 800, 100.00%
 ----------------------------------------
@@ -1923,66 +1924,66 @@ 
Model selection using WAIC

 Sampling ... 
 Batch: 200 of 800, 25.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     36.0        0.33488
-    1   alpha[1]    44.0        0.40093
-    2   beta[1]     24.0        0.25821
-    2   alpha[1]    48.0        0.23364
-    3   beta[1]     32.0        0.25310
-    3   alpha[1]    40.0        0.13081
-    4   beta[1]     40.0        0.32175
-    4   alpha[1]    36.0        0.61020
-    5   beta[1]     36.0        0.14749
-    5   alpha[1]    56.0        0.20722
-    6   beta[1]     44.0        0.19129
-    6   alpha[1]    28.0        0.31538
-    1   kappa       48.0        2.00743
-    2   kappa       52.0        0.54709
-    3   kappa       28.0        0.29430
-    4   kappa       32.0        4.12412
-    5   kappa       40.0        3.11694
-    6   kappa       44.0        3.30968
+    1   beta[1]     44.0        0.32825
+    1   alpha[1]    48.0        0.38521
+    2   beta[1]     52.0        0.26875
+    2   alpha[1]    40.0        0.24318
+    3   beta[1]     36.0        0.24318
+    3   alpha[1]    40.0        0.13890
+    4   beta[1]     44.0        0.30914
+    4   alpha[1]    40.0        0.63510
+    5   beta[1]     36.0        0.14457
+    5   alpha[1]    36.0        0.23836
+    6   beta[1]     36.0        0.18379
+    6   alpha[1]    36.0        0.30302
+    1   kappa       40.0        2.13156
+    2   kappa       48.0        0.62930
+    3   kappa       44.0        0.33183
+    4   kappa       32.0        3.88395
+    5   kappa       44.0        3.24415
+    6   kappa       48.0        3.37654
 -------------------------------------------------
 Batch: 400 of 800, 50.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     28.0        0.30914
-    1   alpha[1]    64.0        0.39299
-    2   beta[1]     32.0        0.27418
-    2   alpha[1]    48.0        0.26343
-    3   beta[1]     48.0        0.26875
-    3   alpha[1]    48.0        0.13346
-    4   beta[1]     44.0        0.31538
-    4   alpha[1]    36.0        0.68800
-    5   beta[1]     48.0        0.15351
-    5   alpha[1]    48.0        0.20722
-    6   beta[1]     56.0        0.17658
-    6   alpha[1]    28.0        0.27418
-    1   kappa       28.0        2.08935
-    2   kappa       48.0        0.59265
-    3   kappa       32.0        0.33183
-    4   kappa       56.0        4.04246
-    5   kappa       36.0        3.17991
-    6   kappa       56.0        3.73166
+    1   beta[1]     48.0        0.33488
+    1   alpha[1]    44.0        0.42572
+    2   beta[1]     56.0        0.27972
+    2   alpha[1]    40.0        0.25310
+    3   beta[1]     52.0        0.26875
+    3   alpha[1]    52.0        0.15047
+    4   beta[1]     32.0        0.32175
+    4   alpha[1]    44.0        0.63510
+    5   beta[1]     48.0        0.13890
+    5   alpha[1]    52.0        0.21568
+    6   beta[1]     36.0        0.19515
+    6   alpha[1]    32.0        0.27972
+    1   kappa       44.0        2.26337
+    2   kappa       36.0        0.58092
+    3   kappa       36.0        0.31881
+    4   kappa       52.0        4.12412
+    5   kappa       48.0        3.17991
+    6   kappa       36.0        3.17991
 -------------------------------------------------
 Batch: 600 of 800, 75.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     52.0        0.32825
-    1   alpha[1]    36.0        0.41729
-    2   beta[1]     44.0        0.27972
-    2   alpha[1]    28.0        0.25310
-    3   beta[1]     48.0        0.25821
-    3   alpha[1]    44.0        0.13081
-    4   beta[1]     32.0        0.32825
-    4   alpha[1]    48.0        0.62252
-    5   beta[1]     56.0        0.14749
-    5   alpha[1]    52.0        0.22448
-    6   beta[1]     40.0        0.16630
-    6   alpha[1]    36.0        0.29701
-    1   kappa       48.0        2.17462
-    2   kappa       52.0        0.58092
-    3   kappa       32.0        0.29430
-    4   kappa       36.0        3.96241
-    5   kappa       48.0        3.17991
-    6   kappa       40.0        3.17991
+    1   beta[1]     36.0        0.29113
+    1   alpha[1]    44.0        0.38521
+    2   beta[1]     36.0        0.26875
+    2   alpha[1]    44.0        0.25310
+    3   beta[1]     56.0        0.25310
+    3   alpha[1]    20.0        0.15047
+    4   beta[1]     48.0        0.32825
+    4   alpha[1]    40.0        0.67437
+    5   beta[1]     36.0        0.13346
+    5   alpha[1]    40.0        0.21141
+    6   beta[1]     36.0        0.18379
+    6   alpha[1]    32.0        0.30914
+    1   kappa       48.0        1.89052
+    2   kappa       36.0        0.60462
+    3   kappa       52.0        0.27716
+    4   kappa       56.0        4.46761
+    5   kappa       44.0        3.24415
+    6   kappa       28.0        3.37654
 -------------------------------------------------
 Batch: 800 of 800, 100.00%
# Negative binomial
@@ -1995,13 +1996,13 @@ 
Model selection using WAIC

 
Currently on species 4 out of 6

 
Currently on species 5 out of 6

 
Currently on species 6 out of 6

-
       elpd        pD      WAIC
-1 -160.5396 17.049287  355.1777
-2 -306.4749 13.083567  639.1169
-3 -640.0901 16.421039 1313.0223
-4 -113.6765  8.339817  244.0327
-5 -269.4515 12.280588  563.4641
-6 -321.0404 13.827084  669.7350

+
       elpd       pD      WAIC
+1 -160.5475 17.51028  356.1156
+2 -306.3167 13.01415  638.6616
+3 -640.0349 16.23307 1312.5360
+4 -113.6294  8.30432  243.8675
+5 -269.5406 12.23824  563.5576
+6 -320.9739 13.81813  669.5841

 
# Poisson
 waicAbund(out.ms, by.sp = TRUE)

 
N.max not specified. Setting upper index of integration of N to 10 plus
@@ -2012,13 +2013,13 @@ 
Model selection using WAIC

 
Currently on species 4 out of 6

 
Currently on species 5 out of 6

 
Currently on species 6 out of 6

-
       elpd         pD      WAIC
-1 -163.6622  14.618585  356.5615
-2 -324.0871  23.530929  695.2361
-3 -951.7242 101.630446 2106.7093
-4 -113.5199   8.560084  244.1600
-5 -269.0655  12.665373  563.4617
-6 -322.3626  13.934981  672.5952

+
       elpd        pD      WAIC
+1 -163.6351 14.762143  356.7944
+2 -324.0913 23.360959  694.9045
+3 -953.3861 99.038745 2104.8498
+4 -113.5324  8.466194  243.9972
+5 -269.1548 12.425717  563.1610
+6 -322.4475 13.752504  672.4000

 
Here we see the negative binomial model provides pretty substantial improvements for the first two species in the data set, while the last four species all have WAIC values that are very similar across the two models.
@@ -2028,7 +2029,7 @@

Prediction

# Look at the resulting object str(out.pred)
List of 3
- $ mu.0.samples: num [1:3000, 1:6, 1:100] 0.00149 0.00184 0.00426 0.00176 0.00566 ...
+ $ mu.0.samples: num [1:3000, 1:6, 1:100] 0.00553 0.00981 0.00387 0.00273 0.00218 ...
  $ N.0.samples : num [1:3000, 1:6, 1:100] 0 0 0 0 0 0 0 0 0 0 ...
  $ call        : language predict.msNMix(object = out.ms.nb, X.0 = X.0, ignore.RE = TRUE)
  - attr(*, "class")= chr "predict.msNMix"
@@ -2047,7 +2048,7 @@

Prediction

theme_bw() + facet_wrap(vars(sp)) + labs(x = 'Covariate', y = 'Expected abundance') -

+

@@ -2104,13 +2105,13 @@

Fitting latent factor multi-species N-mixture models with lfMsNMix()lambda.inits[lower.tri(lambda.inits)] <- rnorm(sum(lower.tri(lambda.inits))) # Check it out. lambda.inits

-
           [,1]       [,2]       [,3]
-[1,]  1.0000000  0.0000000  0.0000000
-[2,] -0.7914937  1.0000000  0.0000000
-[3,]  0.4053592 -1.0237382  1.0000000
-[4,] -0.3204730 -0.2138603 -0.3888896
-[5,]  0.1519413 -0.4604921  0.3511065
-[6,]  2.6696774 -1.3026908  2.1361886
+
            [,1]       [,2]        [,3]
+[1,]  1.00000000  0.0000000  0.00000000
+[2,]  0.11177665  1.0000000  0.00000000
+[3,] -1.02873242  0.1356658  1.00000000
+[4,]  0.02480722 -0.6285882  0.74969371
+[5,] -1.08207620 -2.4490989  0.29699355
+[6,]  0.86191993 -0.4828777 -0.05669138
w.inits <- matrix(0, n.factors, ncol(dataNMixSim$y))
 ms.inits <- list(alpha.comm = 0,
                  beta.comm = 0,
@@ -2174,48 +2175,48 @@ 
Fitting latent factor multi-species N-mixture models with lfMsNMix()Fitting latent factor multi-species N-mixture models with lfMsNMix()Fitting latent factor multi-species N-mixture models with lfMsNMix()
summary(out.lf.ms)
@@ -2332,105 +2333,105 @@

Fitting latent factor multi-species N-mixture models with lfMsNMix() +(Intercept)-sp1 0.0218 0.3493 -0.6960 0.0404 0.6147 1.0310 311 +(Intercept)-sp2 -0.2829 0.2576 -0.8546 -0.2628 0.1658 1.0741 217 +(Intercept)-sp3 0.3081 0.1477 0.0120 0.3172 0.5786 1.0315 301 +(Intercept)-sp4 -0.2440 0.3440 -0.9588 -0.2401 0.4135 1.0072 999 +(Intercept)-sp5 -0.8802 0.3312 -1.6155 -0.8389 -0.3173 1.0607 81 +(Intercept)-sp6 1.4879 0.2096 1.0645 1.4883 1.8896 1.0136 843 +scale(det.cov.1)-sp1 0.1331 0.2204 -0.3182 0.1370 0.5599 1.0004 1123 +scale(det.cov.1)-sp2 0.5285 0.1367 0.2629 0.5266 0.8067 1.0136 1030 +scale(det.cov.1)-sp3 1.0672 0.0961 0.8788 1.0657 1.2545 1.0223 502 +scale(det.cov.1)-sp4 0.1745 0.3128 -0.4494 0.1841 0.7472 1.0018 1894 +scale(det.cov.1)-sp5 0.6639 0.1322 0.4191 0.6589 0.9385 1.0196 518 +scale(det.cov.1)-sp6 0.6316 0.1383 0.3686 0.6331 0.9038 1.0020 1801 +scale(det.cov.2)-sp1 0.6063 0.2159 0.1951 0.5996 1.0466 1.0012 1107 +scale(det.cov.2)-sp2 1.1275 0.1894 0.7531 1.1227 1.5006 1.0044 838 +scale(det.cov.2)-sp3 0.8181 0.0892 0.6399 0.8179 0.9976 1.0046 877 +scale(det.cov.2)-sp4 2.9374 0.6604 1.8464 2.8763 4.4199 1.0018 1541 +scale(det.cov.2)-sp5 -0.1030 0.1068 -0.3160 -0.1034 0.1048 1.0038 1298 +scale(det.cov.2)-sp6 -0.8084 0.1587 -1.1315 -0.8037 -0.5044 1.0172 1355

We see adequate convergence of all model parameters shown in the summary output. Note that the factor loadings themselves are not shown in the summary() output, but are available in the lambda.samples portion of the model list object.

# Rhats for lambda (the factor loadings)
 out.lf.ms$rhat$lambda.lower.tri
-
 [1] 1.265012 1.289588 1.004099 1.089482 1.093478 1.344807 1.058168 1.246524
- [9] 1.024881 1.086373 1.073848 1.075239
+
 [1] 1.143639 1.112604 1.007551 1.031422 1.004691 1.200399 1.033574 1.126968
+ [9] 1.011803 1.019237 1.042185 1.002496
# ESS for lambda
 out.lf.ms$ESS$lambda
    sp1-1     sp2-1     sp3-1     sp4-1     sp5-1     sp6-1     sp1-2     sp2-2 
-  0.00000 100.69735  96.91901 510.25489 301.24250 228.63224   0.00000   0.00000 
+  0.00000  76.66690  82.99070 549.53733 268.76485 250.95721   0.00000   0.00000 
     sp3-2     sp4-2     sp5-2     sp6-2     sp1-3     sp2-3     sp3-3     sp4-3 
- 89.58583 431.02178 153.86812 227.50797   0.00000   0.00000   0.00000 365.22931 
+ 75.15241 338.30477 182.35198 213.65427   0.00000   0.00000   0.00000 335.12783 
     sp5-3     sp6-3 
-167.14141 155.84477
+145.69290 239.37097 
# Posterior quantiles for the latent factor loadings
 summary(out.lf.ms$lambda.samples)$quantiles
-
             2.5%         25%         50%         75%        97.5%
-sp1-1  1.00000000  1.00000000  1.00000000  1.00000000  1.000000000
-sp2-1 -0.86342142 -0.58785077 -0.43155941 -0.25910723  0.174385751
-sp3-1  0.12790428  0.41067096  0.54155475  0.68421763  0.978136593
-sp4-1 -0.58063009 -0.22712171 -0.05204231  0.11036285  0.411968395
-sp5-1 -0.08074725  0.12379780  0.22286782  0.32865332  0.554909887
-sp6-1 -0.48228897 -0.25911451 -0.13348675 -0.01876527  0.203846993
-sp1-2  0.00000000  0.00000000  0.00000000  0.00000000  0.000000000
-sp2-2  1.00000000  1.00000000  1.00000000  1.00000000  1.000000000
-sp3-2  0.32508345  0.64582800  0.77419890  0.90160205  1.114518636
-sp4-2 -0.46835862 -0.12795533  0.06820461  0.26581595  0.622402620
-sp5-2 -0.59866727 -0.40798739 -0.30316176 -0.18990846  0.059981443
-sp6-2 -0.43539788 -0.24153294 -0.15196622 -0.05894365  0.117152810
-sp1-3  0.00000000  0.00000000  0.00000000  0.00000000  0.000000000
-sp2-3  0.00000000  0.00000000  0.00000000  0.00000000  0.000000000
-sp3-3  1.00000000  1.00000000  1.00000000  1.00000000  1.000000000
-sp4-3 -0.40749538 -0.01542192  0.18358990  0.36575029  0.725477022
-sp5-3 -0.36544065 -0.14500001 -0.03791512  0.06384060  0.287054785
-sp6-3 -0.64769066 -0.47056542 -0.37397316 -0.26771448 -0.002994657
+
             2.5%         25%         50%         75%         97.5%
+sp1-1  1.00000000  1.00000000  1.00000000  1.00000000  1.0000000000
+sp2-1 -0.99894551 -0.67232826 -0.48424672 -0.27557071  0.2263571560
+sp3-1  0.03485244  0.34791763  0.50311626  0.65505983  0.9179470059
+sp4-1 -0.61745625 -0.26821195 -0.09145247  0.08096110  0.3996860281
+sp5-1 -0.10921399  0.11165803  0.22404852  0.33234963  0.5427660045
+sp6-1 -0.46493884 -0.26058390 -0.14275230 -0.02629113  0.1790912706
+sp1-2  0.00000000  0.00000000  0.00000000  0.00000000  0.0000000000
+sp2-2  1.00000000  1.00000000  1.00000000  1.00000000  1.0000000000
+sp3-2  0.30514619  0.64608911  0.78153539  0.90877211  1.1420370455
+sp4-2 -0.48402277 -0.14412586  0.03404940  0.22520681  0.5895873842
+sp5-2 -0.62316223 -0.41838344 -0.30800474 -0.19843323  0.0650580641
+sp6-2 -0.47049149 -0.26567432 -0.16410072 -0.06581812  0.1162878023
+sp1-3  0.00000000  0.00000000  0.00000000  0.00000000  0.0000000000
+sp2-3  0.00000000  0.00000000  0.00000000  0.00000000  0.0000000000
+sp3-3  1.00000000  1.00000000  1.00000000  1.00000000  1.0000000000
+sp4-3 -0.31621784  0.02969278  0.20836047  0.37201835  0.7213465937
+sp5-3 -0.35200296 -0.12637166 -0.01190123  0.10939469  0.3155836871
+sp6-3 -0.62821825 -0.45469430 -0.36458875 -0.24954831 -0.0001262112

Notice the Rhat values are only reported for the elements of the factor loadings matrix that are estimated and not fixed at any specific value, while the ESS values are reported for all elements of the factor loadings matrix, and take value 0 for those parameters that are fixed for identifiability purposes. We can inspect the latent factor loadings as well as the latent factors (stored in out.lf.ms$w.samples) to provide information on any groupings that arise from the species in the modeled community. See Appendix S1 in Doser, Finley, and Banerjee (2023) for a discussion on using factor models as a model-based ordination technique.

@@ -2458,17 +2459,17 @@

Posterior predictive checks

---------------------------------------- Community Level ---------------------------------------- -Bayesian p-value: 0.5726 +Bayesian p-value: 0.5507 ---------------------------------------- Species Level ---------------------------------------- -sp1 Bayesian p-value: 0.4277 -sp2 Bayesian p-value: 0.6153 -sp3 Bayesian p-value: 0.5283 -sp4 Bayesian p-value: 0.637 -sp5 Bayesian p-value: 0.6593 -sp6 Bayesian p-value: 0.5677 +sp1 Bayesian p-value: 0.4493 +sp2 Bayesian p-value: 0.5867 +sp3 Bayesian p-value: 0.5217 +sp4 Bayesian p-value: 0.5647 +sp5 Bayesian p-value: 0.6213 +sp6 Bayesian p-value: 0.5607 Fit statistic: chi-squared
@@ -2485,12 +2486,12 @@

Model selection using WAIC

Currently on species 5 out of 6
Currently on species 6 out of 6
       elpd       pD      WAIC
-1 -134.0612 27.59438  323.3112
-2 -243.4556 48.82173  584.5547
-3 -464.5515 84.42150 1097.9459
-4 -109.0633 14.45701  247.0406
-5 -254.7556 23.57167  556.6546
-6 -297.7195 28.75067  652.9404
+1 -133.9715 27.58378 323.1107 +2 -242.5925 48.85054 582.8862 +3 -464.9209 84.97852 1099.7988 +4 -108.9876 14.28511 246.5454 +5 -254.4298 24.49645 557.8526 +6 -297.9585 28.65566 653.2282 
# Without latent factors
 waicAbund(out.ms, by.sp = TRUE)
N.max not specified. Setting upper index of integration of N to 10 plus
@@ -2501,13 +2502,13 @@ 
Model selection using WAIC

 
Currently on species 4 out of 6

 
Currently on species 5 out of 6

 
Currently on species 6 out of 6

-
       elpd         pD      WAIC
-1 -163.6622  14.618585  356.5615
-2 -324.0871  23.530929  695.2361
-3 -951.7242 101.630446 2106.7093
-4 -113.5199   8.560084  244.1600
-5 -269.0655  12.665373  563.4617
-6 -322.3626  13.934981  672.5952

+
       elpd        pD      WAIC
+1 -163.6351 14.762143  356.7944
+2 -324.0913 23.360959  694.9045
+3 -953.3861 99.038745 2104.8498
+4 -113.5324  8.466194  243.9972
+5 -269.1548 12.425717  563.1610
+6 -322.4475 13.752504  672.4000

Prediction

@@ -2516,7 +2517,7 @@

Prediction

# Look at the resulting object str(out.pred)
List of 3
- $ mu.0.samples: num [1:3000, 1:6, 1:100] 0.00971 0.00252 0.00303 0.00788 0.00218 ...
+ $ mu.0.samples: num [1:3000, 1:6, 1:100] 0.0003 0.000769 0.003598 0.003017 0.003083 ...
  $ N.0.samples : int [1:3000, 1:6, 1:100] 0 0 0 0 0 0 0 0 0 0 ...
  $ call        : language predict.lfMsNMix(object = out.lf.ms, X.0 = X.0, ignore.RE = TRUE)
  - attr(*, "class")= chr "predict.lfMsNMix"
@@ -2534,7 +2535,7 @@

Prediction

theme_bw() + facet_wrap(vars(sp)) + labs(x = 'Covariate', y = 'Expected abundance') -

+

@@ -2597,13 +2598,13 @@

Fitting spatial factor multi-species N-mixture models with sfMsNMix()< lambda.inits[lower.tri(lambda.inits)] <- rnorm(sum(lower.tri(lambda.inits))) # Check it out. lambda.inits

-
           [,1]        [,2]       [,3]
-[1,]  1.0000000  0.00000000  0.0000000
-[2,] -0.6554955  1.00000000  0.0000000
-[3,]  1.2301792  0.35756069  1.0000000
-[4,] -1.0778473 -1.68998630 -0.6272787
-[5,]  1.0211612 -0.69993175 -2.8331440
-[6,]  1.8368089  0.03369004  0.2451326
+
           [,1]      [,2]        [,3]
+[1,]  1.0000000 0.0000000  0.00000000
+[2,]  1.1415164 1.0000000  0.00000000
+[3,] -0.4009632 0.5965179  1.00000000
+[4,]  0.9566143 1.6568582  0.03883005
+[5,]  1.2062924 0.6397995  0.12112501
+[6,] -0.3879460 0.8318398 -1.05219916
w.inits <- matrix(0, n.factors, ncol(dataNMixSim$y))
 # Pair-wise distances between all sites
 dist.mat <- dist(dataNMixSim$coords)
@@ -2686,57 +2687,57 @@ 
Fitting spatial factor multi-species N-mixture models with sfMsNMix()<
 Sampling ... 
 Batch: 200 of 800, 25.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     48.0        0.26875
-    1   alpha[1]    44.0        0.40903
-    2   beta[1]     48.0        0.16630
-    2   alpha[1]    48.0        0.24318
+    1   beta[1]     44.0        0.26343
+    1   alpha[1]    44.0        0.44309
+    2   beta[1]     28.0        0.17308
+    2   alpha[1]    32.0        0.22901
     3   beta[1]     36.0        0.08425
-    3   alpha[1]    48.0        0.13346
-    4   beta[1]     32.0        0.30914
-    4   alpha[1]    48.0        0.64793
-    5   beta[1]     44.0        0.12320
-    5   alpha[1]    36.0        0.19910
-    6   beta[1]     32.0        0.17308
-    6   alpha[1]    48.0        0.29701
-    1   phi     8.0     1.19346
-    2   phi     28.0        0.70953
-    3   phi     36.0        0.47561
+    3   alpha[1]    56.0        0.13346
+    4   beta[1]     44.0        0.32175
+    4   alpha[1]    40.0        0.64793
+    5   beta[1]     40.0        0.13081
+    5   alpha[1]    32.0        0.22448
+    6   beta[1]     44.0        0.17308
+    6   alpha[1]    40.0        0.30302
+    1   phi     48.0        1.12395
+    2   phi     52.0        0.72387
+    3   phi     24.0        0.50503
 -------------------------------------------------
 Batch: 400 of 800, 50.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     28.0        0.24809
-    1   alpha[1]    28.0        0.35559
-    2   beta[1]     32.0        0.14749
-    2   alpha[1]    60.0        0.25310
-    3   beta[1]     32.0        0.08258
-    3   alpha[1]    28.0        0.12822
-    4   beta[1]     24.0        0.32825
-    4   alpha[1]    32.0        0.62252
-    5   beta[1]     28.0        0.12569
-    5   alpha[1]    36.0        0.22003
-    6   beta[1]     36.0        0.16301
-    6   alpha[1]    48.0        0.29701
-    1   phi     52.0        1.34562
-    2   phi     20.0        0.58092
-    3   phi     16.0        0.44792
+    1   beta[1]     44.0        0.24809
+    1   alpha[1]    60.0        0.37010
+    2   beta[1]     48.0        0.14171
+    2   alpha[1]    16.0        0.24809
+    3   beta[1]     28.0        0.08425
+    3   alpha[1]    32.0        0.14171
+    4   beta[1]     36.0        0.31538
+    4   alpha[1]    32.0        0.63510
+    5   beta[1]     40.0        0.12822
+    5   alpha[1]    40.0        0.19910
+    6   beta[1]     44.0        0.17308
+    6   alpha[1]    44.0        0.30302
+    1   phi     24.0        1.42883
+    2   phi     52.0        0.73849
+    3   phi     56.0        0.47561
 -------------------------------------------------
 Batch: 600 of 800, 75.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     60.0        0.25310
-    1   alpha[1]    48.0        0.36277
-    2   beta[1]     24.0        0.15661
-    2   alpha[1]    60.0        0.25821
-    3   beta[1]     52.0        0.08595
-    3   alpha[1]    40.0        0.14171
-    4   beta[1]     48.0        0.37010
-    4   alpha[1]    36.0        0.66102
-    5   beta[1]     40.0        0.11837
-    5   alpha[1]    44.0        0.20722
-    6   beta[1]     36.0        0.16301
-    6   alpha[1]    52.0        0.30914
-    1   phi     36.0        1.42883
-    2   phi     32.0        0.69548
-    3   phi     96.0        0.54709
+    1   beta[1]     52.0        0.26875
+    1   alpha[1]    24.0        0.37010
+    2   beta[1]     32.0        0.15661
+    2   alpha[1]    40.0        0.24318
+    3   beta[1]     40.0        0.07623
+    3   alpha[1]    60.0        0.15351
+    4   beta[1]     40.0        0.34165
+    4   alpha[1]    28.0        0.63510
+    5   beta[1]     44.0        0.13346
+    5   alpha[1]    56.0        0.20722
+    6   beta[1]     64.0        0.16966
+    6   alpha[1]    52.0        0.30302
+    1   phi     68.0        1.29285
+    2   phi     28.0        0.61684
+    3   phi     44.0        0.52564
 -------------------------------------------------
 Batch: 800 of 800, 100.00%
 ----------------------------------------
@@ -2745,57 +2746,57 @@ 
Fitting spatial factor multi-species N-mixture models with sfMsNMix()<
 Sampling ... 
 Batch: 200 of 800, 25.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     28.0        0.26343
-    1   alpha[1]    36.0        0.37758
-    2   beta[1]     44.0        0.16301
-    2   alpha[1]    44.0        0.23364
-    3   beta[1]     40.0        0.08258
-    3   alpha[1]    44.0        0.14457
-    4   beta[1]     56.0        0.32175
-    4   alpha[1]    56.0        0.64793
-    5   beta[1]     12.0        0.13346
-    5   alpha[1]    36.0        0.21141
-    6   beta[1]     36.0        0.17658
-    6   alpha[1]    24.0        0.30302
-    1   phi     60.0        1.40053
-    2   phi     40.0        0.76863
-    3   phi     64.0        0.47561
+    1   beta[1]     36.0        0.27418
+    1   alpha[1]    48.0        0.39299
+    2   beta[1]     24.0        0.15351
+    2   alpha[1]    64.0        0.22901
+    3   beta[1]     28.0        0.08769
+    3   alpha[1]    44.0        0.13615
+    4   beta[1]     40.0        0.31538
+    4   alpha[1]    44.0        0.59811
+    5   beta[1]     40.0        0.14749
+    5   alpha[1]    60.0        0.22448
+    6   beta[1]     44.0        0.15978
+    6   alpha[1]    32.0        0.27418
+    1   phi     56.0        1.24216
+    2   phi     56.0        1.10170
+    3   phi     36.0        0.52564
 -------------------------------------------------
 Batch: 400 of 800, 50.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     28.0        0.25821
-    1   alpha[1]    40.0        0.38521
-    2   beta[1]     68.0        0.15661
-    2   alpha[1]    40.0        0.26343
-    3   beta[1]     52.0        0.08595
-    3   alpha[1]    40.0        0.13890
-    4   beta[1]     24.0        0.32175
-    4   alpha[1]    48.0        0.63510
-    5   beta[1]     32.0        0.11602
-    5   alpha[1]    40.0        0.20722
-    6   beta[1]     44.0        0.16630
-    6   alpha[1]    44.0        0.30914
-    1   phi     28.0        1.61100
-    2   phi     44.0        0.66821
-    3   phi     44.0        0.50503
+    1   beta[1]     32.0        0.25310
+    1   alpha[1]    48.0        0.37758
+    2   beta[1]     48.0        0.16966
+    2   alpha[1]    40.0        0.23836
+    3   beta[1]     40.0        0.08095
+    3   alpha[1]    12.0        0.13346
+    4   beta[1]     36.0        0.34165
+    4   alpha[1]    44.0        0.62252
+    5   beta[1]     36.0        0.12569
+    5   alpha[1]    40.0        0.21568
+    6   beta[1]     24.0        0.17308
+    6   alpha[1]    32.0        0.29113
+    1   phi     36.0        1.48714
+    2   phi     36.0        0.84947
+    3   phi     40.0        0.65498
 -------------------------------------------------
 Batch: 600 of 800, 75.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     56.0        0.26343
-    1   alpha[1]    44.0        0.39299
-    2   beta[1]     48.0        0.16630
-    2   alpha[1]    36.0        0.23364
-    3   beta[1]     48.0        0.07934
-    3   alpha[1]    52.0        0.13615
-    4   beta[1]     36.0        0.34165
-    4   alpha[1]    36.0        0.64793
-    5   beta[1]     36.0        0.12822
-    5   alpha[1]    56.0        0.22003
-    6   beta[1]     48.0        0.18750
-    6   alpha[1]    24.0        0.30914
-    1   phi     20.0        1.71061
-    2   phi     52.0        0.69548
-    3   phi     40.0        0.48522
+    1   beta[1]     36.0        0.26875
+    1   alpha[1]    32.0        0.39299
+    2   beta[1]     32.0        0.14749
+    2   alpha[1]    32.0        0.22901
+    3   beta[1]     60.0        0.07777
+    3   alpha[1]    44.0        0.13615
+    4   beta[1]     20.0        0.34165
+    4   alpha[1]    44.0        0.64793
+    5   beta[1]     52.0        0.13346
+    5   alpha[1]    48.0        0.22003
+    6   beta[1]     48.0        0.16630
+    6   alpha[1]    40.0        0.28537
+    1   phi     16.0        1.61100
+    2   phi     80.0        0.62930
+    3   phi     36.0        0.49502
 -------------------------------------------------
 Batch: 800 of 800, 100.00%
 ----------------------------------------
@@ -2804,57 +2805,57 @@ 
Fitting spatial factor multi-species N-mixture models with sfMsNMix()<
 Sampling ... 
 Batch: 200 of 800, 25.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     32.0        0.25310
-    1   alpha[1]    44.0        0.39299
-    2   beta[1]     56.0        0.15047
-    2   alpha[1]    44.0        0.24809
-    3   beta[1]     40.0        0.07777
-    3   alpha[1]    44.0        0.12822
-    4   beta[1]     56.0        0.31538
-    4   alpha[1]    32.0        0.67437
-    5   beta[1]     48.0        0.12569
-    5   alpha[1]    36.0        0.20312
-    6   beta[1]     32.0        0.16630
-    6   alpha[1]    64.0        0.29701
-    1   phi     72.0        1.48714
-    2   phi     36.0        0.81616
-    3   phi     32.0        0.46620
+    1   beta[1]     40.0        0.23364
+    1   alpha[1]    40.0        0.38521
+    2   beta[1]     28.0        0.15047
+    2   alpha[1]    36.0        0.21141
+    3   beta[1]     44.0        0.08095
+    3   alpha[1]    40.0        0.13890
+    4   beta[1]     20.0        0.31538
+    4   alpha[1]    56.0        0.66102
+    5   beta[1]     32.0        0.11602
+    5   alpha[1]    32.0        0.22003
+    6   beta[1]     40.0        0.16301
+    6   alpha[1]    44.0        0.32825
+    1   phi     48.0        1.34562
+    2   phi     32.0        0.68171
+    3   phi     60.0        0.53625
 -------------------------------------------------
 Batch: 400 of 800, 50.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     44.0        0.24318
-    1   alpha[1]    48.0        0.39299
-    2   beta[1]     52.0        0.16630
-    2   alpha[1]    36.0        0.24318
-    3   beta[1]     40.0        0.08095
-    3   alpha[1]    72.0        0.13346
-    4   beta[1]     16.0        0.34165
-    4   alpha[1]    24.0        0.64793
-    5   beta[1]     44.0        0.13081
-    5   alpha[1]    52.0        0.19910
-    6   beta[1]     32.0        0.17658
-    6   alpha[1]    36.0        0.31538
-    1   phi     68.0        1.54783
-    2   phi     64.0        0.72387
-    3   phi     44.0        0.50503
+    1   beta[1]     36.0        0.24318
+    1   alpha[1]    44.0        0.41729
+    2   beta[1]     36.0        0.13890
+    2   alpha[1]    48.0        0.22901
+    3   beta[1]     40.0        0.07934
+    3   alpha[1]    24.0        0.13081
+    4   beta[1]     40.0        0.33488
+    4   alpha[1]    48.0        0.61020
+    5   beta[1]     32.0        0.12822
+    5   alpha[1]    36.0        0.22901
+    6   beta[1]     36.0        0.16966
+    6   alpha[1]    36.0        0.29113
+    1   phi     52.0        1.74517
+    2   phi     36.0        0.62930
+    3   phi     48.0        0.48522
 -------------------------------------------------
 Batch: 600 of 800, 75.00%
     Species Parameter   Acceptance  Tuning
-    1   beta[1]     44.0        0.24318
-    1   alpha[1]    40.0        0.37758
-    2   beta[1]     40.0        0.13890
-    2   alpha[1]    44.0        0.24809
-    3   beta[1]     44.0        0.08095
-    3   alpha[1]    44.0        0.13615
-    4   beta[1]     60.0        0.34165
-    4   alpha[1]    32.0        0.67437
-    5   beta[1]     40.0        0.13081
-    5   alpha[1]    52.0        0.22003
-    6   beta[1]     24.0        0.16301
-    6   alpha[1]    36.0        0.29701
-    1   phi     44.0        1.37280
-    2   phi     28.0        0.81616
-    3   phi     32.0        0.51523
+    1   beta[1]     44.0        0.25310
+    1   alpha[1]    32.0        0.38521
+    2   beta[1]     32.0        0.16301
+    2   alpha[1]    40.0        0.24318
+    3   beta[1]     32.0        0.07623
+    3   alpha[1]    36.0        0.14171
+    4   beta[1]     36.0        0.36277
+    4   alpha[1]    36.0        0.66102
+    5   beta[1]     60.0        0.13890
+    5   alpha[1]    56.0        0.21141
+    6   beta[1]     40.0        0.17308
+    6   alpha[1]    56.0        0.30302
+    1   phi     56.0        1.48714
+    2   phi     56.0        0.81616
+    3   phi     32.0        0.46620
 -------------------------------------------------
 Batch: 800 of 800, 100.00%
summary(out.sf.ms)
@@ -2871,82 +2872,82 @@

Fitting spatial factor multi-species N-mixture models with sfMsNMix()< Thinning Rate: 10 Number of Chains: 3 Total Posterior Samples: 3000 -Run Time (min): 13.7853 +Run Time (min): 16.2479 ---------------------------------------- Community Level ---------------------------------------- Abundance Means (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) -0.8823 0.4745 -1.7705 -0.9051 0.1510 1.0123 735 -scale(abund.cov.1) 0.3698 0.3760 -0.3954 0.3668 1.1082 1.0045 2783 +(Intercept) -0.8258 0.4845 -1.7462 -0.8470 0.1382 1.0554 306 +scale(abund.cov.1) 0.3580 0.3766 -0.4288 0.3738 1.1058 1.0000 3000 Abundance Variances (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) 1.3554 2.5078 0.2336 0.8860 4.8125 1.1411 2609 -scale(abund.cov.1) 0.9482 0.8778 0.2204 0.6783 3.4457 1.0108 3000 +(Intercept) 1.3720 1.3989 0.2249 0.9873 4.9608 1.0738 340 +scale(abund.cov.1) 0.9417 0.8606 0.2310 0.6878 3.4506 1.0039 2671 Abundance Random Effect Variances (log scale): - Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept)-abund.factor.1 0.6873 0.1775 0.4101 0.662 1.1034 1.0029 1145 + Mean SD 2.5% 50% 97.5% Rhat ESS +(Intercept)-abund.factor.1 0.7047 0.1902 0.4133 0.6759 1.164 1.0036 1102 Detection Means (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) 0.0767 0.4115 -0.7193 0.0883 0.9009 1.0023 2412 -scale(det.cov.1) 0.5255 0.2367 0.0434 0.5315 0.9862 1.0095 2426 -scale(det.cov.2) 0.6479 0.5891 -0.5326 0.6520 1.8142 1.0131 3000 +(Intercept) 0.0833 0.4115 -0.7712 0.0915 0.8977 1.0072 1907 +scale(det.cov.1) 0.5301 0.2387 0.0253 0.5375 0.9893 1.0015 2791 +scale(det.cov.2) 0.6599 0.6022 -0.5454 0.6579 1.8484 1.0005 3000 Detection Variances (logit scale): - Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept) 1.1506 1.5302 0.2337 0.7984 3.9874 1.0526 1420 -scale(det.cov.1) 0.3106 0.4214 0.0573 0.2079 1.1461 1.0181 3000 -scale(det.cov.2) 2.5933 3.2349 0.5205 1.7670 9.3077 1.0123 3000 + Mean SD 2.5% 50% 97.5% Rhat ESS +(Intercept) 1.0537 1.2432 0.2302 0.7407 3.8911 1.0267 1802 +scale(det.cov.1) 0.3011 0.3085 0.0585 0.2107 1.1383 1.0110 3000 +scale(det.cov.2) 2.6583 3.6453 0.4861 1.7887 10.0362 1.0347 2802 ---------------------------------------- Species Level ---------------------------------------- Abundance (log scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept)-sp1 -2.0800 0.3571 -2.7829 -2.0727 -1.3886 1.0092 369 -(Intercept)-sp2 -0.7772 0.4640 -1.6435 -0.8039 0.2170 1.0956 50 -(Intercept)-sp3 0.2165 0.4707 -0.6972 0.1880 1.2281 1.6901 23 -(Intercept)-sp4 -1.6503 0.3589 -2.3861 -1.6404 -0.9570 1.0116 467 -(Intercept)-sp5 -0.5581 0.3301 -1.1826 -0.5714 0.1620 1.3383 82 -(Intercept)-sp6 -0.7974 0.3060 -1.4087 -0.8072 -0.2131 1.1878 116 -scale(abund.cov.1)-sp1 1.1530 0.1905 0.7840 1.1496 1.5368 1.0285 530 -scale(abund.cov.1)-sp2 -0.7754 0.1190 -1.0141 -0.7744 -0.5488 1.0293 700 -scale(abund.cov.1)-sp3 0.3255 0.0840 0.1623 0.3249 0.4969 1.0060 615 -scale(abund.cov.1)-sp4 0.2873 0.1707 -0.0379 0.2864 0.6134 1.0075 1412 -scale(abund.cov.1)-sp5 1.2779 0.1064 1.0830 1.2726 1.4908 1.0876 307 -scale(abund.cov.1)-sp6 0.0228 0.0858 -0.1432 0.0214 0.1894 1.0002 1855 +(Intercept)-sp1 -2.0528 0.3654 -2.7473 -2.0567 -1.3097 1.0301 324 +(Intercept)-sp2 -0.6454 0.4523 -1.5929 -0.6324 0.2478 1.4014 61 +(Intercept)-sp3 0.3995 0.4926 -0.7762 0.4380 1.2709 1.3368 24 +(Intercept)-sp4 -1.6657 0.3618 -2.3901 -1.6657 -0.9735 1.0283 427 +(Intercept)-sp5 -0.6269 0.3304 -1.2975 -0.6188 -0.0326 1.0335 111 +(Intercept)-sp6 -0.7827 0.3131 -1.4222 -0.7785 -0.1775 1.0010 155 +scale(abund.cov.1)-sp1 1.1448 0.1940 0.7758 1.1391 1.5424 1.0181 604 +scale(abund.cov.1)-sp2 -0.7797 0.1148 -1.0063 -0.7790 -0.5572 1.0019 469 +scale(abund.cov.1)-sp3 0.3258 0.0864 0.1577 0.3267 0.4961 1.0408 417 +scale(abund.cov.1)-sp4 0.2902 0.1681 -0.0377 0.2872 0.6158 1.0007 1622 +scale(abund.cov.1)-sp5 1.2884 0.1000 1.0825 1.2909 1.4813 1.0120 389 +scale(abund.cov.1)-sp6 0.0241 0.0865 -0.1450 0.0235 0.1966 1.0145 1768 Detection (logit scale): Mean SD 2.5% 50% 97.5% Rhat ESS -(Intercept)-sp1 0.0471 0.3332 -0.6643 0.0662 0.6274 1.0770 297 -(Intercept)-sp2 -0.2172 0.2318 -0.6945 -0.2062 0.2211 1.0154 286 -(Intercept)-sp3 0.2798 0.1406 0.0066 0.2792 0.5490 1.0071 318 -(Intercept)-sp4 -0.2212 0.3399 -0.8838 -0.2243 0.4278 1.0065 1226 -(Intercept)-sp5 -0.9047 0.3239 -1.6269 -0.8807 -0.3448 1.1350 88 -(Intercept)-sp6 1.4842 0.2016 1.0866 1.4821 1.8854 1.0027 925 -scale(det.cov.1)-sp1 0.1148 0.2227 -0.3393 0.1142 0.5415 1.0107 1208 -scale(det.cov.1)-sp2 0.5622 0.1332 0.3098 0.5599 0.8305 1.0089 1263 -scale(det.cov.1)-sp3 1.0892 0.0962 0.9050 1.0877 1.2859 1.0050 694 -scale(det.cov.1)-sp4 0.1606 0.3117 -0.4668 0.1742 0.7421 1.0024 1948 -scale(det.cov.1)-sp5 0.6332 0.1318 0.3924 0.6281 0.9049 1.0493 507 -scale(det.cov.1)-sp6 0.6270 0.1373 0.3599 0.6255 0.9038 1.0030 1826 -scale(det.cov.2)-sp1 0.6032 0.2100 0.2251 0.5921 1.0275 1.0205 1183 -scale(det.cov.2)-sp2 1.0852 0.1791 0.7498 1.0760 1.4557 1.0042 1150 -scale(det.cov.2)-sp3 0.8212 0.0889 0.6509 0.8214 0.9982 1.0006 843 -scale(det.cov.2)-sp4 2.9490 0.6783 1.8234 2.8774 4.4670 1.0003 1848 -scale(det.cov.2)-sp5 -0.1116 0.1086 -0.3309 -0.1095 0.0895 1.0582 1094 -scale(det.cov.2)-sp6 -0.8087 0.1544 -1.1211 -0.8050 -0.5192 1.0023 1631 +(Intercept)-sp1 0.0586 0.3263 -0.6158 0.0653 0.6666 1.0326 373 +(Intercept)-sp2 -0.1989 0.2239 -0.6492 -0.1932 0.2210 1.1185 356 +(Intercept)-sp3 0.2868 0.1471 -0.0066 0.2908 0.5669 1.0203 287 +(Intercept)-sp4 -0.2196 0.3291 -0.8781 -0.2160 0.3957 1.0043 1381 +(Intercept)-sp5 -0.8304 0.3096 -1.5580 -0.8072 -0.3172 1.1500 91 +(Intercept)-sp6 1.4897 0.1999 1.1041 1.4863 1.8887 1.0014 960 +scale(det.cov.1)-sp1 0.1175 0.2213 -0.3367 0.1281 0.5412 1.0007 1436 +scale(det.cov.1)-sp2 0.5730 0.1361 0.3056 0.5730 0.8404 1.0004 1392 +scale(det.cov.1)-sp3 1.0918 0.0996 0.8962 1.0928 1.2956 1.0111 607 +scale(det.cov.1)-sp4 0.1795 0.3008 -0.4378 0.1857 0.7404 1.0012 2338 +scale(det.cov.1)-sp5 0.6554 0.1307 0.4091 0.6534 0.9226 1.0279 575 +scale(det.cov.1)-sp6 0.6261 0.1398 0.3600 0.6239 0.9044 1.0048 1985 +scale(det.cov.2)-sp1 0.6069 0.2098 0.2226 0.5969 1.0346 1.0216 1110 +scale(det.cov.2)-sp2 1.1016 0.1784 0.7548 1.0931 1.4652 1.0329 1238 +scale(det.cov.2)-sp3 0.8277 0.0913 0.6515 0.8294 1.0069 1.0249 870 +scale(det.cov.2)-sp4 2.9203 0.6506 1.7648 2.8804 4.3084 1.0171 1442 +scale(det.cov.2)-sp5 -0.1113 0.1068 -0.3244 -0.1093 0.0949 1.0092 1414 +scale(det.cov.2)-sp6 -0.8111 0.1534 -1.1200 -0.8119 -0.5181 1.0053 1500 ---------------------------------------- Spatial Covariance ---------------------------------------- - Mean SD 2.5% 50% 97.5% Rhat ESS -phi-1 28.1327 8.6746 10.9327 28.8698 41.3981 1.0174 372 -phi-2 4.0653 1.6130 2.2102 3.6851 8.1291 1.0731 96 -phi-3 7.1629 2.9746 2.4202 6.8891 13.5079 1.2337 58 + Mean SD 2.5% 50% 97.5% Rhat ESS +phi-1 27.7902 8.8918 9.3575 28.4294 41.2575 1.0348 350 +phi-2 3.9349 1.4117 2.2058 3.6387 7.4255 1.0078 178 +phi-3 7.7100 2.8441 2.8847 7.4518 13.8732 1.0098 130

Posterior predictive checks

@@ -2973,17 +2974,17 @@

Posterior predictive checks

---------------------------------------- Community Level ---------------------------------------- -Bayesian p-value: 0.5083 +Bayesian p-value: 0.4917 ---------------------------------------- Species Level ---------------------------------------- -sp1 Bayesian p-value: 0.5177 -sp2 Bayesian p-value: 0.4647 -sp3 Bayesian p-value: 0.5887 -sp4 Bayesian p-value: 0.4333 -sp5 Bayesian p-value: 0.7383 -sp6 Bayesian p-value: 0.307 +sp1 Bayesian p-value: 0.5417 +sp2 Bayesian p-value: 0.4717 +sp3 Bayesian p-value: 0.532 +sp4 Bayesian p-value: 0.394 +sp5 Bayesian p-value: 0.6267 +sp6 Bayesian p-value: 0.384 Fit statistic: freeman-tukey
@@ -2999,8 +3000,8 @@

Model selection using WAIC

Currently on species 4 out of 6
Currently on species 5 out of 6
Currently on species 6 out of 6
-
     elpd        pD      WAIC 
--1503.607   227.617  3462.448
+
      elpd         pD       WAIC 
+-1502.8609   228.8501  3463.4219 
# Spatial factor model
 waicAbund(out.sf.ms, by.sp = FALSE)
N.max not specified. Setting upper index of integration of N to 10 plus
@@ -3012,7 +3013,7 @@ 
Model selection using WAIC

 
Currently on species 5 out of 6

 
Currently on species 6 out of 6

 
      elpd         pD       WAIC 
--1514.6710   200.4276  3430.1972 

+-1515.3108   200.6227  3431.8671

Prediction

@@ -3026,7 +3027,7 @@

Prediction

out.pred <- predict(out.sf.ms, X.0, ignore.RE = TRUE, include.sp = FALSE)
 str(out.pred)
List of 4
- $ mu.0.samples: num [1:3000, 1:6, 1:100] 0.000462 0.00159 0.002161 0.010415 0.027684 ...
+ $ mu.0.samples: num [1:3000, 1:6, 1:100] 0.003421 0.002539 0.002348 0.00193 0.000919 ...
  $ N.0.samples : int [1:3000, 1:6, 1:100] 0 0 0 0 0 0 0 0 0 0 ...
  $ call        : language predict.sfMsNMix(object = out.sf.ms, X.0 = X.0, ignore.RE = TRUE, include.sp = FALSE)
  $ object.class: chr "sfMsNMix"
@@ -3045,7 +3046,7 @@ 
Prediction

   theme_bw() +
   facet_wrap(vars(sp)) +
   labs(x = 'Covariate', y = 'Expected abundance')
-

+