simCountData {MigConnectivity}

Simulates Breeding Bird Survey-style count data

Description

Recently updated (version 0.4.3) to more properly match current BBS models. modelCountDataJAGS has not been updated yet

Usage

simCountData(
  nStrata,
  routesPerStrata,
  nYears,
  alphaStrat,
  beta = 0,
  eta = 0,
  sdRoute = 0,
  sdYear = 0,
  sdObs = 0,
  sdCount = 0,
  model = c("S", "Sh", "D", "Dh"),
  obsSurvival = 1,
  fixedyear = round(nYears/2),
  nuCount = 1
)

Arguments

nStrata

Number of populations/regions/strata

routesPerStrata

Vector of length 1 or nStrata containing the number of routes (i.e. counts) per stratum. If length(routesPerStrata) == 1, number of routes is identical for each population

nYears

Number of years surveys were conducted

alphaStrat

Vector of length 1 or nStrata containing the log expected number of individuals counted at each route for each population. If length(alphaStrat) == 1, expected counts are identical for each population

beta

Coefficient of linear year effect (default 0)

eta

Coefficient of first time run effect (default 0)

sdRoute

Standard deviation of random route-level variation. Default is 0, and if you're setting sdObs, it's probably best to keep it that way

sdYear

Standard deviation of random year-level variation (default 0)

sdObs

Standard deviation of random observer variation (default 0)

sdCount

Standard deviation of random count-level variation (default 0)

model

One of "S" (default), "Sh", "D", and "Dh". See Link et al. (2020) for descriptions of these models

obsSurvival

Annual probability that the observer that ran a route one year will run it the next. Default 1 (each route has only one observer)

fixedyear

The year within nYears that alphaStrat applies directly to (default halfway through)

nuCount

For the "h" models, parameter for extra variation in counts (default 1)

Value

simCountData returns a list containing:

nStrata

Number of populations/regions.

nRoutes

Total number of routes.

nYears

Number of years.

routesPerStrata

Number of routes per population.

year

Vector of length nYears with standardized year values.

strat

Vector of length nRoutes indicating the population/region in which each route is located.

alphaStrat

log expected count for each populations.

epsRoute

realized deviation from alphaStrat for each route.

epsYear

realized deviation from alphaStrat for each year.

beta

linear year effect.

sdRoute

standard deviation of random route-level variation.

sdYear

standard deviation of random year-level variation.

expectedCount

nRoutes by nYears matrix containing deterministic expected counts.

C

nRoutes by nYears matrix containing observed counts.

References

Cohen, E. B., J. A. Hostetler, M. T. Hallworth, C. S. Rushing, T. S. Sillett, and P. P. Marra. 2018. Quantifying the strength of migratory connectivity. Methods in Ecology and Evolution 9: 513-524. doi:10.1111/2041-210X.12916

Link, W. A., J. R. Sauer, and D. K. Niven. 2020. Model selection for the North American Breeding Bird Survey. Ecological Applications 30: e02137. doi:10.1002/eap.2137

Examples

set.seed(150)

### Set parameters for simulation ----

# Number of populations
nStrata. <- 4
# Number of routes w/i each population (assumed to be balanced)
routePerStrat. <- 30 # reduced from 90 for example speed
# Number of years
nYears. <- 5 # reduced from 10 for example speed
# log(Expected number of birds counted at each route)
alphaStrat. <- 1.95
# standard deviation of normal distribution assumed for route/observer random
# effects
sdRoute. <- 0.6
# standard deviation of normal distribution assumed for year random effects
sdYear. <- 0.18


# Number of simulated datasets to create and model
nsims <- 50 # reduced from 100 for example speed
# Number of MCMC iterations
ni. <- 1000 # reduced from 15000 for example speed
# Number of iterations to thin from posterior (reduced from 5)
nt. <- 1
# Number of iterations to discard as burn-in
nb. <- 500 # reduced from 5000 for example speed
# Number of MCMC chains
nc. <- 1 # reduced from 3 for example speed

### Create empty matrix to store model output ---
sim_in <- vector("list", nsims)
sim_out <- vector("list", nsims)


# Simulation ---

system.time(for(s in 1:nsims){
  cat("Simulation",s,"of",nsims,"\n")

  # Simulate data
  sim_data <- simCountData(nStrata = nStrata., routesPerStrata = routePerStrat.,
                           nYears = nYears., alphaStrat = alphaStrat.,
                           sdRoute = sdRoute., sdYear = sdYear.)
  sim_in[[s]] <- sim_data


  # Estimate population-level abundance
  out_mcmc <- modelCountDataJAGS(count_data = sim_data, ni = ni., nt = nt.,
                                 nb = nb., nc = nc.)

  # Store model output
  sim_out[[s]] <- out_mcmc
  remove(out_mcmc)

})


### Check that relative abundance is, on average, equal for each population
prop.table(sapply(sim_in, function(x) return(rowsum(colSums(x$C), x$strat))), 2)

rel_names <- paste0('relN[', 1:nStrata., ']')
rel_abund1 <- data.frame(sim=1:nsims,
                         ra1.mean=NA, ra2.mean=NA, ra3.mean=NA, ra4.mean=NA,
                         ra1.low=NA, ra2.low=NA, ra3.low=NA, ra4.low=NA,
                         ra1.high=NA, ra2.high=NA, ra3.high=NA, ra4.high=NA,
                         ra1.cover=0, ra2.cover=0, ra3.cover=0, ra4.cover=0)
for (s in 1:nsims) {
  rel_abund1[s, 2:5] <- summary(sim_out[[s]])$statistics[rel_names, "Mean"]
  rel_abund1[s, 6:9] <- summary(sim_out[[s]])$quantiles[rel_names, 1]
  rel_abund1[s, 10:13] <- summary(sim_out[[s]])$quantiles[rel_names, 5]
}
rel_abund1 <- transform(rel_abund1,
                        ra1.cover = (ra1.low<=0.25 & ra1.high>=0.25),
                        ra2.cover = (ra2.low<=0.25 & ra2.high>=0.25),
                        ra3.cover = (ra3.low<=0.25 & ra3.high>=0.25),
                        ra4.cover = (ra4.low<=0.25 & ra4.high>=0.25))

summary(rel_abund1)