Package 'spAbundance'

Description

Function for fitting univariate abundance generalized linear (mixed) models

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, ...)

Arguments

Value

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

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

Author(s)

Jeffrey W. Doser [email protected],
Andrew O. Finley [email protected],

References

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.

Examples

set.seed(1010)
J.x <- 15
J.y <- 15
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))
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')

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])

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)

tuning <- list(kappa = 0.2, beta = 0.1, beta.star = 0.2)
n.batch <- 5
batch.length <- 25
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) 
summary(out)

Description

Count data on 6 warblers in Pennsylvania, USA in 2018. Data come from the North American Breeding Bird Survey. Data indicate the total numer of individuals for each of 6 species counted at 50 stops along a 40km route (95 routes in the data set). The six species included in the data set are: (1) American Redstart (AMRE); (2) Blackburnian Warbler (BLBW); (3) Black-throated Blue Warbler; (4) Black-throated Green Warbler; (5) Hooded Warbler; and (6) Magnolia Warbler. Covariate data include three bioclimatic variables derived from PRISM. Two landcover variables (forest cover and developed land cover) come from USGS Earth Resources Observation and Science (EROS) Center.

Usage

data(bbsData)

Format

bbsData is a list with three elements:

y: a two-dimensional matrix of count data with rows corresponding to species (6) and columns corresponding to sites (95).

covs: a data frame with 95 rows and 8 columns consisting of covariates for use in modeling relative abundance.

coords: a numeric matrix with 95 rows and two columns containing the site coordinates. The proj4string is "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

Source

U.S. Geological Survey. Downloaded from https://www.sciencebase.gov/catalog/item/52b1dfa8e4b0d9b325230cd9 on August 25, 2023.

References

Daly, C., Halbleib, M., Smith, J. I., Gibson, W. P., Doggett, M. K., Taylor, G. H., Curtis, J., and Pasteris, P. P. (2008). Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous united states. International Journal of Climatology: a Journal of the Royal Meteorological Society, 28(15):2031–2064

Ziolkowski Jr., D.J., Lutmerding, M., English, W.B., Aponte, V.I., and Hudson, M-A.R., 2023, North American Breeding Bird Survey Dataset 1966 - 2022: U.S. Geological Survey data release, https://doi.org/10.5066/P9GS9K64.

Sohl, T., Reker, R., Bouchard, M., Sayler, K., Dornbierer, J., Wika, S., ... & Friesz, A. (2016). Modeled historical land use and land cover for the conterminous United States. Journal of Land Use Science, 11(4), 476-499.

Description

Bioclimatic and land cover variables extracted at a 12km resolution across the state of Pennsylvania, USA for use in predicting relative abundance of six warbler species across the state. Land cover data come from USGS EROS, while climate data come from PRISM.

Usage

data(bbsPredData)

Format

bbsPredData is a data frame with seven columns:

bio2: bioclim variable 2.

bio8: bioclim variable 8.

bio18: bioclim variable 18.

forest: proportion of forest cover within a 5km radius.

devel: proportion of developed land cover within a 5km radius.

x: the x coordinate of the point. The proj4strig is "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs".

y: the y coordinate of the point. The proj4string is "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs".

References

Daly, C., Halbleib, M., Smith, J. I., Gibson, W. P., Doggett, M. K., Taylor, G. H., Curtis, J., and Pasteris, P. P. (2008). Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous united states. International Journal of Climatology: a Journal of the Royal Meteorological Society, 28(15):2031–2064

Sohl, T., Reker, R., Bouchard, M., Sayler, K., Dornbierer, J., Wika, S., ... & Friesz, A. (2016). Modeled historical land use and land cover for the conterminous United States. Journal of Land Use Science, 11(4), 476-499.

Description

A simulated data set of repeated count data for 6 species across 225 sites with a maximum of 3 replicate surveys performed at a given site.

Usage

data(dataNMixSim)

Format

dataNMixSim is a list with four elements:

y: a three-dimensional array of count data with dimensions of species (6), sites (225) and replicates (3).

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.

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

coords: a numeric matrix with 225 rows and two columns containing the 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).

Examples

set.seed(6)
J.x <- 15
J.y <- 15
J <- J.x * J.y
n.rep <- sample(3, size = J, replace = TRUE)
# n.rep <- rep(5, J) 
n.sp <- 6
# Community-level covariate effects
# Occurrence
beta.mean <- c(-1, 0.5)
p.abund <- length(beta.mean)
tau.sq.beta <- c(0.4, 1.2)
# Detection
alpha.mean <- c(0, 0.5, 0.8)
tau.sq.alpha <- c(0.2, 1, 1.5)
p.det <- length(alpha.mean)
# Random effects
mu.RE <- list()
mu.RE <- list(levels = c(10),
	       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)
alpha <- matrix(NA, nrow = n.sp, ncol = p.det)
for (i in 1:p.abund) {
  beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i]))
}
for (i in 1:p.det) {
  alpha[, i] <- rnorm(n.sp, alpha.mean[i], sqrt(tau.sq.alpha[i]))
}
alpha.true <- alpha
sp <- TRUE 
n.factors <- 3
factor.model <- TRUE
phi <- runif(n.factors, 3/1, 3 / .2)
kappa <- runif(n.sp, 0.1, 1)
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)

y <- dat$y
X <- dat$X
X.p <- dat$X.p
X.re <- dat$X.re
X.p.re <- dat$X.p.re
coords <- dat$coords
dimnames(coords)[[2]] <- c('X', 'Y')

# Package all data into a list
abund.covs <- cbind(X, X.re)
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]) 
dataNMixSim <- list(y = y, 
		    abund.covs = abund.covs, 
		    det.covs = det.covs, 
                    coords = coords)

Description

Function for fitting single-sepcies hierarchical distance sampling models

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, ...)

Arguments

Value

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

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

Author(s)

Jeffrey W. Doser [email protected],

References

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.

Royle, J. A., Dawson, D. K., & Bates, S. (2004). Modeling abundance effects in distance sampling. Ecology, 85(6), 1591-1597.

Examples

set.seed(123)
J.x <- 10
J.y <- 10
J <- J.x * J.y
# Number of distance bins from which to simulate data.
n.bins <- 5
# Length of each bin. This should be of length n.bins
bin.width <- c(.10, .10, .20, .3, .1)
# Abundance coefficients
beta <- c(1.0, 0.2, 0.3, -0.2)
p.abund <- length(beta)
# Detection coefficients
alpha <- c(-1.0, -0.3)
p.det <- length(alpha)
# Detection decay function
det.func <- 'halfnormal'
mu.RE <- list()
p.RE <- list()
sp <- FALSE
family <- 'NB'
kappa <- 0.1
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,
             offset = offset, transect = transect)

y <- dat$y
X <- dat$X
X.re <- dat$X.re
X.p <- dat$X.p
X.p.re <- dat$X.p.re
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')

data.list <- list(y = y,
                  covs = covs,
                  dist.breaks = dist.breaks,
                  offset = offset)

# Priors
prior.list <- list(beta.normal = list(mean = 0, var = 10),
                   alpha.normal = list(mean = 0,
                                       var = 10),
                   kappa.unif = c(0, 100))
# Starting values
inits.list <- list(alpha = 0,
                   beta = 0,
                   kappa = 1)
# Tuning values
tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1,
               kappa = 0.2)

out <- DS(abund.formula = ~ abund.cov.1 + abund.cov.2 + abund.cov.3,
          det.formula = ~ det.cov.1,
          data = data.list,
          n.batch = 10,
          batch.length = 25,
          inits = inits.list,
          family = 'NB',
          det.func = 'halfnormal',
          transect = 'point',
          tuning = tuning,
          priors = prior.list,
          accept.rate = 0.43,
          n.omp.threads = 1,
          verbose = TRUE,
          n.report = 100,
          n.burn = 100,
          n.thin = 1,
          n.chains = 1)
summary(out)

Description

Method for extracting model fitted values from a fitted GLMM (abund).

Usage

## S3 method for class 'abund'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values for fitted model objects of class abund.

Value

A three-dimensional numeric array of fitted values for use in Goodness of Fit assessments. Array dimensions correspond to MCMC samples, sites, and replicates

Description

Method for extracting model fitted values and cell-specific detection probabilities from a hierarchical distance sampling (DS) model.

Usage

## S3 method for class 'DS'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probabilities for fitted model objects of class DS.

Value

A list comprised of:

Description

Method for extracting model fitted values from a fitted latent factor multivariate GLMM (lfMsAbund).

Usage

## S3 method for class 'lfMsAbund'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values for fitted model objects of class lfMsAbund.

Value

A four-dimensional numeric array of fitted values for use in Goodness of Fit assessments. Array dimensions correspond to MCMC samples, species, sites, and replicates.

Description

Method for extracting model fitted values and cell-specific detection probabilities from a latent factor multi-species hierarchical distance sampling (lfMsDS) model.

Usage

## S3 method for class 'lfMsDS'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probabilities for fitted model objects of class lfMsDS.

Value

A list comprised of:

Description

Method for extracting model fitted values and detection probability values from a fitted latent factor multi-species N-mixture (lfMsNMix) model.

Usage

## S3 method for class 'lfMsNMix'
fitted(object, type = 'marginal', ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probability values for fitted model objects of class lfMsNMix.

Value

A list comprised of:

Description

Method for extracting model fitted values from a fitted multivariate GLMM (msAbund).

Usage

## S3 method for class 'msAbund'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values for fitted model objects of class msAbund.

Value

A four-dimensional numeric array of fitted values for use in Goodness of Fit assessments. Array dimensions correspond to MCMC samples, species, sites, and replicates.

Description

Method for extracting model fitted values and cell-specific detection probabilities from a multi-species hierarchical distance sampling (msDS) model.

Usage

## S3 method for class 'msDS'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probabilities for fitted model objects of class msDS.

Value

A list comprised of:

Description

Method for extracting model fitted values and detection probability values from a fitted multi-species N-mixture (msNMix) model.

Usage

## S3 method for class 'msNMix'
fitted(object, type = 'marginal', ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probability values for fitted model objects of class msNMix.

Value

A list comprised of:

Description

Method for extracting model fitted values and detection probabilities from a fitted N-mixture (NMix) model.

Usage

## S3 method for class 'NMix'
fitted(object, type = 'marginal', ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probabilities for fitted model objects of class NMix.

Value

A list comprised of:

Description

Method for extracting model fitted values from a fitted spatial factor multivariate GLMM (sfMsAbund).

Usage

## S3 method for class 'sfMsAbund'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values for fitted model objects of class sfMsAbund.

Value

A four-dimensional numeric array of fitted values for use in Goodness of Fit assessments. Array dimensions correspond to MCMC samples, species, sites, and replicates.

Description

Method for extracting model fitted values and cell-specific detection probabilities from a spatial factor multi-species hierarchical distance sampling (sfMsDS) model.

Usage

## S3 method for class 'sfMsDS'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probabilities for fitted model objects of class sfMsDS.

Value

A list comprised of:

Description

Method for extracting model fitted values and detection probability values from a fitted spatial factor multi-species N-mixture (sfMsNMix) model.

Usage

## S3 method for class 'sfMsNMix'
fitted(object, type = 'marginal', ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probability values for fitted model objects of class sfMsNMix.

Value

A list comprised of:

Description

Method for extracting model fitted values from a fitted spatial GLMM (spAbund).

Usage

## S3 method for class 'spAbund'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values for fitted model objects of class spAbund.

Value

A three-dimensional numeric array of fitted values for use in Goodness of Fit assessments. Array dimensions correspond to MCMC samples, sites, and replicates

Description

Method for extracting model fitted values and cell-specific detection probabilities from a spatial hierarchical distance sampling (spDS) model.

Usage

## S3 method for class 'spDS'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probabilities for fitted model objects of class spDS.

Value

A list comprised of:

Description

Method for extracting model fitted values and detection probabilities from a fitted spatial N-mixture (spNMix) model.

Usage

## S3 method for class 'spNMix'
fitted(object, type = 'marginal', ...)

Arguments

Details

A method to the generic fitted function to extract fitted values and detection probabilities for fitted model objects of class spNMix.

Value

A list comprised of:

Description

Method for extracting model fitted values from a fitted spatially-varying coefficient GLMM (svcAbund).

Usage

## S3 method for class 'svcAbund'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values for fitted model objects of class svcAbund.

Value

A three-dimensional numeric array of fitted values for use in Goodness of Fit assessments. Array dimensions correspond to MCMC samples, sites, and replicates

Description

Method for extracting model fitted values from a fitted multivatiate spatially-varying coefficient GLMM (svcMsAbund).

Usage

## S3 method for class 'svcMsAbund'
fitted(object, ...)

Arguments

Details

A method to the generic fitted function to extract fitted values for fitted model objects of class svcMsAbund.

Value

A four-dimensional numeric array of fitted values for use in Goodness of Fit assessments. Array dimensions correspond to MCMC samples, species, sites, and replicates.

Description

Repeated count data of 12 foliage gleaning bird species in 2015 in the Hubbard Brook Experimental Forest (HBEF) in New Hampshire, USA. Data were collected at 373 sites over three replicate point counts each of 10 minutes in length, with a detection radius of 100m. Some sites were not visited for all three replicates. The 12 species included in the data set are as follows: (1) AMRE: American Redstart; (2) BAWW: Black-and-white Warbler; (3) BHVI: Blue-headed Vireo; (4) BLBW: Blackburnian Warbler; (5) BLPW: Blackpoll Warbler; (6) BTBW: Black-throated Blue Warbler; (7) BTNW: BLack-throated Green Warbler; (8) CAWA: Canada Warbler; (9) MAWA: Magnolia Warbler; (10) NAWA: Nashville Warbler; (11) OVEN: Ovenbird; (12) REVI: Red-eyed Vireo.

Usage

data(hbefCount2015)

Format

hbefCount2015 is a list with four elements:

y: a three-dimensional array of count data with dimensions of species (12), sites (373) and replicates (3).

abund.covs: a data frame with 373 rows and one column consisting of the elevation at each site.

det.covs: a list of two numeric matrices with 373 rows and 3 columns. The first element is the day of year when the survey was conducted for a given site and replicate. The second element is the time of day when the survey was conducted.

coords: a numeric matrix with 373 rows and two columns containing the site coordinates (Easting and Northing) in UTM Zone 19. The proj4string is "+proj=utm +zone=19 +units=m +datum=NAD83".

Source

Rodenhouse, N. and S. Sillett. 2019. Valleywide Bird Survey, Hubbard Brook Experimental Forest, 1999-2016 (ongoing) ver 3. Environmental Data Initiative. doi:10.6073/pasta/faca2b2cf2db9d415c39b695cc7fc217 (Accessed 2021-09-07)

References

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

Description

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.

Usage

lfMsAbund(formula, data, inits, priors, tuning, 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

Value

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

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

Author(s)

Jeffrey W. Doser [email protected],
Andrew O. Finley [email protected],

References

Roberts, G.O. and Rosenthal J.S. (2009) Examples of adaptive MCMC. Journal of Computational and Graphical Statistics, 18(2):349-367.

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.

Examples

set.seed(408)
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(3, size = J, replace = TRUE)
n.sp <- 6
# Community-level covariate effects
beta.mean <- c(-2, 0.5)
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))
# 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 <- FALSE
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')

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)
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)
tuning.list <- list(kappa = 0.3, beta = 0.1, beta.star = 0.1, 
                    lambda = 0.5, w = 0.5)

# Small
n.batch <- 2
batch.length <- 25
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)
summary(out)

Description

Function for fitting latent factor multi-species hierarchical distance sampling models.

Usage

lfMsDS(abund.formula, det.formula, data, inits, priors,  
       tuning, n.factors, 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

Value

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

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

Author(s)

Jeffrey W. Doser [email protected],

References

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.

Royle, J. A. (2004). N‐mixture models for estimating population size from spatially replicated counts. Biometrics, 60(1), 108-115.

Sollmann, R., Gardner, B., Williams, K. A., Gilbert, A. T., & Veit, R. R. (2016). A hierarchical distance sampling model to estimate abundance and covariate associations of species and communities. Methods in Ecology and Evolution, 7(5), 529-537.

Examples

set.seed(210)
J.x <- 10
J.y <- 10 
J <- J.x * J.y
# Number of distance bins from which to simulate data. 
n.bins <- 5
# Length of each bin. This should be of length n.bins
bin.width <- c(.10, .10, .20, .3, .1)
# Number of species
n.sp <- 5
# Community-level abundance coefficients
beta.mean <- c(-1, 0.2, 0.3, -0.2)
p.abund <- length(beta.mean)
tau.sq.beta <- c(0.2, 0.3, 0.5, 0.4)
# Detection coefficients
alpha.mean <- c(-1.0, -0.3)
p.det <- length(alpha.mean)
tau.sq.alpha <- c(0.1, 0.2)
# Detection decay function
det.func <- 'halfnormal'
mu.RE <- list()
p.RE <- list()
# Draw species-level effects from community means.
beta <- matrix(NA, nrow = n.sp, ncol = p.abund)
alpha <- matrix(NA, nrow = n.sp, ncol = p.det)
for (i in 1:p.abund) {
  beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i]))
}
for (i in 1:p.det) {
  alpha[, i] <- rnorm(n.sp, alpha.mean[i], sqrt(tau.sq.alpha[i]))
}
sp <- FALSE 
family <- 'Poisson'
kappa <- runif(n.sp, 0.3, 3) 
offset <- pi * .8^2
transect <- 'line'
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, 
               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, 
               n.factors = n.factors)

y <- dat$y
X <- dat$X
X.p <- dat$X.p
coords <- dat$coords
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')

data.list <- list(y = y, 
                  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), 
                   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)

tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1, 
               kappa = 0.8, lambda = 1, w = 1) 

n.batch <- 4 
batch.length <- 25
n.burn <- 0
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) 
summary(out, level = 'community')

Description

Function for fitting multi-species N-mixture models with species correlations (i.e., an abundance-based joint species distribution model with imperfect detection). We use a factor modeling approach for dimension reduction.

Usage

lfMsNMix(abund.formula, det.formula, data, inits, priors,
         tuning, 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, ...)

Arguments

Value

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

The return object will include additional objects used for subsequent prediction and/or model fit evaluation. Note that detection probability estimated values are not included in the model object, but can be extracted using fitted().

Author(s)

Jeffrey W. Doser [email protected],

References

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.

Royle, J. A. (2004). N‐mixture models for estimating population size from spatially replicated counts. Biometrics, 60(1), 108-115.

Yamaura, Y., Royle, J. A., Shimada, N., Asanuma, S., Sato, T., Taki, H., & Makino, S. I. (2012). Biodiversity of man-made open habitats in an underused country: a class of multispecies abundance models for count data. Biodiversity and Conservation, 21(6), 1365-1380.

Examples

set.seed(408)
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(5, size = J, replace = TRUE)
n.sp <- 6
# Community-level covariate effects
# Abundance
beta.mean <- c(0, 0.5)
p.abund <- length(beta.mean)
tau.sq.beta <- c(0.2, 1.2)
# Detection
alpha.mean <- c(0, 0.5, 0.8)
tau.sq.alpha <- c(0.2, 1, 1.5)
p.det <- length(alpha.mean)
# Random effects
mu.RE <- list()
p.RE <- list()
# Draw species-level effects from community means.
beta <- matrix(NA, nrow = n.sp, ncol = p.abund)
alpha <- matrix(NA, nrow = n.sp, ncol = p.det)
for (i in 1:p.abund) {
  beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i]))
}
for (i in 1:p.det) {
  alpha[, i] <- rnorm(n.sp, alpha.mean[i], sqrt(tau.sq.alpha[i]))
}
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,
                 mu.RE = mu.RE, p.RE = p.RE, sp = FALSE, family = 'Poisson',
                 factor.model = TRUE, n.factors = n.factors)

y <- dat$y
X <- dat$X
X.p <- dat$X.p
X.re <- dat$X.re
X.p.re <- dat$X.p.re
coords <- dat$coords

# Package all data into a list
abund.covs <- X
colnames(abund.covs) <- c('int', 'abund.cov.1')
det.covs <- list(det.cov.1 = as.data.frame(X.p[, , 2]),
                 det.cov.2 = as.data.frame(X.p[, , 3]))
data.list <- list(y = y,
                  abund.covs = abund.covs,
                  det.covs = det.covs,
                  coords = coords)
prior.list <- list(beta.comm.normal = list(mean = rep(0, p.abund),
                                           var = rep(100, p.abund)),
                   alpha.comm.normal = list(mean = rep(0, p.det),
                                            var = rep(2.72, p.det)),
                   tau.sq.beta.ig = list(a = 0.1, b = 0.1),
                   tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
inits.list <- list(beta.comm = 0, alpha.comm = 0,
                   beta = 0, alpha = 0,
                   tau.sq.beta = 0.5, tau.sq.alpha = 0.5,
                   N = apply(y, c(1, 2), max, na.rm = TRUE))
tuning.list <- list(beta = 0.5, alpha = 0.5, lambda = 0.5, w = 0.5)

n.batch <- 4
batch.length <- 25
n.burn <- 0
n.thin <- 1
n.chains <- 1

out <- lfMsNMix(abund.formula = ~ abund.cov.1,
                det.formula = ~ det.cov.1 + det.cov.2,
                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 = n.factors,
                verbose = TRUE,
                n.report = 1,
                n.burn = n.burn,
                n.thin = n.thin,
                n.chains = n.chains)

summary(out, level = 'community')

Description

The function 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.burn = round(.10 * n.batch * batch.length), n.thin = 1, n.chains = 1,
        save.fitted = TRUE, ...)

Arguments

Value

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

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

Author(s)

Jeffrey W. Doser [email protected],
Andrew O. Finley [email protected],

References

Roberts, G.O. and Rosenthal J.S. (2009) Examples of adaptive MCMC. Journal of Computational and Graphical Statistics, 18(2):349-367.

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.

Examples

set.seed(408)
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(3, size = J, replace = TRUE)
n.sp <- 6
# Community-level covariate effects
beta.mean <- c(-2, 0.5)
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))
# 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 <- 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')

y <- dat$y
X <- dat$X
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)
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)
tuning.list <- list(kappa = 0.3, beta = 0.1, beta.star = 0.1)

# Small
n.batch <- 2
batch.length <- 25
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)
summary(out)

Description

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, ...)

Arguments

Value

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

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

Author(s)

Jeffrey W. Doser [email protected],

References

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.

Royle, J. A. (2004). N‐mixture models for estimating population size from spatially replicated counts. Biometrics, 60(1), 108-115.

Sollmann, R., Gardner, B., Williams, K. A., Gilbert, A. T., & Veit, R. R. (2016). A hierarchical distance sampling model to estimate abundance and covariate associations of species and communities. Methods in Ecology and Evolution, 7(5), 529-537.

Examples

set.seed(210)
J.x <- 10
J.y <- 10 
J <- J.x * J.y
# Number of distance bins from which to simulate data. 
n.bins <- 5
# Length of each bin. This should be of length n.bins
bin.width <- c(.10, .10, .20, .3, .1)
# Number of species
n.sp <- 5
# Community-level abundance coefficients
beta.mean <- c(-1, 0.2, 0.3, -0.2)
p.abund <- length(beta.mean)
tau.sq.beta <- c(0.2, 0.3, 0.5, 0.4)
# Detection coefficients
alpha.mean <- c(-1.0, -0.3)
p.det <- length(alpha.mean)
tau.sq.alpha <- c(0.1, 0.2)
# Detection decay function
det.func <- 'halfnormal'
mu.RE <- list()
p.RE <- list()
# Draw species-level effects from community means.
beta <- matrix(NA, nrow = n.sp, ncol = p.abund)
alpha <- matrix(NA, nrow = n.sp, ncol = p.det)
for (i in 1:p.abund) {
  beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i]))
}
for (i in 1:p.det) {
  alpha[, i] <- rnorm(n.sp, alpha.mean[i], sqrt(tau.sq.alpha[i]))
}
sp <- FALSE 
family <- 'Poisson'
kappa <- runif(n.sp, 0.3, 3) 
offset <- pi * .8^2
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, 
               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)

y <- dat$y
X <- dat$X
X.re <- dat$X.re
X.p <- dat$X.p
X.p.re <- dat$X.p.re
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')

data.list <- list(y = y, 
                  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), 
                   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)

tuning <- list(beta = 0.1, alpha = 0.1, beta.star = 0.3, alpha.star = 0.1, 
               kappa = 0.8) 

n.batch <- 4 
batch.length <- 25
n.burn <- 0
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) 
summary(out, level = 'community')

Description

Function for fitting multi-species N-mixture models.

Usage

msNMix(abund.formula, det.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, ...)

Arguments

Value

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

The return object will include additional objects used for subsequent prediction and/or model fit evaluation. Note that detection probability estimated values are not included in the model object, but can be extracted using fitted().

Author(s)

Jeffrey W. Doser [email protected],

References

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.

Royle, J. A. (2004). N‐mixture models for estimating population size from spatially replicated counts. Biometrics, 60(1), 108-115.

Yamaura, Y., Royle, J. A., Shimada, N., Asanuma, S., Sato, T., Taki, H., & Makino, S. I. (2012). Biodiversity of man-made open habitats in an underused country: a class of multispecies abundance models for count data. Biodiversity and Conservation, 21(6), 1365-1380.

Examples

set.seed(408)
J.x <- 8
J.y <- 8
J <- J.x * J.y
n.rep <- sample(5, size = J, replace = TRUE)
n.sp <- 6
# Community-level covariate effects
# Abundance
beta.mean <- c(0, 0.5)
p.abund <- length(beta.mean)
tau.sq.beta <- c(0.2, 1.2)
# Detection
alpha.mean <- c(0, 0.5, 0.8)
tau.sq.alpha <- c(0.2, 1, 1.5)
p.det <- length(alpha.mean)
# Random effects
mu.RE <- list()
p.RE <- list()
# Draw species-level effects from community means.
beta <- matrix(NA, nrow = n.sp, ncol = p.abund)
alpha <- matrix(NA, nrow = n.sp, ncol = p.det)
for (i in 1:p.abund) {
  beta[, i] <- rnorm(n.sp, beta.mean[i], sqrt(tau.sq.beta[i]))
}
for (i in 1:p.det) {
  alpha[, i] <- rnorm(n.sp, alpha.mean[i], sqrt(tau.sq.alpha[i]))
}

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 = FALSE, family = 'Poisson')

y <- dat$y
X <- dat$X
X.p <- dat$X.p
X.re <- dat$X.re
X.p.re <- dat$X.p.re

# Package all data into a list
abund.covs <- X
colnames(abund.covs) <- c('int', 'abund.cov.1')
det.covs <- list(det.cov.1 = as.data.frame(X.p[, , 2]),
                 det.cov.2 = as.data.frame(X.p[, , 3]))
data.list <- list(y = y,
                  abund.covs = abund.covs,
                  det.covs = det.covs)
prior.list <- list(beta.comm.normal = list(mean = rep(0, p.abund),
                                           var = rep(100, p.abund)),
                   alpha.comm.normal = list(mean = rep(0, p.det),
                                            var = rep(2.72, p.det)),
                   tau.sq.beta.ig = list(a = 0.1, b = 0.1),
                   tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
inits.list <- list(beta.comm = 0, alpha.comm = 0,
                   beta = 0, alpha = 0,
                   tau.sq.beta = 0.5, tau.sq.alpha = 0.5,
                   N = apply(y, c(1, 2), max, na.rm = TRUE))
tuning.list <- list(beta = 0.5, alpha = 0.5)

n.batch <- 4
batch.length <- 25
n.burn <- 0
n.thin <- 1
n.chains <- 1

out <- msNMix(abund.formula = ~ abund.cov.1,
              det.formula = ~ det.cov.1 + det.cov.2,
              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)

summary(out, level = 'community')

Description

Distance sampling data of 16 bird species in 2018 in the Disney Wilderness Preserve in Florida, USA. These data were collected as part of the National Ecological Observatory Network (NEON). Data were collected at 90 sites where observers recorded the number of all bird species observed during a six minute, unlimited radius point count survey once 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) 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

data(neonDWP)

Format

neonDWP is a list with five elements:

y: a three-dimensional array of distance sampling data with dimensions of species (16), sites (90) and distance bin (4).

covs: a data frame with 90 rows and four columns consisting of covariates for use in modeling abundance and/or detection.

dist.breaks: a vector of five values indicating the break points of the four distance bands in the data.

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.

coords: a numeric matrix with 373 rows and two columns containing the site coordinates (Easting and Northing) in UTM Zone 17N. The EPSG is 32617.

Note

These data were updated on November 28, 2023 (v0.1.1), and thus any analyses run with the data from v0.1.0 will not give the same results. The data were updated after NEON discovered an error in how the bird count data was associated with the spatial locations of each site. These data have been updated following the procedures outlined at https://www.neonscience.org/impact/observatory-blog/bird-point-ids-within-grids-were-transposed-resulting-inaccurate-point.

Source

NEON (National Ecological Observatory Network). Breeding landbird point counts (DP1.10003.001), RELEASE-2023. https://doi.org/10.48443/00pg-vm19. Dataset accessed from https://data.neonscience.org on May 25, 2023

References

Barnett, D. T., Duffy, P. A., Schimel, D. S., Krauss, R. E., Irvine, K. M., Davis, F. W.,Gross, J. E., Azuaje, E. I., Thorpe, A. S., Gudex-Cross, D., et al. (2019). The terrestrial organism and biogeochemistry spatial sampling design for the national ecological observatory network. Ecosphere, 10(2):e02540.

Description

Land cover covariates (forest cover and grassland cover) extracted at a 1km resolution across the Disney Wildnerss Preserve for use in predicting density across the park. Land cover data come from USGS EROS.

Usage

data(neonPredData)

Format

neonPredData is a data frame with four columns:

forest: proportion of forest cover within 1km radius.

grass: proportion of grassland cover within 1km radius.

Easting: the x coordinate of the point. The EPSG is 32617 (UTM Zone 17N).

Northing: the y coordinate of the point. The EPSG is 32617 (UTM Zone 17N).

Source

USGS Earth Resources Observation and Science Ceneter https://www.usgs.gov/centers/eros

Description

Function for fitting single-species N-mixture models.

Usage

NMix(abund.formula, det.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, ...)

Arguments

Value

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

The return object will include additional objects used for subsequent prediction and/or model fit evaluation. Note that detection probability estimated values are not included in the model object, but can be extracted using fitted().

Author(s)

Jeffrey W. Doser [email protected],

References

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.

Roberts, G.O. and Rosenthal J.S. (2009) Examples of adaptive MCMC. Journal of Computational and Graphical Statistics, 18(2):349-367.

Royle, J. A. (2004). N‐mixture models for estimating population size from spatially replicated counts. Biometrics, 60(1), 108-115. doi:10.1111/j.0006-341X.2004.00142.x.

Examples

set.seed(1010)
J.x <- 10
J.y <- 10
J <- J.x * J.y
n.rep <- sample(3, J, replace = TRUE)
beta <- c(0, -1.5)
p.abund <- length(beta)
alpha <- c(0.5, 1.2, -0.5)
p.det <- length(alpha)
mu.RE <- list()
p.RE <- list()
phi <- 3/.6
sigma.sq <- 2
kappa <- 0.3
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, 
               phi = phi, sigma.sq = sigma.sq, cov.model = cov.model, 
               family = 'NB')

y <- dat$y
X <- dat$X
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])

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)), 
                   kappa.unif = c(0, 10)) 
# Starting values
inits.list <- list(alpha = 0,
                   beta = 0,
                   kappa = kappa, 
                   N = apply(y, 1, max, na.rm = TRUE))

n.batch <- 4
batch.length <- 25
n.burn <- 50
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)

Description

Function for performing posterior predictive checks on spAbundance model objects.

Usage

ppcAbund(object, fit.stat, group, type = 'marginal', ...)

Arguments

Value

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

The return object will include additional objects used for standard extractor functions.

Note

ppcAbund will return an error for Gaussian or zero-inflated Gaussian models. For Gaussian models, standard residual diagnostics can be used to assess model fit.

Author(s)

Jeffrey W. Doser [email protected],

Examples

set.seed(1010)
J.x <- 10
J.y <- 10
J <- J.x * J.y
n.rep <- sample(3, J, replace = TRUE)
beta <- c(0, -1.5)
p.abund <- length(beta)
alpha <- c(0.5, 1.2, -0.5)
p.det <- length(alpha)
mu.RE <- list()
p.RE <- list()
phi <- 3/.6
sigma.sq <- 2
kappa <- 0.3
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, 
               phi = phi, sigma.sq = sigma.sq, cov.model = cov.model, 
               family = 'NB')

y <- dat$y
X <- dat$X
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])

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)), 
                   kappa.unif = c(0, 10)) 
# Starting values
inits.list <- list(alpha = 0, beta = 0, kappa = kappa, 
                   N = apply(y, 1, max, na.rm = TRUE))


tuning <- 0.5
n.batch <- 4
batch.length <- 25
n.burn <- 50
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) 

# Posterior predictive check
ppc.out <- ppcAbund(out, fit.stat = 'chi-squared', group = 0)
summary(ppc.out)

Description

The function predict collects posterior predictive samples for a set of new locations given an object of class 'abund'.

Usage

## S3 method for class 'abund'
predict(object, X.0, ignore.RE = FALSE, z.0.samples, ...)

Arguments

Value

A list object of class predict.abund. The list consists of:

The return object will include additional objects used for standard extractor functions.

Note

When ignore.RE = FALSE, both sampled levels and non-sampled levels of random effects are supported for prediction. For sampled levels, the posterior distribution for the random effect corresponding to that level of the random effect will be used in the prediction. For non-sampled levels, random values are drawn from a normal distribution using the posterior samples of the random effect variance, which results in fully propagated uncertainty in predictions with models that incorporate random effects.

Author(s)

Jeffrey W. Doser [email protected],

Examples

set.seed(1010)
J.x <- 15
J.y <- 15
J <- J.x * J.y
n.rep <- sample(1, J, replace = TRUE)
beta <- c(0, -1.5, 0.3, -0.8)
p.abund <- length(beta)
mu.RE <- list()
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')

# Split into fitting and prediction data set
pred.indx <- sample(1:J, round(J * .25), replace = FALSE)
y <- dat$y[-pred.indx, ]
# Abundance covariates
X <- dat$X[-pred.indx, , , drop = FALSE]
# Prediction covariates
X.0 <- dat$X[pred.indx, , ]
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])

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)

n.batch <- 5
batch.length <- 25
n.burn <- 0
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) 

# Predict at new locations ------------------------------------------------
colnames(X.0) <- c('intercept', 'abund.cov', 'abund.cov.2', 'abund.cov.3')
out.pred <- predict(out, X.0)
mu.0.quants <- apply(out.pred$mu.0.samples, 2, quantile, c(0.025, 0.5, 0.975))
plot(dat$mu[pred.indx], mu.0.quants[2, ], pch = 19, xlab = 'True', 
     ylab = 'Fitted', ylim = c(min(mu.0.quants), max(mu.0.quants)))
segments(dat$mu[pred.indx], mu.0.quants[1, ], dat$mu[pred.indx], mu.0.quants[3, ])
lines(dat$mu[pred.indx], dat$mu[pred.indx])