Willow Tit
Golas:
- perform a logistic regression in JAGS
- incorporate detection probability in occupancy models
We will use data from presence/absence of “Willow Tit” (Poecile montanus) in sampling points from Switzerland (Royle & Dorazio (2008). Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations and communities. Academic Press).
Here’s the data:
y.1 <- c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1,
0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1,
0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0,
1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0,
1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
y.2 <- c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1,
0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1,
0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, NA, 1, 1, 0, NA, 0, 0, 0, 1, 0, 0,
0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,
1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
y.3 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,
0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, NA, 1, 1, NA, 0, 1, 0, 0, 0, 1, 0, 1,
0, 1, 0, 1, NA, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, NA,
NA, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, NA, 1, 1, 0, NA, 0, 0, 0, 1,
0, 0, 1, 0, 1, 0, 0, 0, NA, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 1, NA, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, NA, 1, 0, NA, 1, 0, 1, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
NA, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, NA, 0, 0, 0, 0, 1, 0, 0,
0, 0, 1, 1, 0, 0, 0, 0, NA, 0, 0, 0, 0, NA, NA, NA, NA, 1, NA, NA, NA, 1,
1, 1, 1, 1, NA, 1, 1, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA)
elev <- c(420, 450, 1050, 1110, 510, 630, 590, 530, 1140, 770, 1220, 460, 1010,
760, 1300, 1270, 380, 550, 390, 1380, 530, 1190, 1490, 920, 620, 540, 820,
1220, 1180, 730, 580, 490, 1270, 930, 510, 1830, 1940, 1090, 2010, 790,
1730, 1120, 620, 590, 1840, 930, 540, 420, 1400, 890, 1460, 1850, 1420,
600, 550, 750, 1790, 680, 990, 520, 760, 370, 2710, 650, 660, 730, 460,
360, 250, 1940, 2420, 1320, 790, 470, 410, 1940, 1570, 1880, 1290, 1340,
880, 710, 630, 450, 2360, 1650, 2050, 790, 730, 490, 840, 450, 1010, 1450,
850, 1240, 570, 1500, 1420, 780, 470, 2080, 2020, 1750, 1630, 640, 530,
660, 1360, 1290, 670, 480, 2000, 800, 1010, 680, 490, 400, 320, 1400, 590,
430, 410, 510, 1930, 1560, 700, 440, 480, 490, 1850, 490, 480, 440, 820,
670, 610, 1700, 900, 1250, 910, 540, 450, 470, 1880, 1400, 470, 830, 1950,
1660, 1210, 380, 1840, 950, 880, 480, 550, 550, 340, 550, 1320, 260, 2010,
1130, 1910, 1630, 1540, 1340, 1180, 830, 540, 1410, 1320, 890, 580, 500,
1220, 710, 1820, 520, 1890, 1650, 1390, 1390, 1290, 1810, 1130, 480, 1040,
410, 2030, 1880, 950, 450, 950, 1070, 1980, 1350, 540, 420, 1420, 2240,
1980, 2330, 2050, 1180, 2010, 2220, 2310, 2190, 1770, 1830, 1870, 1840,
2420, 1820, 1850, 2250, 2250, 2450, 2350, 2650, 2750, 2450, 1950, 2150,
1850, 2050, 1850, 2050, 2150, 2050, 1950, 1950, 2250, 2250, 2350)
forest <- c(3, 21, 32, 35, 2, 60, 5, 13, 50, 57, 84, 17, 58, 26, 32, 66, 45,
31, 8, 78, 37, 18, 54, 6, 3, 27, 56, 66, 43, 40, 5, 3, 52, 38, 49, 15, 62,
82, 0, 5, 58, 98, 33, 53, 17, 16, 60, 2, 52, 58, 50, 12, 29, 16, 83, 79,
51, 8, 45, 5, 64, 14, 0, 0, 30, 17, 65, 30, 0, 5, 0, 66, 8, 75, 68, 72,
95, 57, 64, 44, 44, 34, 21, 23, 1, 43, 30, 26, 7, 87, 49, 39, 93, 48, 11,
73, 12, 42, 26, 13, 2, 15, 54, 74, 45, 61, 33, 54, 69, 47, 5, 36, 60, 72,
29, 7, 7, 13, 2, 16, 62, 9, 37, 71, 41, 19, 24, 10, 0, 70, 45, 18, 20, 3,
20, 36, 19, 35, 87, 51, 21, 11, 29, 55, 32, 32, 3, 58, 36, 86, 75, 23, 0,
71, 69, 6, 67, 67, 42, 72, 79, 47, 26, 56, 18, 33, 21, 39, 69, 42, 98, 52,
90, 50, 54, 9, 51, 35, 57, 18, 12, 34, 68, 63, 48, 33, 90, 33, 37, 8, 36,
66, 56, 0, 41, 63, 0, 35, 14, 1, 27, 0, 6, 0, 1, 58, 9, 14, 2, 30, 93, 82,
81, 82, 0, 69, 53, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 15, 0, 0, 2, 0, 0,
0, 0)There are 237 sites visited (up to tree times). The vector y.1 holds data for the first visit (NA if the site was not visited that time, 1 if the species was detected, and 0 if not).
y.2 is for the second visit and y.3 for the third. Vector elev holds the elevation value (meters above sea level) for every site and forest holds the percent forest cover for every site.
The goal is to relate the presence of this species to forest cover and elevation. When we have predictors with continous values and with different scales or units it is useful to center and standardize them before the analysis. Below we do that and we also create a new variable with the square of elevation so that we can test for possible optimal elevation for this species. Whe then count how many times each site was visited and we combine the data from the different visits into a single variable:
alti <- as.vector(scale(elev, center = TRUE))
bosq <- as.vector(scale(forest, center = TRUE))
alti2 <- alti * alti
ym <- as.matrix(cbind(y.1, y.2, y.3)) # combine all visits into a matrix 'ym'
J <- rowSums(!is.na(ym)) # Count the number of visits to each site
M <- nrow(ym) # number of visited sites
y <- rowSums(ym, na.rm = TRUE) # number of times that the species was register at each site.Now we can fit a logistic regression in JAGS. First, we have to write the model in bugslanguage:
cat(file = "WillowBin.bug",
"
model {
for(i in 1:M){
logit(z[i]) <- b0 + b1*alti[i] + b2*alti2[i] + b3*bosq[i]
y[i] ~ dbin(z[i],J[i])
}
b0 ~ dnorm(0,.001)
b1 ~ dnorm(0,.001)
b2 ~ dnorm(0,.001)
b3 ~ dnorm(0,.001)
logit(psi0) <- b0 # aquí guardamos la pr de ocurrencia para condiciones promedio
}
")Now we have to make a list of the data that we need to pass to JAGS, a function to generate appropriate random starting points for the Markov Chains.
data <- list("y", "J", "M", "bosq", "alti", "alti2")
inits <- function() {
list(b0 = rnorm(1), b1 = rnorm(1), b2 = rnorm(1), b3 = rnorm(1))
}
parameters <- c("b0", "b1", "b2", "b3", "psi0")
ni <- 5000 # número de iteraciones
nt <- 4 # descartamos 'thin' algunos valores
nb <- 1000 # cuantas iteraciones usamos de 'burn in'
nc <- 3 # y cuantas cadenas corremosNow we can call JAGS…
library(jagsUI)
wb.sim <- jags(data, inits, parameters, "willowBin.bug", n.chains = nc, n.thin = nt,
n.iter = ni, n.burnin = nb)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 1993
##
## Initializing model
##
## Adaptive phase, 100 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 4000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.As usual, we have to check that all the chains have converged, and that we have reasonable effective sample sizes for all posteriors. To see a summary of the results we write:
print(wb.sim)## JAGS output for model 'willowBin.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## burn-in = 1000 iterations and thin rate = 4,
## yielding 3000 total samples from the joint posterior.
## MCMC ran for 0.146 minutes at time 2016-02-18 22:12:36.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## b0 -0.874 0.179 -1.239 -0.869 -0.548 FALSE 1 1.030 73
## b1 2.128 0.185 1.770 2.123 2.520 FALSE 1 1.084 35
## b2 -0.968 0.182 -1.320 -0.972 -0.625 FALSE 1 1.018 134
## b3 0.864 0.127 0.611 0.864 1.113 FALSE 1 1.008 270
## psi0 0.296 0.037 0.225 0.295 0.366 FALSE 1 1.029 74
## deviance 456.296 2.738 452.963 455.640 463.072 FALSE 1 1.013 190
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 3.7 and DIC = 460.007
## DIC is an estimate of expected predictive error (lower is better).We see that both elevation and cover have coefficients larger than zero while elevation squared has a negative coefficient (implying that there is a maximum probability of occurrence at some intermediate elevation). We can also see that the occurrence probability for average environmental conditions is 0.297. Lets plot the changes in the probability of occurrence based on the posterior estimates
op <- par(mfrow = c(1, 2), mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.3,
font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(elev, plogis(wb.sim$mean$b0 + wb.sim$mean$b1 * alti + wb.sim$mean$b2 *
alti2), xlab = "Elevation", ylab = "Probability of Occurrence")
plot(forest, plogis(wb.sim$mean$b0 + wb.sim$mean$b3 * bosq), xlab = "Forest Cover",
ylab = "Probability of Occurrence")
par(op)Imperfect Detection
If we consider that it is possible that we do not detect the presence of the species in a site even if it is actually there, we could model separately the “true” presence/absence from the “observed” one. To do this, we define a hidden variable z as a stochastic variable following a Bernoulli distribution (a Binomial with N=1) where the probability of success (presence) depends on covariates as in the previous case. Then we connect this hidden variable z with the observations through a detection probability p. In the bugs model below the trick is to multiply z times p and to use this as the probability of success for the Binomial that models the data. When z is \(0\), the probability of success is \(0\) and we have a “structural zero”, but when z is equal to 1, then the probability of success in the Binomial is given by the detection probability p and we could have false negatives.
cat(file = "willow.bug",
"
model {
for(i in 1:M){
z[i] ~ dbin(psi[i],1)
logit(psi[i]) <- b0 + b1 * alti[i] + b2 * alti2[i] + b3 * bosq[i]
tmp[i] <- z[i] * p
y[i] ~ dbin(tmp[i], J[i])
}
p ~ dunif(0,1)
b0 ~ dnorm(0,.001)
b1 ~ dnorm(0,.001)
b2 ~ dnorm(0,.001)
b3 ~ dnorm(0,.001)
logit(psi0)<-b0
}
")We use the same data that we used before but in the inits function we include initial values for z and p.
data <- list("y", "J", "M", "bosq", "alti", "alti2")
inits <- function() {
list(z = rbinom(M, 1, 1), b0 = 4, p = runif(1), b1 = rnorm(1), b2 = rnorm(1),
b3 = rnorm(1))
}
parameters <- c("p", "b0", "b1", "b2", "b3", "psi0")
ni <- 5000 # número de iteraciones
nt <- 4 # descartamos 'thin' algunos valores
nb <- 1000 # cuantas iteraciones usamos de 'burn in'
nc <- 3 # y cuantas cadenas corremos
w.sim <- jags(data, inits, parameters, "willow.bug", n.chains = nc, n.thin = nt,
n.iter = ni, n.burnin = nb)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 2469
##
## Initializing model
##
## Adaptive phase, 100 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 4000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.print(w.sim)## JAGS output for model 'willow.bug', generated by jagsUI.
## Estimates based on 3 chains of 5000 iterations,
## burn-in = 1000 iterations and thin rate = 4,
## yielding 3000 total samples from the joint posterior.
## MCMC ran for 0.251 minutes at time 2016-02-18 22:12:45.
##
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## p 0.787 0.028 0.732 0.788 0.841 FALSE 1.000 1.000 3000
## b0 -0.215 0.286 -0.790 -0.212 0.329 TRUE 0.774 1.001 1583
## b1 2.088 0.326 1.487 2.078 2.771 FALSE 1.000 1.001 1830
## b2 -1.149 0.282 -1.723 -1.134 -0.618 FALSE 1.000 1.000 3000
## b3 0.890 0.240 0.443 0.879 1.390 FALSE 1.000 1.002 963
## psi0 0.448 0.069 0.312 0.447 0.582 FALSE 1.000 1.001 1545
## deviance 173.146 10.128 160.633 171.028 197.356 FALSE 1.000 1.001 1544
##
## Successful convergence based on Rhat values (all < 1.1).
## Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## For each parameter, n.eff is a crude measure of effective sample size.
##
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
##
## DIC info: (pD = var(deviance)/2)
## pD = 51.3 and DIC = 224.403
## DIC is an estimate of expected predictive error (lower is better).We can see that the detection probability is relatively high and that the coefficients for the environmental variables have similar values to those of the previos model. Also, the probability of occurrence for the average environmental conditions increased. We can plot the posterior for this probability writing:
plot(density(w.sim$sims.list$psi0), main = "")
Lets compare the different model predictions (model with imperfect detection in red).
op <- par(mfrow = c(1, 2), mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.3,
font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
plot(elev, plogis(wb.sim$mean$b0 + wb.sim$mean$b1 * alti + wb.sim$mean$b2 *
alti2), xlab = "Elevation", ylab = "Probability of Occurrence", ylim = c(0,
1))
points(elev, plogis(w.sim$mean$b0 + w.sim$mean$b1 * alti + w.sim$mean$b2 * alti2),
xlab = "Elevation", ylab = "Probability of Occurrence", col = 2)
plot(forest, plogis(wb.sim$mean$b0 + wb.sim$mean$b3 * bosq), xlab = "Forest Cover",
ylab = "Probability of Occurrence", ylim = c(0, 1))
points(forest, plogis(w.sim$mean$b0 + w.sim$mean$b3 * bosq), xlab = "Forest Cover",
ylab = "Probability of Occurrence", col = 2)
par(op)Juan Manuel Morales. 2016-02-18