Function for Fitting Integrated Multi-Species Occupancy Models Using Polya-Gamma Latent Variables

Function for fitting integrated multi-species occupancy models using Polya-Gamma latent variables.

Usage

intMsPGOcc(occ.formula, det.formula, data, inits, priors, n.samples,
           n.omp.threads = 1, verbose = TRUE, n.report = 100, 
           n.burn = round(.10 * n.samples), n.thin = 1, n.chains = 1,
           ...)

Arguments

occ.formula: a symbolic description of the model to be fit for the occurrence portion of the model using R's model syntax. Only right-hand side of formula is specified. See example below. Random intercepts are allowed using lme4 syntax (Bates et al. 2015).
det.formula: a list of symbolic descriptions of the models to be fit for the detection portion of the model using R's model syntax for each data set. Each element in the list is a formula for the detection model of a given data set. Only right-hand side of formula is specified. Random effects are not currently supported. See example below.
data: a list containing data necessary for model fitting. Valid tags are y, occ.covs, det.covs, sites, and species. y is a list of three-dimensional arrays. Each element of the list has first dimension equal to the number of species observed in that data source, second dimension equal to the number of sites observed in that data source, and thir dimension equal to the maximum number of replicates at a given site. occ.covs is a matrix or data frame containing the variables used in the occurrence portion of the model, with the number of rows being the number of sites with at least one data source for each column (variable). det.covs is a list of variables included in the detection portion of the model for each data source. det.covs should have the same number of elements as y, where each element is itself a list. Each element of the list for a given data source is a different detection covariate, which can be site-level or observational-level. Site-level covariates are specified as a vector with length equal to the number of observed sites of that data source, while observational-level covariates are specified as a matrix or data frame with the number of rows equal to the number of observed sites of that data source and number of columns equal to the maximum number of replicates at a given site. sites is a list of site indices with number of elements equal to the number of data sources being modeled. Each element contains a vector of length equal to the number of sites that specific data source contains. Each value in the vector indicates the row in occ.covs that corresponds with the specific row of the detection-nondetection data for the data source. This is used to properly link sites across data sets. species is a list with number of data sources being modeled. Each element of the list is a vector of codes (these can be numeric or character) that indicate the species modeled in the specific data set.
inits: a list with each tag corresponding to a parameter name. Valid tags are alpha.comm, beta.comm, beta, alpha, tau.sq.beta, tau.sq.alpha, sigma.sq.psi, and z. The value portion of each tag is the parameter's initial value. See priors description for definition of each parameter name. Additionally, the tag fix can be set to TRUE to fix the starting values across all chains. If fix is not specified (the default), starting values are varied randomly across chains.
priors: a list with each tag corresponding to a parameter name. Valid tags are beta.comm.normal, alpha.comm.normal, tau.sq.beta.ig, tau.sq.alpha.ig, sigma.sq.psi.ig, and sigma.sq.p.ig. Community-level occurrence (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 corresponding to the mean and variance of the normal distribution, which are each specified as vectors of length equal to the number of coefficients to be estimated or of length one if priors are the same for all coefficients. If not specified, prior means are set to 0 and prior variances set to 2.72. For the community-level detection means (alpha.comm), the mean and variance hyperparameters are themselves passed in as lists, with each element of the list corresponding to the specific hyperparameters for the detection parameters in a given data source. If not specified, prior means are set to 0 and prior variances set to 2.72. Community-level variance parameters for occurrence (tau.sq.beta) and detection (tau.sq.alpha) are assumed to follow an inverse Gamma distribution. For the occurrence parameters, the hyperparameters of the inverse gamma distribution are passed as a list of length two with the first and second elements corresponding to the shape and scale parameters, which are each specified as vectors of length equal to the number of coefficients to be estimated or a single value if all parameters are assigned the same prior. If not specified, prior shape and scale parameters are set to 0.1. For the detection community-level variance parameters (tau.sq.alpha), the shape and scale parameters are passed in as lists, with each element of the list corresponding to the specific hyperparameters for the detection variances in a given data source. sigma.sq.psi and are the random effect variances for any occurrence random effects, respectively, 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 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.
n.samples: the number of posterior samples to collect in each chain.
n.omp.threads: a positive integer indicating the number of threads to use for SMP parallel processing within chains. This will have no impact on model run times for non-spatial models. The package must be compiled for OpenMP support. For most Intel-based machines, we recommend setting n.omp.threads up to the number of hypterthreaded cores. Note, n.omp.threads > 1 might not work on some systems. Currently only relevant for spatial models.
verbose: if TRUE, messages about data preparation, model specification, and progress of the sampler are printed to the screen. Otherwise, no messages are printed.
n.report: the interval to report MCMC progress.
n.burn: the number of samples out of the total n.samples to discard as burn-in for each chain. By default, the first 10% of samples is discarded.
n.thin: the thinning interval for collection of MCMC samples. The thinning occurs after the n.burn samples are discarded. Default value is set to 1.
n.chains: the number of chains to run in sequence.
...: currently no additional arguments

References

Polson, N.G., J.G. Scott, and J. Windle. (2013) Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables. Journal of the American Statistical Association, 108:1339-1349.

Bates, Douglas, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01 .

Dorazio, R. M., and Royle, J. A. (2005). Estimating size and composition of biological communities by modeling the occurrence of species. Journal of the American Statistical Association, 100(470), 389-398.

Doser, J. W., Leuenberger, W., Sillett, T. S., Hallworth, M. T. & Zipkin, E. F. (2022). Integrated community occupancy models: A framework to assess occurrence and biodiversity dynamics using multiple data sources. Methods in Ecology and Evolution, 00, 1-14. doi:10.1111/2041-210X.13811

Author

Jeffrey W. Doser doserjef@msu.edu,

Value

An object of class intMsPGOcc that is a list comprised of:

beta.comm.samples: a coda object of posterior samples for the community level occurrence regression coefficients.
alpha.comm.samples: a coda object of posterior samples for the community level detection regression coefficients for all data sources.
tau.sq.beta.samples: a coda object of posterior samples for the occurrence community variance parameters.
tau.sq.alpha.samples: a coda object of posterior samples for the detection community variance parameters for all data sources.
beta.samples: a coda object of posterior samples for the species level occurrence regression coefficients.
alpha.samples: a coda object of posterior samples for the species level detection regression coefficients for all data sources.
z.samples: a three-dimensional array of posterior samples for the latent occurrence values for each species.
psi.samples: a three-dimensional array of posterior samples for the latent occurrence probability values for each species.
sigma.sq.psi.samples: a coda object of posterior samples for variances of random intercepts included in the occurrence portion of the model. Only included if random intercepts are specified in occ.formula.
beta.star.samples: a coda object of posterior samples for the occurrence random effects. Only included if random intercepts are specified in occ.formula.
like.samples: a three-dimensional array of posterior samples for the likelihood value associated with each site and species. Used for calculating WAIC.
rhat: a list of Gelman-Rubin diagnostic values for some of the model parameters.
ESS: a list of effective sample sizes for some of the model parameters.
run.time: MCMC sampler execution time reported using proc.time().

The return object will include additional objects used for subsequent prediction and/or model fit evaluation.

Examples

set.seed(91)
J.x <- 10
J.y <- 10
# Total number of data sources across the study region
J.all <- J.x * J.y
# Number of data sources.
n.data <- 2
# Sites for each data source.
J.obs <- sample(ceiling(0.2 * J.all):ceiling(0.5 * J.all), n.data, replace = TRUE)
n.rep <- list()
n.rep[[1]] <- rep(3, J.obs[1])
n.rep[[2]] <- rep(4, J.obs[2])

# Number of species observed in each data source
N <- c(8, 3)

# Community-level covariate effects
# Occurrence
beta.mean <- c(0.2, 0.5)
p.occ <- length(beta.mean)
tau.sq.beta <- c(0.4, 0.3)
# Detection
# Detection covariates
alpha.mean <- list()
tau.sq.alpha <- list()
# Number of detection parameters in each data source
p.det.long <- c(4, 3)
for (i in 1:n.data) {
  alpha.mean[[i]] <- runif(p.det.long[i], -1, 1)
  tau.sq.alpha[[i]] <- runif(p.det.long[i], 0.1, 1)
}
# Random effects
psi.RE <- list()
p.RE <- list()
beta <- matrix(NA, nrow = max(N), ncol = p.occ)
for (i in 1:p.occ) {
  beta[, i] <- rnorm(max(N), beta.mean[i], sqrt(tau.sq.beta[i]))
}
alpha <- list()
for (i in 1:n.data) {
  alpha[[i]] <- matrix(NA, nrow = N[i], ncol = p.det.long[i])
  for (t in 1:p.det.long[i]) {
    alpha[[i]][, t] <- rnorm(N[i], alpha.mean[[i]][t], sqrt(tau.sq.alpha[[i]])[t])
  }
}
sp <- FALSE
factor.model <- FALSE
# Simulate occupancy data
dat <- simIntMsOcc(n.data = n.data, J.x = J.x, J.y = J.y,
                   J.obs = J.obs, n.rep = n.rep, N = N, beta = beta, alpha = alpha,
                   psi.RE = psi.RE, p.RE = p.RE, sp = sp, factor.model = factor.model,
                   n.factors = n.factors)
J <- nrow(dat$coords.obs)
y <- dat$y
X <- dat$X.obs
X.p <- dat$X.p
X.re <- dat$X.re.obs
X.p.re <- dat$X.p.re
sites <- dat$sites
species <- dat$species

# Package all data into a list
occ.covs <- cbind(X)
colnames(occ.covs) <- c('int', 'occ.cov.1')
#colnames(occ.covs) <- c('occ.cov')
det.covs <- list()
# Add covariates one by one
det.covs[[1]] <- list(det.cov.1.1 = X.p[[1]][, , 2], 
                      det.cov.1.2 = X.p[[1]][, , 3], 
                      det.cov.1.3 = X.p[[1]][, , 4])
det.covs[[2]] <- list(det.cov.2.1 = X.p[[2]][, , 2], 
                      det.cov.2.2 = X.p[[2]][, , 3]) 

data.list <- list(y = y, 
                  occ.covs = occ.covs, 
                  det.covs = det.covs, 
                  sites = sites, 
                  species = species)
# Take a look at the data.list structure for integrated multi-species
# occupancy models.
# Priors 
prior.list <- list(beta.comm.normal = list(mean = 0,var = 2.73),
                   alpha.comm.normal = list(mean = list(0, 0),
                                            var = list(2.72, 2.72)), 
                   tau.sq.beta.ig = list(a = 0.1, b = 0.1), 
                   tau.sq.alpha.ig = list(a = list(0.1, 0.1), 
                                          b = list(0.1, 0.1)))
inits.list <- list(alpha.comm = list(0, 0), 
                   beta.comm = 0, 
                   tau.sq.beta = 1, 
                   tau.sq.alpha = list(1, 1), 
                   alpha = list(a = matrix(rnorm(p.det.long[1] * N[1]), N[1], p.det.long[1]), 
                                b = matrix(rnorm(p.det.long[2] * N[2]), N[2], p.det.long[2])),
                   beta = 0)

# Fit the model. 
# Note that this is just a test case and more iterations/chains may need to 
# be run to ensure convergence.
out <- intMsPGOcc(occ.formula = ~ occ.cov.1,
                  det.formula = list(f.1 = ~ det.cov.1.1 + det.cov.1.2 + det.cov.1.3,
                                     f.2 = ~ det.cov.2.1 + det.cov.2.2),
                  inits = inits.list,
                  priors = prior.list,
                  data = data.list, 
                  n.samples = 100, 
                  n.omp.threads = 1, 
                  verbose = TRUE, 
                  n.report = 10, 
                  n.burn = 50, 
                  n.thin = 1, 
                  n.chains = 1) 
#> ----------------------------------------
#> 	Preparing to run the model
#> ----------------------------------------
#> z is not specified in initial values.
#> Setting initial values based on observed data
#> ----------------------------------------
#> 	Model description
#> ----------------------------------------
#> Integrated Multispecies Occupancy Model with Polya-Gamma latent
#> variable fit with 57 sites and 8 species.
#> 
#> Integrating 2 occupancy data sets.
#> 
#> Samples per Chain: 100 
#> Burn-in: 50 
#> Thinning Rate: 1 
#> Number of Chains: 1 
#> Total Posterior Samples: 50 
#> 
#> Source compiled with OpenMP support and model fit using 1 thread(s).
#> 
#> ----------------------------------------
#> 	Chain 1
#> ----------------------------------------
#> Sampling ... 
#> Sampled: 10 of 100, 10.00%
#> -------------------------------------------------
#> Sampled: 20 of 100, 20.00%
#> -------------------------------------------------
#> Sampled: 30 of 100, 30.00%
#> -------------------------------------------------
#> Sampled: 40 of 100, 40.00%
#> -------------------------------------------------
#> Sampled: 50 of 100, 50.00%
#> -------------------------------------------------
#> Sampled: 60 of 100, 60.00%
#> -------------------------------------------------
#> Sampled: 70 of 100, 70.00%
#> -------------------------------------------------
#> Sampled: 80 of 100, 80.00%
#> -------------------------------------------------
#> Sampled: 90 of 100, 90.00%
#> -------------------------------------------------
#> Sampled: 100 of 100, 100.00%
summary(out, level = 'community')
#> 
#> Call:
#> intMsPGOcc(occ.formula = ~occ.cov.1, det.formula = list(f.1 = ~det.cov.1.1 + 
#>     det.cov.1.2 + det.cov.1.3, f.2 = ~det.cov.2.1 + det.cov.2.2), 
#>     data = data.list, inits = inits.list, priors = prior.list, 
#>     n.samples = 100, n.omp.threads = 1, verbose = TRUE, n.report = 10, 
#>     n.burn = 50, n.thin = 1, n.chains = 1)
#> 
#> Samples per Chain: 100
#> Burn-in: 50
#> Thinning Rate: 1
#> Number of Chains: 1
#> Total Posterior Samples: 50
#> Run Time (min): 0.0014
#> 
#> ----------------------------------------
#> 	Community Level
#> ----------------------------------------
#> Occurrence Means (logit scale): 
#>                Mean     SD    2.5%     50%  97.5% Rhat ESS
#> (Intercept) -0.3109 0.2390 -0.7323 -0.3295 0.2033   NA  22
#> occ.cov.1    0.6516 0.2269  0.1991  0.6525 1.1044   NA  32
#> 
#> Occurrence Variances (logit scale): 
#>               Mean     SD   2.5%    50%  97.5% Rhat ESS
#> (Intercept) 0.2679 0.2573 0.0238 0.2026 0.9775   NA  17
#> occ.cov.1   0.3124 0.2167 0.0836 0.2630 0.6550   NA  50
#> 
#> -----------------------------
#> 	Data source 1
#> -----------------------------
#> Detection Means (logit scale): 
#>                Mean     SD    2.5%     50%   97.5% Rhat ESS
#> (Intercept)  0.0946 0.3880 -0.6708  0.1426  0.9458   NA  50
#> det.cov.1.1 -0.0395 0.2637 -0.6586  0.0049  0.3893   NA  27
#> det.cov.1.2  1.2311 0.3277  0.7509  1.1615  1.8435   NA  22
#> det.cov.1.3 -0.8574 0.4296 -1.5549 -0.8380 -0.1302   NA  26
#> 
#> Detection Variances (logit scale): 
#>               Mean     SD   2.5%    50%  97.5% Rhat ESS
#> (Intercept) 1.6672 1.1805 0.1461 1.4077 4.3266   NA  11
#> det.cov.1.1 0.5857 0.7484 0.1050 0.2943 2.6104   NA  20
#> det.cov.1.2 0.5243 0.5657 0.0661 0.3538 2.5072   NA  33
#> det.cov.1.3 1.0231 0.9256 0.2351 0.7942 3.1089   NA  15
#> 
#> -----------------------------
#> 	Data source 2
#> -----------------------------
#> Detection Means (logit scale): 
#>                Mean     SD    2.5%     50%  97.5% Rhat ESS
#> (Intercept) -0.0430 0.7845 -1.4123 -0.1250 1.4883   NA  50
#> det.cov.2.1 -0.8899 0.9956 -2.4205 -1.0049 1.3380   NA  50
#> det.cov.2.2  0.4634 0.4806 -0.3250  0.4960 1.2488   NA  50
#> 
#> Detection Variances (logit scale): 
#>               Mean      SD   2.5%    50%   97.5% Rhat ESS
#> (Intercept) 2.9028  4.9285 0.1832 1.1772 16.0598   NA  50
#> det.cov.2.1 9.6785 16.3462 0.2914 3.8987 42.6686   NA  32
#> det.cov.2.2 1.7419  5.9932 0.1117 0.5554  6.5815   NA  50
#>