# Functions for the book Applied Hierarchical Modeling in Ecology (AHM)
# Marc Kery & Andy Royle, Academic Press, 2016.
# data.fn - AHM1 section 4.3 p135
# Function to simulate data for binomial mixture model
# (generally: spatially and temporally replicated counts)
# (introduced in AHM1 Section 4.3)
data.fn <- function(M = 267, J = 3, mean.lambda = 2, beta1 = -2, beta2 = 2, beta3 = 1, mean.detection = 0.3, alpha1 = 1, alpha2 = -3, alpha3 = 0, show.plot = TRUE){
#
# Function to simulate point counts replicated at M sites during J occasions.
# Population closure is assumed for each site.
# Expected abundance may be affected by elevation (elev),
# forest cover (forest) and their interaction.
# Expected detection probability may be affected by elevation,
# wind speed (wind) and their interaction.
# Function arguments:
# M: Number of spatial replicates (sites)
# J: Number of temporal replicates (occasions)
# mean.lambda: Mean abundance at value 0 of abundance covariates
# beta1: Main effect of elevation on abundance
# beta2: Main effect of forest cover on abundance
# beta3: Interaction effect on abundance of elevation and forest cover
# mean.detection: Mean detection prob. at value 0 of detection covariates
# alpha1: Main effect of elevation on detection probability
# alpha2: Main effect of wind speed on detection probability
# alpha3: Interaction effect on detection of elevation and wind speed
# show.plot: if TRUE, plots of the data will be displayed;
# set to FALSE if you are running simulations.
if(FALSE) x <- NULL # deals with R CMD check issues with 'curve'
logit <- plogis # so 'logit' is displayed on y axis
# Create covariates
elev <- runif(n = M, -1, 1) # Scaled elevation
forest <- runif(n = M, -1, 1) # Scaled forest cover
wind <- array(runif(n = M*J, -1, 1), dim = c(M, J)) # Scaled wind speed
# Model for abundance
beta0 <- log(mean.lambda) # Mean abundance on link scale
lambda <- exp(beta0 + beta1*elev + beta2*forest + beta3*elev*forest)
N <- rpois(n = M, lambda = lambda) # Realised abundance
Ntotal <- sum(N) # Total abundance (all sites)
psi.true <- mean(N>0) # True occupancy in sample
# Model for observations
alpha0 <- qlogis(mean.detection) # mean detection on link scale
p <- plogis(alpha0 + alpha1*elev + alpha2*wind + alpha3*elev*wind)
C <- matrix(NA, nrow = M, ncol = J) # Prepare matrix for counts
for (i in 1:J){ # Generate counts by survey
C[,i] <- rbinom(n = M, size = N, prob = p[,i])
}
summaxC <- sum(apply(C,1,max)) # Sum of max counts (all sites)
psi.obs <- mean(apply(C,1,max)>0) # Observed occupancy in sample
# All the plots
if(show.plot){
oldpar <- par(mfrow = c(2, 2), cex.main = 1)
oldAsk <- devAskNewPage(ask = dev.interactive(orNone = TRUE))
on.exit({par(oldpar) ; devAskNewPage(oldAsk)})
tryPlot <- try( {
# Page 1: lambda relationships
curve(exp(beta0 + beta1*x), -1, 1, col = "red", main = "Relationship lambda-elevation \nat average forest cover", frame.plot = FALSE, xlab = "Scaled elevation")
plot(elev, lambda, xlab = "Scaled elevation", main = "Relationship lambda-elevation \nat observed forest cover", frame.plot = FALSE)
curve(exp(beta0 + beta2*x), -1, 1, col = "red", main = "Relationship lambda-forest \ncover at average elevation", xlab = "Scaled forest cover", frame.plot = FALSE)
plot(forest, lambda, xlab = "Scaled forest cover", main = "Relationship lambda-forest cover \nat observed elevation", frame.plot = FALSE)
# Page 2: p relationships
par(mfrow = c(2, 2))
curve(logit(alpha0 + alpha1*x), -1, 1, col = "red", main = "Relationship p-elevation \nat average wind speed", xlab = "Scaled elevation", frame.plot = FALSE)
matplot(elev, p, xlab = "Scaled elevation", main = "Relationship p-elevation\n at observed wind speed", pch = "*", frame.plot = FALSE)
curve(logit(alpha0 + alpha2*x), -1, 1, col = "red", main = "Relationship p-wind speed \n at average elevation", xlab = "Scaled wind speed", frame.plot = FALSE)
matplot(wind, p, xlab = "Scaled wind speed", main = "Relationship p-wind speed \nat observed elevation", pch = "*", frame.plot = FALSE)
# Page 3: counts
matplot(elev, C, xlab = "Scaled elevation", main = "Relationship counts and elevation", pch = "*", frame.plot = FALSE)
matplot(forest, C, xlab = "Scaled forest cover", main = "Relationship counts and forest cover", pch = "*", frame.plot = FALSE)
matplot(wind, C, xlab = "Scaled wind speed", main = "Relationship counts and wind speed", pch = "*", frame.plot = FALSE)
desc <- paste('Counts at', M, 'sites during', J, 'surveys')
# hist(C, main = desc, breaks = 50, col = "grey")
histCount(C, NULL, color = "grey", main = desc)
}, silent = TRUE )
if(inherits(tryPlot, "try-error"))
tryPlotError(tryPlot)
}
# Output
return(list(M = M, J = J, mean.lambda = mean.lambda, beta0 = beta0, beta1 = beta1, beta2 = beta2, beta3 = beta3, mean.detection = mean.detection, alpha0 = alpha0, alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, elev = elev, forest = forest, wind = wind, lambda = lambda, N = N, p = p, C = C, Ntotal = Ntotal, psi.true = psi.true, summaxC = summaxC, psi.obs = psi.obs))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.