############################################################################# # # Bayesian population analysis using WinBUGS - a hierarchical perspective # Marc Kéry & Michael Schaub # ############################################################################# # # Code for the solution of the exercises # # last up-date: 6.2.2012 # ############################################################################# ############################################################################# # # Chapter 1 # ############################################################################# ############## # Exercise 2 ############## p <- 0.4 N <- 1:100 counts <- array(NA, dim = c(100, 100)) for (i in 1:100){ counts[i,] <- rbinom(n = 100, size = N[i], prob = p) } obs.variance <- apply(counts, 1, var) theo.variance <- N*p*(1-p) plot(N, obs.variance, col = "blue", type = "l", xlab = "Population size (N)", ylab = "Variance of counts", las = 1, lwd = 3, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.2, main = "Detection probability p = 0.4") lines(N, theo.variance, col = "red", type = "l", lwd = 3) legend(5, 25, c('Variance (theoretical)', 'Variance (observed)'), col=c("red", "blue"), lty = c(1,1), lwd = 2, cex = 1.2) p <- seq(0,1,0.01) N <- 16 counts <- array(NA, dim = c(101, 100)) for (i in 1:101){ counts[i,] <- rbinom(n = 100, size = N, prob = p[i]) } obs.variance <- apply(counts, 1, var) theo.variance <- N*p*(1-p) plot(p, obs.variance, col = "blue", type = "l", xlab = "Detection probability (p)", ylab = "Variance of counts", las = 1, lwd = 3, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.2, main = "Population size N = 16") lines(p, theo.variance, col = "red", type = "l", lwd = 3) legend(0.25, 1, c('Variance (theoretical)', 'Variance (observed)'), col=c("red", "blue"), lty = c(1,1), lwd = 2, cex = 1.2) ############################################################################# # # Chapter 2 # ############################################################################# ############## # No exercises ############## ############################################################################# # # Chapter 3 # ############################################################################# ############## # Exercise 1 ############## data.fn <- function(n = 40, alpha = 4.32671, beta1 = 1.22677, beta2 = 0.05552, beta3 = -0.23758){ # n: Number of years # alpha, beta1, beta2, beta3: coefficients of a # cubic polynomial of count on year (standardised to have mean 0 # and sd 1) # Generate values of time covariate and scale it year <- 1:n year <- as.numeric(scale(year)) # Signal: Build up systematic part of the GLM log.expected.count <- alpha + beta1 * year + beta2 * year^2 + beta3 * year^3 expected.count <- exp(log.expected.count) # Noise: generate random part of the GLM: Poisson noise around expected counts C <- rpois(n = n, lambda = expected.count) # Plot simulated data plot(year, C, type = "b", lwd = 2, col = "black", main = "", las = 1, ylab = "Population size", xlab = "Year", cex.lab = 1.2, cex.axis = 1.2) lines(year, expected.count, type = "l", lwd = 3, col = "red") return(list(n = n, alpha = alpha, beta1 = beta1, beta2 = beta2, beta3 = beta3, year = year, expected.count = expected.count, C = C)) } data <- data.fn() fm <- glm(C ~ year + I(year^2) + I(year^3), family = poisson, data = data) summary(fm) bugs.dir <- "c:/Program files/WinBUGS14/" # Specify model in BUGS language sink("GLM_Poisson.txt") cat(" model { # Priors alpha ~ dunif(-20, 20) beta1 ~ dunif(-10, 10) beta2 ~ dunif(-10, 10) beta3 ~ dunif(-10, 10) # Likelihood: Note key components of a GLM on one line each for (i in 1:n){ C[i] ~ dpois(lambda[i]) # 1. Distribution for random part log(lambda[i]) <- log.lambda[i] # 2. Link function log.lambda[i] <- alpha + beta1 * year[i] + beta2 * pow(year[i],2) + beta3 * pow(year[i],3) # 3. Linear predictor } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = data$C, n = length(data$C), year = data$year) # Initial values inits <- function() list(alpha = runif(1, -2, 2), beta1 = runif(1, -3, 3)) # Parameters monitored params <- c("alpha", "beta1", "beta2", "beta3", "lambda") # MCMC settings ni <- 2000 nt <- 2 nb <- 1000 nc <- 3 # Call WinBUGS from R out <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = "GLM_Poisson.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out$summary[1:4,], dig = 3) print(fm$coef, dig = 4) plot(1:40, data$C, type = "b", lwd = 2, col = "black", main = "", las = 1, ylab = "Population size", xlab = "Year") R.predictions <- predict(glm(C ~ year + I(year^2) + I(year^3), family = poisson, data = data), type = "response") lines(1:40, R.predictions, type = "l", lwd = 3, col = "green") WinBUGS.predictions <- out$mean$lambda lines(1:40, WinBUGS.predictions, type = "l", lwd = 3, col = "blue", lty = 2) ############## # Exercise 2 ############## n <- c(22, 8, 10, 7, 10, 6, 11) r <- c(20, 7, 10, 6, 0, 1, 2) X <- c(0, 3, 1, 4, 5, 8, 10) plot(X, r/n, type = "p", xlab = "Covariate X", ylab = "Ratio r/n", las = 1, lwd = 3, cex.main = 1.5, cex.lab = 1.5, cex.axis = 1.2, main = "") # Specify model in BUGS language sink("logistic_regression.txt") cat(" model { # Priors alpha ~ dnorm(0, 0.001) beta ~ dnorm(0, 0.001) # Likelihood for (i in 1:sample.size){ r[i] ~ dbin(p[i], n[i]) logit(p[i]) <- alpha + beta * X[i] } } ",fill = TRUE) sink() # Bundle data win.data <- list(r = r, n = n, sample.size = length(X), X = X) # Initial values inits <- function() list(alpha = runif(1, -1, 1), beta = runif(1, -1, 1)) # Parameters monitored params <- c("alpha", "beta", "p") # MCMC settings ni <- 2500 ; nt <- 2 ; nb <- 500 ; nc <- 3 # Call WinBUGS from R (BRT < 1 min) out <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = "logistic_regression.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Plot predictions (don’t get confused with order!) predictions <- out$mean$p lines(sort(X), predictions[order(X)], lwd = 3, col = "blue", lty = 2) ############## # Exercise 3 ############## data.generation <- function(){ n <- 200 # Number of sites X <- sort(runif(n, -1, 1)) # Random uniform on range -1 to 1 a <- -2 # Intercept ... b <- 5 # ... and slope of logit-linear relationship p <- plogis(a + b * X) # Success probability y <- rbinom(n = n, size = 1, prob = p) # Observed binary data plot(X, p, xlim = c(-1, 1), ylim = c(0,1), type = "l", lwd = 3, col = "red") points(X, y) return(list(X = X, y = y)) } data <- data.generation() attach(data) # Define the model in the BUGS language and write a text file sink("model.txt") cat(" model { # Priors alpha ~ dunif(-20, 20) beta ~ dunif(-20, 20) # Likelihood: Note key components of a GLM in one line each for (i in 1:nobs){ y[i] ~ dbern(p[i]) # 1. Distribution for random part logit(p[i]) <- lp[i] # 2. Link function lp[i] <- alpha + beta * X[i] # 3. Linear predictor } # Form predictions using pred.x for (i in 1:nobs.pred){ logit(pred.p[i]) <- alpha + beta * pred.x[i] } } # end model ",fill=TRUE) sink() # Bundle data pred.x = seq(-1, 1, length.out = 100) win.data <- list(y = y, X = X, nobs = length(y), pred.x = pred.x, nobs.pred = length(pred.x)) # Initial values (not required for all) inits <- function() list(alpha = runif(1, -2, 2)) # Define parameters to be monitored params <- c("alpha", "beta", "pred.p") # MCMC settings ni <- 600 ; nt <- 2 ; nb <- 100 ; nc <- 3 # Call WinBUGS from R out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = FALSE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posterior distributions print(out, dig = 3) # Plot model predictions into the same plot points(pred.x, out$mean$pred.p, col = "blue", lwd = 3, type = "l") ############## # Exercise 4 ############## # Specify Binomial GLM in BUGS language sink("GLM_Binomial.txt") cat(" model { # Priors alpha ~ dnorm(0, 0.001) beta1 ~ dnorm(0, 0.001) beta2 ~ dnorm(0, 0.001) # Likelihood for (i in 1:nyears){ C[i] ~ dbin(p[i], N[i]) # 1. Distribution for random part logit(p[i]) <- alpha + beta1 * year[i] + beta2 * pow(year[i],2) # link function and linear predictor } } ",fill = TRUE) sink() # Bundle data win.data <- list(C = R.Pairs, N = Pairs, nyears = length(Year), year = (Year-1984)/20) # Initial values inits <- function() list(alpha = runif(1, -1, 1), beta1 = runif(1, -1, 1), beta2 = runif(1, -1, 1)) # Parameters monitored params <- c("alpha", "beta1", "beta2", "p") # MCMC settings ni <- 2500 nt <- 2 nb <- 500 nc <- 3 # Call WinBUGS from R (BRT < 1 min) out1 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = "GLM_Binomial.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors and plot estimates of proportion of successful pairs print(out1, dig = 3) plot(Year, R.Pairs/Pairs, type = "b", lwd = 2, col = "black", main = "", las = 1, ylab = "Proportion successful pairs", xlab = "Year", ylim = c(0,1)) lines(Year, out1$mean$p, type = "l", lwd = 3, col = "blue") C <- R.Pairs N <- Pairs year <- (Year-1984)/20 fm1 <- glm(cbind(C, N-C) ~ year + I(year^2), family = binomial) summary(fm1) print(out1$summary[1:3,], dig = 3) print(summary(fm1)$coef[,1:2]) # Specify Poisson GLM in BUGS language sink("GLM_Poisson.txt") cat(" model { # Priors alpha ~ dnorm(0, 0.001) beta1 ~ dnorm(0, 0.001) beta2 ~ dnorm(0, 0.001) # Likelihood for (i in 1:nyears){ C[i] ~ dpois(p[i]) # 1. Distribution for random part log(p[i]) <- alpha + beta1 * year[i] + beta2 * pow(year[i],2) # link function and linear predictor } } ",fill = TRUE) sink() # Bundle data win.data <- list(C = R.Pairs, nyears = length(Year), year = (Year-1984)/20) # Initial values inits <- function() list(alpha = runif(1, -1, 1), beta1 = runif(1, -1, 1), beta2 = runif(1, -1, 1)) # Parameters monitored params <- c("alpha", "beta1", "beta2", "p") # MCMC settings ni <- 2500 nt <- 2 nb <- 500 nc <- 3 # Call WinBUGS from R (BRT < 1 min) out2 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = "GLM_Poisson.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors and plot estimates of counts of successful pairs print(out2, dig = 3) plot(Year, R.Pairs, type = "b", lwd = 2, col = "black", main = "", las = 1, ylab = "Number of successful pairs", xlab = "Year", ylim = c(0,150)) lines(Year, out2$mean$p, type = "l", lwd = 3, col = "blue") C <- R.Pairs year <- (Year-1984)/20 fm2 <- glm(C ~ year + I(year^2), family = poisson) summary(fm2) print(out2$summary[1:3,], dig = 3) print(summary(fm2)$coef[,1:2]) # Specify Normal GLM in BUGS language sink("GLM_Normal.txt") cat(" model { # Priors alpha ~ dnorm(0, 0.00001) beta1 ~ dnorm(0, 0.00001) beta2 ~ dnorm(0, 0.00001) tau <- pow(sigma, -2) sigma ~ dunif(0, 100) # Likelihood for (i in 1:nyears){ C[i] ~ dnorm(p[i], tau) # 1. Distribution for random part p[i] <- alpha + beta1 * year[i] + beta2 * pow(year[i],2) # link function and linear predictor } } ",fill = TRUE) sink() # Bundle data win.data <- list(C = R.Pairs, nyears = length(Year), year = (Year-1984)/20) # Initial values inits <- function() list(alpha = runif(1, -1, 1), beta1 = runif(1, -1, 1), beta2 = runif(1, -1, 1)) # Parameters monitored params <- c("alpha", "beta1", "beta2", "p", "sigma") # MCMC settings ni <- 2500 nt <- 2 nb <- 500 nc <- 3 # Call WinBUGS from R (BRT < 1 min) out3 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = "GLM_Normal.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors and plot estimates of counts of successful pairs print(out3, dig = 3) plot(Year, R.Pairs, type = "b", lwd = 2, col = "black", main = "", las = 1, ylab = "Number of successful pairs", xlab = "Year", ylim = c(0,150)) lines(Year, out3$mean$p, type = "l", lwd = 3, col = "blue") C <- R.Pairs year <- (Year-1984)/20 fm3 <- glm(C ~ year + I(year^2), family = gaussian) summary(fm3) print(out3$summary[1:3,], dig = 3) print(summary(fm3)$coef[,1:2]) plot(Year, R.Pairs, type = "b", lwd = 2, col = "black", main = "", las = 1, ylab = "Number of successful pairs", xlab = "Year", ylim = c(0,150)) lines(Year, Pairs*out1$mean$p, type = "l", lwd = 3, col = "blue") lines(Year, out2$mean$p, type = "l", lwd = 3, col = "green") lines(Year, out3$mean$p, type = "l", lwd = 3, col = "brown") legend(1965, 150, c('Binomial GLM', 'Poisson GLM', 'Normal GLM'), col = c("blue", "green", "brown"), lty = 1, lwd = 2, cex = 1.2) ############################################################################# # # Chapter 4 # ############################################################################# ############## # Exercise 1 ############## data <- data.fn() # Specify model in BUGS language sink("GLM_Poisson.txt") cat(" model { # Priors alpha ~ dunif(-20, 20) beta1 ~ dunif(-10, 10) beta2 ~ dunif(-10, 10) beta3 ~ dunif(-10, 10) # Likelihood for (i in 1:n){ C[i] ~ dpois(lambda[i]) log(lambda[i]) <- log.lambda[i] log.lambda[i] <- alpha + beta1 * year[i] + beta2 * pow(year[i],2) + beta3 * pow(year[i],3) } } ",fill = TRUE) sink() # Bundle data (note standardization of covariate) win.data <- list(C = data$C, n = length(data$C), year = (data$year-20)/20) # Initial values inits <- function() list(alpha = runif(1, -2, 2), beta1 = runif(1, -3, 3)) # Parameters monitored params <- c("alpha", "beta1", "beta2", "beta3", "lambda") # MCMC settings ni <- 2000 nt <- 2 nb <- 1000 nc <- 3 # Call WinBUGS from R out1 <- bugs(data = win.data, inits = inits, parameters.to.save = params, model.file = "GLM_Poisson.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out1, dig = 3) # Specify model in BUGS language sink("GLMM_Poisson.txt") cat(" model { # Priors alpha ~ dunif(-20, 20) beta1 ~ dunif(-10, 10) beta2 ~ dunif(-10, 10) beta3 ~ dunif(-10, 10) tau <- 1 / (sd*sd) sd ~ dunif(0, 5) # Likelihood for (i in 1:n){ C[i] ~ dpois(lambda[i]) log(lambda[i]) <- log.lambda[i] log.lambda[i] <- alpha + beta1 * year[i] + beta2 * pow(year[i],2) + beta3 * pow(year[i],3) + eps[i] eps[i] ~ dnorm(0, tau) } } ",fill = TRUE) sink() # Bundle data (note standardization of covariate) win.data <- list(C = data$C, n = length(data$C), year = (data$year-20)/20) # Initial values inits <- function() list(alpha = runif(1, -2, 2), beta1 = runif(1, -3, 3), sd = runif(1, 0,1)) # Parameters monitored params <- c("alpha", "beta1", "beta2", "beta3", "lambda", "sd", "eps") # MCMC settings ni <- 30000 nt <- 10 nb <- 20000 nc <- 3 # Call WinBUGS from R (BRT <1 min) out2 <- bugs(win.data, inits, params, "GLMM_Poisson.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out2, dig = 2) GLM.estimate <- out1$summary[1:4,1:2] GLMM.estimate <- out2$summary[1:4,1:2] cbind(GLM.estimate, GLMM.estimate) mean(1 - (GLM.estimate[,2] / GLMM.estimate[,2])) plot(data$year, data$expected.count, type = "l", col = "red", lwd = 2, main = "", las = 1, ylab = "Population size", xlab = "Year", ylim = c(0, 320), cex.axis = 1.2, cex.lab = 1.2, las = 1) lines(1:40-0.2, out1$mean$lambda, type = "b", col = "blue", lwd = 2, lty = 2) segments(1:40-0.2, out1$summary[5:44,3], 1:40-0.2, out1$summary[5:44,7], col = "blue", lwd = 2) lines(1:40+0.2, out2$mean$lambda, type = "b", col = "green", lwd = 2, lty = 2) segments(1:40+0.2, out2$summary[5:44,3], 1:40+0.2, out2$summary[5:44,7], col = "green", lwd = 2) legend(0, 325, c('True expected counts', 'Poisson GLM (with 95% CRI)', 'Poisson GLMM (with 95% CRI)'), col=c("red", "blue", "green"), lty = 1, lwd = 2, cex = 1.2) ############## # Exercise 2 ############## # Bundle new data without first year win.data <- list(C = t(C[,-1]), nsite = nrow(C), nyear = ncol(C)-1, first = t(first[,-1])) # Specify model in BUGS language sink("GLMM3.txt") cat(" model { # Priors mu ~ dnorm(0, 0.01) # Overall mean beta2 ~ dnorm(0, 0.01) # First-year observer effect for (j in 1:nsite){ alpha[j] ~ dnorm(0, tau.alpha) # Random site effects } tau.alpha <- 1/ (sd.alpha * sd.alpha) sd.alpha ~ dunif(0, 5) for (i in 1:nyear){ eps[i] ~ dnorm(0, tau.eps) # Random year effects } tau.eps <- 1/ (sd.eps * sd.eps) sd.eps ~ dunif(0, 5) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- mu + beta2 * first[i,j] + alpha[j] + eps[i] } #j } #i } ",fill = TRUE) sink() # Initial values inits <- function() list(mu = runif(1, 0, 4), beta2 = runif(1, -1, 1), alpha = runif(235, -2, 2), eps = runif(8, -1, 1)) # Parameters monitored params <- c("mu", "beta2", "alpha", "eps", "sd.alpha", "sd.eps") # MCMC settings ni <- 6000 nt <- 5 nb <- 1000 nc <- 3 # Call WinBUGS from R (BRT 3 min) out5.minus1 <- bugs(win.data, inits, params, "GLMM3.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posterior of first-year observer effect in table and figure print(out5.minus1$summary[1:2,], dig = 3) hist(out5.minus1$sims.list$beta2, breaks = 100, col = "grey") abline(v = 0, col = "red", lwd = 3) ############## # Exercise 3 ############## # Specify model in BUGS language: this is the model as in the book sink("GLMM2A.txt") cat(" model { # Priors mu ~ dnorm(0, 0.01) # Grand mean for (j in 1:nsite){ alpha[j] ~ dnorm(0, tau.alpha) # Random site effects } tau.alpha <- 1/ (sd.alpha * sd.alpha) sd.alpha ~ dunif(0, 5) for (i in 1:nyear){ eps[i] ~ dnorm(0, tau.eps) # Random year effects } tau.eps <- 1/ (sd.eps * sd.eps) sd.eps ~ dunif(0, 3) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- mu + alpha[j] + eps[i] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C)) # Initial values (not required for all) inits <- function() list(mu = runif(1, 0, 4), alpha = runif(235, -2, 2), eps = runif(9, -1, 1)) # Parameters monitored params <- c("mu", "sd.alpha", "sd.eps", "alpha", "eps") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT 3 min) out4A <- bugs(win.data, inits, params, "GLMM2A.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors of hyperparameters print(out4A$summary[c(1, 246:247),], dig = 4) # Specify model in BUGS language sink("GLMM2B.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.alpha) # Random site effects with grand mean } mu ~ dnorm(0, 0.01) tau.alpha <- 1/ (sd.alpha * sd.alpha) sd.alpha ~ dunif(0, 5) for (i in 1:nyear){ eps[i] ~ dnorm(0, tau.eps) # Random year effects } tau.eps <- 1/ (sd.eps * sd.eps) sd.eps ~ dunif(0, 3) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + eps[i] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C)) # Initial values (not required for all) inits <- function() list(mu = runif(1, 0, 4), alpha = runif(235, -2, 2), eps = runif(9, -1, 1)) # Parameters monitored params <- c("mu", "sd.alpha", "sd.eps", "alpha", "eps") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT 3 min) out4B <- bugs(win.data, inits, params, "GLMM2B.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors of hyperparameters print(out4B$summary[c(1, 246:247),], dig = 4) # Produce traceplots for hyperparameters from both analyses par(mfrow = c(3,2)) matplot(out4A$sims.array[,,1], type = "l", col = c("red", "blue", "green"), main = "Parameterization A") matplot(out4B$sims.array[,,1], type = "l", col = c("red", "blue", "green"), main = "Parameterization B") matplot(out4A$sims.array[,,246], col = c("red", "blue", "green"), type = "l") matplot(out4B$sims.array[,,246], col = c("red", "blue", "green"), type = "l") matplot(out4A$sims.array[,,247], col = c("red", "blue", "green"), type = "l") matplot(out4B$sims.array[,,247], col = c("red", "blue", "green"), type = "l") ############## # Exercise 4 ############## # Specify model 1 in BUGS language sink("model1.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.site) # Random site effects } mu ~ dnorm(0, 0.01) tau.site <- pow(sd.site, -2) sd.site ~ dunif(0, 5) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C)) # Initial values inits <- function() list(mu = runif(1, 2, 3)) # Parameters monitored params <- c("mu", "sd.site", "alpha") # MCMC settings ni <- 1200 nt <- 2 nb <- 200 nc <- 3 # Call WinBUGS from R (BRT 1 min) out.model1 <- bugs(win.data, inits, params, "model1.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model1, dig = 2) # Summarize posteriors for hyperparams only print(out.model1$summary[236:237,], dig = 4) # Specify model 2 in BUGS language sink("model2.txt") cat(" model { # Priors for (i in 1:nyear){ beta[i] ~ dnorm(mu, tau.year) # Random year effects } mu ~ dnorm(0, 0.01) tau.year <- pow(sd.year, -2) sd.year ~ dunif(0, 3) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- beta[i] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C)) # Initial values inits <- function() list(mu = runif(1, 2, 3)) # Parameters monitored params <- c("mu", "sd.year", "beta") # MCMC settings ni <- 1200 nt <- 2 nb <- 200 nc <- 3 # Call WinBUGS from R (BRT 1 min) out.model2 <- bugs(win.data, inits, params, "model2.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model2, dig = 2) # Summarize posteriors for hyperparams only print(out.model2$summary[10:11,], dig = 4) # Specify model 3 in BUGS language sink("model3.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.site) # Random site effects } mu ~ dnorm(0, 0.01) tau.site <- pow(sd.site, -2) sd.site ~ dunif(0, 5) for (i in 1:nyear){ beta[i] ~ dnorm(0, tau.year) # Random year effects } tau.year <- pow(sd.year, -2) sd.year ~ dunif(0, 3) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + beta[i] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C)) # Initial values (A) inits <- function() list(mu = runif(1, 0, 4)) # Initial values (B) inits <- function() list(mu = runif(1, 0, 4), alpha = runif(235, -2, 2), beta = runif(9, -1, 1)) # Parameters monitored params <- c("mu", "sd.year", "sd.site", "alpha", "beta") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT 6 min) out.model3 <- bugs(win.data, inits, params, "model3.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model3, dig = 2) # Summarize posteriors for hyperparams only print(out.model3$summary[1:3,], dig = 4) # Specify model 4 in BUGS language sink("model4.txt") cat(" model { # Priors for (i in 1:nyear){ for (j in 1:nsite){ alpha[i,j] ~ dnorm(mu, tau) # Random site-year effects } #j } #i mu ~ dnorm(0, 0.01) tau <- pow(sd, -2) sd ~ dunif(0, 5) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[i, j] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C)) # Initial values (A) inits <- function() list(mu = runif(1, 0, 4)) # Initial values (B) inits <- function() list(mu = runif(1, 0, 4), alpha = runif(9*235, -2, 2)) # Parameters monitored params <- c("mu", "sd", "alpha") # MCMC settings (so chosen that we don’t have to save too many samples) ni <- 10000 nt <- 16 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT 9 min) out.model4 <- bugs(win.data, inits, params, "model4.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model4, dig = 2) # Summarize posteriors for hyperparams only print(out.model4$summary[1:2,], dig = 4) ############## # Exercise 5 ############## # Specify model 1R in BUGS language sink("model1R.txt") cat(" model { # Priors for (i in 1:nyear){ beta[i] ~ dnorm(mu, tau.year) # Random year effects } mu ~ dnorm(0, 0.01) tau.year <- pow(sd.year, -2) sd.year ~ dunif(0, 3) for (k in 1:nobs){ gamma[k] ~ dnorm(0, tau.obs) # Random observer effects } tau.obs <- pow(sd.obs, -2) sd.obs ~ dunif(0, 3) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- beta[i] + gamma[newobs[i,j]] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), newobs = t(newobs), nsite = nrow(C), nyear = ncol(C), nobs = length(unique(c(newobs)))) # Initial values (A) inits <- function() list(mu = runif(1, 1, 3)) # Initial values (B) inits <- function() list(mu = runif(1, 1, 3), beta = runif(9, -1, 1), gamma = runif(272, -1, 1)) # Parameters monitored params <- c("mu", "sd.year", "sd.obs", "beta", "gamma") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # MCMC test settings ni <- 120 nt <- 2 nb <- 20 nc <- 3 # Call WinBUGS from R (BRT 6 min) out.model1R <- bugs(win.data, inits, params, "model1R.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model1R, dig = 2) # Summarize posteriors for hyperparams only print(out.model1R$summary[1:3,], dig = 4) # Specify model 1R in BUGS language sink("model1F.txt") cat(" model { # Priors for (i in 1:nyear){ beta[i] ~ dnorm(0, 0.01) # Fixed year effects } for (k in 1:(nobs-1)){ # gamma[k] ~ dnorm(0, 0.1) # Fixed observer effects gamma[k] ~ dunif(-10, 10) # Fixed observer effects } gamma[nobs] <- 0 # Set effect of last observer to zero # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- beta[i] + gamma[newobs[i,j]] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), newobs = t(newobs), nsite = nrow(C), nyear = ncol(C), nobs = length(unique(c(newobs)))) # Initial values inits <- function() list(beta = runif(9, -1, 1)) # Parameters monitored params <- c("beta", "gamma") # MCMC settings ni <- 6000 nt <- 5 nb <- 1000 nc <- 3 # Call WinBUGS from R (BRT 3 min) out.model1F <- bugs(win.data, inits, params, "model1F.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model1F, dig = 2) # Summarize posteriors for year effects only print(out.model1F$summary[1:9,], dig = 3) par(mfrow = c(3, 3)) for (i in 1:271){ hist(out.model1F$sims.list$gamma[,i], xlim = c(-10, 10), col = "grey", main = paste("Observer Nr.", i), xlab = "") abline(v = c(-10, 10), col = "red", lwd = 3) browser() } # Specify model 2R in BUGS language sink("model2R.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.site) # Random site effects } mu ~ dnorm(0, 0.01) tau.site <- pow(sd.site, -2) sd.site ~ dunif(0, 5) for (i in 1:nyear){ beta[i] ~ dnorm(0, tau.year) # Random year effects } tau.year <- pow(sd.year, -2) sd.year ~ dunif(0, 3) for (k in 1:nobs){ gamma[k] ~ dnorm(0, tau.obs) # Random observer effects } tau.obs <- pow(sd.obs, -2) sd.obs ~ dunif(0, 3) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + beta[i] + gamma[newobs[i,j]] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), newobs = t(newobs), nsite = nrow(C), nyear = ncol(C), nobs = length(unique(c(newobs)))) # Initial values (B) inits <- function() list(mu = runif(1, 0, 4), alpha = runif(235, -2, 2), beta = runif(9, -1, 1)) # Parameters monitored params <- c("mu", "sd.year", "sd.site", "sd.obs", "alpha", "beta", "gamma") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT 11 min) out.model2R <- bugs(win.data, inits, params, "model2R.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model2R, dig = 2) # Summarize posteriors for hyperparams only print(out.model2R$summary[1:4,], dig = 3) # Specify model 2F in BUGS language sink("model2F.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(0, 0.01) # Fixed site effects } for (i in 2:nyear){ beta[i] ~ dnorm(0, 0.1) # Fixed year effects } beta[1] <- 1 # Effect of first year set to zero for (k in 1:(nobs-1)){ # gamma[k] ~ dunif(-10, 10) # Fixed observer effects gamma[k] ~ dnorm(0, 0.1) # Fixed observer effects } gamma[nobs] <- 0 # Effect of last observer set to zero # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + beta[i] + gamma[newobs[i,j]] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), newobs = t(newobs), nsite = nrow(C), nyear = ncol(C), nobs = length(unique(c(newobs)))) # Initial values (NOTE: NA inits for params fixed at zero) inits <- function() list(alpha = runif(235, 1, 3), beta = c(NA, runif(8, -1, 1)), gamma = c(runif(271, -1, 1), NA)) # Parameters monitored params <- c("alpha", "beta", "gamma") # MCMC settings ni <- 6000 nt <- 5 nb <- 1000 nc <- 3 # Call WinBUGS from R (BRT 7 min) out.model2F <- bugs(win.data, inits, params, "model2F.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors for all quantities print(out.model2F, dig = 2) length(which(out.model2F$summary[,8] > 1.2)) ############## # Exercise 6 ############## # Specify model 0 in BUGS language sink("model0.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.site) # Random site effects with grand mean } mu ~ dnorm(0, 0.01) tau.site <- 1/ (sd.site * sd.site) sd.site ~ dunif(0, 5) for (i in 1:nyear){ beta[i] ~ dnorm(0, tau.year) # Random year effects } tau.year <- 1/ (sd.year * sd.year) sd.year ~ dunif(0, 3) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + beta[i] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C)) # Initial values inits <- function() list(mu = runif(1, 0, 4), alpha = runif(235, -2, 2), beta = runif(9, -1, 1)) # Parameters monitored params <- c("mu", "sd.site", "sd.year", "alpha", "beta") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT 6 min) out0 <- bugs(win.data, inits, params, "model0.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors of hyperparameters print(out0$summary[1:3,], dig = 4) # Bundle data (including covariates) forst <- as.numeric(scale(tits$forest)) # forest_standardised forst[is.na(forst)] <- 0 win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C), forst = forst) # Specify model 1 in BUGS language sink("model1.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.site) # Random site effects with grand mean } mu ~ dnorm(0, 0.01) tau.site <- 1/ (sd.site * sd.site) sd.site ~ dunif(0, 5) for (i in 1:nyear){ beta[i] ~ dnorm(0, tau.year) # Random year effects } tau.year <- 1/ (sd.year * sd.year) sd.year ~ dunif(0, 3) gamma1 ~ dnorm(0, 0.01) # Linear effect of forest gamma2 ~ dnorm(0, 0.01) # Quadratic effect of forest # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + beta[i] + gamma1 * forst[j] + gamma2 * pow(forst[j],2) } #j } #i } ",fill = TRUE) sink() # Initial values inits <- function() list(mu = runif(1, 1, 2), alpha = runif(235, -1, 1), beta = runif(9, -1, 1), gamma1 = runif(1), gamma2 = runif(1)) # Parameters monitored params <- c("mu", "sd.site", "sd.year", "alpha", "beta", "gamma1", "gamma2") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # MCMC settings ni <- 120 ; nt <- 2 ; nb <- 20 ; nc <- 3 # Call WinBUGS from R (BRT 7 min) out1 <- bugs(win.data, inits, params, "model1.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors of hyperparameters print(out1$summary[c(1:3, 248:249),], dig = 3) forest.pred.original <- 1:100 forest.pred.st <- (forest.pred.original - mean(tits$forest, na.rm = TRUE)) / sd(tits$forest, na.rm = TRUE) pred.count <- exp(out1$mean$mu + out1$mean$gamma1 * forest.pred.st + out1$mean$gamma2 * forest.pred.st^2) par(mar = c(5,6,3,2), cex.main = 1.2, cex.lab = 1.5, cex.axis = 1.2) plot(forest.pred.original, pred.count, main = "Relationship forest cover and coal tit counts", xlab = "Forest cover (%)", ylab = "Expected count \n(whatever that means biologically)", las = 1, type = "l", col = "blue", lwd = 3, ylim = c(0,25)) (out0$mean$sd.site^2 - out1$mean$sd.site^2) / out0$mean$sd.site^2 # Bundle data (including covariates) elest <- as.numeric(scale(tits$elevation)) # elevation_standardised win.data <- list(C = t(C), nsite = nrow(C), nyear = ncol(C), elest = elest) # Specify model 2 in BUGS language sink("model2.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.site) # Random site effects with grand mean } mu ~ dnorm(0, 0.01) tau.site <- 1/ (sd.site * sd.site) sd.site ~ dunif(0, 5) for (i in 1:nyear){ beta[i] ~ dnorm(0, tau.year) # Random year effects } tau.year <- 1/ (sd.year * sd.year) sd.year ~ dunif(0, 3) gamma3 ~ dnorm(0, 0.01) # Linear effect of elevation gamma4 ~ dnorm(0, 0.01) # Quadratic effect of elevation # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + beta[i] + gamma3 * elest[j] + gamma3 * pow(elest[j],2) } #j } #i } ",fill = TRUE) sink() # Initial values inits <- function() list(mu = runif(1, 1, 2), alpha = runif(235, -1, 1), beta = runif(9, -1, 1), gamma3 = runif(1), gamma4 = runif(1)) # Parameters monitored params <- c("mu", "sd.site", "sd.year", "alpha", "beta", "gamma3", "gamma4") # MCMC settings ni <- 10000 nt <- 2 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT 7 min) out2 <- bugs(win.data, inits, params, "model2.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors of hyperparameters print(out2$summary[c(1:3, 248:249),], dig = 3) elevation.pred.original <- 250:2000 elevation.pred.st <- (elevation.pred.original - mean(tits$elevation)) / sd(tits$elevation) pred.count <- exp(out2$mean$mu + out2$mean$gamma3 * elevation.pred.st) par(mar = c(5,6,3,2), cex.main = 1.2, cex.lab = 1.5, cex.axis = 1.2) plot(elevation.pred.original, pred.count, main = "Relationship elevation and coal tit counts", xlab = "Elevation (m a.s.l.)", ylab = "Expected count", las = 1, type = "l", col = "blue", lwd = 3, ylim = c(0, 25)) (out0$mean$sd.site^2 - out2$mean$sd.site^2) / out0$mean$sd.site^2 ############## # Exercise 7 ############## data <- data.fn(nsite = 10, nyear = 40, sd.site = 0.3, sd.year = 0.2) # Specify model in BUGS language sink("GLMM_Poisson_a.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.alpha) } mu ~ dnorm(0, 0.01) tau.alpha <- 1 / (sd.alpha*sd.alpha) sd.alpha ~ dunif(0, 2) beta ~ dnorm(0, 0.01) tau.year <- 1 / (sd.year*sd.year) sd.year ~ dunif(0, 5) # Likelihood for (i in 1:nyear){ eps[i] ~ dnorm(0, tau.year) for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + beta * year[i] + eps[i] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = data$C, nsite = ncol(data$C), nyear = nrow(data$C), year = (data$year-20) / 20) # Initial values inits <- function() list(mu = runif(1, 0, 2), alpha = runif(data$nsite, -1, 1), beta = runif(1, -1, 1), sd.alpha = runif(1, 0, 0.1), sd.year = runif(1, 0, 0.1)) # Parameters monitored params <- c("mu", "alpha", "beta", "sd.alpha", "sd.year") # MCMC settings ni <- 25000 nt <- 5 nb <- 15000 nc <- 3 # Call WinBUGS from R (BRT 4 min) out.a <- bugs(win.data, inits, params, "GLMM_Poisson_a.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out.a, dig = 2) # Specify model in BUGS language sink("GLM_Poisson_b.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(0, 0.01) } for (i in 2:nyear){ eps[i] ~ dnorm(0, 0.01) } eps[1] <- 0 # This is the corner constraint to avoid overparamaterization # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dpois(lambda[i,j]) lambda[i,j] <- exp(log.lambda[i,j]) log.lambda[i,j] <- alpha[j] + eps[i] } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = data$C, nsite = ncol(data$C), nyear = nrow(data$C)) # Initial values inits <- function() list(alpha = runif(data$nsite, -1, 1), eps = c(NA, runif(39, 0, 1))) # Parameters monitored params <- c("alpha", "eps") # MCMC settings ni <- 25000 nt <- 5 nb <- 15000 nc <- 3 # Call WinBUGS from R (BRT 2 min) out.b <- bugs(win.data, inits, params, "GLM_Poisson_b.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out.b, dig = 2) ############## # Exercise 8 ############## # Specify model in BUGS language sink("GLMM_Normal.txt") cat(" model { # Priors for (j in 1:nsite){ alpha[j] ~ dnorm(mu, tau.alpha) # 4. Random site effects } mu ~ dnorm(0, 0.01) # Hyperparameter 1 tau.alpha <- 1 / (sd.alpha*sd.alpha) # Hyperparameter 2 sd.alpha ~ dunif(0, 500) for (p in 1:3){ beta[p] ~ dnorm(0, 0.01) } tau.res <- pow(sd.res, -2) sd.res ~ dunif(0, 100) # Likelihood for (i in 1:nyear){ for (j in 1:nsite){ C[i,j] ~ dnorm(lambda[i,j], tau.res)# 1. Distribution for random part lambda[i,j] <- alpha[j] + beta[1] * year[i] + beta[2] * pow(year[i],2) + beta[3] * pow(year[i],3) # 3. Linear predictor includes random site effects, link is identity } #j } #i } ",fill = TRUE) sink() # Bundle data win.data <- list(C = data$C, nsite = ncol(data$C), nyear = nrow(data$C), year = (data$year-20) / 20) # Initial values inits <- function() list(mu = runif(1, 0, 2), alpha = runif(data$nsite, -1, 1), beta = runif(3, -1, 1), sd.alpha = runif(1, 0, 1), sd.res = runif(1, 0, 1)) # Parameters monitored (may want to add "lambda") params <- c("mu", "alpha", "beta", "sd.alpha", "sd.res") # MCMC settings (for normal GLM shorter chains suffice) ni <- 2500 nt <- 2 nb <- 500 nc <- 3 # Call WinBUGS from R (BRT <1 min) out.b <- bugs(win.data, inits, params, "GLMM_Normal.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out.b, dig = 2) ############################################################################# # # Chapter 5 # ############################################################################# ############## # Exercise 1 ############## # Data simulation # A) 25 years # Simulate the development of the population n.years <- 25 # Number of years N <- rep(50, n.years) # Simulate detection probability and counts p <- runif(n.years, 0.3, 0.7) y <- numeric() for (t in 1:n.years){ y[t] <- rbinom(1, N[t],p[t]) } # Data analysis # Specify model in BUGS language sink("ssm.bug") cat(" model { # Priors and constraints N.est[1] ~ dunif(0, 500) # Prior for initial population size mean.lambda ~ dunif(0, 10) sigma.proc ~ dunif(0, 10) tau.proc <- pow(sigma.proc, -2) sigma2.proc <- pow(sigma.proc, 2) sigma.obs ~ dunif(0, 100) tau.obs <- pow(sigma.obs, -2) sigma2.obs <- pow(sigma.obs, 2) # State process for (t in 1:(T-1)){ lambda[t] ~ dnorm(mean.lambda, tau.proc) N.est[t+1] <- N.est[t] * lambda[t] } # Observation process for (t in 1:T) { y[t] ~ dnorm(N.est[t], tau.obs) } } ",fill=TRUE) sink() # Bundle data bugs.data <- list(y = y, T = n.years) # Initial values inits <- function(){list(sigma.proc = runif(1, 0, 1), mean.lambda = runif(1, 0.1, 2), sigma.obs = runif(1, 5, 50), N.est = c(runif(1, 20, 40), rep(NA, (n.years-1))))} # Parameters monitored parameters <- c("lambda", "mean.lambda", "sigma2.obs", "sigma2.proc", "N.est") # MCMC settings niter <- 25000 nthin <- 3 nburn <- 10000 nchains <- 3 # Call WinBUGS from R (BRT < 1) ssmex <- bugs(bugs.data, inits, parameters, "ssm.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir) # Define function for visualisation of results graph.ssm <- function(ssm, N, y){ fitted <- lower <- upper <- numeric() n.years <- length(y) for (i in 1:n.years){ fitted[i] <- mean(ssm$sims.list$N.est[,i]) lower[i] <- quantile(ssm$sims.list$N.est[,i], 0.025) upper[i] <- quantile(ssm$sims.list$N.est[,i], 0.975)} m1 <- min(c(y, fitted, N, lower)) m2 <- max(c(y, fitted, N, upper)) par(mar = c(4.5, 4, 1, 1)) plot(0, 0, ylim = c(m1, m2), xlim = c(1, n.years), ylab = "Population size", xlab = "Time", las = 1, col = "black", type = "l", lwd = 2) polygon(x=c(1:n.years,n.years:1), y = c(lower, upper[n.years:1]), col = "grey90", border = "grey90") points(N, type = "l", col = "blue", lwd = 2) points(y, type = "l", col = "red", lwd = 2) points(fitted, type = "l", col = "grey30", lwd = 2) legend(x = 1, y = m2, legend = c("True", "Observed", "Estimated"), lty = c(1, 1, 1),lwd = c(2, 2, 2), col = c("blue", "red", "grey30"), bty = "n", cex = 1.5) } # Apply graph function graph.ssm(ssmex, N, y) # B) 50 years n.years <- 50 # Number of years ############## # Exercise 2 ############## # Read in data set # Load data: House martin population from Magden hm <- c(271, 261, 309, 318, 231, 216, 208, 226, 195, 226, 233, 209, 226, 192, 191, 225, 245, 205, 191, 174) year <- 1990:2009 # Data analysis # Specify model in BUGS language sink("ssm.bug") cat(" model { # Priors and constraints logN.est[1] ~ dnorm(5.6, 0.01) # Prior for initial population size mean.r ~ dnorm(1, 0.001) # Prior for mean growth rate sigma.proc ~ dunif(0, 1) # Prior for sd of state process sigma2.proc <- pow(sigma.proc, 2) tau.proc <- pow(sigma.proc, -2) sigma.obs1 ~ dunif(0, 1) # Prior for sd of obs. process period 1 sigma2.obs1 <- pow(sigma.obs1, 2) tau.obs1 <- pow(sigma.obs1, -2) sigma.obs2 ~ dunif(0, 1) # Prior for sd of obs. process period 2 sigma2.obs2 <- pow(sigma.obs2, 2) tau.obs2 <- pow(sigma.obs2, -2) # State process for (t in 1:(T-1)){ r[t] ~ dnorm(mean.r, tau.proc) logN.est[t+1] <- logN.est[t] + r[t] } # Observation process: the observation error changes for (t in 1:8) { y[t] ~ dnorm(logN.est[t], tau.obs1) } for (t in 9:T) { y[t] ~ dnorm(logN.est[t], tau.obs2) } # Population sizes on real scale for (t in 1:T) { N.est[t] <- exp(logN.est[t]) } } ",fill=TRUE) sink() # Bundle data bugs.data <- list(y = log(hm), T = length(year)) # Initial values inits <- function(){list(sigma.proc = runif(1, 0, 1), mean.r = rnorm(1), sigma.obs1 = runif(1, 0, 1), sigma.obs2 = runif(1, 0, 1), logN.est = c(rnorm(1, 5.6, 0.1), rep(NA, (length(year)-1))))} # Parameters monitored parameters <- c("r", "mean.r", "sigma2.obs1", "sigma2.obs2", "sigma2.proc", "N.est") # MCMC settings niter <- 50000 nthin <- 3 nburn <- 25000 nchains <- 3 # Call WinBUGS from R (BRT < 1) hm.562 <- bugs(bugs.data, inits, parameters, "ssm.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(hm.562, digits = 3) # Here is the second solution with the GLM formulation: # Specify model in BUGS language sink("ssm.bug") cat(" model { # Priors and constraints logN.est[1] ~ dnorm(5.6, 0.01) # Prior for initial population size mean.r ~ dnorm(1, 0.001) # Prior for mean growth rate sigma.proc ~ dunif(0, 1) # Prior for sd of state process sigma2.proc <- pow(sigma.proc, 2) tau.proc <- pow(sigma.proc, -2) for (i in 1:2){ sigma.obs[i] ~ dunif(0, 100) # Priors for sd of obs proccesses tau.obs[i] <- pow(sigma.obs[i], -2) sigma2.obs[i] <- pow(sigma.obs[i], 2) } # State process for (t in 1:(T-1)){ r[t] ~ dnorm(mean.r, tau.proc) logN.est[t+1] <- logN.est[t] + r[t] } # Observation process: the observation error changes for (t in 1:T) { y[t] ~ dnorm(logN.est[t], tau.obs[period[t]]) } # Population sizes on real scale for (t in 1:T) { N.est[t] <- exp(logN.est[t]) } } ",fill=TRUE) sink() # Bundle data bugs.data <- list(y = log(hm), T = length(year), period = c(rep(1, 8), rep(2, 12))) # Initial values inits <- function(){list(sigma.proc = runif(1, 0, 1), mean.r = rnorm(1), sigma.obs = runif(2, 0, 1), logN.est = c(rnorm(1, 5.6, 0.1), rep(NA, (length(year)-1))))} # Parameters monitored parameters <- c("r", "mean.r", "sigma2.obs", "sigma2.proc", "N.est") # MCMC settings niter <- 50000 nthin <- 3 nburn <- 25000 nchains <- 3 # Call WinBUGS from R (BRT 2 min) hm.562alt <- bugs(bugs.data, inits, parameters, "ssm.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(hm.562alt, digits = 3) ############## # Exercise 3 ############## # Read in the data # Load data peregrine <- read.table("falcons.txt", header = TRUE) # Data analysis # Specify model in BUGS language sink("ssm.bug") cat(" model { # Priors and constraints N.est[1] ~ dunif(0, 200) # Prior for initial population size mean.lambda ~ dunif(0, 2) # Prior for mean growth rate sigma.proc ~ dunif(0, 5) # Prior sd of state process sigma2.proc <- pow(sigma.proc, 2) tau.proc <- pow(sigma.proc, -2) sigma.obs ~ dunif(0, 20) # Prior sd of observation process sigma2.obs <- pow(sigma.obs, 2) tau.obs <- pow(sigma.obs, -2) # State process for (t in 1:(T-1)){ lambda[t] ~ dnorm(mean.lambda, tau.proc) N.est[t+1] <- N.est[t] * lambda[t] } # Observation process for (t in 1:T) { y[t] ~ dnorm(N.est[t], tau.obs) } } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = peregrine$Pairs, T = length(peregrine$Pairs)) # Initial values inits <- function(){list(sigma.proc = runif(1, 0, 1), mean.lambda = runif(1, 0.9, 1.1), sigma.obs = runif(1, 0.5, 5), N.est = c(runif(1, 20, 40), rep(NA, (length(peregrine$Pairs)-1))))} # Parameters monitored parameters <- c("lambda", "mean.lambda", "sigma2.obs", "sigma2.proc", "N.est") # MCMC settings niter <- 100000 nthin <- 10 nburn <- 50000 nchains <- 3 # Call WinBUGS from R (BRT 8 min) per1 <- bugs(bugs.data, inits, parameters, "ssm.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir) print(per1, 3) # Specify model in BUGS language sink("ssm.bug") cat(" model { # Priors and constraints logN[1] ~ dnorm(0, 0.01) # Prior for initial population size mean.r ~ dnorm(0, 0.01) # Prior for mean stochastic growth rate sigma.proc ~ dunif(0, 2) # Prior sd of state process sigma2.proc <- pow(sigma.proc, 2) tau.proc <- pow(sigma.proc, -2) sigma.obs ~ dunif(0, 10) # Prior sd of observation process sigma2.obs <- pow(sigma.obs, 2) tau.obs <- pow(sigma.obs, -2) # State process for (t in 1:(T-1)){ r[t] ~ dnorm(mean.r, tau.proc) logN[t+1] <- logN[t] + r[t] } # Observation process for (t in 1:T) { y[t] ~ dnorm(logN[t], tau.obs) N.est[t] <- exp(logN[t]) # Backtransformation from log-scale } } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = log(peregrine$Pairs), T = length(peregrine$Pairs)) # Initial values inits <- function(){list(sigma.proc = runif(1, 0, 1), mean.r = rnorm(1), sigma.obs = runif(1, 0.5, 5), logN = c(runif(1, -1, 5), rep(NA, (length(peregrine$Pairs)-1))))} # Parameters monitored parameters <- c("r", "mean.r", "sigma2.obs", "sigma2.proc", "N.est") # MCMC settings niter <- 100000 nthin <- 10 nburn <- 50000 nchains <- 3 # Call WinBUGS from R (BRT 8 min) per2 <- bugs(bugs.data, inits, parameters, "ssm.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir) print(per2, 3) plot(peregrine$Pairs, type = "l", ylab = "Population size", xlab = "Year", lwd = 2, las = 1) points(per1$mean$N.est, type = "l", col = "blue", lwd = 2) points(per2$mean$N.est, type = "l", col = "blue", lwd = 2, lty = 2) legend(y = 200, x = 1, legend = c("Counts", "SSM", "SSM on log scale"), lty = c(1, 1, 2), col = c("black", "blue", "blue"), bty = "n", lwd = rep(2, 3)) # Specify model in BUGS language sink("ssm.bug") cat(" model { # Priors and constraints logN[1] ~ dnorm(0, 0.01) # Prior for initial population size alpha ~ dnorm(0, 0.01) # Prior for mean stochastic growth rate beta ~ dnorm(0, 0.01) # Prior for mean stochastic growth rate sigma.proc ~ dunif(0, 2) # Prior sd of state process sigma2.proc <- pow(sigma.proc, 2) tau.proc <- pow(sigma.proc, -2) sigma.obs ~ dunif(0, 10) # Prior sd of observation process sigma2.obs <- pow(sigma.obs, 2) tau.obs <- pow(sigma.obs, -2) # State process for (t in 1:(T-1)){ r[t] ~ dnorm(mu.r[t], tau.proc) mu.r[t] <- alpha + beta * t logN[t+1] <- logN[t] + r[t] } # Observation process for (t in 1:T) { y[t] ~ dnorm(logN[t], tau.obs) N.est[t] <- exp(logN[t]) # Backtransformation from log-scale } } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = log(peregrine$Pairs), T = length(peregrine$Pairs)) # Initial values inits <- function(){list(alpha = rnorm(1), beta = rnorm(1), sigma.obs = runif(1, 0.5, 5), sigma.proc = runif(1, 0.5, 2), logN = c(runif(1, -1, 5), rep(NA, (length(peregrine$Pairs)-1))))} # Parameters monitored parameters <- c("r", "alpha", "beta", "sigma2.obs", "sigma2.proc", "N.est") # MCMC settings niter <- 100000 nthin <- 10 nburn <- 50000 nchains <- 3 # Call WinBUGS from R (BRT 8 min) per3 <- bugs(bugs.data, inits, parameters, "ssm.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir) print(per3, 3) # Specify model in BUGS language sink("ssmTrendError.bug") cat(" model { # Priors and constraints logN[1] ~ dnorm(3.5, 100) # Prior for initial population size mean.r ~ dnorm(0, 0.01)I(-20,20) beta ~ dnorm(0, 0.01)I(-20,20) # Prior for slope mean.err ~ dnorm(0, 0.01) I(-20,20) # Prior for mean log(obs error) sigma.proc ~ dunif(0, 5) # Prior sd of state process sigma2.proc <- pow(sigma.proc, 2) tau.proc <- pow(sigma.proc, -2) # State process for (t in 1:(T-1)){ r[t] ~ dnorm(mean.r, tau.proc) logN[t+1] <- logN[t] + r[t] } # Observation process for (t in 1:T) { y[t] ~ dnorm(logN[t], tau.obs[t]) tau.obs[t] <- 1/sigma2.obs[t] log(sigma2.obs[t]) <- mean.err + beta * period[t] N.est[t] <- exp(logN[t]) # Backtransformation from log-scale } } ",fill=TRUE) sink() # Bundle data bugs.data <- list(y = log(peregrine$Pairs), period = scale(1: length(peregrine$Pairs))[,1], T = length(peregrine$Pairs)) # Initial values inits <- function(){list(sigma.proc = runif(1, 0, 1), mean.r = rnorm(1), beta = runif(1, -1, 1), mean.err = rnorm(1), logN = c(runif(1, 3, 4), rep(NA, (length(peregrine$Pairs)-1))))} # Parameters monitored parameters <- c("r", "mean.r", "mean.err", "sigma2.obs", "sigma2.proc", "N.est", "beta") # MCMC settings niter <- 100000 nthin <- 10 nburn <- 50000 nchains <- 3 # Call WinBUGS from R (BRT 8 min) per4 <- bugs(bugs.data, inits, parameters, "ssmTrendError.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir) print(per4, 3) ############################################################################# # # Chapter 6 # ############################################################################# ############## # Exercise 1 ############## data <- data.fn(N = 100, mean.p = 0.4, T = 5, sd = 1) attach(data) # Augment data set nz <- 100 yaug <- rbind(yobs, array(0, dim=c(nz, T))) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors omega ~ dunif(0, 1) mean.lp <- logit(mean.p) mean.p ~ dunif(0, 1) tau <- 1 / (sd * sd) sd ~ dunif(0, 5) # Likelihood for (i in 1:M){ z[i] ~ dbern(omega) logit(p[i]) <- eps[i] eps[i] ~ dnorm(mean.lp, tau)I(-16, 16) for(j in 1:T){ p.eff[i, j] <- z[i] * p[i] y[i, j] ~ dbern(p.eff[i, j]) } #j } #i # Derived quantities N <- sum(z[]) } ",fill = TRUE) sink() # Bundle data win.data <- list(y = yaug, M = nrow(yaug), T = ncol(yaug)) # Initial values inits <- function() list(z = rep(1, nrow(yaug)), sd = runif(1, 0.1, 0.9)) # Parameters monitored params <- c("N", "mean.p", "sd", "omega") # MCMC settings ni <- 25000 nt <- 2 nb <- 5000 nc <- 3 # Call WinBUGS from R (BRT 4 min) out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out, dig = 3) hist(out$sims.list$N, nclass = 50, col = "gray", main = "", xlab = "Population size", las = 1, xlim = c(80, 200)) abline(v = data$C, col = "black", lwd = 3) ############## # Exercise 2 ############## # Construct the array X, which indicates capture ever before # of individual i at occasion j (solution due to Tomas Telensky) X <- as.matrix(y) for (i in 1:nrow(y)){ seen.before <- 0 for (j in 1:ncol(y)){ X[i,j] <- seen.before if (y[i,j] == 1) seen.before <- 1 } } # Check whether X correct (it is) head(y) head(X) # Bundle data, including array X win.data <- list(y = as.matrix(y), M = nrow(y), T = ncol(y), X = as.matrix(X)) # Specify model in BUGS language sink("M_tbh.txt") cat(" model { # Priors omega ~ dunif(0, 1) for (j in 1:T){ alpha[j] <- log(mean.p[j] / (1-mean.p[j])) # Define logit mean.p[j] ~ dunif(0, 1) # Detection intercepts } gamma ~ dnorm(0, 0.01) tau <- 1 / (sd * sd) sd ~ dunif(0, 5) # Likelihood for (i in 1:M){ z[i] ~ dbern(omega) eps[i] ~ dnorm(0, tau)I(-16, 16) # First occasion: no term for recapture (gamma) y[i,1] ~ dbern(p.eff[i,1]) p.eff[i,1] <- z[i] * p[i,1] p[i,1] <- 1 / (1 + exp(-lp[i,1])) lp[i,1] <- alpha[1] + eps[i] # All subsequent occasions: includes recapture term (gamma) for (j in 2:T){ y[i,j] ~ dbern(p.eff[i,j]) p.eff[i,j] <- z[i] * p[i,j] p[i,j] <- 1 / (1 + exp(-lp[i,j])) lp[i,j] <- alpha[j] + eps[i] + gamma * X[i,j] ## Only change } #j } #i # Derived quantities N <- sum(z[]) } ",fill = TRUE) sink() # Initial values inits <- function() list(z = rep(1, nrow(y)), sd = runif(1, 0.1, 1)) # Parameters monitored params <- c("N", "mean.p", "gamma", "sd", "omega") # MCMC settings ni <- 50000 nt <- 4 nb <- 10000 nc <- 3 # Call WinBUGS from R (BRT 33 min) out <- bugs(win.data, inits, params, "M_tbh.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors and plot posteriors of N and gamma print(out, dig = 2) par(mfrow = c(2,1)) hist(out$sims.list$N, breaks = 100, col = "gray", main = "", xlab = "Community size (N)", las = 1, xlim = c(30, 100), freq = FALSE) abline(v = C, col = "black", lwd = 3) hist(out$sims.list$gamma, breaks = 100, col = "gray", main = "", xlab = "Behavioural response effect (gamma)", las = 1, xlim = c(-4, 4), freq = FALSE) abline(v = 0, col = "black", lwd = 3) mean(out$sims.list$gamma < 0) ############## # Exercise 3 ############## # New definition of data simulation function for model Mh data.fn <- function(N = 100, mean.p = 0.4, T = 5, sd = 1){ yfull <- yobs <- array(NA, dim = c(N, T)) mean.lp <- log(mean.p / (1-mean.p)) p.vec <- plogis(mean.lp+ rnorm(N, 0, sd)) for (i in 1:N){ yfull[i,] <- rbinom(n = T, size = 1, prob = p.vec[i]) } ever.detected <- apply(yfull, 1, max) C <- sum(ever.detected) yobs <- yfull[ever.detected == 1,] yfreq <- sort(apply(yobs, 1, sum), decreasing = TRUE) cat(C, "out of", N, "animals present were detected.\n") hist(p.vec, xlim = c(0,1), nclass = 20, col = "gray", main = "", xlab = "Detection probability", las = 1) return(list(N = N, p.vec = p.vec, mean.lp = mean.lp, C = C, T = T, yfull = yfull, yobs = yobs, yfreq = yfreq)) } # Define a function to do data augmentation, run analysis using Mh and return results all at once model.fn <- function(nz = 200, ni = 25000, nt = 2, nb = 5000, nc = 3, data.file = data, debg = FALSE){ # Function arguments: # nz -- number of fake DA individuals # ni/nt/nb/nc -- MCMC settings # data.file -- name of the object with the simulated data # debg -- setting of DEBUG argument in bugs() # Do data augmentation and bundle data yaug <- c(data.file$yfreq, rep(0, nz)) win.data <- list(y = yaug, M = length(yaug), T = data.file$T) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors omega ~ dunif(0, 1) mean.lp <- logit(mean.p) mean.p ~ dunif(0, 1) tau <- pow(sd, -2) sd ~ dunif(0, 5) # Might have to be be adapted depending on your data set # Likelihood for (i in 1:M){ z[i] ~ dbern(omega) logit(p[i]) <- eps[i] eps[i] ~ dnorm(mean.lp, tau)I(-16, 16) p.eff[i] <- z[i] * p[i] y[i] ~ dbin(p.eff[i], T) } # Derived quantities N <- sum(z[]) } ",fill = TRUE) sink() # Initial values inits <- function() list(z = rep(1, length(yaug)), sd = runif(1, 0.1, 0.9)) # Parameters monitored params <- c("N", "mean.p", "sd", "omega") # Call WinBUGS from R (BRT 6 min) out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = debg, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out, dig = 3) hist(out$sims.list$N, nclass = 50, col = "gray", main = "", xlab = "Population size", las = 1, xlim = c(0, 1.5*nz)) abline(v = data$C, col = "black", lwd = 3) return(post.estimates = out$summary) } data <- data.fn(N = 20, mean.p = 0.5, T = 5, sd = 1) estimates <- model.fn(nz = 150, ni = 250, nt = 2, nb = 50, nc = 3, data.file = data, debg = TRUE) # Set up data structures to hold the results simreps <- 100 data.sets <- array(NA, dim = c(20, simreps)) solutions <- array(NA, dim = c(5, 9, simreps)) rownames(solutions) <- rownames(estimates) colnames(solutions) <- colnames(estimates) # Run data generation/data analysis cycle simrep times for (i in 1:simreps){ cat(paste("\n\n*** SimRep", i, "***\n")) data <- data.fn(N = 20, mean.p = 0.5, T = 5, sd = 1) data.sets[1:data$C,i] <- data$yfreq estimates <- model.fn(nz = 50, ni = 25000, nt = 5, nb = 10000, nc = 3, data.file = data, debg = FALSE) solutions[,,i] <- estimates } # Get the number of observed individuals as a simple estimator of N Nobs <- array(NA, dim = simreps) for (i in 1:simreps){ Nobs[i] <- sum(!is.na(data.sets[,i])) } # Summarize simulation results par(mfrow = c(3, 3)) hist(Nobs, breaks = 25, col = "grey", xlab = "Estimates of N", main = "N: # Observed individuals", xlim = c(0, 50), las = 1) abline(v = 20, col = "red", lwd = 3) abline(v = mean(Nobs), col = "blue", lwd = 3) hist(solutions[1,1,], breaks = 25, col = "grey", xlab = "Estimates of N", main = "N: posterior mean", xlim = c(0, 50), las = 1) abline(v = 20, col = "red", lwd = 3) abline(v = mean(solutions[1,1,]), col = "blue", lwd = 3) hist(solutions[1,5,], breaks = 25, col = "grey", xlab = "Estimates of N", main = "N: posterior median", xlim = c(0, 50), las = 1) abline(v = 20, col = "red", lwd = 3) abline(v = mean(solutions[1,5,]), col = "blue", lwd = 3) hist(solutions[5,1,], breaks = 25, col = "grey", xlab = "Estimates of mean.p", main = "mean.p: not available", xlim = c(0, 1), las = 1) hist(solutions[2,1,], breaks = 25, col = "grey", xlab = "Estimates of mean.p", main = "mean.p: posterior mean", xlim = c(0, 1), las = 1) abline(v = 0.5, col = "red", lwd = 3) abline(v = mean(solutions[2,1,]), col = "blue", lwd = 3) hist(solutions[2,5,], breaks = 25, col = "grey", xlab = "Estimates of mean.p", main = "mean.p: posterior median", xlim = c(0, 1), las = 1) abline(v = 0.5, col = "red", lwd = 3) abline(v = mean(solutions[2,5,]), col = "blue", lwd = 3) hist(solutions[5,1,], breaks = 25, col = "grey", xlab = "Estimates of sd", main = "sd: not available", xlim = c(0, 1), las = 1) hist(solutions[3,1,], breaks = 25, col = "grey", xlab = "Estimates of sd", main = "sd: posterior mean", xlim = c(0, 5), las = 1) abline(v = 1, col = "red", lwd = 3) abline(v = mean(solutions[3,1,]), col = "blue", lwd = 3) hist(solutions[3,5,], breaks = 25, col = "grey", xlab = "Estimates of sd", main = "sd: posterior median", xlim = c(0, 5), las = 1) abline(v = 1, col = "red", lwd = 3) abline(v = mean(solutions[3,5,]), col = "blue", lwd = 3) ############## # Exercise 4 ############## par(mfrow = c(2,2)) data1 <- data.fn(N = 1000, mean.p = 0.2, T = 3, sd = 1) data2 <- data.fn(N = 1000, mean.p = 0.6, T = 3, sd = 1) data3 <- data.fn(N = 1000, mean.p = 0.2, T = 3, sd = 5) data4 <- data.fn(N = 1000, mean.p = 0.6, T = 3, sd = 5) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors omega ~ dunif(0, 1) p ~ dunif(0, 1) # Likelihood for (i in 1:M){ z[i] ~ dbern(omega) # Inclusion indicators for (j in 1:T){ yaug[i,j] ~ dbern(p.eff[i,j]) p.eff[i,j] <- z[i] * p # Can only be detected if z=1 } #j } #i # Derived quantities N <- sum(z[]) } ",fill = TRUE) sink() # Initial values inits <- function() list(z = rep(1, nrow(yaug)), p = runif(1, 0, 1)) # Parameters monitored params <- c("N", "p", "omega") # MCMC settings ni <- 2500 nt <- 2 nb <- 500 nc <- 3 # Augment data set and bundle data nz <- 1000 yaug <- rbind(data1$yobs, array(0, dim = c(nz, data1$T))) win.data <- list(yaug = yaug, M = nrow(yaug), T = ncol(yaug)) # Call WinBUGS from R (BRT <1 min) out1 <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out1, dig = 2) # Augment data set and bundle data nz <- 1000 yaug <- rbind(data2$yobs, array(0, dim = c(nz, data2$T))) win.data <- list(yaug = yaug, M = nrow(yaug), T = ncol(yaug)) # Call WinBUGS from R (BRT <1 min) out2 <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out2, dig = 2) # Augment data set and bundle data nz <- 1000 yaug <- rbind(data3$yobs, array(0, dim = c(nz, data3$T))) win.data <- list(yaug = yaug, M = nrow(yaug), T = ncol(yaug)) # Call WinBUGS from R (BRT <1 min) out3 <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out3, dig = 2) # Augment data set and bundle data nz <- 1000 yaug <- rbind(data4$yobs, array(0, dim = c(nz, data4$T))) win.data <- list(yaug = yaug, M = nrow(yaug), T = ncol(yaug)) # Call WinBUGS from R (BRT <1 min) out4 <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out4, dig = 2) # Compare posteriors for N for 4 data sets par(mfrow = c(2, 2)) hist(out1$sims.list$N, nclass = 50, col = "gray", main = "mean.p = 0.2, sd = 1", xlab = "Population size (N)", las = 1, xlim = c(500, 1000)) abline(v = data1$C, lwd = 2, col = "green") abline(v = 1000, lwd = 2, col = "red") hist(out2$sims.list$N, nclass = 50, col = "gray", main = "mean.p = 0.6, sd = 1", xlab = "Population size (N)", las = 1, xlim = c(500, 1000)) abline(v = data2$C, lwd = 2, col = "green") abline(v = 1000, lwd = 2, col = "red") hist(out3$sims.list$N, nclass = 50, col = "gray", main = "mean.p = 0.2, sd = 5", xlab = "Population size (N)", las = 1, xlim = c(500, 1000)) abline(v = data3$C, lwd = 2, col = "green") abline(v = 1000, lwd = 2, col = "red") hist(out4$sims.list$N, nclass = 50, col = "gray", main = "mean.p = 0.6, sd = 5", xlab = "Population size (N)", las = 1, xlim = c(500, 1000)) abline(v = data4$C, lwd = 2, col = "green") abline(v = 1000, lwd = 2, col = "red") plot(1:4, c(out1$mean$p, out2$mean$p, out3$mean$p, out4$mean$p), ylim = c(0, 1), col = c("blue", "brown", "blue", "brown"), pch = 16, cex = 1.2, xlab = "Data set number", ylab = "Detection probability", cex.lab = 1.2, cex.axis = 1.2) abline(h = c(0.2, 0.6), col = c("blue", "brown")) segments(1, out1$summary[2,3], 1, out1$summary[2,7], col = "blue") segments(2, out2$summary[2,3], 2, out2$summary[2,7], col = "brown") segments(3, out3$summary[2,3], 3, out3$summary[2,7], col = "blue") segments(4, out4$summary[2,3], 4, out4$summary[2,7], col = "brown") ############## # Exercise 5 ############## # Define function to simulate data under Mbt with T = 2 fixed data.fn <- function(N = 100, p1 = 0.3, c = 0.2, deltaT2 = 0.4){ yfull <- yobs <- array(NA, dim = c(N, 2) ) # First capture occasion yfull[,1] <- rbinom(n = N, size = 1, prob = p1) # Second capture occasions p2 <- p1 + deltaT2 + yfull[,1]*c yfull[,2] <- rbinom(n = N, size = 1, prob = p2) ever.detected <- apply(yfull, 1, max) C <- sum(ever.detected) yobs <- yfull[ever.detected == 1,] cat(C, "out of", N, "animals present were detected.\n") return(list(N = N, p1 = p1, p2 = p2, c = c, deltaT2 = deltaT2, C = C, T = 2, yfull = yfull, yobs = yobs)) } str(data <- data.fn(N = 200, p1 = 0.3, c = 0.2, deltaT2 = 0.4)) # Specify model M_bt in BUGS language sink("model.txt") cat(" model { # Priors omega ~ dunif(0, 1) p1 ~ dunif(0, 1) # Cap prob during 1st occasion c ~ dunif(-1, 1) # Diff. in cap.prob. during 2nd occasion # if captured during 1st deltaT2 ~ dunif(-1, 1)# Diff. in cap.prob. during 2nd occasion # Likelihood for (i in 1:M){ z[i] ~ dbern(omega) # First occasion yaug[i,1] ~ dbern(p1.eff[i,1]) p1.eff[i,1] <- z[i] * p1 # Second occasion yaug[i,2] ~ dbern(p2.eff[i,2]) p2.eff[i,2] <- z[i] * (p1 + deltaT2 + yaug[i,1] * c) } #i # Derived quantities N <- sum(z[]) } # end model ",fill = TRUE) sink() # Initial values (chose random starts for z, but fix for params) inits <- function() list(z = round(runif(nrow(yaug), 0, 1)), p1 = 0.5, c = 0.1, deltaT2 = 0.2) # Parameters monitored params <- c("N", "p1", "c", "deltaT2", "omega") # MCMC settings ni <- 4000 nt <- 1 nb <- 3000 nc <- 3 data <- data.fn(N = 1000, p1 = 0.3, c = 0.2, deltaT2 = 0.4) # Augment data set and bundle data nz <- 500 yaug <- rbind(data$yobs, array(0, dim = c(nz, data$T))) win.data <- list(yaug = yaug, M = nrow(yaug)) # TRY initials right on spot (chose random starts for z) inits <- function() list(z = round(runif(nrow(yaug), 0, 1)), p1 = 0.4, c = 0.1, deltaT2 = 0.2) # MCMC settings ni <- 43000 nt <- 4 nb <- 3000 nc <- 3 # Call WinBUGS from R (BRT 1 min) out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out, dig = 3) ############## # Exercise 6 ############## data <- data.fn(N = 10000, T = 2, p = 0.3, c = 0.4) data <- data.fn(N = 10000, T = 2, p = 0.6, c = 0.2) data <- data.fn(N = 10000, T = 2, p = 0.1, c = 0.9) ############## # Exercise 7 ############## # Define function to simulate data under Mt with random time effects data.fn <- function(N = 100, mean.p = 0.5, T = 3, sd.lp = 1){ yfull <- yobs <- array(NA, dim = c(N, T)) mean.lp <- log(mean.p / (1 - mean.p)) p.vec <- plogis(rnorm(T, mean.lp, sd.lp)) # Draw p for each T for (j in 1:T){ yfull[,j] <- rbinom(n = N, size = 1, prob = p.vec[j]) } ever.detected <- apply(yfull, 1, max) C <- sum(ever.detected) yobs <- yfull[ever.detected == 1,] cat(C, "out of", N, "animals present were detected.\n") return(list(N = N, p.vec = p.vec, C = C, T = T, yfull = yfull, yobs = yobs, mean.p = mean.p, mean.lp = mean.lp, sd.lp = sd.lp)) } str(data <- data.fn(N = 100, mean.p = 0.3, T = 3, sd.lp = 1)) str(data <- data.fn(N = 100, mean.p = 0.1, T = 10, sd.lp = 0.5)) # Augment data set nz <- 150 yaug <- rbind(data$yobs, array(0, dim = c(nz, data$T))) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors omega ~ dunif(0, 1) mean.lp ~ dnorm(0, 0.001) tau.lp <- pow(sd.lp, -2) sd.lp ~ dunif(0, 3) # Random effects distribution taken outside of loop below, to avoid # multiple definitions of p[j] for (j in 1:T){ logit(p[j]) <- lp[j] lp[j] ~ dnorm(mean.lp, tau.lp)I(-16, 16) } #j # Likelihood for (i in 1:M){ z[i] ~ dbern(omega) for (j in 1:T){ yaug[i,j] ~ dbern(p.eff[i,j]) p.eff[i,j] <- z[i] * p[j] } #j } #i # Derived quantities N <- sum(z[]) } # end model ",fill = TRUE) sink() # Bundle data win.data <- list(yaug = yaug, M = nrow(yaug), T = ncol(yaug)) # Initial values inits <- function() list(z = rep(1, nrow(yaug)), mean.lp = 0, sd.lp = runif(1, 0, 1)) # Parameters monitored params <- c("N", "p", "mean.lp", "sd.lp", "omega") # MCMC settings ni <- 2500 nt <- 2 nb <- 500 nc <- 3 # Call WinBUGS from R (BRT 1 min) out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors in table and graphs print(out, dig = 2) # Observed value of N and estimate hist(out$sims.list$N, nclass = 40, col = "gray", main = "", xlab = "Population size", las = 1, xlim = c(70, 150)) abline(v = data$C, col = "black", lwd = 3) cbind(data$p.vec, out$summary[2:11,c(1:3, 7)]) ############## # Exercise 8 ############## # Data generation from Exercise 7 data.fn <- function(N = 100, mean.p = 0.5, T = 3, sd.lp = 1){ yfull <- yobs <- array(NA, dim = c(N, T)) mean.lp <- log(mean.p / (1 - mean.p)) p.vec <- plogis(rnorm(T, mean.lp, sd.lp)) # Draw p for each T for (j in 1:T){ yfull[,j] <- rbinom(n = N, size = 1, prob = p.vec[j]) } ever.detected <- apply(yfull, 1, max) C <- sum(ever.detected) yobs <- yfull[ever.detected == 1,] cat(C, "out of", N, "animals present were detected.\n") return(list(N = N, p.vec = p.vec, C = C, T = T, yfull = yfull, yobs = yobs, mean.p = mean.p, mean.lp = mean.lp, sd.lp = sd.lp)) } str(data <- data.fn(N = 10000, mean.p = 0.2, T = 5, sd.lp = 1)) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors omega ~ dunif(0, 1) p ~ dunif(0, 1) # Likelihood for (i in 1:M){ z[i] ~ dbern(omega) for (j in 1:T){ yaug[i,j] ~ dbern(p.eff[i,j]) p.eff[i,j] <- z[i] * p } #j } #i # Derived quantities N <- sum(z[]) } ",fill = TRUE) sink() # Initial values inits <- function() list(z = rep(1, nrow(yaug)), p = runif(1, 0, 1)) # Parameters monitored params <- c("N", "p", "omega") # MCMC settings ni <- 2500 nt <- 2 nb <- 500 nc <- 3 # Augment data set and bundle data nz <- 3000 yaug <- rbind(data$yobs, array(0, dim = c(nz, data$T))) win.data <- list(yaug = yaug, M = nrow(yaug), T = ncol(yaug)) # Call WinBUGS from R (BRT 15 min) out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors in table and sketch them in graphs print(out, dig = 2) par(mfrow = c(2, 1)) hist(out$sims.list$N, nclass = 50, col = "gray", main = "", xlab = "Population size (N)", las = 1, xlim = c(8000, 12000)) abline(v = data$C, lwd = 2, col = "green") abline(v = 10000, lwd = 2, col = "red") hist(out$sims.list$p, nclass = 50, col = "gray", main = "", xlab = "Detection probability (p)", las = 1, xlim = c(0.2, 0.3)) #abline(v = data$mean.p, lwd = 2, col = "green") abline(v = mean(data$p.vec), lwd = 2, col = "red") str(data <- data.fn(N = 1000, mean.p = 0.4, T = 3, sd.lp = 1)) # MCMC settings ni <- 1200 nt <- 1 nb <- 200 nc <- 3 # Augment data set and bundle data nz <- 1500 yaug <- rbind(data$yobs, array(0, dim = c(nz, data$T))) win.data <- list(yaug = yaug, M = nrow(yaug), T = ncol(yaug)) # Call WinBUGS from R (BRT 15 min) out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors in table and sketch them in graphs print(out, dig = 2) par(mfrow = c(2, 1)) hist(out$sims.list$N, nclass = 50, col = "gray", main = "", xlab = "Population size (N)", las = 1, xlim = c(500, 1700)) abline(v = data$C, lwd = 2, col = "green") abline(v = 1000, lwd = 2, col = "red") hist(out$sims.list$p, nclass = 50, col = "gray", main = "", xlab = "Detection probability (p)", las = 1, xlim = c(0.0, 0.3)) #abline(v = data$mean.p, lwd = 2, col = "green") abline(v = mean(data$p.vec), lwd = 2, col = "red") ############## # Exercise 9 ############## p610 <- read.table("p610.txt", header = TRUE) y <- as.matrix(p610[,5:9]) # Grab counts and convert to matrix y[y > 1] <- 1 # Convert to det-nondetections dimnames(y) <- NULL # Specify model in BUGS language sink("M_t+X.txt") cat(" model { # Priors omega ~ dunif(0, 1) for (j in 1:T){ alpha[j] <- log(mean.p[j] / (1-mean.p[j])) mean.p[j] ~ dunif(0, 1) } beta ~ dnorm(0, 0.01) # Likelihood for (i in 1:Nspec){ # Loop over individuals z[i] ~ dbern(omega) for (j in 1:T){ # Loop over occasions y[i,j] ~ dbern(p.eff[i,j]) p.eff[i,j] <- z[i] * p[i,j] p[i,j] <- 1 / (1 + exp(-lp[i,j])) lp[i,j] <- alpha[j] + beta * size[i] } #j } #i # Derived quantities N <- sum(z[]) } ",fill = TRUE) sink() # Bundle data win.data <- list(y = y, size = log(p610$bm)-mean(log(p610$bm)), Nspec = nrow(y), T = ncol(y)) # Initial values inits <- function() list(z = rep(1, nrow(y)), beta = runif(1, 0, 1)) # Parameters monitored params <- c("N", "mean.p", "beta", "omega") # MCMC settings ni <- 5000 nt <- 2 nb <- 1000 nc <- 3 # Call WinBUGS from R (BRT 1 min) outX2 <- bugs(win.data, inits, params, "M_t+X.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors and plot posterior for N print(outX2, dig = 3) par(mfrow = c(1, 2)) hist(outX2$sims.list$N, breaks = 100, col = "gray", main = "", xlab = "Community size", las = 1, xlim = c(30, 100), freq = FALSE) abline(v = 31, col = "black", lwd = 3) pred.wt <- seq(5, 2000, length.out = 100) # Cov. vals for prediction pred.wt.st <- log(pred.wt)- mean(log(p610$bm)) # Transform them in the same was as in the analysis pred.p<- plogis(log(mean(outX2$mean$mean.p)/(1- mean(outX2$mean$mean.p))) + outX2$mean$beta * pred.wt.st) # Compute predicted response plot(pred.wt, pred.p, type = "l", lwd = 3, col = "blue", las = 1, frame.plot = FALSE, ylim = c(0, 0.5)) ############################################################################# # # Chapter 7 # ############################################################################# ############## # Exercise 1 ############## # Specify model in BUGS language sink("cjs-c-c.bug") cat(" model { # Priors and constraints mean.phi ~ dunif(0, 1) # Prior for mean survival mean.p ~ dunif(0, 1) # Prior for mean recapture # Likelihood for (i in 1:nind){ # Define latent state at first capture z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- mean.phi * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- mean.p * z[i,t] } #t } #i } ",fill = TRUE) sink() ############## # Exercise 2 ############## # Data simulation # Define the parameters n.occasions <- 15 # Number of capture occasions marked <- rep(30, n.occasions-1) # Number of newly marked individuals phi.m <- rep(0.6, n.occasions-1) phi.f <- rep(0.5, n.occasions-1) p <- rep(0.4, n.occasions-1) # Define a matrix with the survival and recapture probabilities PHI.M <- matrix(rep(phi.m, sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE) PHI.F <- matrix(rep(phi.f, sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE) P <- matrix(rep(p, sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE) # Apply simulation function CH.M <- simul.cjs(PHI.M, P, marked) CH.F <- simul.cjs(PHI.F, P, marked) # Create the m-arrays from the capture-histories marr.m <- marray(CH.M) marr.f <- marray(CH.F) # Data analysis # Specify model in BUGS language sink("cjs-mnl-g.bug") cat(" model { # Priors and constraints for (t in 1:(n.occasions-1)){ phi.m[t] <- mean.phim phi.f[t] <- mean.phif p[t] <- mean.p } mean.phim ~ dunif(0, 1) # Prior for mean survival mean.phif ~ dunif(0, 1) # Prior for mean survival mean.p ~ dunif(0, 1) # Prior for mean recapture # Define the multinomial likelihood for (t in 1:(n.occasions-1)){ marr.m[t,1:n.occasions] ~ dmulti(pr.m[t, ], r.m[t]) marr.f[t,1:n.occasions] ~ dmulti(pr.f[t, ], r.f[t]) } # Calculate the number of birds released each year for (t in 1:(n.occasions-1)){ r.m[t] <- sum(marr.m[t, ]) r.f[t] <- sum(marr.f[t, ]) } # Define the cell probabilities of the m-array # Main diagonal for (t in 1:(n.occasions-1)){ q[t] <- 1-p[t] # Probability of non-recapture pr.m[t,t] <- phi.m[t]*p[t] pr.f[t,t] <- phi.f[t]*p[t] # Above main diagonal for (j in (t+1):(n.occasions-1)){ pr.m[t,j] <- prod(phi.m[t:j])*prod(q[t:(j-1)])*p[j] pr.f[t,j] <- prod(phi.f[t:j])*prod(q[t:(j-1)])*p[j] } #j # Below main diagonal for (j in 1:(t-1)){ pr.m[t,j] <- 0 pr.f[t,j] <- 0 } #j } #t # Last column: probability of non-recapture for (t in 1:(n.occasions-1)){ pr.m[t,n.occasions] <- 1-sum(pr.m[t,1:(n.occasions-1)]) pr.f[t,n.occasions] <- 1-sum(pr.f[t,1:(n.occasions-1)]) } #t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr.m = marr.m, marr.f = marr.f, n.occasions = dim(marr.m)[2]) # Initial values inits <- function(){list(mean.phim = runif(1, 0, 1), mean.phif = runif(1, 0, 1), mean.p = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.phim", "mean.phif", "mean.p") # MCMC settings niter <- 10000 nthin <- 3 nburn <- 5000 nchains <- 3 # Call WinBUGS from R (BRT 0.8 min) cjs <- bugs(bugs.data, inits, parameters, "cjs-mnl-g.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(cjs, 3) ############## # Exercise 3 ############## # Data simulation # Define the parameters n.occasions <- 10 # Number of recapture occasions marked.j <- rep(30, n.occasions-1) # Annual number of newly marked juveniles marked.a <- rep(20, n.occasions-1) # Annual number of newly marked adults phi.juvm <- 0.3 # Juvenile annual survival of males phi.juvf <- 0.2 # Juvenile annual survival of females phi.ad <- 0.65 # Adult annual survival p.m <- c(0.5, 0.6, 0.4, 0.4, 0.7, 0.5, 0.8, 0.3, 0.8) # Recapture males diff <- 0.3 p.f <- plogis(qlogis(p.m)+diff) phi.jm <- c(phi.juvm, rep(phi.ad, n.occasions-2)) phi.jf <- c(phi.juvf, rep(phi.ad, n.occasions-2)) phi.a <- rep(phi.ad, n.occasions-1) # Define matrices with the survival and recapture probabilities PHI.JM <- matrix(0, ncol = n.occasions-1, nrow = sum(marked.j)) for (i in 1:(length(marked.j)-1)){ PHI.JM[(sum(marked.j[1:i])-marked.j[i]+1):sum(marked.j[1:i]),i:(n.occasions-1)] <- matrix(rep(phi.jm[1:(n.occasions-i)],marked.j[i]), ncol = n.occasions-i, byrow = TRUE) } PHI.JF <- matrix(0, ncol = n.occasions-1, nrow = sum(marked.j)) for (i in 1:(length(marked.j)-1)){ PHI.JF[(sum(marked.j[1:i])-marked.j[i]+1):sum(marked.j[1:i]),i:(n.occasions-1)] <- matrix(rep(phi.jf[1:(n.occasions-i)],marked.j[i]), ncol = n.occasions-i, byrow = TRUE) } PHI.A <- matrix(rep(phi.a, sum(marked.a)), ncol = n.occasions-1, nrow = sum(marked.a), byrow = TRUE) P.JM <- matrix(rep(p.m, sum(marked.j)), ncol = n.occasions-1, nrow = sum(marked.j), byrow = TRUE) P.AM <- matrix(rep(p.m, sum(marked.a)), ncol = n.occasions-1, nrow = sum(marked.a), byrow = TRUE) P.JF <- matrix(rep(p.f, n.occasions*sum(marked.j)), ncol = n.occasions-1, nrow = sum(marked.j), byrow = TRUE) P.AF <- matrix(rep(p.f, n.occasions*sum(marked.a)), ncol = n.occasions-1, nrow = sum(marked.a), byrow = TRUE) # Apply simulation function CH.JM <- simul.cjs(PHI.JM, P.JM, marked.j) CH.AM <- simul.cjs(PHI.A, P.AM, marked.a) CH.JF <- simul.cjs(PHI.JF, P.JF, marked.j) CH.AF <- simul.cjs(PHI.A, P.AF, marked.a) # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) f.jm <- apply(CH.JM, 1, get.first) f.am <- apply(CH.AM, 1, get.first) f.jf <- apply(CH.JF, 1, get.first) f.af <- apply(CH.AF, 1, get.first) # Create matrices X indicating the age class x.jm <- matrix(NA, ncol = dim(CH.JM)[2]-1, nrow = dim(CH.JM)[1]) x.am <- matrix(NA, ncol = dim(CH.AM)[2]-1, nrow = dim(CH.AM)[1]) x.jf <- matrix(NA, ncol = dim(CH.JF)[2]-1, nrow = dim(CH.JF)[1]) x.af <- matrix(NA, ncol = dim(CH.AF)[2]-1, nrow = dim(CH.AF)[1]) for (i in 1:dim(CH.JM)[1]){ for (t in f.jm[i]:(dim(CH.JM)[2]-1)){ x.jm[i,t] <- 2 x.jm[i,f.jm[i]] <- 1 } #t } #i for (i in 1:dim(CH.AM)[1]){ for (t in f.am[i]:(dim(CH.AM)[2]-1)){ x.am[i,t] <- 2 } #t } #i for (i in 1:dim(CH.JF)[1]){ for (t in f.jf[i]:(dim(CH.JF)[2]-1)){ x.jf[i,t] <- 2 x.jf[i,f.jf[i]] <- 1 } #t } #i for (i in 1:dim(CH.AF)[1]){ for (t in f.af[i]:(dim(CH.AF)[2]-1)){ x.af[i,t] <- 2 } #t } #i # Combine the data, and create group variable CH <- rbind(CH.JM, CH.AM, CH.JF, CH.AF) f <- c(f.jm, f.am, f.jf, f.af) x <- rbind(x.jm, x.am, x.jf, x.af) group <- c(rep(1, dim(CH.JM)[1]), rep(1, dim(CH.AM)[1]), rep(2, dim(CH.JF)[1]), rep(2, dim(CH.AF)[1])) # Data analysis # Specify model in BUGS language sink("cjs-age2.bug") cat(" model { # Priors and constraints for (i in 1:nind){ for (t in f[i]:(n.occasions-1)){ phi[i,t] <- beta[x[i,t],group[i]] logit(p[i,t]) <- gamma[t] + delta[group[i]] } #t } #i beta[1,1] ~ dunif(0, 1) # Prior for juv survival of male beta[1,2] ~ dunif(0, 1) # Prior for juv survival of female beta[2,1] ~ dunif(0, 1) # Prior for ad survival of male beta[2,2] <- beta[2,1] # Ad survival of female identical to male for (t in 1:(n.occasions-1)){ gamma[t] ~ dnorm(0,0.001)I(-15,15) # Prior for recapture of males } delta[1] <- 0 delta[2] ~ dnorm(0, 0.001)I(-15,15) # Prior for recapture (diff. between sexes) # Back-transformations for (t in 1:(n.occasions-1)){ p.m[t] <- 1/(1+exp(-gamma[t])) p.f[t] <- 1/(1+exp(-gamma[t]-delta[2])) } # Likelihood for (i in 1:nind){ # Ensures that individuals enter the sample with probability 1 z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- phi[i,t-1] * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- p[i,t-1] * z[i,t] } # t } # i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = CH, f = f, nind = dim(CH)[1], n.occasions = dim(CH)[2], z = known.state.cjs(CH), x = x, group = group) # Initial values inits <- function(){list(z = cjs.init.z(CH, f), beta = matrix(c(runif(3, 0, 1), NA), ncol = 2, byrow = TRUE), gamma = rnorm(dim(CH)[2]-1), delta = c(NA, rnorm(1)))} # Parameters monitored parameters <- c("beta", "p.m", "p.f") # MCMC settings niter <- 5000 nthin <- 3 nburn <- 3000 nchains <- 3 # Call WinBUGS from R (BRT 5 min) cjs.age <- bugs(bugs.data, inits, parameters, "cjs-age2.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(cjs.age, 3) ############## # Exercise 4 ############## # Specify model in BUGS language sink("cjs-c-c.bug") cat(" model { # Priors and constraints for (i in 1:nind){ for (t in f[i]:(n.occasions-1)){ phi[i,t] <- mean.phi p[i,t] <- mean.p } #t } #i mean.phi ~ dunif(0, 1) # Prior for mean survival mean.p ~ dunif(0, 1) # Prior for mean recapture # Define the likelihood for (i in 1:nind){ # Ensures that individuals enter the sample with probability 1 z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- phi[i,t-1] * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- p[i,t-1] * z[i,t] } # t } # i } ",fill = TRUE) sink() # MCMC settings niter <- 2000 nthin <- 1 nburn <- 1000 nchains <- 1 # Define the simulation parameters and data structures to store the output nsim <- 500 phi.est <- p.est <- numeric() # Start loop for the simulation for (sim in 1:nsim){ # Define the parameters n.occasions <- 6 # Number of capture occasions marked <- rep(30, n.occasions-1) phi <- rep(0.65, n.occasions-1) p <- rep(0.4, n.occasions-1) # Define a matrix with the survival and recapture probabilities PHI <- matrix(rep(phi, sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE) P <- matrix(rep(p, sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE) # Apply simulation function CH <- simul.cjs(PHI, P, marked) # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) f <- apply(CH, 1, get.first) # Bundle data bugs.data <- list(y = CH, f = f, nind = dim(CH)[1], n.occasions = dim(CH)[2], z = known.state.cjs(CH)) # Initial values inits <- function(){list(z = cjs.init.z(CH, f), mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.phi", "mean.p") # Call WinBUGS from R out <- bugs(bugs.data, inits, parameters, "cjs-c-c.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = FALSE, bugs.directory = bugs.dir, working.directory = getwd()) # Store results phi.est[sim] <- out$mean$mean.phi p.est[sim] <- out$mean$mean.p } #sim par(mfrow = c(1, 2)) hist(phi.est, nclass = 15, col = "lightgrey", las = 1, ylab = "Frequency", xlab = expression(phi), main = "") abline(v = 0.65, col = "red", lwd = 2) hist(p.est, nclass = 15, col = "lightgrey", las = 1, ylab = "Frequency", xlab = "p", main = "") abline(v = 0.4, col = "red", lwd = 2) mean(phi.est)-phi[1] mean(p.est)-p[1] sd(phi.est) sd(p.est) # Define the parameters n.occasions <- 6 # Number of recapture occasions marked <- rep(10000, n.occasions-1) # Annual number of newly marked individuals phi <- rep(0.65, n.occasions-1) p <- rep(0.4, n.occasions-1) # Define a matrix with the survival and recapture probabilities PHI <- matrix(rep(phi, sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = T) P <- matrix(rep(p, sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = T) # Simulate capture-recapture data CH <- simul.cjs(PHI, P, marked) # Specify model in BUGS language sink("cjs-mnl.bug") cat(" model { # Priors and constraints for (t in 1:(n.occasions-1)){ phi[t] <- mean.phi p[t] <- mean.p } mean.phi ~ dunif(0, 1) # Prior for mean survival mean.p ~ dunif(0, 1) # Prior for mean recapture # Define the multinomial likelihood for (t in 1:(n.occasions-1)){ marr[t,1:n.occasions] ~ dmulti(pr[t, ], r[t]) } # Calculate the number of birds released each year for (t in 1:(n.occasions-1)){ r[t] <- sum(marr[t, ]) } # Define the cell probabilities of the m-array # Main diagonal for (t in 1:(n.occasions-1)){ q[t] <- 1-p[t] # Probability of non-recapture pr[t,t] <- phi[t]*p[t] # Above main diagonal for (j in (t+1):(n.occasions-1)){ pr[t,j] <- prod(phi[t:j])*prod(q[t:(j-1)])*p[j] } #j # Below main diagonal for (j in 1:(t-1)){ pr[t,j] <- 0 } #j } #t # Last column: probability of non-recapture for (t in 1:(n.occasions-1)){ pr[t,n.occasions] <- 1-sum(pr[t,1:(n.occasions-1)]) } #t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr = marray(CH), n.occasions = dim(CH)[2]) # Initial values inits <- function(){list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.phi", "mean.p") # MCMC settings niter <- 2000 nthin <- 3 nburn <- 1000 nchains <- 3 # Call WinBUGS from R (BRT 0.03 min) cjs.sim <- bugs(bugs.data, inits, parameters, "cjs-mnl.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = T, bugs.directory = bugs.dir, working.directory = getwd()) print(cjs.sim, 3) ############## # Exercise 5 ############## # Data simulation # Define the parameters n.occasions <- 10 # Number of recapture occasions marked <- rep(200, n.occasions-1) # Annual number of newly marked juv. phi.juv <- 0.3 # Juvenile annual survival phi.ad <- 0.65 # Adult annual survival p <- rep(0.5, n.occasions-1) # Recapture phi.j <- c(phi.juv, rep(phi.ad, n.occasions-2)) # Define matrices with the survival and recapture probabilities PHI.J <- matrix(0, ncol = n.occasions-1, nrow = sum(marked)) for (i in 1:(length(marked)-1)){ PHI.J[(sum(marked[1:i])-marked[i]+1):sum(marked[1:i]),i:(n.occasions-1)] <- matrix(rep(phi.j[1:(n.occasions-i)],marked[i]), ncol = n.occasions-i, byrow = TRUE) } P <- matrix(rep(p, n.occasions*sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE) # Apply simulation function CH <- simul.cjs(PHI.J, P, marked) # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) f <- apply(CH, 1, get.first) # Create matrix X indicating the age class x <- matrix(NA, ncol = dim(CH)[2]-1, nrow = dim(CH)[1]) for (i in 1:dim(CH)[1]){ for (t in f[i]:(dim(CH)[2]-1)){ x[i,t] <- t-f[i]+1 } #t } #i # Data analysis # Specify model in BUGS language sink("cjs-age2.bug") cat(" model { # Priors and constraints for (i in 1:nind){ for (t in f[i]:(n.occasions-1)){ logit(phi[i,t]) <- beta[x[i,t]] p[i,t] <- mean.p } #t } #i beta[1] ~ dnorm(0, 0.001) # Prior for first year survival for (u in 2:(n.occasions-1)){ beta[u] <- mu + gamma*(u-1) # Linear model for age > 1y } mu ~ dnorm(0, 0.001)I(-15,15) # Prior for intercept of linear model gamma ~ dnorm(0, 0.001)I(-15,15) # Prior for slope of linear model mean.p ~ dunif(0, 1) # Prior for mean recapture # Back-transformations for (t in 1:(n.occasions-1)){ phi.a[t] <- 1/(1+exp(-beta[t])) } # Define the likelihood for (i in 1:nind){ # Ensures that individuals enter the sample with probability 1 z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- phi[i,t-1] * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- p[i,t-1] * z[i,t] } # t } # i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = CH, f = f, nind = dim(CH)[1], n.occasions = dim(CH)[2], z = known.state.cjs(CH), x = x) # Initial values inits <- function(){list(z = cjs.init.z(CH, f), beta = c(rnorm(1), rep(NA, n.occasions-2)), mu = rnorm(1), gamma = rnorm(1), mean.p = runif(1, 0, 1))} # Parameters monitored parameters <- c("phi.a", "mu", "gamma", "mean.p") # MCMC settings niter <- 2000 nthin <- 3 nburn <- 1000 nchains <- 3 # Call WinBUGS from R (BRT 6 min) cjs.age <- bugs(bugs.data, inits, parameters, "cjs-age2.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(cjs.age, 3) ############## # Exercise 6 ############## # Data simulation # Define the parameters n.occasions <- 15 # Number of recapture occasions marked <- rep(100, n.occasions-1) # Number of newly marked juveniles phi.juv <- 0.4 # Juvenile annual survival var.phi <- 0.3 phi.ad <- 0.8 # Adult annual survival p <- rep(0.6, n.occasions-1) # Recapture phi.juv.t <- plogis(rnorm(n.occasions-1, qlogis(phi.juv), var.phi^0.5)) phi.j <- matrix(0, nrow = n.occasions-1, ncol = n.occasions-1) for (t in 1 :(n.occasions-1)){ phi.j[t,t:(n.occasions-1)] <- c(phi.juv.t[t], rep(phi.ad, n.occasions-1-t)) } # Define matrices with the survival and recapture probabilities PHI.J <- matrix(0, ncol = n.occasions-1, nrow = sum(marked)) for (i in 1:length(marked)){ PHI.J[(sum(marked[1:i])-marked[i]+1):sum(marked[1:i]),i:(n.occasions-1)] <- matrix(rep(phi.j[i,i:(n.occasions-1)],marked[i]), ncol = n.occasions-i, byrow = TRUE) } P <- matrix(rep(p, n.occasions*sum(marked)), ncol = n.occasions-1, nrow = sum(marked), byrow = TRUE) # Apply simulation function CH <- simul.cjs(PHI.J, P, marked) # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) f <- apply(CH, 1, get.first) # Create matrices X indicating the age class x <- matrix(NA, ncol = dim(CH)[2]-1, nrow = dim(CH)[1]) for (i in 1:dim(CH)[1]){ for (t in f[i]:(dim(CH)[2]-1)){ x[i,t] <- 2 x[i,f[i]] <- 1 } #t } #i # Data analysis # a) State-space likelihood # Specify model in BUGS language sink("cjs-age2.bug") cat(" model { # Priors and constraints for (i in 1:nind){ for (t in f[i]:(n.occasions-1)){ logit(phi[i,t]) <- beta[x[i,t],t] p[i,t] <- mean.p } #t } #i for (t in 1:(n.occasions-1)){ beta[1,t] <- mu + epsilon[t] # Juvenile survival epsilon[t] ~ dnorm(0, tau) # Temporal variation of juv survival phi.j[t] <- 1/(1+exp(-beta[1,t])) beta[2,t] <- lphi.ad # Constrain ad survival to be constant } mu <- log(mean.phij / (1-mean.phij)) mean.phij ~ dunif(0, 1) # Prior for mean juv survival sigma ~ dunif(0, 10) # Prior on sd of temp. var tau <- pow(sigma, -2) sigma2 <- pow(sigma, 2) lphi.ad <- log(mean.phiad / (1-mean.phiad)) mean.phiad ~ dunif(0, 1) # Prior for mean ad survival mean.p ~ dunif(0, 1) # Prior for mean recapture # Define the likelihood for (i in 1:nind){ # Ensures that individuals enter the sample with probability 1 z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- phi[i,t-1] * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- p[i,t-1] * z[i,t] } # t } # i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = CH, f = f, nind = dim(CH)[1], n.occasions = dim(CH)[2], z = known.state.cjs(CH), x = x) # Initial values inits <- function(){list(z = cjs.init.z(CH, f), mean.phij = runif(1, 0, 1), mean.phiad = runif(1, 0, 1), sigma = runif(1, 0, 5), mean.p = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.phij", "phi.j", "sigma2", "mean.phiad", "mean.p") # MCMC settings niter <- 2000 nthin <- 3 nburn <- 1000 nchains <- 3 # Call WinBUGS from R (BRT 9 min) cjs.1 <- bugs(bugs.data, inits, parameters, "cjs-age2.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(cjs.1, digits = 3) # b) Multinomial likelihood # Create m-arrays cap <- apply(CH, 1, sum) ind <- which(cap >= 2) CH.R <- CH[ind,] # Juvenile CH recaptured at least once CH.N <- CH[-ind,] # Juvenile CH never recaptured # Remove first capture first <- numeric() for (i in 1:dim(CH.R)[1]){ first[i] <- min(which(CH.R[i,]==1)) } CH.R1 <- CH.R for (i in 1:dim(CH.R)[1]){ CH.R1[i,first[i]] <- 0 } # Create m-array of those recaptured at least once CH.A.marray <- marray(CH.R1) # Create CH matrix for juveniles, ignoring subsequent recaptures second <- numeric() for (i in 1:dim(CH.R1)[1]){ second[i] <- min(which(CH.R1[i,]==1)) } CH.R2 <- matrix(0, nrow = dim(CH.R)[1], ncol = dim(CH.R)[2]) for (i in 1:dim(CH.R)[1]){ CH.R2[i,first[i]] <- 1 CH.R2[i,second[i]] <- 1 } # Create m-array for these CH.R.marray <- marray(CH.R2) # The last column ought to show the number of juveniles not recaptured again and should all be zeros, since all of them are released as adults CH.R.marray[,dim(CH)[2]] <- 0 # Create the m-array for juveniles never recaptured and add it to the previous m-array CH.N.marray <- marray(CH.N) CH.J.marray <- CH.R.marray + CH.N.marray # Specify model in BUGS language sink("cjs-mnl-2age.bug") cat(" model { # Priors and constraints for (t in 1:(n.occasions-1)){ logit(phi.juv[t]) <- mu + epsilon[t] epsilon[t] ~ dnorm(0, tau)I(-15,15) # Range restriction phi.j[t] <- 1/(1+exp(-mu-epsilon[t])) phi.ad[t] <- mean.phiad p[t] <- mean.p } mu <- log(mean.phij / (1-mean.phij)) mean.phij ~ dunif(0, 1) # Prior for mean juv survival sigma ~ dunif(0, 10) # Prior on sd of temp. var tau <- pow(sigma, -2) sigma2 <- pow(sigma, 2) mean.phiad ~ dunif(0, 1) # Prior for mean ad survival mean.p ~ dunif(0, 1) # Prior for mean recapture # Define the multinomial likelihood for (t in 1:(n.occasions-1)){ marrj[t,1:n.occasions] ~ dmulti(prj[t,], rj[t]) marra[t,1:n.occasions] ~ dmulti(pra[t,], ra[t]) } # Calculate the number of birds released each year for (t in 1:(n.occasions-1)){ rj[t] <- sum(marrj[t,]) ra[t] <- sum(marra[t,]) } # Define the cell probabilities of the m-arrays # Main diagonal for (t in 1:(n.occasions-1)){ q[t] <- 1-p[t] # Probability of non-recapture prj[t,t] <- phi.juv[t]*p[t] pra[t,t] <- phi.ad[t]*p[t] # Above main diagonal for (j in (t+1):(n.occasions-1)){ prj[t,j] <- phi.juv[t]*prod(phi.ad[(t+1):j])*prod(q[t:(j-1)])*p[j] pra[t,j] <- prod(phi.ad[t:j])*prod(q[t:(j-1)])*p[j] } # j # Below main diagonal for (j in 1:(t-1)){ prj[t,j] <- 0 pra[t,j] <- 0 } # j } # t # Last column: probability of non-recapture for (t in 1:(n.occasions-1)){ prj[t,n.occasions] <- 1-sum(prj[t,1:(n.occasions-1)]) pra[t,n.occasions] <- 1-sum(pra[t,1:(n.occasions-1)]) } # t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marrj = CH.J.marray, marra = CH.A.marray, n.occasions = dim(CH)[2]) # Initial values inits <- function(){list(mean.phij = runif(1, 0, 1), mean.phiad = runif(1, 0, 1), sigma = runif(1, 0, 5), mean.p = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.phij", "phi.j", "sigma2", "mean.phiad", "mean.p") # MCMC settings niter <- 5000 nthin <- 6 nburn <- 2500 nchains <- 3 # Call WinBUGS from R (BRT 1 min) cjs.2 <- bugs(bugs.data, inits, parameters, "cjs-mnl-2age.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = T, bugs.directory = bugs.dir, working.directory = getwd()) print(cjs.2, 3) ############################################################################# # # Chapter 8 # ############################################################################# ############## # Exercise 1 ############## # Data simulation # Define the parameters n.occasions <- 10 # Number of occasions marked <- rep(50, n.occasions) # Annual number of newly marked individuals s <- rep(0.5, n.occasions) r1 <- rep(0.1, n.occasions) r2 <- rep(0.2, n.occasions) # Define matrices with the survival and recovery probabilities S <- matrix(rep(s, sum(marked)), ncol = n.occasions, nrow = sum(marked), byrow = TRUE) R1 <- matrix(rep(r1, sum(marked)), ncol = n.occasions, nrow = sum(marked), byrow = TRUE) R2 <- matrix(rep(r2, sum(marked)), ncol = n.occasions, nrow = sum(marked), byrow = TRUE) # Apply function MR1 <- simul.mr(S, R1, marked) MR2 <- simul.mr(S, R2, marked) # Merge capture-histories MR <- rbind(MR1, MR2) # Create group variable group <- c(rep(1, dim(MR1)[1]), rep(2, dim(MR2)[1])) # Create vector with occasion of marking get.first <- function(x) min(which(x!=0)) f <- apply(MR, 1, get.first) # Create m-arrays marr1 <- marray.dead(MR1) marr2 <- marray.dead(MR2) # Data analyses # a) Multinomial likelihood # Specify model in BUGS language sink("mr-mnl.bug") cat(" model { # Priors and constraints for (t in 1:n.occasions){ s[t] <- mean.s r1[t] <- mean.r1 r2[t] <- mean.r2 } mean.s ~ dunif(0, 1) mean.r1 ~ dunif(0, 1) mean.r2 ~ dunif(0, 1) # Define the multinomial likelihoods for (t in 1:n.occasions){ marr1[t,1:(n.occasions+1)] ~ dmulti(pr1[t,], rel1[t]) marr2[t,1:(n.occasions+1)] ~ dmulti(pr2[t,], rel2[t]) } # Calculate the number of birds released each year for (t in 1:n.occasions){ rel1[t] <- sum(marr1[t,]) rel2[t] <- sum(marr2[t,]) } # Define the cell probabilities of the m-array # Main diagonal for (t in 1:n.occasions){ pr1[t,t] <- (1-s[t])*r1[t] pr2[t,t] <- (1-s[t])*r2[t] # Above main diagonal for (j in (t+1):n.occasions){ pr1[t,j] <- prod(s[t:(j-1)])*(1-s[j])*r1[j] pr2[t,j] <- prod(s[t:(j-1)])*(1-s[j])*r2[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr1[t,j] <- 0 pr2[t,j] <- 0 } # j } # t # Last column: probability of non-recovery for (t in 1:n.occasions){ pr1[t,n.occasions+1] <- 1-sum(pr1[t,1:n.occasions]) pr2[t,n.occasions+1] <- 1-sum(pr2[t,1:n.occasions]) } # t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr1 = marr1, marr2 = marr2, n.occasions = dim(marr1)[2]-1) # Initial values inits <- function(){list(mean.s = runif(1, 0, 1), mean.r1 = runif(1, 0, 1), mean.r2 = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.s", "mean.r1", "mean.r2") # MCMC settings niter <- 5000 nthin <- 6 nburn <- 2000 nchains <- 3 # Call WinBUGS from R (BRT 0.1 min) mr.age <- bugs(bugs.data, inits, parameters, "mr-mnl.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir) # Inspect results print(mr.age, 3) # b) State-space likelihood # Specify model in BUGS language sink("mr.ss.bug") cat(" model { # Priors and constraints for (i in 1:nind){ for (t in 1:n.occasions){ s[i,t] <- mean.s r[i,t] <- mean.r[group[i]] } #t } #i mean.s ~ dunif(0, 1) for (u in 1:g){ mean.r[u] ~ dunif(0, 1) } # Likelihood for (i in 1:nind){ # Define latent state at first capture z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- s[i,t-1] * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- r[i,t-1] * (z[i,t-1] - z[i,t]) } #t } #i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = MR, f = f, group = group, g = length(unique(group)), nind = dim(MR)[1], n.occasions = dim(MR)[2], z = known.state.mr(MR)) # Initial values inits <- function(){list(z = mr.init.z(MR), mean.s = runif(1, 0, 1), mean.r = runif(2, 0, 1))} # Parameters monitored parameters <- c("mean.s", "mean.r") # MCMC settings niter <- 7000 nthin <- 3 nburn <- 5000 nchains <- 3 # Call WinBUGS from R (BRT 11 min) mr <- bugs(bugs.data, inits, parameters, "mr.ss.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, working.directory = getwd(), bugs.directory = bugs.dir) # Inspect results print(mr, digits = 3) ############## # Exercise 2 ############## # Data simulation n.occasions <- 15 # Number of occasions marked.j <- rep(200, n.occasions) # Annual number of newly marked young sjuv <- 0.3 # First year survival probability sad <- 0.7 # Adult survival probability rjuv <- 0.25 # First year recovery probability rad <- 0.15 # Adult recovery probability sj <- c(sjuv, rep(sad, n.occasions-1)) rj <- c(rjuv, rep(rad, n.occasions-1)) sa <- rep(sad, n.occasions) ra <- rep(rad, n.occasions) # Define matrices with the survival and recovery probabilities S <- matrix(0, ncol = n.occasions, nrow = sum(marked.j)) for (i in 1:length(marked.j)){ S[(sum(marked.j[1:i])-marked.j[i]+1):sum(marked.j[1:i]),i:n.occasions] <- matrix(rep(sj[1:(n.occasions-i+1)], marked.j[i]), ncol = n.occasions-i+1, byrow = TRUE) } R <- matrix(0, ncol = n.occasions, nrow = sum(marked.j)) for (i in 1:length(marked.j)){ R[(sum(marked.j[1:i])-marked.j[i]+1):sum(marked.j[1:i]),i:n.occasions] <- matrix(rep(rj[1:(n.occasions-i+1)], marked.j[i]), ncol = n.occasions-i+1, byrow = TRUE) } # Apply simulation function MR <- simul.mr(S, R, marked.j) # Create m-arrays marr <- marray.dead(MR) # Data analysis # a) Using the data generating model {sa2, ra2} # Specify model in BUGS language sink("mr-mnl-age1.bug") cat(" model { # Priors and constraints for (t in 1:n.occasions){ sj[t] <- mean.sj sa[t] <- mean.sa rj[t] <- mean.rj ra[t] <- mean.ra } mean.sj ~ dunif(0, 1) mean.sa ~ dunif(0, 1) mean.rj ~ dunif(0, 1) mean.ra ~ dunif(0, 1) # Define the multinomial likelihoods for (t in 1:n.occasions){ marr[t,1:(n.occasions+1)] ~ dmulti(pr[t,], rel[t]) } # Calculate the number of birds released each year for (t in 1:n.occasions){ rel[t] <- sum(marr[t,]) } # Define the cell probabilities of the m-array # Main diagonal for (t in 1:n.occasions){ pr[t,t] <- (1-sj[t])*rj[t] # Further above main diagonal for (j in (t+2):n.occasions){ pr[t,j] <- sj[t]*prod(sa[(t+1):(j-1)])*(1-sa[j])*ra[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr[t,j] <- 0 } # j } # t for (t in 1:(n.occasions-1)){ # One above main diagonal pr[t,t+1] <- sj[t]*(1-sa[t+1])*ra[t+1] } # t # Last column: probability of non-recovery for (t in 1:n.occasions){ pr[t,n.occasions+1] <- 1-sum(pr[t,1:n.occasions]) } # t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr = marr, n.occasions = dim(marr)[2]-1) # Initial values inits <- function(){list(mean.sj = runif(1, 0, 1), mean.sa = runif(1, 0, 1), mean.rj = runif(1, 0, 1), mean.ra = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.sj", "mean.sa", "mean.rj", "mean.ra") # MCMC settings niter <- 20000 nthin <- 6 nburn <- 10000 nchains <- 3 # Call WinBUGS from R (BRT 2 min) mr.age1 <- bugs(bugs.data, inits, parameters, "mr-mnl-age1.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, working.directory = getwd(), bugs.directory = bugs.dir) # Inspect results print(mr.age1, 3) # Plot posterior distribution of the parameters par(mfrow = c(2, 2)) plot(density(mr.age1$sims.list$mean.sj), xlab = "Juvenile survival", main = "") segments(0, 1, 1, 1, lty = 2) plot(density(mr.age1$sims.list$mean.sa), xlab = "Adult survival", main = "") segments(0.6, 1, 0.85, 1, lty = 2) plot(density(mr.age1$sims.list$mean.rj), xlab = "Juvenile recovery", main = "") segments(0, 1, 1, 1, lty = 2) plot(density(mr.age1$sims.list$mean.ra), xlab = "Adult recovery", main = "") segments(0, 1, 1, 1, lty = 2) # b) Using the model {sa2, r} # Specify model in BUGS language sink("mr-mnl-age2.bug") cat(" model { # Priors and constraints for (t in 1:n.occasions){ sj[t] <- mean.sj sa[t] <- mean.sa r[t] <- mean.r } mean.sj ~ dunif(0, 1) mean.sa ~ dunif(0, 1) mean.r ~ dunif(0, 1) # Define the multinomial likelihoods for (t in 1:n.occasions){ marr[t,1:(n.occasions+1)] ~ dmulti(pr[t,], rel[t]) } # Calculate the number of birds released each year for (t in 1:n.occasions){ rel[t] <- sum(marr[t,]) } # Define the cell probabilities of the m-array # Main diagonal for (t in 1:n.occasions){ pr[t,t] <- (1-sj[t])*r[t] # Further above main diagonal for (j in (t+2):n.occasions){ pr[t,j] <- sj[t]*prod(sa[(t+1):(j-1)])*(1-sa[j])*r[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr[t,j] <- 0 } # j } # t for (t in 1:(n.occasions-1)){ # One above main diagonal pr[t,t+1] <- sj[t]*(1-sa[t+1])*r[t+1] } # t # Last column: probability of non-recovery for (t in 1:n.occasions){ pr[t,n.occasions+1] <- 1-sum(pr[t,1:n.occasions]) } # t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr = marr, n.occasions = dim(marr)[2]-1) # Initial values inits <- function(){list(mean.sj = runif(1, 0, 1), mean.sa = runif(1, 0, 1), mean.r = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.sj", "mean.sa", "mean.r") # MCMC settings niter <- 20000 nthin <- 6 nburn <- 10000 nchains <- 3 # Call WinBUGS from R (BRT 2 min) mr.age2 <- bugs(bugs.data, inits, parameters, "mr-mnl-age2.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, working.directory = getwd(), bugs.directory = bugs.dir) # Inspect results print(mr.age2, 3) # Plot posterior distribution of the parameters par(mfrow = c(2, 2)) plot(density(mr.age2$sims.list$mean.sj), xlab = "Juvenile survival", main = "") segments(0.15, 1, 0.32, 1, lty = 2) plot(density(mr.age2$sims.list$mean.sa), xlab = "Adult survival", main = "") segments(0.5, 1, 0.86, 1, lty = 2) plot(density(mr.age2$sims.list$mean.r), xlab = "Recovery", main = "") segments(0.1, 1, 0.27, 1, lty = 2) ############## # Exercise 3 ############## # Data simulation n.occasions <- 20 # Number of occasions marked <- rep(500, n.occasions) # Annual number of newly marked young s <- seq(0.8, 0.6, length.out = n.occasions) # Survival probability r <- rep(0.05, n.occasions) # Recovery probability # Define matrices with the survival and recovery probabilities S <- matrix(rep(s, sum(marked)), ncol = n.occasions, nrow = sum(marked), byrow = TRUE) R <- matrix(rep(r, sum(marked)), ncol = n.occasions, nrow = sum(marked), byrow = TRUE) # Apply simulation function MR <- simul.mr(S, R, marked) # Create m-arrays marr <- marray.dead(MR) # Data analysis # Specify model in BUGS language sink("mr-trend.bug") cat(" model { # Priors and constraints for (t in 1:n.occasions){ logit(s[t]) <- mu + slope*(t-n.occasions/2) # standardise trend r[t] <- mean.r } mu ~ dnorm(0, 0.001) slope ~ dnorm(0, 0.001) mean.r ~ dunif(0, 1) # Define the multinomial likelihoods for (t in 1:n.occasions){ marr[t,1:(n.occasions+1)] ~ dmulti(pr[t,], rel[t]) } # Calculate the number of birds released each year for (t in 1:n.occasions){ rel[t] <- sum(marr[t,]) } # Define the cell probabilities of the m-array # Main diagonal for (t in 1:n.occasions){ pr[t,t] <- (1-s[t])*r[t] # Above main diagonal for (j in (t+1):n.occasions){ pr[t,j] <- prod(s[t:(j-1)])*(1-s[j])*r[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr[t,j] <- 0 } # j } # t # Last column: probability of non-recovery for (t in 1:n.occasions){ pr[t,n.occasions+1] <- 1-sum(pr[t,1:n.occasions]) } # t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr = marr, n.occasions = dim(marr)[2]-1) # Initial values inits <- function(){list(mu = rnorm(1), slope = rnorm(1), mean.r = runif(1, 0, 1))} # Parameters monitored parameters <- c("s", "mu", "slope", "mean.r") # MCMC settings niter <- 10000 nthin <- 6 nburn <- 5000 nchains <- 3 # Call WinBUGS from R (BRT 2 min) mr.trend <- bugs(bugs.data, inits, parameters, "mr-trend.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, working.directory = getwd(), bugs.directory = bugs.dir) # Inspect results print(mr.trend, digits=3) ############## # Exercise 4 ############## # Data simulation n.occasions <- 10 # Number of occasions marked <- rep(100, n.occasions) # Annual number of newly marked young s <- rep(0.7, n.occasions) # Survival probability mean.r <- 0.2 v.ind <- 0.7 r <- plogis(rnorm(sum(marked), qlogis(mean.r), v.ind^0.5)) # Define matrices with the survival and recovery probabilities S <- matrix(rep(s, sum(marked)), ncol = n.occasions, nrow = sum(marked), byrow = TRUE) R <- matrix(rep(r, n.occasions), ncol = n.occasions, nrow = sum(marked), byrow = FALSE) # Apply simulation function MR <- simul.mr(S, R, marked) # Compute vector with occasion of first capture get.first <- function(x) min(which(x!=0)) f <- apply(MR, 1, get.first) # Data analysis # a) Model with individual variation in recovery probability # Specify model in BUGS language sink("mr.indvar.bug") cat(" model { # Priors and constraints for (i in 1:nind){ for (t in 1:n.occasions){ s[i,t] <- mean.s logit(r[i,t]) <- mu + epsilon[i] } #t epsilon[i] ~ dnorm(0, tau)I(-15,15) } #i mean.s ~ dunif(0, 1) mu <- log(mean.r / (1-mean.r)) # Logit transformation mean.r ~ dunif(0, 1) tau <- pow(sigma, -2) sigma ~ dunif(0, 5) # Prior on standard deviation sigma2 <- pow(sigma, 2) # Likelihood for (i in 1:nind){ # Define latent state at first capture z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- s[i,t-1] * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- r[i,t-1] * (z[i,t-1] - z[i,t]) } #t } #i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = MR, f = f, nind = dim(MR)[1], n.occasions = dim(MR)[2], z = known.state.mr(MR)) # Initial values inits <- function(){list(z = mr.init.z(MR), mean.s = runif(1, 0, 1), mean.r = runif(1, 0, 1), sigma = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.s", "mean.r", "sigma2") # MCMC settings niter <- 50000 nthin <- 6 nburn <- 30000 nchains <- 3 # Call WinBUGS from R (BRT 115 min) mr.indvar <- bugs(bugs.data, inits, parameters, "mr.indvar.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, working.directory = getwd(), bugs.directory = bugs.dir) # Inspect results print(mr.indvar, digits = 3) # b) Model without individual variation in recovery probability # Specify model in BUGS language sink("mr.bug") cat(" model { # Priors and constraints for (i in 1:nind){ for (t in 1:n.occasions){ s[i,t] <- mean.s r[i,t] <- mean.r } #t } #i mean.s ~ dunif(0, 1) mean.r ~ dunif(0, 1) # Likelihood for (i in 1:nind){ # Define latent state at first capture z[i,f[i]] <- 1 for (t in (f[i]+1):n.occasions){ # State process z[i,t] ~ dbern(mu1[i,t]) mu1[i,t] <- s[i,t-1] * z[i,t-1] # Observation process y[i,t] ~ dbern(mu2[i,t]) mu2[i,t] <- r[i,t-1] * (z[i,t-1] - z[i,t]) } #t } #i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = MR, f = f, nind = dim(MR)[1], n.occasions = dim(MR)[2], z = known.state.mr(MR)) # Initial values inits <- function(){list(z = mr.init.z(MR), mean.s = runif(1, 0, 1), mean.r = runif(1, 0, 1))} # Parameters monitored parameters <- c("mean.s", "mean.r") # MCMC settings niter <- 20000 nthin <- 3 nburn <- 10000 nchains <- 3 # Call WinBUGS from R (BRT 25 min) mr <- bugs(bugs.data, inits, parameters, "mr.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, working.directory = getwd(), bugs.directory = bugs.dir) # Inspect results print(mr, digits = 3) ############################################################################# # # Chapter 9 # ############################################################################# ############## # Exercise 1 ############## # Data simulation # Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals phiAm <- 0.5 phiBm <- 0.6 pAm <- 0.3 pBm <- 0.7 phiAf <- 0.7 phiBf <- 0.6 pAf <- 0.4 pBf <- 0.8 psiAB <- 0.2 psiBA <- 0.5 n.occasions <- 6 n.states <- 3 n.obs <- 3 marked <- matrix(NA, ncol = n.states, nrow = n.occasions) marked[,1] <- rep(20, n.occasions) marked[,2] <- rep(20, n.occasions) marked[,3] <- rep(0, n.occasions) # Simulate male data # Define arrays with survival, transition and recapture probabilities # These are 4-dimensional arrays, with # Dimension 1: state of departure # Dimension 2: state of arrival # Dimension 3: individual # Dimension 4: time # 1. State process array totrel <- sum(marked) PSI.STATE <- array(NA, dim = c(n.states, n.states, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.STATE[,,i,t] <- matrix(c( phiAm*(1-psiAB), phiAm*psiAB, 1-phiAm, phiBm*psiBA, phiBm*(1-psiBA), 1-phiBm, 0, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # 2. Observation array PSI.OBS <- array(NA, dim = c(n.states, n.obs, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.OBS[,,i,t] <- matrix(c( pAm, 0, 1-pAm, 0, pBm, 1-pBm, 0, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # Execute simulation function sim <- simul.ms(PSI.STATE, PSI.OBS, marked) CHm <- sim$CH # Simulate female data # Define arrays with survival, transition and recapture probabilities # 1. State process array totrel <- sum(marked) PSI.STATE <- array(NA, dim = c(n.states, n.states, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.STATE[,,i,t] <- matrix(c( phiAf*(1-psiAB), phiAf*psiAB, 1-phiAf, phiBf*psiBA, phiBf*(1-psiBA), 1-phiBf, 0, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # 2. Observation array PSI.OBS <- array(NA, dim = c(n.states, n.obs, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.OBS[,,i,t] <- matrix(c( pAf, 0, 1-pAf, 0, pBf, 1-pBf, 0, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # Execute simulation function sim <- simul.ms(PSI.STATE, PSI.OBS, marked) CHf <- sim$CH # Compute vector with occasion of first capture get.first <- function(x) min(which(x!=0)) fm <- apply(CHm, 1, get.first) ff <- apply(CHf, 1, get.first) # Recode CH matrix: note, a 0 is not allowed! # 1 = seen alive in A, 2 = seen alive in B, 3 = not seen rCHm <- CHm # recoded CH rCHm[rCHm==0] <- 3 rCHf <- CHf # recoded CH rCHf[rCHf==0] <- 3 # Combine data sets rCH <- rbind(rCHm, rCHf) # Create group variable group <- c(rep(1, nrow(rCHm)), rep(2, nrow(rCHf))) # Compute vector with occasion of first capture get.first <- function(x) min(which(x!=3)) f <- apply(rCH, 1, get.first) # Data analysis # Specify model in BUGS language sink("ms.bug") cat(" model { #################################################### # Parameters: # phiA: survival probability at site A # phiB: survival probability at site B # psiAB: movement probability from site A to site B # psiBA: movement probability from site B to site A # pA: recapture probability at site A # pB: recapture probability at site B ################################################### # States (S) # 1 alive at A # 2 alive at B # 3 dead # Observations (O) # 1 seen at A # 2 seen at B # 3 not seen ################################################### # Priors and constraints for (i in 1: nind){ for (t in 1:(n.occasions-1)){ phiA[i,t] <- mean.phiA[group[i]] phiB[i,t] <- mean.phiB[group[i]] psiAB[i,t] <- mean.psi[1] psiBA[i,t] <- mean.psi[2] pA[i,t] <- mean.pA[group[i]] pB[i,t] <- mean.pB[group[i]] } #t } #i for (u in 1:2){ mean.phiA[u] ~ dunif(0, 1) # Priors for mean state-spec. survival (at A) mean.phiB[u] ~ dunif(0, 1) # Priors for mean state-spec. survival (at B) mean.psi[u] ~ dunif(0, 1) # Priors for mean transitions mean.pA[u] ~ dunif(0, 1) # Priors for mean state-spec. recapture (at A) mean.pB[u] ~ dunif(0, 1) # Priors for mean state-spec. recapture (at B) } # Define state and observation matrices for (i in 1:nind){ # Define probabilities of state S(t+1) given S(t) for (t in f[i]:(n.occasions-1)){ ps[1,i,t,1] <- phiA[i,t] * (1-psiAB[i,t]) ps[1,i,t,2] <- phiA[i,t] * psiAB[i,t] ps[1,i,t,3] <- 1-phiA[i,t] ps[2,i,t,1] <- phiB[i,t] * psiBA[i,t] ps[2,i,t,2] <- phiB[i,t] * (1-psiBA[i,t]) ps[2,i,t,3] <- 1-phiB[i,t] ps[3,i,t,1] <- 0 ps[3,i,t,2] <- 0 ps[3,i,t,3] <- 1 # Define probabilities of O(t) given S(t) po[1,i,t,1] <- pA[i,t] po[1,i,t,2] <- 0 po[1,i,t,3] <- 1-pA[i,t] po[2,i,t,1] <- 0 po[2,i,t,2] <- pB[i,t] po[2,i,t,3] <- 1-pB[i,t] po[3,i,t,1] <- 0 po[3,i,t,2] <- 0 po[3,i,t,3] <- 1 } #t } #i # Likelihood for (i in 1:nind){ # Define latent state at first capture z[i,f[i]] <- y[i,f[i]] for (t in (f[i]+1):n.occasions){ # State process: draw S(t) given S(t-1) z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,]) # Observation process: draw O(t) given S(t) y[i,t] ~ dcat(po[z[i,t], i, t-1,]) } #t } #i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = rCH, group = group, f = f, n.occasions = dim(rCH)[2], nind = dim(rCH)[1], z = known.state.ms(rCH, 3)) # Initial values inits <- function(){list(mean.phiA = runif(2, 0, 1), mean.phiB = runif(2, 0, 1), mean.psi = runif(2, 0, 1), mean.pA = runif(2, 0, 1), mean.pB = runif(2, 0, 1), z = ms.init.z(rCH, f))} # Parameters monitored parameters <- c("mean.phiA", "mean.phiB", "mean.psi", "mean.pA", "mean.pB") # MCMC settings ni <- 10000 nt <- 1 nb <- 5000 nc <- 3 # Call WinBUGS from R (BRT 7 min) ms <- bugs(bugs.data, inits, parameters, "ms.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(ms, 3) ############## # Exercise 2 ############## # Data simulation # Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals phiA <- c(0.5, 0.6, 0.3, 0.7, 0.5, 0.65, 0.55) phiB <- 0.6 pA <- 0.3 pB <- 0.7 psiAB <- 0.2 psiBA <- 0.5 n.occasions <- 8 n.states <- 3 n.obs <- 3 marked <- matrix(NA, ncol = n.states, nrow = n.occasions) marked[,1] <- rep(20, n.occasions) marked[,2] <- rep(20, n.occasions) marked[,3] <- rep(0, n.occasions) # Define arrays with survival, transition and recapture probabilities # These are 4-dimensional arrays, with # Dimension 1: state of departure # Dimension 2: state of arrival # Dimension 3: individual # Dimension 4: time # 1. State process array totrel <- sum(marked) PSI.STATE <- array(NA, dim = c(n.states, n.states, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.STATE[,,i,t] <- matrix(c( phiA[t]*(1-psiAB), phiA[t]*psiAB, 1-phiA[t], phiB*psiBA, phiB*(1-psiBA), 1-phiB, 0, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # 2. Observation array PSI.OBS <- array(NA, dim = c(n.states, n.obs, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.OBS[,,i,t] <- matrix(c( pA, 0, 1-pA, 0, pB, 1-pB, 0, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # Execute simulation function sim <- simul.ms(PSI.STATE, PSI.OBS, marked) CH <- sim$CH # Compute vector with occasion of first capture get.first <- function(x) min(which(x!=0)) f <- apply(CH, 1, get.first) # Recode CH matrix: note, a 0 is not allowed! # 1 = seen alive in A, 2 = seen alive in B, 3 = not seen rCH <- CH # recoded CH rCH[rCH==0] <- 3 # Data analysis # a) Fixed time effects # Specify model in BUGS language sink("ms.bug") cat(" model { #################################################### # Parameters: # phiA: survival probability at site A # phiB: survival probability at site B # psiAB: movement probability from site A to site B # psiBA: movement probability from site B to site A # pA: recapture probability at site A # pB: recapture probability at site B ################################################### # States (S) # 1 alive at A # 2 alive at B # 3 dead # Observations (O) # 1 seen at A # 2 seen at B # 3 not seen ################################################### # Priors and constraints for (t in 1:(n.occasions-1)){ phiA[t] ~ dunif(0, 1) phiB[t] <- mean.phiB psiAB[t] <- mean.psi[1] psiBA[t] <- mean.psi[2] pA[t] <- mean.p[1] pB[t] <- mean.p[2] } mean.phiB ~ dunif(0, 1) for (u in 1:2){ mean.psi[u] ~ dunif(0, 1) mean.p[u] ~ dunif(0, 1) } # Define parameters for (i in 1:nind){ # Define probabilities of state S(t+1) given S(t) for (t in f[i]:(n.occasions-1)){ # loop over time # First index = states at time t-1, last index = states at time t ps[1,i,t,1] <- phiA[t] * (1-psiAB[t]) ps[1,i,t,2] <- phiA[t] * psiAB[t] ps[1,i,t,3] <- 1-phiA[t] ps[2,i,t,1] <- phiB[t] * psiBA[t] ps[2,i,t,2] <- phiB[t] * (1-psiBA[t]) ps[2,i,t,3] <- 1-phiB[t] ps[3,i,t,1] <- 0 ps[3,i,t,2] <- 0 ps[3,i,t,3] <- 1 # Define probabilities of O(t) given S(t) # First index = states at time t, last index = observations at time t po[1,i,t,1] <- pA[t] po[1,i,t,2] <- 0 po[1,i,t,3] <- 1-pA[t] po[2,i,t,1] <- 0 po[2,i,t,2] <- pB[t] po[2,i,t,3] <- 1-pB[t] po[3,i,t,1] <- 0 po[3,i,t,2] <- 0 po[3,i,t,3] <- 1 } #t } #i # State-space model likelihood for (i in 1:nind){ z[i,f[i]] <- y[i,f[i]] for (t in (f[i]+1):n.occasions){ # loop over time # State process: draw S(t) given S(t-1) z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,]) # Observation process: draw O(t) given S(t) y[i,t] ~ dcat(po[z[i,t], i, t-1,]) } #t } # i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = rCH, f = f, n.occasions = dim(rCH)[2], nind = dim(rCH)[1], z = known.state.ms(rCH, 3)) # Initial values inits <- function(){list(phiA = runif(n.occasions-1, 0, 1), mean.phiB = runif(1, 0, 1), mean.psi = runif(2, 0, 1), mean.p = runif(2, 0, 1), z = ms.init.z(rCH, f))} # Parameters monitored parameters <- c("phiA", "mean.phiB", "mean.psi", "mean.p") # MCMC settings ni <- 2000 nt <- 3 nb <- 1000 nc <- 3 # Call WinBUGS from R (BRT 1 min) msf <- bugs(bugs.data, inits, parameters, "ms.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Inspect results print(msf, 3) # b) Random time effects # Specify model in BUGS language sink("msrand.bug") cat(" model { #################################################### # Parameters: # phiA: survival probability at site A # phiB: survival probability at site B # psiAB: movement probability from site A to site B # psiBA: movement probability from site B to site A # pA: recapture probability at site A # pB: recapture probability at site B ################################################### # States (S) # 1 alive at A # 2 alive at B # 3 dead # Observations (O) # 1 seen at A # 2 seen at B # 3 not seen ################################################### # Priors and constraints for (t in 1:(n.occasions-1)){ logit(phiA[t]) <- mu + epsilon[t] epsilon[t] ~ dnorm(0, tau)I(-15,15) phiB[t] <- mean.phiB psiAB[t] <- mean.psi[1] psiBA[t] <- mean.psi[2] pA[t] <- mean.p[1] pB[t] <- mean.p[2] } mean.phiB ~ dunif(0, 1) for (u in 1:2){ mean.psi[u] ~ dunif(0, 1) mean.p[u] ~ dunif(0, 1) } mu <- log(mean.phiA / (1- mean.phiA )) # Logit transformation mean.phiA ~ dunif(0, 1) # Prior for mean survival tau <- pow(sigma, -2) sigma ~ dunif(0, 10) # Prior on standard deviation sigma2 <- pow(sigma, 2) # Temporal variance # Define parameters for (i in 1:nind){ # Define probabilities of state S(t+1) given S(t) for (t in f[i]:(n.occasions-1)){ # loop over time # First index = states at time t-1, last index = states at time t ps[1,i,t,1] <- phiA[t] * (1-psiAB[t]) ps[1,i,t,2] <- phiA[t] * psiAB[t] ps[1,i,t,3] <- 1-phiA[t] ps[2,i,t,1] <- phiB[t] * psiBA[t] ps[2,i,t,2] <- phiB[t] * (1-psiBA[t]) ps[2,i,t,3] <- 1-phiB[t] ps[3,i,t,1] <- 0 ps[3,i,t,2] <- 0 ps[3,i,t,3] <- 1 # Define probabilities of O(t) given S(t) # First index = states at time t, last index = observations at time t po[1,i,t,1] <- pA[t] po[1,i,t,2] <- 0 po[1,i,t,3] <- 1-pA[t] po[2,i,t,1] <- 0 po[2,i,t,2] <- pB[t] po[2,i,t,3] <- 1-pB[t] po[3,i,t,1] <- 0 po[3,i,t,2] <- 0 po[3,i,t,3] <- 1 } #t } #i # State-space model likelihood for (i in 1:nind){ z[i,f[i]] <- y[i,f[i]] for (t in (f[i]+1):n.occasions){ # loop over time # State process: draw S(t) given S(t-1) z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,]) # Observation process: draw O(t) given S(t) y[i,t] ~ dcat(po[z[i,t], i, t-1,]) } #t } # i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = rCH, f = f, n.occasions = dim(rCH)[2], nind = dim(rCH)[1], z = known.state.ms(rCH, 3)) # Initial values inits <- function(){list(mean.phiA = runif(1, 0, 1), sigma = runif(1, 0, 1), mean.phiB = runif(1, 0, 1), mean.psi = runif(2, 0, 1), mean.p = runif(2, 0, 1), z = ms.init.z(rCH, f))} # Parameters monitored parameters <- c("phiA", "mean.phiA", "sigma2", "mean.phiB", "mean.psi", "mean.p") # MCMC settings ni <- 10000 nt <- 3 nb <- 5000 nc <- 3 # Call WinBUGS from R (BRT 8 min) ms <- bugs(bugs.data, inits, parameters, "msrand.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Inspect results print(ms, 3) plot(msf$mean$phiA, type = "b", las = 1, ylab = "Survival at site A", xlab = "Time", lwd = 2) points(ms$mean$phiA, type="b", col="red", lwd = 2) abline(h = ms$mean$mean.phiA, lty = 2, col = "red") legend(x = 4.5, y = 0.45, legend = c("Time fixed effect", "Time random effect", "Mean"), lty = c(1,1,2), col = c("black", "red", "red"), bty = "n", lwd = c(2, 2, 1)) ############## # Exercise 3 ############## # Data simulation # Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals phi <- 0.7 psiIO <- 0.4 psiOI <- 0.8 mean.p <- 0.5 v.ind <- 0.4 n.occasions <- 10 n.states <- 3 n.obs <- 2 marked <- matrix(NA, ncol = n.states, nrow = n.occasions) marked[,1] <- rep(100, n.occasions) # present marked[,2] <- rep(0, n.occasions) # absent marked[,3] <- rep(0, n.occasions) # dead # Draw individual recapture probabilities logit.p <- rnorm(sum(marked), qlogis(mean.p), v.ind^0.5) p <- plogis(logit.p) # Define arrays with survival, transition and recapture probabilities # These are 4-dimensional arrays, with # Dimension 1: state of departure # Dimension 2: state of arrival # Dimension 3: individual # Dimension 4: time # 1. State process array totrel <- sum(marked) PSI.STATE <- array(NA, dim = c(n.states, n.states, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.STATE[,,i,t] <- matrix(c( phi*(1-psiIO), phi*psiIO, 1-phi, phi*psiOI, phi*(1-psiOI), 1-phi, 0, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # 2.Observation array PSI.OBS <- array(NA, dim = c(n.states, n.obs, totrel, n.occasions-1)) for (i in 1:totrel){ for (t in 1:(n.occasions-1)){ PSI.OBS[,,i,t] <- matrix(c( p[i], 1-p[i], 0, 1, 0, 1 ), nrow = n.states, byrow = TRUE) } #t } #i # Execute simulation function sim <- simul.ms(PSI.STATE, PSI.OBS, marked) CH <- sim$CH # Compute vector with occasion of first capture get.first <- function(x) min(which(x!=0)) f <- apply(CH, 1, get.first) # Recode CH matrix: note, a 0 is not allowed! # 1 = seen alive, 2 = not seen rCH <- CH # recoded CH rCH[rCH==0] <- 2 # Data analysis # Specify model in BUGS language sink("tempemi.bug") cat(" model { #################################################### # Parameters: # phi: survival probability # psiIO: probability to emigrate # psiOI: probability to immigrate # p: recapture probability ################################## # States (S) # 1 alive and present # 2 alive and absent # 3 dead # Observations (O) # 1 seen # 2 not seen ################################## # Priors and constraints for (t in 1:(n.occasions-1)){ phi[t] <- mean.phi psiIO[t] <- mean.psiIO psiOI[t] <- mean.psiOI } mean.phi ~ dunif(0, 1) mean.psiIO ~ dunif(0, 1) mean.psiOI ~ dunif(0, 1) for (i in 1:nind){ for (t in 1:(n.occasions-1)){ logit(p[i,t]) <- mu + epsilon[i] } #t epsilon[i] ~ dnorm(0, tau)I(-15,15) } # i mu <- log(mean.p / (1-mean.p)) # Logit transformation mean.p ~ dunif(0, 1) # Prior for mean recapture tau <- pow(sigma, -2) sigma ~ dunif(0, 5) # Prior for standard deviation sigma2 <-pow(sigma, 2) # Define parameters for (i in 1:nind){ # Define probabilities of state S(t+1) given S(t) for (t in f[i]:(n.occasions-1)){ # loop over time ps[1,i,t,1] <- phi[t] * (1-psiIO[t]) ps[1,i,t,2] <- phi[t] * psiIO[t] ps[1,i,t,3] <- 1-phi[t] ps[2,i,t,1] <- phi[t] * psiOI[t] ps[2,i,t,2] <- phi[t] * (1-psiOI[t]) ps[2,i,t,3] <- 1-phi[t] ps[3,i,t,1] <- 0 ps[3,i,t,2] <- 0 ps[3,i,t,3] <- 1 # Define probabilities of O(t) given S(t) po[1,i,t,1] <- p[i,t] po[1,i,t,2] <- 1-p[i,t] po[2,i,t,1] <- 0 po[2,i,t,2] <- 1 po[3,i,t,1] <- 0 po[3,i,t,2] <- 1 } #t } #i # State-space model likelihood for (i in 1:nind){ z[i,f[i]] <- y[i,f[i]] for (t in (f[i]+1):n.occasions){ # loop over time # State process: draw S(t) given S(t-1) z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,]) # Observation process: draw O(t) given S(t) y[i,t] ~ dcat(po[z[i,t], i, t-1,]) } #t } # i } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = rCH, f = f, n.occasions = dim(rCH)[2], nind = dim(rCH)[1], z = known.state.ms(rCH, 2)) # Initial values inits <- function(){list(mean.phi = runif(1, 0, 1), mean.psiIO = runif(1, 0, 1), mean.psiOI = runif(1, 0, 1), mean.p = runif(1, 0, 1), sigma = runif(1, 0, 1), z = ms.init.z(rCH, f))} # Parameters monitored parameters <- c("mean.phi", "mean.psiIO", "mean.psiOI", "mean.p", "sigma2") # MCMC settings ni <- 30000 nt <- 3 nb <- 20000 nc <- 3 # Call WinBUGS from R (BRT 113 min) tempemi <- bugs(bugs.data, inits, parameters, "tempemi.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(tempemi, 3) ############################################################################# # # Chapter 10 # ############################################################################# ############## # Exercise 1 ############## # Data simulation # Define the parameters n.occasions <- 8 N <- 300 phi.m <- rep(0.75, n.occasions-1) phi.f <- rep(0.5, n.occasions-1) b <- c(0.3, rep(0.1, n.occasions-1)) p <- rep(0.4, n.occasions) PHI.M <- matrix(rep(phi.m, (n.occasions-1)*N), ncol = n.occasions-1, nrow = N, byrow = T) PHI.F <- matrix(rep(phi.f, (n.occasions-1)*N), ncol = n.occasions-1, nrow = N, byrow = T) P <- matrix(rep(p, n.occasions*N), ncol = n.occasions, nrow = N, byrow = T) # Apply simulation function sim.m <- simul.js(PHI.M, P, b, N) sim.f <- simul.js(PHI.F, P, b, N) CH.m <- sim.m$CH CH.f <- sim.f$CH # Data analysis # Specify model in BUGS language sink("js-rest.occ.bug") cat(" model { # Priors and constraints for (i in 1:M){ for (t in 1:(n.occasions-1)){ phi[i,t] <- beta[group[i]] } # t for (t in 1:n.occasions){ p[i,t] <- mean.p } #t } #i for (i in 1:2){ beta[i] ~ dunif(0, 1) } mean.p ~ dunif(0, 1) for (t in 1:n.occasions){ gamma[t] ~ dunif(0, 1) } #t # Define the likelihoods for (i in 1:M){ # First occasion # State process z[i,1] ~ dbern(gamma[1]) mu1[i] <- z[i,1] * p[i,1] # Observation process y[i,1] ~ dbern(mu1[i]) # Subsequent occasions for (t in 2:n.occasions){ # State process q[i,t-1] <- 1-z[i,t-1] mu2[i,t] <- phi[i,t-1]*z[i,t-1] + gamma[t]*prod(q[i,1:(t-1)]) z[i,t] ~ dbern(mu2[i,t]) # Observation process mu3[i,t] <- z[i,t] * p[i,t] y[i,t] ~ dbern(mu3[i,t]) } # t } # i # Calculate derived population parameters for (t in 1:n.occasions){ qgamma[t] <- 1-gamma[t] } cprob[1] <- gamma[1] for (t in 2:n.occasions){ cprob[t] <- gamma[t] * prod(qgamma[1:(t-1)]) } # t psi <- sum(cprob[]) # Inclusion probability for (t in 1:n.occasions){ b[t] <- cprob[t] / psi # Entry probability } # t for (i in 1:M){ recruit[i,1] <- z[i,1] for (t in 2:n.occasions){ recruit[i,t] <- (1-z[i,t-1]) * z[i,t] } # t } # i for (t in 1:n.occasions){ Nm[t] <- sum(z[1:mm,t]) # Actual population size of males Nf[t] <- sum(z[(mm+1):M,t]) # Actual population size of females Bm[t] <- sum(recruit[1:mm,t]) # Number of entries of males Bf[t] <- sum(recruit[(mm+1):M,t]) # Number of entries of females } # t for (i in 1:M){ Nind[i] <- sum(z[i,1:n.occasions]) Nalive[i] <- 1-equals(Nind[i], 0) } # i Nsuperm <- sum(Nalive[1:mm]) # Size of superpopulation of males Nsuperf <- sum(Nalive[(mm+1):M]) # Size of superpopulation of females } ",fill=TRUE) sink() # Augment the capture-histories by pseudo-individuals nz <- 500 CHm.aug <- rbind(CH.m, matrix(0, ncol = dim(CH.m)[2], nrow = nz)) m <- rep(1, dim(CHm.aug)[1]) CHf.aug <- rbind(CH.f, matrix(0, ncol = dim(CH.f)[2], nrow = nz)) f <- rep(2, dim(CHf.aug)[1]) group <- c(m, f) y <- rbind(CHm.aug, CHf.aug) # Bundle data bugs.data <- list(y = y, n.occasions = dim(y)[2], M = dim(y)[1], group = group, mm = length(m)) # Initial values inits <- function(){list(beta = runif(2, 0, 1), mean.p = runif(1, 0, 1), z = y)} # Parameters monitored parameters <- c("psi", "mean.p", "beta", "b", "Nsuperm", "Nsuperf", "Nm", "Nf", "Bm", "Bf") # MCMC settings niter <- 5000 nthin <- 3 nburn <- 3000 nchains <- 3 # Call WinBUGS from R (BRT 32 min) js.occ <- bugs(bugs.data, inits, parameters, "js-rest.occ.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(js.occ, 3) ############## # Exercise 2 ############## # Data simulation # Define the parameters n.occasions <- 7 N <- 500 phi <- rep(0.5, n.occasions-1) b <- c(0.4, rep(0.1, n.occasions-1)) p <- rep(0.5, n.occasions) PHI <- matrix(rep(phi, (n.occasions-1)*N), ncol = n.occasions-1, nrow = N, byrow = T) P <- matrix(rep(p, n.occasions*N), ncol = n.occasions, nrow = N, byrow = T) # Apply simulation function sim <- simul.js(PHI, P, b, N) CH <- sim$CH # Data analysis # Specify model in BUGS language sink("js-super.bug") cat(" model { # Priors and constraints for (i in 1:M){ for (t in 1:(n.occasions-1)){ phi[i,t] <- mean.phi } # t for (t in 1:n.occasions){ p[i,t] <- mean.p } #t } #i mean.phi ~ dunif(0, 1) # Prior for mean survival mean.p ~ dunif(0, 1) # Prior for mean capture psi ~ dunif(0, 1) # Prior for inclusion probability # Choose priors for the first two free b lb[1] ~ dunif(0, 1) lb[2] ~ dunif(0, 1) for (t in 3:n.occasions){ lb[t] <- lb[2] } # Weigh the lb such that they sum to 1 for (t in 1:n.occasions){ b[t] <- lb[t] / sum(lb[]) } # Convert entry probs to conditional entry probs nu[1] <- b[1] for (t in 2:n.occasions){ nu[t] <- b[t] / (1-sum(b[1:(t-1)])) } # t # Define the likelihood for (i in 1:M){ # First occasion # State process w[i] ~ dbern(psi) # Draw latent inclusion variable z[i,1] ~ dbern(nu[1]) # Observation process mu1[i] <- z[i,1] * p[i,1] * w[i] y[i,1] ~ dbern(mu1[i]) # Subsequent occasions for (t in 2:n.occasions){ # State process q[i,t-1] <- 1-z[i,t-1] mu2[i,t] <- phi[i,t-1] * z[i,t-1] + nu[t] * prod(q[i,1:(t-1)]) z[i,t] ~ dbern(mu2[i,t]) # Observation process mu3[i,t] <- z[i,t] * p[i,t] * w[i] y[i,t] ~ dbern(mu3[i,t]) } #t } #i # Calculate derived population parameters for (i in 1:M){ for (t in 1:n.occasions){ u[i,t] <- z[i,t]*w[i] # Deflated latent state (u) } } for (i in 1:M){ recruit[i,1] <- u[i,1] for (t in 2:n.occasions){ recruit[i,t] <- (1-u[i,t-1]) * u[i,t] } #t } #i for (t in 1:n.occasions){ N[t] <- sum(u[1:M,t]) # Actual population size B[t] <- sum(recruit[1:M,t]) # Number of entries } #t for (i in 1:M){ Nind[i] <- sum(u[i,1:n.occasions]) Nalive[i] <- 1-equals(Nind[i], 0) } #i Nsuper <- sum(Nalive[]) # Size of superpopulation } ",fill=TRUE) sink() # Augment the capture-histories by nz pseudo-individuals nz <- 500 CH.aug <- rbind(CH, matrix(0, ncol = dim(CH)[2], nrow = nz)) # Bundle data bugs.data <- list(y = CH.aug, n.occasions = dim(CH.aug)[2], M = dim(CH.aug)[1]) # Initial values inits <- function(){list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), psi = runif(1, 0, 1), lb = c(runif(2, 0, 0.1), rep(NA, n.occasions-2)), z = CH.aug)} # Parameters monitored parameters <- c("psi", "mean.p", "mean.phi", "b", "Nsuper", "N", "B") # MCMC settings niter <- 5000 nthin <- 3 nburn <- 3000 nchains <- 3 # Call WinBUGS from R (BRT 24 min) js.super <- bugs(bugs.data, inits, parameters, "js-super.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(js.super, digits = 3) ############## # Exercise 3 ############## # Data simulation # Function to simulate capture-recapture data for JS analysis that also includes dead recoveries # Dead recoveries are coded with a 2 in the resulting capture histories simul.jsrecov <- function(PHI, P, R, b, N){ B <- rmultinom(1, N, b) # Generate no. of entering ind. per occasion n.occasions <- dim(PHI)[2] + 1 CH.sur <- CH.p <- CH.r <- matrix(0, ncol = n.occasions, nrow = N) # Define a vector with the occasion of entering the population ent.occ <- numeric() for (t in 1:n.occasions){ ent.occ <- c(ent.occ, rep(t, B[t])) } # Modeling survival for (i in 1:N){ CH.sur[i, ent.occ[i]] <- 1 # Write 1 when ind. enters the pop. if (ent.occ[i] == n.occasions) next for (t in (ent.occ[i]+1):n.occasions){ # Bernoulli trial: has individual survived occasion? sur <- rbinom(1, 1, PHI[i,t-1]) if (sur==1) CH.sur[i,t] <- 1 if (sur==0) CH.sur[i,t] <- 2 if (sur==0) break } #t } #i # Modeling capture for (i in 1:N){ CH.p[i,] <- rbinom(n.occasions, 1, P[i,]) } #i # Modeling recovery for (i in 1:N){ CH.r[i,] <- rbinom(n.occasions, 1, R[i,]) } #i # Full capture-recapture matrix CH <- CH.sur * CH.p CH.2 <- CH.sur * CH.r CH[CH==2] <- 0 CH[CH.2==2] <- 2 # Remove individuals never captured cap.sum <- rowSums(CH) never <- which(cap.sum == 0) CH <- CH[-never,] # Remove individuals that were found dead, but never marked cap.sum <- rowSums(CH) two <- which(cap.sum==2) CH.help <- CH[two,] CH.help[CH.help==0] <- 1 w <- apply(CH.help, 1, prod) rem <- two[which(w==2)] CH <- CH[-rem,] # Actual population size CH.pop <- CH.sur CH.pop[CH.pop==2] <- 0 Nt <- colSums(CH.pop) return(list(CH=CH, B=B, N=Nt)) } # Define the parameters n.occasions <- 10 N <- 500 phi <- 0.5 b <- rep(0.1, 10) p <- 0.6 r <- 0.2 PHI <- matrix(rep(phi, (n.occasions-1)*N), ncol = n.occasions-1, nrow = N, byrow = T) P <- matrix(rep(p, n.occasions*N), ncol = n.occasions, nrow = N, byrow = T) R <- matrix(rep(r, n.occasions*N), ncol = n.occasions, nrow = N, byrow = T) # Apply simulation function sim <- simul.jsrecov(PHI, P, R, b, N) CH <- sim$CH # Specify model in BUGS language sink("js-ms.bug") cat(" model { #-------------------------------------- # Parameters: # phi: survival probability # gamma: removal entry probability # p: capture probability #-------------------------------------- # States (S): # 1 not yet entered # 2 alive # 3 recently dead and recovered # 4 dead or recently dead, but not recovered # Observations (O): # 1 seen # 2 recovered # 3 neither seen nor recovered #-------------------------------------- # Priors and constraints for (t in 1:(n.occasions-1)){ phi[t] <- mean.phi gamma[t] ~ dunif(0, 1) # Prior for entry probabilities p[t] <- mean.p r[t] <- mean.r } mean.phi ~ dunif(0, 1) # Prior for mean survival mean.p ~ dunif(0, 1) # Prior for mean capture mean.r ~ dunif(0, 1) # Prior for mean recovery # Define state-transition and observation matrices for (i in 1:M){ # Define probabilities of state S(t+1) given S(t) for (t in 1:(n.occasions-1)){ ps[1,i,t,1] <- 1-gamma[t] ps[1,i,t,2] <- gamma[t] ps[1,i,t,3] <- 0 ps[1,i,t,4] <- 0 ps[2,i,t,1] <- 0 ps[2,i,t,2] <- phi[t] ps[2,i,t,3] <- (1-phi[t])*r[t] ps[2,i,t,4] <- (1-phi[t])*(1-r[t]) ps[3,i,t,1] <- 0 ps[3,i,t,2] <- 0 ps[3,i,t,3] <- 0 ps[3,i,t,4] <- 1 ps[4,i,t,1] <- 0 ps[4,i,t,2] <- 0 ps[4,i,t,3] <- 0 ps[4,i,t,4] <- 1 # Define probabilities of O(t) given S(t) po[1,i,t,1] <- 0 po[1,i,t,2] <- 0 po[1,i,t,3] <- 1 po[2,i,t,1] <- p[t] po[2,i,t,2] <- 0 po[2,i,t,3] <- 1-p[t] po[3,i,t,1] <- 0 po[3,i,t,2] <- 1 po[3,i,t,3] <- 0 po[4,i,t,1] <- 0 po[4,i,t,2] <- 0 po[4,i,t,3] <- 1 } #t } #i # Likelihood for (i in 1:M){ # Define latent state at first occasion z[i,1] <- 1 # Make sure that all M individuals are in state 1 at t=1 for (t in 2:n.occasions){ # State process: draw S(t) given S(t-1) z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,]) # Observation process: draw O(t) given S(t) y[i,t] ~ dcat(po[z[i,t], i, t-1,]) } #t } #i # Calculate derived population parameters for (t in 1:(n.occasions-1)){ qgamma[t] <- 1-gamma[t] } cprob[1] <- gamma[1] for (t in 2:(n.occasions-1)){ cprob[t] <- gamma[t] * prod(qgamma[1:(t-1)]) } #t psi <- sum(cprob[]) # Inclusion probability for (t in 1:(n.occasions-1)){ b[t] <- cprob[t] / psi # Entry probability } #t for (i in 1:M){ for (t in 2:n.occasions){ al[i,t-1] <- equals(z[i,t], 2) } #t for (t in 1:(n.occasions-1)){ d[i,t] <- equals(z[i,t]-al[i,t],0) } #t alive[i] <- sum(al[i,]) } #i for (t in 1:(n.occasions-1)){ N[t] <- sum(al[,t]) # Actual population size B[t] <- sum(d[,t]) # Number of entries } #t for (i in 1:M){ w[i] <- 1-equals(alive[i],0) } #i Nsuper <- sum(w[]) # Superpopulation size } ",fill = TRUE) sink() # Add dummy occasion CH.du <- cbind(rep(0, dim(CH)[1]), CH) # Augment data nz <- 500 CH.ms <- rbind(CH.du, matrix(0, ncol = dim(CH.du)[2], nrow = nz)) # Recode CH matrix: a 0 is not allowed in WinBUGS! CH.ms[CH.ms==0] <- 3 # Not seen = 3, seen = 1, recovered = 3 # Bundle data bugs.data <- list(y = CH.ms, n.occasions = dim(CH.ms)[2], M = dim(CH.ms)[1]) # Initial values inits <- function(){list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), mean.r = runif(1, 0, 0.5), z = cbind(rep(NA, dim(CH.ms)[1]), CH.ms[,-1]))} # Parameters monitored parameters <- c("mean.p", "mean.r", "mean.phi", "b", "Nsuper", "N", "B") # MCMC settings ni <- 50000 nt <- 3 nb <- 25000 nc <- 3 # Call WinBUGS from R (BRT 203 min) js.ms <- bugs(bugs.data, inits, parameters, "js-ms.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(js.ms, digits = 3) # Remove recoveries CH.msalt <- CH.ms CH.msalt[CH.msalt==2] <- 3 CH.msalt[CH.msalt==3] <- 2 # 1: seen alive, 2: not seen # Specify model in BUGS language sink("js-ms.bug") cat(" model { #-------------------------------------- # Parameters: # phi: survival probability # gamma: removal entry probability # p: capture probability #-------------------------------------- # States (S): # 1 not yet entered # 2 alive # 3 dead # Observations (O): # 1 seen # 2 not seen #-------------------------------------- # Priors and constraints for (t in 1:(n.occasions-1)){ phi[t] <- mean.phi gamma[t] ~ dunif(0, 1) # Prior for entry probabilities p[t] <- mean.p } mean.phi ~ dunif(0, 1) # Prior for mean survival mean.p ~ dunif(0, 1) # Prior for mean capture # Define state-transition and observation matrices for (i in 1:M){ # Define probabilities of state S(t+1) given S(t) for (t in 1:(n.occasions-1)){ ps[1,i,t,1] <- 1-gamma[t] ps[1,i,t,2] <- gamma[t] ps[1,i,t,3] <- 0 ps[2,i,t,1] <- 0 ps[2,i,t,2] <- phi[t] ps[2,i,t,3] <- 1-phi[t] ps[3,i,t,1] <- 0 ps[3,i,t,2] <- 0 ps[3,i,t,3] <- 1 # Define probabilities of O(t) given S(t) po[1,i,t,1] <- 0 po[1,i,t,2] <- 1 po[2,i,t,1] <- p[t] po[2,i,t,2] <- 1-p[t] po[3,i,t,1] <- 0 po[3,i,t,2] <- 1 } #t } #i # Likelihood for (i in 1:M){ # Define latent state at first occasion z[i,1] <- 1 # Make sure that all M individuals are in state 1 at t=1 for (t in 2:n.occasions){ # State process: draw S(t) given S(t-1) z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,]) # Observation process: draw O(t) given S(t) y[i,t] ~ dcat(po[z[i,t], i, t-1,]) } #t } #i # Calculate derived population parameters for (t in 1:(n.occasions-1)){ qgamma[t] <- 1-gamma[t] } cprob[1] <- gamma[1] for (t in 2:(n.occasions-1)){ cprob[t] <- gamma[t] * prod(qgamma[1:(t-1)]) } #t psi <- sum(cprob[]) # Inclusion probability for (t in 1:(n.occasions-1)){ b[t] <- cprob[t] / psi # Entry probability } #t for (i in 1:M){ for (t in 2:n.occasions){ al[i,t-1] <- equals(z[i,t], 2) } #t for (t in 1:(n.occasions-1)){ d[i,t] <- equals(z[i,t]-al[i,t],0) } #t alive[i] <- sum(al[i,]) } #i for (t in 1:(n.occasions-1)){ N[t] <- sum(al[,t]) # Actual population size B[t] <- sum(d[,t]) # Number of entries } #t for (i in 1:M){ w[i] <- 1-equals(alive[i],0) } #i Nsuper <- sum(w[]) # Superpopulation size } ",fill = TRUE) sink() # Bundle data bugs.data <- list(y = CH.msalt, n.occasions = dim(CH.msalt)[2], M = dim(CH.msalt)[1]) # Initial values inits <- function(){list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = cbind(rep(NA, dim(CH.msalt)[1]), CH.msalt[,-1]))} # Parameters monitored parameters <- c("mean.p", "mean.phi", "b", "Nsuper", "N", "B") # MCMC settings ni <- 50000 nt <- 3 nb <- 25000 nc <- 3 # Call WinBUGS from R (BRT 240 min) js.msalt <- bugs(bugs.data, inits, parameters, "js-ms.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(js.msalt, 3) ############################################################################# # # Chapter 11 # ############################################################################# ############## # Exercise 1 ############## # Load data nyears <- 9 # Number of years # Capture recapture data: m-array of juveniles and adults (these are males and females together) marray.j <- matrix (c(15, 3, 0, 0, 0, 0, 0, 0, 198, 0, 34, 9, 1, 0, 0, 0, 0, 287, 0, 0, 56, 8, 1, 0, 0, 0, 455, 0, 0, 0, 48, 3, 1, 0, 0, 518, 0, 0, 0, 0, 45, 13, 2, 0, 463, 0, 0, 0, 0, 0, 27, 7, 0, 493, 0, 0, 0, 0, 0, 0, 37, 3, 434, 0, 0, 0, 0, 0, 0, 0, 39, 405), nrow = 8, ncol = 9, byrow = TRUE) marray.a <- matrix(c(14, 2, 0, 0, 0, 0, 0, 0, 43, 0, 22, 4, 0, 0, 0, 0, 0, 44, 0, 0, 34, 2, 0, 0, 0, 0, 79, 0, 0, 0, 51, 3, 0, 0, 0, 94, 0, 0, 0, 0, 45, 3, 0, 0, 118, 0, 0, 0, 0, 0, 44, 3, 0, 113, 0, 0, 0, 0, 0, 0, 48, 2, 99, 0, 0, 0, 0, 0, 0, 0, 51, 90), nrow = 8, ncol = 9, byrow = TRUE) # Population population count data popcount <- c(32, 42, 64, 85, 82, 78, 73, 69, 79) # Reproductive success J <- c(189, 274, 398, 538, 520, 476, 463, 438, 507) # number offspring R <- c(28, 36, 57, 77, 81, 83, 77, 72, 85) # number of surveyed broods # Number of years with predictions t.pred <- 3 # Data analysis # Specify model in BUGS language sink("ipm.hoopoe.bug") cat(" model { ########################################################################## # Integrated population model # - Age structured model with 2 age classes: # 1-year old and at least 2 years old # - Age at first breeding = 1 year # - Pre-breeding census, female-based # - All vital rates are assumed to be time-dependent (random) # - Explicit estimate of immigration ########################################################################## ############################################################ # 1. Define the priors for the parameters ############################################################ # Initial population sizes N1[1] ~ dnorm(10, 0.001)I(0,) # 1-year old individuals NadSurv[1] ~ dnorm(10, 0.001)I(0,) # Adults >= 2 years Nadimm[1] ~ dnorm(10, 0.001)I(0,) # Immigrants # Mean demographic parameters mphij ~ dunif(0, 1) mphia ~ dunif(0, 1) mfec ~ dunif(0, 15) mim ~ dunif(0, 3) mp ~ dunif(0, 1) # Precision of standard deviations of temporal variability sig.phij ~ dunif(0, 10) tau.phij <- pow(sig.phij, -2) sig.phia ~ dunif(0, 10) tau.phia <- pow(sig.phia, -2) sig.fec ~ dunif(0, 10) tau.fec <- pow(sig.fec, -2) sig.im ~ dunif(0, 10) tau.im <- pow(sig.im, -2) # Distribution of error terms (Bounded to help with convergence) for (t in 1:(nyears-1+t.pred)){ epsilon.phij[t] ~ dnorm(0, tau.phij)I(-15,15) epsilon.phia[t] ~ dnorm(0, tau.phia)I(-15,15) epsilon.fec[t] ~ dnorm(0, tau.fec)I(-15,15) epsilon.im[t] ~ dnorm(0, tau.im)I(-15,15) } ############################################################ # 2. Constrain parameters ############################################################ for (t in 1:(nyears-1+t.pred)){ logit(phij[t]) <- l.mphij + epsilon.phij[t] # Juv. apparent survival logit(phia[t]) <- l.mphia + epsilon.phia[t] # Adult apparent survival log(f[t]) <- l.mfec + epsilon.fec[t] # Fecundity log(omega[t]) <- l.mim + epsilon.im[t] # Immigration p[t] <- mp # Recapture probability } ############################################################ # 3. Derived parameters ############################################################ l.mphij <- log(mphij / (1-mphij)) # Logit mean juv. survival l.mphia <- log(mphia / (1-mphia)) # Logit mean adult survival l.mfec <- log(mfec) # Log mean fecundity l.mim <- log(mim) # Log mean immigration rate # Population growth rate for (t in 1:(nyears-1+t.pred)){ lambda[t] <- Ntot[t+1] / Ntot[t] logla[t] <- log(lambda[t]) } mlam <- exp((1/(nyears-1))*sum(logla[1:(nyears-1)])) # Geometric mean ###################################################################### # 4. The likelihoods of the single data sets ###################################################################### ################################################################### # 4.1. Likelihood for population count data (state-space model) ################################################################### ############################# # 4.1.1 System process ############################# for (t in 2:nyears+t.pred){ mean1[t] <- 0.5 * f[t-1] * phij[t-1] * Ntot[t-1] N1[t] ~ dpois(mean1[t]) NadSurv[t] ~ dbin(phia[t-1], Ntot[t-1]) mpo[t] <- Ntot[t-1] * omega[t-1] Nadimm[t] ~ dpois(mpo[t]) } for (t in 1:nyears+t.pred){ Ntot[t] <- NadSurv[t] + Nadimm[t] + N1[t] } ############################# # 4.1.2 Observation process ############################# for (t in 1:nyears){ y[t] ~ dpois(Ntot[t]) } ###################################################################### # 4.2 Likelihood for capture-recapture data: CJS model (2 age classes) ###################################################################### # Multinomial likelihood for (t in 1:(nyears-1)){ marray.j[t,1:nyears] ~ dmulti(pr.j[t,], r.j[t]) marray.a[t,1:nyears] ~ dmulti(pr.a[t,], r.a[t]) } # Calculate number of released individuals for (t in 1:(nyears-1)){ r.j[t] <- sum(marray.j[t,]) r.a[t] <- sum(marray.a[t,]) } # m-array cell probabilities for juveniles for (t in 1:(nyears-1)){ q[t] <- 1-p[t] # Main diagonal pr.j[t,t] <- phij[t]*p[t] # Above main diagonal for (j in (t+1):(nyears-1)){ pr.j[t,j] <- phij[t]*prod(phia[(t+1):j])*prod(q[t:(j-1)])*p[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr.j[t,j] <- 0 } # j # Last column pr.j[t,nyears] <- 1-sum(pr.j[t,1:(nyears-1)]) } # t # m-array cell probabilities for adults for (t in 1:(nyears-1)){ # Main diagonal pr.a[t,t] <- phia[t]*p[t] # above main diagonal for (j in (t+1):(nyears-1)){ pr.a[t,j] <- prod(phia[t:j])*prod(q[t:(j-1)])*p[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr.a[t,j] <- 0 } # j # Last column pr.a[t,nyears] <- 1-sum(pr.a[t,1:(nyears-1)]) } # t ############################################################ # 4.3. Likelihood for reproductive data: Poisson regression ############################################################ for (t in 1:nyears){ J[t] ~ dpois(rho[t]) rho[t] <- R[t] * f[t] } } # End Model ",fill = TRUE) sink() # Bundle data bugs.data <- list(nyears = nyears, marray.j = marray.j, marray.a = marray.a, y = popcount, J = J, R = R, t.pred = t.pred) # Initial values inits <- function(){list(mphij = runif(1, 0.05, 0.2), mphia = runif(1, 0.2, 0.5), mfec = runif(1, 4, 6), mim = runif(1, 0, 0.3), mp = runif(1, 0.5, 0.8), sig.phij = runif(1, 0.1, 1), sig.phia = runif(1, 0.1, 1), sig.fec = runif(1, 0.1, 1), sig.im = runif(1, 0.1, 1), N1 = round(runif(nyears+t.pred, 10, 20), 0), NadSurv = round(runif(nyears+t.pred, 10, 30), 0), Nadimm = round(runif(nyears+t.pred, 1, 20), 0))} # Parameters monitored parameters <- c("phij", "phia", "f", "omega", "p", "lambda", "mphij", "mphia", "mfec", "mim", "mlam", "sig.phij", "sig.phia", "sig.fec", "sig.im", "N1", "NadSurv", "Nadimm", "Ntot") # MCMC settings niter <- 10000 nthin <- 3 nburn <- 5000 nchains <- 3 # Call WinBUGS from R (BRT 4 min) ipm.hoopoe.1 <- bugs(bugs.data, inits, parameters, "ipm.hoopoe.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Inspect results print(ipm.hoopoe1, 3) # Compute extinction probability mean(ipm.hoopoe.1$sims.list$Ntot[,12]<5) ############## # Exercise 2 ############## # Load data nyears <- 9 # Number of years # Capture recapture data: m-array of juveniles and adults (these are males and females together) marray.j <- matrix (c(15, 3, 0, 0, 0, 0, 0, 0, 198, 0, 34, 9, 1, 0, 0, 0, 0, 287, 0, 0, 56, 8, 1, 0, 0, 0, 455, 0, 0, 0, 48, 3, 1, 0, 0, 518, 0, 0, 0, 0, 45, 13, 2, 0, 463, 0, 0, 0, 0, 0, 27, 7, 0, 493, 0, 0, 0, 0, 0, 0, 37, 3, 434, 0, 0, 0, 0, 0, 0, 0, 39, 405), nrow = 8, ncol = 9, byrow = TRUE) marray.a <- matrix(c(14, 2, 0, 0, 0, 0, 0, 0, 43, 0, 22, 4, 0, 0, 0, 0, 0, 44, 0, 0, 34, 2, 0, 0, 0, 0, 79, 0, 0, 0, 51, 3, 0, 0, 0, 94, 0, 0, 0, 0, 45, 3, 0, 0, 118, 0, 0, 0, 0, 0, 44, 3, 0, 113, 0, 0, 0, 0, 0, 0, 48, 2, 99, 0, 0, 0, 0, 0, 0, 0, 51, 90), nrow = 8, ncol = 9, byrow = TRUE) # Population population count data popcount <- c(32, 42, NA, 85, NA, 78, 73, 69, 79) # Reproductive success J <- c(189, 274, 398, 538, 520, 476, 463, 438, 507) # number offspring R <- c(28, 36, 57, 77, 81, 83, 77, 72, 85) # number of surveyed broods # Data analysis # Specify model in BUGS language sink("ipm.hoopoe.bug") cat(" model { ########################################################################## # Integrated population model # - Age structured model with 2 age classes: # 1-year old and at least 2 years old # - Age at first breeding = 1 year # - Pre-breeding census, female-based # - All vital rates are assumed to be time-dependent (random) # - Explicit estimate of immigration ########################################################################## ############################################################ # 1. Define the priors for the parameters ############################################################ # Initial population sizes N1[1] ~ dnorm(10, 0.001)I(0,) # 1-year old individuals NadSurv[1] ~ dnorm(10, 0.001)I(0,) # Adults >= 2 years Nadimm[1] ~ dnorm(10, 0.001)I(0,) # Immigrants # Mean demographic parameters mphij ~ dunif(0, 1) mphia ~ dunif(0, 1) mfec ~ dunif(0, 15) mim ~ dunif(0, 3) mp ~ dunif(0, 1) # Precision of standard deviations of temporal variability sig.phij ~ dunif(0, 10) tau.phij <- pow(sig.phij, -2) sig.phia ~ dunif(0, 10) tau.phia <- pow(sig.phia, -2) sig.fec ~ dunif(0, 10) tau.fec <- pow(sig.fec, -2) sig.im ~ dunif(0, 10) tau.im <- pow(sig.im, -2) # Distribution of error terms (Bounded to help with convergence) for (t in 1:(nyears-1)){ epsilon.phij[t] ~ dnorm(0, tau.phij)I(-15,15) epsilon.phia[t] ~ dnorm(0, tau.phia)I(-15,15) epsilon.im[t] ~ dnorm(0, tau.im)I(-15,15) } for (t in 1:nyears){ epsilon.fec[t] ~ dnorm(0, tau.fec)I(-15,15) } ############################################################ # 2. Constrain parameters ############################################################ for (t in 1:(nyears-1)){ logit(phij[t]) <- l.mphij + epsilon.phij[t] # Juv. apparent survival logit(phia[t]) <- l.mphia + epsilon.phia[t] # Adult apparent survival log(omega[t]) <- l.mim + epsilon.im[t] # Immigration p[t] <- mp # Recapture probability } # Fecundity: note data from an additional year compared to the other vital rates is available for (t in 1:nyears){ log(f[t]) <- l.mfec + epsilon.fec[t] } ############################################################ # 3. Derived parameters ############################################################ l.mphij <- log(mphij / (1-mphij)) # Logit mean juv. survival l.mphia <- log(mphia / (1-mphia)) # Logit mean adult survival l.mfec <- log(mfec) # Log mean fecundity l.mim <- log(mim) # Log mean immigration rate # Population growth rate for (t in 1:(nyears-1)){ lambda[t] <- Ntot[t+1] / Ntot[t] logla[t] <- log(lambda[t]) } mlam <- exp((1/(nyears-1))*sum(logla[1:(nyears-1)])) # Geometric mean ###################################################################### # 4. The likelihoods of the single data sets ###################################################################### ################################################################### # 4.1. Likelihood for population count data (state-space model) ################################################################### ############################# # 4.1.1 System process ############################# for (t in 2:nyears){ mean1[t] <- 0.5 * f[t-1] * phij[t-1] * Ntot[t-1] N1[t] ~ dpois(mean1[t]) NadSurv[t] ~ dbin(phia[t-1], Ntot[t-1]) mpo[t] <- Ntot[t-1] * omega[t-1] Nadimm[t] ~ dpois(mpo[t]) } ############################# # 4.1.2 Observation process ############################# for (t in 1:nyears){ Ntot[t] <- NadSurv[t] + Nadimm[t] + N1[t] y[t] ~ dpois(Ntot[t]) } ###################################################################### # 4.2 Likelihood for capture-recapture data: CJS model (2 age classes) ###################################################################### # Multinomial likelihood for (t in 1:(nyears-1)){ marray.j[t,1:nyears] ~ dmulti(pr.j[t,], r.j[t]) marray.a[t,1:nyears] ~ dmulti(pr.a[t,], r.a[t]) } # Calculate number of released individuals for (t in 1:(nyears-1)){ r.j[t] <- sum(marray.j[t,]) r.a[t] <- sum(marray.a[t,]) } # m-array cell probabilities for juveniles for (t in 1:(nyears-1)){ q[t] <- 1-p[t] # Main diagonal pr.j[t,t] <- phij[t]*p[t] # Above main diagonal for (j in (t+1):(nyears-1)){ pr.j[t,j] <- phij[t]*prod(phia[(t+1):j])*prod(q[t:(j-1)])*p[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr.j[t,j] <- 0 } # j # Last column pr.j[t,nyears] <- 1-sum(pr.j[t,1:(nyears-1)]) } # t # m-array cell probabilities for adults for (t in 1:(nyears-1)){ # Main diagonal pr.a[t,t] <- phia[t]*p[t] # above main diagonal for (j in (t+1):(nyears-1)){ pr.a[t,j] <- prod(phia[t:j])*prod(q[t:(j-1)])*p[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr.a[t,j] <- 0 } # j # Last column pr.a[t,nyears] <- 1-sum(pr.a[t,1:(nyears-1)]) } # t ############################################################ # 4.3. Likelihood for reproductive data: Poisson regression ############################################################ for (t in 1:nyears){ J[t] ~ dpois(rho[t]) rho[t] <- R[t] * f[t] } } # End Model ",fill = TRUE) sink() # Bundle data bugs.data <- list(nyears = nyears, marray.j = marray.j, marray.a = marray.a, y = popcount, J = J, R = R) # Initial values inits <- function(){list(mphij = runif(1, 0.05, 0.2), mphia = runif(1, 0.2, 0.5), mfec = runif(1, 4, 6), mim = runif(1, 0, 0.3), mp = runif(1, 0.5, 0.8), sig.phij = runif(1, 0.1, 1), sig.phia = runif(1, 0.1, 1), sig.fec = runif(1, 0.1, 1), sig.im = runif(1, 0.1, 1), N1 = round(runif(nyears, 10, 20), 0), NadSurv = round(runif(nyears, 10, 30), 0), Nadimm = round(runif(nyears, 1, 20), 0))} # Parameters monitored parameters <- c("phij", "phia", "f", "omega", "p", "lambda", "mphij", "mphia", "mfec", "mim", "mlam", "sig.phij", "sig.phia", "sig.fec", "sig.im", "N1", "NadSurv", "Nadimm", "Ntot") # MCMC settings niter <- 10000 nthin <- 3 nburn <- 5000 nchains <- 3 # Call WinBUGS from R (BRT 4 min) ipm.hoopoe.2 <- bugs(bugs.data, inits, parameters, "ipm.hoopoe.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Inspect results print(imp.hoopoe.2, 3) ############## # Exercise 3 ############## # Load data # Population count data y <- c(45, 48, 44, 59, 62, 62, 55, 51, 46, 42) # Capture-recapture data (in m-array format) m <- matrix(c(11, 0, 0, 0, 0, 0, 0, 0, 0, 70, 0, 12, 0, 1, 0, 0, 0, 0, 0, 52, 0, 0, 15, 5, 1, 0, 0, 0, 0, 42, 0, 0, 0, 8, 3, 0, 0, 0, 0, 51, 0, 0, 0, 0, 4, 3, 0, 0, 0, 61, 0, 0, 0, 0, 0, 12, 2, 3, 0, 66, 0, 0, 0, 0, 0, 0, 16, 5, 0, 44, 0, 0, 0, 0, 0, 0, 0, 12, 0, 46, 0, 0, 0, 0, 0, 0, 0, 0, 11, 71, 10, 2, 0, 0, 0, 0, 0, 0, 0, 13, 0, 7, 0, 1, 0, 0, 0, 0, 0, 27, 0, 0, 13, 2, 1, 1, 0, 0, 0, 14, 0, 0, 0, 12, 2, 0, 0, 0, 0, 20, 0, 0, 0, 0, 10, 2, 0, 0, 0, 21, 0, 0, 0, 0, 0, 11, 2, 1, 1, 14, 0, 0, 0, 0, 0, 0, 12, 0, 0, 18, 0, 0, 0, 0, 0, 0, 0, 11, 1, 21, 0, 0, 0, 0, 0, 0, 0, 0, 10, 26), ncol = 10, byrow = TRUE) # Data on productivity J <- c(64, 132, 86, 154, 156, 134, 116, 106, 110, 144) R <- c(21, 28, 26, 38, 35, 33, 31, 30, 33, 34) # Data analysis # Specify model in BUGS language sink("ipm.bug") cat(" model { ########################################################################## # Integrated population model # - Age structured model with 2 age classes: # 1-year old and at least 2 years old # - Age at first breeding = 1 year # - Pre-breeding census, female-based # - All vital rates assumed to be constant ########################################################################## ############################################################ # 1. Define the priors for the parameters ############################################################ # Observation error tauy <- pow(sigma.y, -2) sigma.y ~ dunif(0, 50) sigma2.y <- pow(sigma.y, 2) # Initial population sizes N1[1] ~ dnorm(100, 0.0001)I(0,) # 1-year Nad[1] ~ dnorm(100, 0.0001)I(0,) # Adults # Survival and recapture probabilities, as well as productivity for (t in 1:(nyears-1)){ sjuv[t] ~ dunif(0, 1) sad[t] ~ dunif(0, 1) p[t] ~ dunif(0, 1) f[t] ~ dunif(0, 20) } ################################################################### # 2. Derived parameters ################################################################### # Population growth rate for (t in 1:(nyears-1)){ lambda[t] <- Ntot[t+1] / Ntot[t] } ###################################################################### # 3. The likelihoods of the single data sets ###################################################################### ################################################################### # 3.1. Likelihood for population count data (state-space model) ################################################################### ############################# # 3.1.1 System process ############################# for (t in 2:nyears){ mean1[t] <- f[t-1] / 2 * sjuv[t-1] * Ntot[t-1] N1[t] ~ dpois(mean1[t]) Nad[t] ~ dbin(sad[t-1], Ntot[t-1]) } for (t in 1:nyears){ Ntot[t] <- Nad[t] + N1[t] } ############################### # 3.1.2 Observation process ############################### for (t in 1:nyears){ y[t] ~ dnorm(Ntot[t], tauy) } ###################################################################### # 3.2 Likelihood for capture-recapture data: CJS model (2 age classes) ###################################################################### # Multinomial likelihood for (t in 1:2*(nyears-1)){ m[t,1:nyears] ~ dmulti(pr[t,], r[t]) } # Calculate the number of released individuals for (t in 1:2*(nyears-1)){ r[t] <- sum(m[t,]) } # m-array cell probabilities for juveniles for (t in 1:(nyears-1)){ # Main diagonal q[t] <- 1-p[t] pr[t,t] <- sjuv[t] * p[t] # Above main diagonal for (j in (t+1):(nyears-1)){ pr[t,j] <- sjuv[t]*prod(sad[(t+1):j])*prod(q[t:(j-1)])*p[j] } # j # Below main diagonal for (j in 1:(t-1)){ pr[t,j] <- 0 } # j # Last column: probability of non-recapture pr[t,nyears] <- 1-sum(pr[t,1:(nyears-1)]) } # t # m-array cell probabilities for adults for (t in 1:(nyears-1)){ # Main diagonal pr[t+nyears-1,t] <- sad[t] * p[t] # Above main diagonal for (j in (t+1):(nyears-1)){ pr[t+nyears-1,j] <- prod(sad[t:j])*prod(q[t:(j-1)])*p[j] } # Below main diagonal for (j in 1:(t-1)){ pr[t+nyears-1,j] <- 0 } # j # Last column pr[t+nyears-1,nyears] <- 1 - sum(pr[t+nyears-1,1:(nyears-1)]) } # t ############################################################ # 3.3. Likelihood for reproductive data: Poisson regression ############################################################ for (t in 1:(nyears-1)){ J[t] ~ dpois(rho[t]) rho[t] <- R[t]*f[t] } } # End Model ",fill = TRUE) sink() # Bundle data bugs.data <- list(m = m, y = y, J = J, R = R, nyears = dim(m)[2]) # Initial values inits <- function(){list(sjuv = runif(dim(m)[2]-1, 0, 1), sad = runif(dim(m)[2]-1, 0, 1), p = runif(dim(m)[2]-1, 0, 1), f = runif(dim(m)[2]-1, 0, 10), N1 = rpois(dim(m)[2], 30), Nad = rpois(dim(m)[2], 30), sigma.y = runif(1 ,0, 10))} # Parameters monitored parameters <- c("sjuv", "sad", "p", "f", "N1", "Nad", "Ntot", "lambda", "sigma2.y") # MCMC settings niter <- 20000 nthin <- 6 nburn <- 5000 nchains <- 3 # Call WinBUGS from R (BRT 3.5 min) ipm.t <- bugs(bugs.data, inits, parameters, "ipm.bug", n.chains = nchains, n.thin = nthin, n.iter = niter, n.burnin = nburn, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Inspect results print(ipm.t, digits = 3) lower <- upper <- lower.t <- upper.t <- numeric() for (i in 1:10){ lower[i] <- quantile(ipm$sims.list$Ntot[,i], 0.025) upper[i] <- quantile(ipm$sims.list$Ntot[,i], 0.975) lower.t[i] <- quantile(ipm.t$sims.list$Ntot[,i], 0.025) upper.t[i] <- quantile(ipm.t$sims.list$Ntot[,i], 0.975) } plot(ipm$mean$Ntot, type = "b", ylim = c(min(c(lower, lower.t)), max(c(upper, upper.t))), ylab = "Population size", xlab = "Year", las = 1, cex = 1.2, pch = 16) segments(1:10, lower, 1:10, upper) points((1:10)+0.2, ipm.t$mean$Ntot, type = "b", cex = 1.2, pch = 19, col = "blue") segments((1:10)+0.2, lower.t, (1:10)+0.2, upper.t, col = "blue") points((1:10)+0.1, y, type = "p", col = "red", pch = 15, cex = 1.5) legend(x = 1, y = 77, legend = c("estimated constant model", "estimated time-dep. model", "observed"), pch = c(16, 19, 15), col = c("black", "blue", "red"), bty = "n") ############## # Exercise 4 ############## # Load data # Lapwing data taken from Brooks et al. 2004 (Animal Biodiversity and Conservation 27.1: 515-529) # National population index (1965-1998) y = c(NA,NA, 1092.23,1100.01, 1234.32, 1460.85, 1570.38, 1819.79,1391.27,1507.60, 1541.44,1631.21,1628.60,1609.33,1801.68,1809.08,1754.74,1779.48,1699.13, 1681.39,1610.46,1918.45,1717.07,1415.69, 1229.02,1082.02,1096.61,1045.84, 1137.03, 981.1, 647.67, 992.65, 968.62, 926.83, 952.96, 865.64) # Mark-recoveries from individuals marked as hatchings (1963-1997) in m-array format dead.recov <- matrix(c( 13, 4, 1, 2, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1124, 0, 16, 4, 3, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1259, 0, 0, 11, 1, 1, 1, 0, 2, 1, 1, 1, 1, 2, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1082, 0, 0, 0, 10, 4, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1595, 0, 0, 0, 0, 11, 1, 5, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1596, 0, 0, 0, 0, 0, 9, 5, 4, 0, 2, 2, 2, 1, 2, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2091, 0, 0, 0, 0, 0, 0, 11, 9, 4, 3, 1, 1, 1, 3, 2, 2, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1964, 0, 0, 0, 0, 0, 0, 0, 8, 4, 2, 0, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1942, 0, 0, 0, 0, 0, 0, 0, 0, 4, 1, 1, 2, 2, 1, 3, 3, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2444, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 2, 2, 2, 6, 1, 5, 2, 1, 3, 1, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3055, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16, 1, 1, 1, 2, 3, 2, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3412, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13, 4, 4, 7, 3, 1, 1, 1, 1, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3907, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 4, 0, 2, 1, 1, 2, 2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2538, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 3, 5, 1, 3, 3, 2, 3, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3270, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12, 5, 0, 5, 4, 2, 1, 2, 3, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 3443, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15, 5, 2, 2, 0, 5, 3, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3132, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 4, 6, 1, 3, 3, 2, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 3275, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13, 8, 1, 2, 4, 5, 3, 0, 1, 2, 0, 0, 1, 0, 0, 0, 0, 0, 3447, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 23, 2, 2, 3, 3, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3902, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 6, 2, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 2860, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 19, 7, 6, 4, 0, 0, 2, 0, 0, 0, 1, 2, 0, 0, 1, 4077, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12, 3, 2, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 4017, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 25, 2, 5, 2, 0, 2, 2, 2, 0, 0, 0, 0, 0, 4827, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14, 4, 3, 4, 4, 2, 2, 1, 0, 2, 0, 1, 4732, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14, 2, 1, 2, 2, 3, 0, 0, 3, 0, 0, 5000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18, 4, 4, 3, 0, 2, 1, 0, 2, 1, 4769, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 4, 2, 4, 2, 2, 3, 1, 1, 3603, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12, 3, 3, 2, 1, 0, 2, 0, 4147, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, 4, 6, 1, 0, 1, 0, 4293, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 18, 3, 1, 2, 0, 1, 3455, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 5, 2, 2, 1, 3673, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12, 4, 6, 0, 3900, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 5, 1, 3578, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 4481, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 4334), ncol= 36, nrow=35, byrow = T) # Normalised number of frost days from 1963-1997 f = c(0.1922,0.3082,0.3082,-0.9676,0.5401,0.3082,1.1995,0.1921,-0.8526,-1.0835,-0.6196,-1.1995,-0.5037,-0.1557,0.0762,2.628,-0.3877,-0.968,1.9318,-0.6196,-0.3877,1.700, 2.2797,0.6561,-0.8516,-1.0835,-1.0835,0.1922,0.1922,-0.1557,-0.5037,-0.8516,0.8880,-0.0398,-1.1995) # Data analysis # a) Age at first breeding is 2 years # Specify model in BUGS language sink("ipm-lapwing1.bug") cat(" model { #------------------------------------------------- # Integrated population model # - Age structured model with 2 age classes: # 1-year old and adult (at least 2 years old) # - Age at first breeding = 2 years # - Prebreeding census, female-based # - All vital rates assumed to be constant #------------------------------------------------- #------------------------------------------------- # 1. Define the priors for the parameters #------------------------------------------------- # Observation error tauy <- pow(sigma.y, -2) sigma.y ~ dunif(0, 50) sigma2.y <- pow(sigma.y, 2) # Initial population sizes N1[1] ~ dnorm(200, 0.0001)I(0,) # 1-year Nad[1] ~ dnorm(1000, 0.0001)I(0,) # Adults # Survival and recapture probabilities, as well as productivity for (t in 1:(n.occasions+1)){ logit(sj[t]) <- mu.sj + ep.sj[t] logit(sa[t]) <- mu.sa + ep.sa[t] logit(rj[t]) <- mu.rj + ep.rj[t] ra[t] <- rj[t] log(fec[t]) <- mu.fec + ep.fec[t] ep.sj[t] ~ dnorm(0, tau.sj)I(-10,10) ep.sa[t] ~ dnorm(0, tau.sa)I(-10,10) ep.rj[t] ~ dnorm(0, tau.rj)I(-10,10) ep.fec[t] ~ dnorm(0, tau.fec)I(-10,10) } mean.sj ~ dunif(0, 1) mu.sj <- log(mean.sj / (1-mean.sj)) mean.sa ~ dunif(0, 1) mu.sa <- log(mean.sa / (1-mean.sa)) mean.rj ~ dunif(0, 1) mu.rj <- log(mean.rj / (1-mean.rj)) mean.fec ~ dunif(0, 5) mu.fec <- log(mean.fec) sigma.sj ~ dunif(0, 10) tau.sj <- pow(sigma.sj, -2) sigma2.sj <- pow(sigma.sj, 2) sigma.sa ~ dunif(0, 10) tau.sa <- pow(sigma.sa, -2) sigma2.sa <- pow(sigma.sa, 2) sigma.rj ~ dunif(0, 10) tau.rj <- pow(sigma.rj, -2) sigma2.rj <- pow(sigma.rj, 2) sigma.fec ~ dunif(0, 10) tau.fec <- pow(sigma.fec, -2) sigma2.fec <- pow(sigma.fec, 2) #------------------------------------------------- # 2. The likelihoods of the single data sets #------------------------------------------------- # 2.1. Likelihood for population count data (state-space model) # 3.1.1 System process for (t in 2:(n.occasions+1)){ mean1[t] <- fec[t-1] / 2 * sj[t-1] * Nad[t-1] N1[t] ~ dpois(mean1[t]) Nad[t] ~ dbin(sa[t-1], Ntot[t-1]) } for (t in 1:(n.occasions+1)){ Ntot[t] <- Nad[t] + N1[t] # only breeding birds are counted } # 3.1.2 Observation process for (t in 3:(n.occasions+1)){ y[t] ~ dnorm(Nad[t], tauy) } # Define the multinomial likelihoods for (t in 1:n.occasions){ marr.j[t,1:(n.occasions+1)] ~ dmulti(pr.j[t,], rel.j[t]) } # Calculate the number of birds released each year for (t in 1:n.occasions){ rel.j[t] <- sum(marr.j[t,]) } # Define the cell probabilities of the juvenile m-array # Main diagonal for (t in 1:n.occasions){ pr.j[t,t] <- (1-sj[t])*rj[t] # Further above main diagonal for (j in (t+2):n.occasions){ pr.j[t,j] <- sj[t]*prod(sa[(t+1):(j-1)])*(1-sa[j])*ra[j] } #j # Below main diagonal for (j in 1:(t-1)){ pr.j[t,j] <- 0 } #j } #t for (t in 1:(n.occasions-1)){ # One above main diagonal pr.j[t,t+1] <- sj[t]*(1-sa[t+1])*ra[t+1] } #t # Last column: probability of non-recovery for (t in 1:n.occasions){ pr.j[t,n.occasions+1] <- 1-sum(pr.j[t,1:n.occasions]) } #t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr.j = dead.recov, y = y, n.occasions = dim(dead.recov)[2]-1) # last year of census not used, because no demographic data available (but, since we have # Initial values inits <- function(){list(mean.sj = runif(1, 0.4, 0.6), mean.sa = runif(1, 0.7, 0.9), mean.rj = runif(1, 0, 0.2), mean.fec = runif(1, 0, 2), sigma.sj = runif(1, 0, 1), sigma.sa = runif(1, 0, 1), sigma.rj = runif(1, 0, 1), sigma.fec = runif(1, 0, 1), N1 = rpois(36, 400), Nad = rpois(36, 1000), sigma.y = runif(1, 0, 10))} # Parameters monitored parameters <- c("sj", "mean.sj", "sigma2.sj", "sa", "mean.sa", "sigma2.sa", "rj", "mean.rj", "sigma2.rj", "fec", "mean.fec", "sigma2.fec", "sigma2.y", "N1", "Nad", "Ntot") # MCMC settings ni <- 10000 nt <- 3 nb <- 5000 nc <- 3 # Call WinBUGS from R (BRT 362 min) ipm.lapwing1 <- bugs(bugs.data, inits, parameters, "ipm-lapwing1.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) save(ipm.lapwing1a, file="ipm.lapwing1a.Rdata") # Inspect results print(ipm.lapwing1, 3) # Produce graph similar to the one in Fig. 11-7 to visualize results par(mfrow = c(2, 2), cex.axis = 1.2, cex.lab = 1.2, mar = c(5, 6, 1.5, 2), las = 1) lower <- upper <- numeric() year <- 1963:1998 nyears <- length(year) for (i in 1:nyears){ lower[i] <- quantile(ipm.lapwing1$sims.list$Nad[,i], 0.025) upper[i] <- quantile(ipm.lapwing1$sims.list$Nad[,i], 0.975)} m1 <- min(c(ipm.lapwing1$mean$Nad, y, lower), na.rm = T) m2 <- max(c(ipm.lapwing1$mean$Nad, y, upper), na.rm = T) plot(0, 0, ylim = c(0, m2), xlim = c(1, nyears), ylab = "Population size", xlab = " ", col = "black", type = "l", axes = F, frame = F) axis(2) axis(1, at = 1:nyears, labels = year) polygon(x = c(1:nyears, nyears:1), y = c(lower, upper[nyears:1]), col = "grey90", border = "grey90") points(y, type = "l", col = "grey30", lwd = 2) points(ipm.lapwing1$mean$Nad, type = "l", col = "blue", lwd = 2) legend(x = 2, y = 500, legend = c("Counts", "Estimates"), lty = c(1, 1),lwd = c(2, 2), col = c("grey30", "blue"), bty = "n", cex = 1) lower <- upper <- numeric() T <- nyears for (t in 1:T){ lower[t] <- quantile(ipm.lapwing1$sims.list$sj[,t], 0.025) upper[t] <- quantile(ipm.lapwing1$sims.list$sj[,t], 0.975)} par(mgp=c(3.8,1,0)) plot(y = ipm.lapwing1$mean$sj, x = (1:T)+0.5, xlim= c(1, T), type = "b", pch = 16, ylim = c(0.2, 1), ylab = "Annual survival probability", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2) axis(2) axis(1, at = 1:(T+1), labels = 1962:1998) segments((1:T)+0.5, lower, (1:T)+0.5, upper, lwd = 2) segments(1, ipm.lapwing1$mean$mean.sj, T+1, ipm.lapwing1$mean$mean.sj, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.sj, 0.025), T+1, quantile(ipm.lapwing1$sims.list$mean.sj, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.sj, 0.975), T+1, quantile(ipm.lapwing1$sims.list$mean.sj, 0.975), lty = 2, col = "red") for (t in 1:T){ lower[t] <- quantile(ipm.lapwing1$sims.list$sa[,t], 0.025) upper[t] <- quantile(ipm.lapwing1$sims.list$sa[,t], 0.975)} points(y=ipm.lapwing1$mean$sa, x = (1:T)+0.5, type = "b", pch = 1, cex = 1.5, lwd = 2) segments((1:T)+0.5, lower, (1:T)+0.5, upper, lwd = 2) segments(1, ipm.lapwing1$mean$mean.sa, T+1, ipm.lapwing1$mean$mean.sa, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.sa, 0.025), T+1, quantile(ipm.lapwing1$sims.list$mean.sa, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.sa, 0.975), T+1, quantile(ipm.lapwing1$sims.list$mean.sa, 0.975), lty = 2, col = "red") legend(x = 4.5, y = 0.4, legend = c("Adults", "Juveniles"), pch = c(1, 16), bty = "n") lower <- upper <- numeric() T <- nyears for (t in 1:T){ lower[t] <- quantile(ipm.lapwing1$sims.list$fec[,t], 0.025) upper[t] <- quantile(ipm.lapwing1$sims.list$fec[,t], 0.975)} plot(y=ipm.lapwing1$mean$fec, x = (1:T), type = "b", pch = 16, ylim = c(0, 2), ylab = "Fecundity (fledgling / female)", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2) axis(2) axis(1, at = 1:T, labels = 1963:1998) segments((1:T), lower, (1:T), upper) segments(1, ipm.lapwing1$mean$mean.fec, T, ipm.lapwing1$mean$mean.fec, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.fec, 0.025), T, quantile(ipm.lapwing1$sims.list$mean.fec, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.fec, 0.975), T, quantile(ipm.lapwing1$sims.list$mean.fec, 0.975), lty = 2, col = "red") lower <- upper <- numeric() T <- nyears for (t in 1:T){ lower[t] <- quantile(ipm.lapwing1$sims.list$rj[,t], 0.025) upper[t] <- quantile(ipm.lapwing1$sims.list$rj[,t], 0.975)} plot(y=ipm.lapwing1$mean$rj, x = (1:T), type = "b", pch = 16, ylim = c(0, 0.04), ylab = "Recovery probability", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2) axis(2) axis(1, at = 1:T, labels = 1963:1998) segments((1:T), lower, (1:T), upper) segments(1, ipm.lapwing1$mean$mean.rj, T, ipm.lapwing1$mean$mean.rj, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.rj, 0.025), T, quantile(ipm.lapwing1$sims.list$mean.rj, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.lapwing1$sims.list$mean.rj, 0.975), T, quantile(ipm.lapwing1$sims.list$mean.rj, 0.975), lty = 2, col = "red") # b) Age at first breeding is 3 years # Specify model in BUGS language sink("ipm-lapwing2.bug") cat(" model { #------------------------------------------------- # Integrated population model # - Age structured model with 2 age classes: # 1-year old and adult (at least 2 years old) # - Age at first breeding = 3 years # - Prebreeding census, female-based # - All vital rates assumed to be constant #------------------------------------------------- #------------------------------------------------- # 1. Define the priors for the parameters #------------------------------------------------- # Observation error tauy <- pow(sigma.y, -2) sigma.y ~ dunif(0, 50) sigma2.y <- pow(sigma.y, 2) # Initial population sizes N1[1] ~ dnorm(200, 0.001)I(0,) # 1-year N2[1] ~ dnorm(200, 0.001)I(0,) # 1-year Nad[1] ~ dnorm(1000, 0.001)I(0,) # Adults # Survival and recapture probabilities, as well as productivity for (t in 1:(n.occasions+1)){ logit(sj[t]) <- mu.sj + ep.sj[t] logit(sa[t]) <- mu.sa + ep.sa[t] logit(rj[t]) <- mu.rj + ep.rj[t] ra[t] <- rj[t] log(fec[t]) <- mu.fec + ep.fec[t] ep.sj[t] ~ dnorm(0, tau.sj)I(-10,10) ep.sa[t] ~ dnorm(0, tau.sa)I(-10,10) ep.rj[t] ~ dnorm(0, tau.rj)I(-10,10) ep.fec[t] ~ dnorm(0, tau.fec)I(-10,10) } mean.sj ~ dunif(0, 1) mu.sj <- log(mean.sj / (1-mean.sj)) mean.sa ~ dunif(0, 1) mu.sa <- log(mean.sa / (1-mean.sa)) mean.rj ~ dunif(0, 1) mu.rj <- log(mean.rj / (1-mean.rj)) mean.fec ~ dunif(0, 5) mu.fec <- log(mean.fec) sigma.sj ~ dunif(0, 10) tau.sj <- pow(sigma.sj, -2) sigma2.sj <- pow(sigma.sj, 2) sigma.sa ~ dunif(0, 10) tau.sa <- pow(sigma.sa, -2) sigma2.sa <- pow(sigma.sa, 2) sigma.rj ~ dunif(0, 10) tau.rj <- pow(sigma.rj, -2) sigma2.rj <- pow(sigma.rj, 2) sigma.fec ~ dunif(0, 10) tau.fec <- pow(sigma.fec, -2) sigma2.fec <- pow(sigma.fec, 2) #------------------------------------------------- # 2. The likelihoods of the single data sets #------------------------------------------------- # 2.1. Likelihood for population count data (state-space model) # 3.1.1 System process for (t in 2:(n.occasions+1)){ mean1[t] <- fec[t-1] / 2 * sj[t-1] * Nad[t-1] N1[t] ~ dpois(mean1[t]) N2[t] ~ dbin(sa[t-1], N1[t-1]) Nad[t] ~ dbin(sa[t-1], Np[t-1]) } for (t in 1:(n.occasions+1)){ Np[t] <- Nad[t] + N2[t] Ntot[t] <- N1[t] + N2[t] + Nad[t] } # 3.1.2 Observation process for (t in 3:(n.occasions+1)){ y[t] ~ dnorm(Nad[t], tauy) } # Define the multinomial likelihoods for (t in 1:n.occasions){ marr.j[t,1:(n.occasions+1)] ~ dmulti(pr.j[t,], rel.j[t]) } # Calculate the number of birds released each year for (t in 1:n.occasions){ rel.j[t] <- sum(marr.j[t,]) } # Define the cell probabilities of the juvenile m-array # Main diagonal for (t in 1:n.occasions){ pr.j[t,t] <- (1-sj[t])*rj[t] # Further above main diagonal for (j in (t+2):n.occasions){ pr.j[t,j] <- sj[t]*prod(sa[(t+1):(j-1)])*(1-sa[j])*ra[j] } #j # Below main diagonal for (j in 1:(t-1)){ pr.j[t,j] <- 0 } #j } #t for (t in 1:(n.occasions-1)){ # One above main diagonal pr.j[t,t+1] <- sj[t]*(1-sa[t+1])*ra[t+1] } #t # Last column: probability of non-recovery for (t in 1:n.occasions){ pr.j[t,n.occasions+1] <- 1-sum(pr.j[t,1:n.occasions]) } #t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr.j = dead.recov, y = y, n.occasions = dim(dead.recov)[2]-1) # Initial values inits <- function(){list(mean.sj = runif(1, 0.4, 0.6), mean.sa = runif(1, 0.7, 0.9), mean.rj = runif(1, 0, 0.2), mean.fec = runif(1, 0, 2), sigma.sj = runif(1, 0, 1), sigma.sa = runif(1, 0, 1), sigma.rj = runif(1, 0, 1), sigma.fec = runif(1, 0, 1), N1 = rpois(36, 200), N2 = rpois(36, 100), Nad = rpois(36, 1000), sigma.y = runif(1, 0, 10))} # Parameters monitored parameters <- c("sj", "mean.sj", "sigma2.sj", "sa", "mean.sa", "sigma2.sa", "rj", "mean.rj", "sigma2.rj", "fec", "mean.fec", "sigma2.fec", "sigma2.y", "N1", "N2", "Nad", "Ntot") # MCMC settings ni <- 10000 nt <- 3 nb <- 5000 nc <- 3 # Call WinBUGS from R (BRT 420 min) ipm.lapwing2 <- bugs(bugs.data, inits, parameters, "ipm-lapwing2.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(ipm.lapwing2, 3) # c) Survival is a function of the number of frost days # Specify model in BUGS language sink("ipm-lapwing3.bug") cat(" model { #------------------------------------------------- # Integrated population model # - Age structured model with 2 age classes: # 1-year old and adult (at least 2 years old) # - Age at first breeding = 2 years # - Prebreeding census, female-based # - All vital rates assumed to be constant #------------------------------------------------- #------------------------------------------------- # 1. Define the priors for the parameters #------------------------------------------------- # Observation error tauy <- pow(sigma.y, -2) sigma.y ~ dunif(0, 50) sigma2.y <- pow(sigma.y, 2) # Initial population sizes N1[1] ~ dnorm(200, 0.0001)I(0,) # 1-year Nad[1] ~ dnorm(1000, 0.0001)I(0,) # Adults # Survival and recapture probabilities, as well as productivity for (t in 1:n.occasions){ logit(sj[t]) <- mu.sj + beta*frost[t] + ep.sj[t] logit(sa[t]) <- mu.sa + beta*frost[t] + ep.sa[t] logit(rj[t]) <- mu.rj + ep.rj[t] ra[t] <- rj[t] log(fec[t]) <- mu.fec + ep.fec[t] ep.sj[t] ~ dnorm(0, tau.sj)I(-10,10) ep.sa[t] ~ dnorm(0, tau.sa)I(-10,10) ep.rj[t] ~ dnorm(0, tau.rj)I(-10,10) ep.fec[t] ~ dnorm(0, tau.fec)I(-10,10) } mean.sj ~ dunif(0, 1) mu.sj <- log(mean.sj / (1-mean.sj)) mean.sa ~ dunif(0, 1) mu.sa <- log(mean.sa / (1-mean.sa)) mean.rj ~ dunif(0, 1) mu.rj <- log(mean.rj / (1-mean.rj)) mean.fec ~ dunif(0, 5) mu.fec <- log(mean.fec) beta ~ dnorm(0, 0.001) sigma.sj ~ dunif(0, 10) tau.sj <- pow(sigma.sj, -2) sigma2.sj <- pow(sigma.sj, 2) sigma.sa ~ dunif(0, 10) tau.sa <- pow(sigma.sa, -2) sigma2.sa <- pow(sigma.sa, 2) sigma.rj ~ dunif(0, 10) tau.rj <- pow(sigma.rj, -2) sigma2.rj <- pow(sigma.rj, 2) sigma.fec ~ dunif(0, 10) tau.fec <- pow(sigma.fec, -2) sigma2.fec <- pow(sigma.fec, 2) #------------------------------------------------- # 2. The likelihoods of the single data sets #------------------------------------------------- # 2.1. Likelihood for population count data (state-space model) # 3.1.1 System process for (t in 2:n.occasions){ mean1[t] <- fec[t-1] / 2 * sj[t-1] * Nad[t-1] N1[t] ~ dpois(mean1[t]) Nad[t] ~ dbin(sa[t-1], Ntot[t-1]) } for (t in 1:n.occasions){ Ntot[t] <- Nad[t] + N1[t] # only breeding birds are counted } # 3.1.2 Observation process for (t in 3:n.occasions){ y[t] ~ dnorm(Nad[t], tauy) } # Define the multinomial likelihoods for (t in 1:n.occasions){ marr.j[t,1:(n.occasions+1)] ~ dmulti(pr.j[t,], rel.j[t]) } # Calculate the number of birds released each year for (t in 1:n.occasions){ rel.j[t] <- sum(marr.j[t,]) } # Define the cell probabilities of the juvenile m-array # Main diagonal for (t in 1:n.occasions){ pr.j[t,t] <- (1-sj[t])*rj[t] # Further above main diagonal for (j in (t+2):n.occasions){ pr.j[t,j] <- sj[t]*prod(sa[(t+1):(j-1)])*(1-sa[j])*ra[j] } #j # Below main diagonal for (j in 1:(t-1)){ pr.j[t,j] <- 0 } #j } #t for (t in 1:(n.occasions-1)){ # One above main diagonal pr.j[t,t+1] <- sj[t]*(1-sa[t+1])*ra[t+1] } #t # Last column: probability of non-recovery for (t in 1:n.occasions){ pr.j[t,n.occasions+1] <- 1-sum(pr.j[t,1:n.occasions]) } #t } ",fill = TRUE) sink() # Bundle data bugs.data <- list(marr.j = dead.recov, y = y[-36], n.occasions = dim(dead.recov)[2]-1, frost = f) # last year of census not used # Initial values inits <- function(){list(mean.sj = runif(1, 0.4, 0.6), mean.sa = runif(1, 0.7, 0.9), mean.rj = runif(1, 0, 0.2), mean.fec = runif(1, 0, 2), sigma.sj = runif(1, 0, 1), sigma.sa = runif(1, 0, 1), sigma.rj = runif(1, 0, 1), sigma.fec = runif(1, 0, 1), N1 = rpois(35, 400), Nad = rpois(35, 1000), sigma.y = runif(1, 0, 10), beta = rnorm(1))} # Parameters monitored parameters <- c("sj", "mean.sj", "sigma2.sj", "sa", "mean.sa", "sigma2.sa", "rj", "mean.rj", "sigma2.rj", "fec", "mean.fec", "sigma2.fec", "sigma2.y", "N1", "Nad", "Ntot", "beta") # MCMC settings ni <- 20 nt <- 1 nb <- 5 nc <- 3 # Call WinBUGS from R (BRT 808 min) ipm.lapwing3 <- bugs(bugs.data, inits, parameters, "ipm-lapwing3.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(ipm.lapwing3, 3) plot(density(ipm.lapwing3$sims.list$beta)) sum(ipm.lapwing3$sims.list$beta<0)/length(ipm.lapwing3$sims.list$beta) ############################################################################# # # Chapter 12 # ############################################################################# ############## # Exercise 1 ############## # Bundle data, including the new sampling covariate X2 y <- data$y X2 <- matrix(rnorm(length(y)), ncol = ncol(y)) win.data <- list(y = y, R = nrow(y), T = ncol(y), X = data$X, X2 = X2) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors alpha0 ~ dunif(-10, 10) alpha1 ~ dunif(-10, 10) beta0 ~ dunif(-10, 10) beta1 ~ dunif(-10, 10) beta2 ~ dunif(-10, 10) # Likelihood # Ecological model for true abundance for (i in 1:R){ N[i] ~ dpois(lambda[i]) log(lambda[i]) <- alpha0 + alpha1 * X[i] # Observation model for replicated counts for (j in 1:T){ y[i,j] ~ dbin(p[i,j], N[i]) p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j])) lp[i,j] <- beta0 + beta1 * X[i] + beta2 * X2[i,j] } #j } #i # Derived quantities totalN <- sum(N[]) } ",fill = TRUE) sink() # Initial values Nst <- apply(y, 1, max) + 1 # Important to give good inits for latent N inits <- function() list(N = Nst, alpha0 = runif(1, -1, 1), alpha1 = runif(1, -1, 1), beta0 = runif(1, -1, 1), beta1 = runif(1, -1, 1), beta2 = runif(1, -1, 1)) # Parameters monitored params <- c("totalN", "alpha0", "alpha1", "beta0", "beta1", "beta2") # MCMC settings ni <- 22000 nt <- 20 nb <- 2000 nc <- 3 # Call WinBUGS from R out1 <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out, dig = 2) library(unmarked) umf <- unmarkedFramePCount(y = y, siteCovs = data.frame(X = data$X), obsCovs = list(X2 = X2)) summary(umf) ############## # Exercise 2 ############## # Specify model in BUGS language sink("Nmix2A.txt") cat(" model { # Priors for (k in 1:7){ alpha.lam[k] ~ dnorm(0, 0.1) beta[k] ~ dnorm(0, 0.1) } # Abundance site and detection site-by-day random effects for (i in 1:R){ eps[i] ~ dnorm(0, tau.lam) } tau.lam <- 1 / (sd.lam * sd.lam) sd.lam ~ dunif(0, 3) tau.p <- 1 / (sd.p * sd.p) sd.p ~ dunif(0, 3) # Likelihood # Ecological model for true abundance for (i in 1:R){ # Loop for R sites (95) for (k in 1:7){ # Loop for days (7) N[i,k] ~ dpois(lambda[i,k]) # Abundance log(lambda[i,k]) <- alpha.lam[k] + eps[i] # Observation model for replicated counts lp[i,k] ~ dnorm(beta[k], tau.p) # random delta def. implicitly p[i,k] <- 1 / (1 + exp(-lp[i,k])) for (j in 1:T){ # Loop for temporal reps (2) y[i,j,k] ~ dbin(p[i,k], N[i,k]) # Detection # Assess model fit using Chi-squared discrepancy # Compute fit statistic for observed data eval[i,j,k] <- p[i,k] * N[i,k] E[i,j,k] <- pow((y[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k]+0.5) # Generate replicate data and compute fit stats for them y.new[i,j,k] ~ dbin(p[i,k], N[i,k]) E.new[i,j,k] <- pow((y.new[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k]+0.5) } #j ik.p[i,k] <- mean(p[i,k]) } #k } #i # Derived and other quantities for (k in 1:7){ totalN[k] <- sum(N[,k]) # Estimate total pop. size across all sites mean.abundance[k] <- mean(lambda[,k]) mean.N[k] <- mean(N[,k]) mean.detection[k] <- mean(ik.p[,k]) } fit <- sum(E[,,]) fit.new <- sum(E.new[,,]) } ",fill = TRUE) sink() # Bundle data R = nrow(y) T = ncol(y) win.data <- list(y = y, R = R, T = T) # Initial values Nst <- apply(y, c(1, 3), max) + 1 Nst[is.na(Nst)] <- 1 inits <- function(){list(N = Nst, alpha.lam = runif(7, -1, 1), beta = runif(7, -1, 1), sd.lam = runif(1, 0, 1), sd.p = runif(1, 0, 1))} # Parameters monitored params <- c("totalN", "alpha.lam", "beta", "sd.lam", "sd.p", "mean.abundance", "mean.N", "mean.detection", "fit", "fit.new") # MCMC settings ni <- 350000 nt <- 300 nb <- 50000 nc <- 3 # Call WinBUGS from R out2A <- bugs(win.data, inits, params, "Nmix2A.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = FALSE, bugs.directory = bugs.dir, working.directory = getwd()) # Evaluation of fit plot(out2A$sims.list$fit, out2A$sims.list$fit.new, main = "", xlab = "Discrepancy actual data", ylab = "Discrepancy replicate data", frame.plot = FALSE, xlim = c(100, 300), ylim = c(100, 300)) abline(0, 1, lwd = 2, col = "black") mean(out2A$sims.list$fit.new > out2A$sims.list$fit) # Summarize posteriors print(out2A, dig = 2) ############## # Exercise 3 ############## # Specify model in BUGS language sink("Nmix2B.txt") cat(" model { # Priors for (k in 1:7){ alpha.lam[k] ~ dnorm(0, 0.1) beta[k] ~ dnorm(0, 0.1) } # Abundance and detection site effects # and detection site-by-day-by-rep random effects for (i in 1:R){ eps.lam[i] ~ dnorm(0, tau.lam) # Random site effects on abundance eps.p[i] ~ dnorm(0, tau.site.p) # Random site effects on detection } tau.lam <- pow(sd.lam, -2) sd.lam ~ dunif(0, 3) tau.site.p <- pow(sd.site.p, -2) sd.site.p ~ dunif(0, 3) tau.p <- pow(sd.p, -2) sd.p ~ dunif(0, 3) # Likelihood # Ecological model for true abundance for (i in 1:R){ # Loop over sites (95) for (k in 1:7){ # Loop over days (7) N[i,k] ~ dpois(lambda[i,k]) # Abundance log(lambda[i,k]) <- alpha.lam[k] + eps.lam[i] # Observation model for replicated counts mu.lp[i,k] <- beta[k] + eps.p[i] for (j in 1:T){ # Loop over reps (2) y[i,j,k] ~ dbin(p[i,j,k], N[i,k]) # Detection p[i,j,k] <- 1 / (1 + exp(-lp[i,j,k])) lp[i,j,k] ~ dnorm(mu.lp[i,k], tau.p) # random delta defined implicitly # Assess model fit using Chi-squared discrepancy # Compute fit statistic for observed data eval[i,j,k] <- p[i,j,k] * N[i,k] E[i,j,k] <- pow((y[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k]+0.5) # Generate replicate data and compute fit stats for them y.new[i,j,k] ~ dbin(p[i,j,k], N[i,k]) E.new[i,j,k] <- pow((y.new[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k]+0.5) } #j ik.p[i,k] <- mean(p[i,,k]) } #k } #i # Derived and other quantities for (k in 1:7){ totalN[k] <- sum(N[,k]) # Estimate total pop. size across all sites mean.abundance[k] <- mean(lambda[,k]) mean.N[k] <- mean(N[,k]) mean.detection[k] <- mean(ik.p[,k]) } fit <- sum(E[,,]) fit.new <- sum(E.new[,,]) } ",fill = TRUE) sink() # Bundle data R = nrow(y) T = ncol(y) win.data <- list(y = y, R = R, T = T) # Initial values Nst <- apply(y, c(1, 3), max) + 1 Nst[is.na(Nst)] <- 1 inits <- function(){list(N = Nst, alpha.lam = runif(7, -3, 3), beta = runif(7, -3, 3), sd.lam = runif(1, 0, 1), sd.p = runif(1, 0, 1))} # Parameters monitored params <- c("totalN", "alpha.lam", "beta", "sd.lam", "sd.site.p", "sd.p", "mean.abundance", "mean.N", "mean.detection", "fit", "fit.new") # MCMC settings ni <- 350000 nt <- 300 nb <- 50000 nc <- 3 # Call WinBUGS from R out2B <- bugs(win.data, inits, params, "Nmix2B.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Evaluation of fit plot(out2B$sims.list$fit, out2B$sims.list$fit.new, main = "", xlab = "Discrepancy actual data", ylab = "Discrepancy replicate data", frame.plot = FALSE, xlim = c(50, 200), ylim = c(50, 200)) abline(0, 1, lwd = 2, col = "black") mean(out2B$sims.list$fit.new > out2B$sims.list$fit) mean(out2B$mean$fit) / mean(out2B$mean$fit.new) # Summarize posteriors print(out2B, dig = 2) max.day.count <- apply(y, c(1, 3), max, na.rm = TRUE) max.day.count[max.day.count == "-Inf"] <- NA mean.max.count <- apply(max.day.count, 2, mean, na.rm = TRUE) mean.max.count jit <- 0.1 # degree of translation of x-axis par(mfrow = c(3, 1)) plot(1:7, mean.max.count, xlab = "Day", ylab = "Mean daily abundance", las = 1, ylim = c(0, 2), type = "b", main = "", frame.plot = FALSE, pch = 16, lwd = 2) plot(1:7-jit, out2$summary[24:30,5], type = "b", pch = 16, col = "blue", lwd = 2, xlab = "Day", ylab = "Mean daily abundance", ylim = c(0, 32), main = "", frame.plot = FALSE) segments(1:7-jit, out2$summary[24:30,3], 1:7-jit, out2$summary[24:30,7], col = "blue") lines(1:7, out2A$summary[24:30,5], type = "b", pch = 16, col = "green", lwd = 2) segments(1:7, out2A$summary[24:30,3], 1:7, out2A$summary[24:30,7], col = "green") lines(1:7+jit, out2B$summary[25:31,5], type = "b", pch = 16, col = "red", lwd = 2) segments(1:7+jit, out2B$summary[25:31,5], 1:7+jit, out2B$summary[24:30,7], col = "red") plot(1:7-jit, out2$summary[38:44,5], xlab = "Day", ylab = "Detection probability ", las = 1, ylim = c(0, 1), type = "b", col = "blue", pch = 16, frame.plot = FALSE, lwd = 2) segments(1:7-jit, out2$summary[38:44,3], 1:7-jit, out2$summary[38:44,7], col = "blue") lines(1:7, out2A$summary[38:44,5], type = "b", col = "green", pch = 16, lwd = 2) segments(1:7, out2A$summary[38:44,3], 1:7, out2A$summary[38:44,7], col = "green") lines(1:7+jit, out2B$summary[39:45,5], type = "b", col = "red", pch = 16, lwd = 2) segments(1:7+jit, out2B$summary[39:45,3], 1:7+jit, out2B$summary[38:44,7], col = "red") par(mfrow = c(2, 1)) hist(out2B$sims.list$sd.site.p, breaks = 100, col = "gray") hist(out2B$sims.list$sd.p, breaks = 100, col = "gray") ############################################################################# # # Chapter 13 # ############################################################################# ############## # Exercise 1 ############## # Generate a ‘seen-before’ covariate sb <- array(NA, dim = dim(y)) for (i in 1:27){ for (j in 1:6){ sb[i,j] <- max(y[i, 1:(j-1)]) } } sb[is.na(y)] <- 0 # Impute ‘irrelevant’ zeroes library(unmarked) bugdata <- unmarkedFrameOccu(y = y, siteCovs = data.frame(edge = edge), obsCovs = list(DATES = DATES, HOURS = HOURS)) fm <- occu(formula = ~ DATES + I(DATES^2) + HOURS + I(HOURS^2) ~ as.factor(edge)-1, data = bugdata) ML.results <- rbind(summary(fm)$state[,1:2], summary(fm)$det[,1:2]) Bayesian.results <- out$summary[c(1:2, 5:9), c(1:3, 7)] print(cbind(ML.results, Bayesian.results)) tapply(apply(y, 1, max, na.rm = TRUE), edge, mean) # Call JAGS from R and get run time library("R2jags") # requires rjags outJAGS <- jags(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb) print(outJAGS, dig = 3) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors alpha.in ~ dnorm(0, 0.01) alpha.edge ~ dnorm(0, 0.01) alpha.p ~ dnorm(0, 0.01) beta1.p ~ dnorm(0, 0.01) beta2.p ~ dnorm(0, 0.01) beta3.p ~ dnorm(0, 0.01) beta4.p ~ dnorm(0, 0.01) # Likelihood # Ecological model for the partially observed true state for (i in 1:R){ z[i] ~ dbern(psi[i]) # True occurrence z at site i psi[i] <- 1 / (1 + exp(-lpsi.lim[i])) lpsi.lim[i] <- min(999, max(-999, lpsi[i])) lpsi[i] <- alpha.in * (1-edge[i]) + alpha.edge * edge[i] # Observation model for the observations for (j in 1:T){ y[i,j] ~ dbern(mu.p[i,j]) # Detection-nondetection at i and j mu.p[i,j] <- z[i] * p[i,j] p[i,j] <- 1 / (1 + exp(-lp.lim[i,j])) lp.lim[i,j] <- min(999, max(-999, lp[i,j])) lp[i,j] <- alpha.p + beta1.p * DATES[i,j] + beta2.p * pow(DATES[i,j], 2) + beta3.p * HOURS[i,j] + beta4.p * pow(HOURS[i,j], 2) } #j } #i # Derived quantities occ.fs <- sum(z[]) # Number of occupied sites mean.p <- exp(alpha.p) / (1 + exp(alpha.p)) # Sort of average detection } ",fill = TRUE) sink() # Bundle data win.data <- list(y = y, R = nrow(y), T = ncol(y), edge = edge, DATES = DATES, HOURS = HOURS) # Initial values zst <- apply(y, 1, max, na.rm = TRUE) # Good starting values crucial inits <- function(){list(z = zst, alpha.p = runif(1, -3, 3))} # Parameters monitored params <- c("alpha.in", "alpha.edge", "mean.p", "occ.fs", "alpha.p", "beta1.p", "beta2.p", "beta3.p", "beta4.p") # MCMC settings ni <- 30000 nt <- 10 nb <- 20000 nc <- 3 # Call WinBUGS from R (BRT < 1 min) out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) print(out, 2) # Specify model in BUGS language sink("model.txt") cat(" model { # Priors alpha.psi ~ dnorm(0, 0.01) # beta.psi ~ dnorm(0, 0.01) # Drop that term alpha.p ~ dnorm(0, 0.01) beta1.p ~ dnorm(0, 0.01) beta2.p ~ dnorm(0, 0.01) beta3.p ~ dnorm(0, 0.01) beta4.p ~ dnorm(0, 0.01) beta5.p ~ dnorm(0, 0.01) # Likelihood # Ecological model for the partially observed true state for (i in 1:R){ z[i] ~ dbern(psi[i]) # True occurrence z at site i psi[i] <- 1 / (1 + exp(-lpsi.lim[i])) lpsi.lim[i] <- min(999, max(-999, lpsi[i])) lpsi[i] <- alpha.psi ### + beta.psi * edge[i] # Drop beta.psi # Observation model for the observations for (j in 1:T){ y[i,j] ~ dbern(mu.p[i,j]) # Detection-nondetection at i and j mu.p[i,j] <- z[i] * p[i,j] p[i,j] <- 1 / (1 + exp(-lp.lim[i,j])) lp.lim[i,j] <- min(999, max(-999, lp[i,j])) lp[i,j] <- alpha.p + beta1.p*DATES[i,j] + beta2.p*pow(DATES[i,j], 2) + beta3.p*HOURS[i,j] + beta4.p*pow(HOURS[i,j], 2) + beta5.p*sb[i,j] } #j } #i # Derived quantities occ.fs <- sum(z[]) # Number of occupied sites mean.p <- exp(alpha.p) / (1 + exp(alpha.p)) # Average first detection } ",fill = TRUE) sink() # Bundle data win.data <- list(y = y, R = nrow(y), T = ncol(y), DATES = DATES, HOURS = HOURS, sb = sb) # Initial values zst <- apply(y, 1, max, na.rm = TRUE) # Good starting values crucial inits <- function(){list(z = zst, alpha.psi=runif(1, -3, 3), alpha.p = runif(1, -3, 3))} # Parameters monitored params <- c("alpha.psi", "mean.p", "occ.fs", "alpha.p", "beta1.p", "beta2.p", "beta3.p", "beta4.p", "beta5.p") # MCMC settings ni <- 25000 nt <- 10 nb <- 5000 nc <- 3 # Call WinBUGS from R out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out, dig = 2) # Plot posterior and prior for alpha.psi in same graph hist(out$sims.list$alpha.psi, breaks = 100, col = "grey", freq = FALSE) lines(dnorm(0:30, mean = 0, sd = sqrt(1 / (0.01))), col = "red", lwd =3) library(unmarked) bugdata <- unmarkedFrameOccu(y = y, siteCovs = data.frame(edge = edge), obsCovs = list(DATES = DATES, HOURS = HOURS, sb = sb)) summary(bugdata) summary(fm <- occu(formula = ~ DATES + I(DATES^2) + HOURS + I(HOURS^2) + sb ~ 1, data = bugdata)) summary(fm <- occu(formula = ~ sb ~ 1, data = bugdata)) summary(fm <- occu(formula = ~ ~ DATES + I(DATES^2) + HOURS + I(HOURS^2) ~ 1, data = bugdata)) ############## # Exercise 2 ############## data <- data.fn(R = 250, J = 3, K = 10, psi1 = 0.6, range.p = c(0.1, 0.9), range.phi = c(0.7, 0.9), range.gamma = c(0.1, 0.5)) attach(data) zobs <- apply(y, c(1, 3), max) # Observed occurrence as inits for z win.data <- list(y = zobs, nsite = dim(zobs)[1], nyear = dim(zobs)[2]) # Specify model in BUGS language sink("Naive.Dynocc.txt") cat(" model { # Specify priors psi1 ~ dunif(0, 1) for (k in 1:(nyear-1)){ phi[k] ~ dunif(0, 1) gamma[k] ~ dunif(0, 1) } # Combined model: No separation of state and observation processes for (i in 1:nsite){ y[i,1] ~ dbern(psi1) for (k in 2:nyear){ muZ[i,k]<- y[i,k-1]*phi[k-1] + (1-y[i,k-1])*gamma[k-1] y[i,k] ~ dbern(muZ[i,k]) } #k } #i # Derived parameters: Population occupancy, growth rate and turnover psi[1] <- psi1 for (k in 2:nyear){ psi[k] <- psi[k-1]*phi[k-1] + (1-psi[k-1])*gamma[k-1] growthr[k] <- psi[k]/psi[k-1] turnover[k-1] <- (1 - psi[k-1]) * gamma[k-1]/psi[k] } } ",fill = TRUE) sink() # Initial values inits <- function(){ list(psi1 = runif(1, 0, 1))} # Parameters monitored params <- c("psi", "phi", "gamma", "growthr", "turnover") # MCMC settings ni <- 1100 ; nt <- 1 ; nb <- 100 ; nc <- 3 # Call WinBUGS from R out <- bugs(win.data, inits, params, "Naive.Dynocc.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = FALSE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out, dig = 2) lines(1:data$K, out$mean$psi, type = "l", col = "green", lwd = 2, lty = "dashed") TO <- GR <- array(dim = 9) for (k in 2:data$K){ GR[k-1] <- data$psi[k] / data$psi[k-1] TO[k-1] <- (1 - data$psi[k-1]) * data$gamma[k-1] / data$psi[k] } # Compare estimates of dynamic parameters with their true values par(mfrow = c(2, 2), mar = c(5, 5, 2, 1)) # Survival probability plot(1:(data$K-1), data$phi, type = "l", xlab = "Yearly Interval", ylab = "Probability", col = "red", xlim = c(0, data$K), ylim = c(0, 1), lwd = 2, lty = 1, frame.plot = FALSE, las = 1, main = "Survival probability") lines(1:(data$K-1), out$mean$phi, type = "l", col = "blue", lwd = 2, lty = 1) segments(1:(data$K-1), out$summary[11:19,3], 1:(data$K-1), out$summary[11:19,7], lwd = 2, col = "blue") # legend(4, 0.15, c('Truth', 'Estimate naïve occupancy model'), col=c("red", "blue"), lty = c(1, 1), lwd = 2, cex = 0.8) # Colonization probability plot(1:(data$K-1), data$gamma, type = "l", xlab = "Yearly Interval", ylab = "Probability", col = "red", xlim = c(0, data$K), ylim = c(0, 1), lwd = 2, lty = 1, frame.plot = FALSE, las = 1, main = "Colonization probability") lines(1:(data$K-1), out$mean$gamma, type = "l", col = "blue", lwd = 2, lty = 1) segments(1:(data$K-1), out$summary[20:28,3], 1:(data$K-1), out$summary[20:28,7], lwd = 2, col = "blue") # legend(0, 0.9, c('Truth', 'Estimate naïve occupancy model'), col=c("red", "blue"), lty = c(1, 1), lwd = 2, cex = 0.8) # Growth rate plot(1:(data$K-1), GR, type = "l", xlab = "Yearly Interval", ylab = "Rate", col = "red", xlim = c(0, data$K), ylim = c(0, 5), lwd = 2, lty = 1, frame.plot = FALSE, las = 1, main = "Occupancy growth rate") lines(1:(data$K-1), out$mean$growthr, type = "l", col = "blue", lwd = 2, lty = 1) segments(1:(data$K-1), out$summary[29:37,3], 1:(data$K-1), out$summary[29:37,7], lwd = 2, col = "blue") abline(h = 1, lwd = 1, col = "black") # legend(0, 5, c('Truth', 'Estimate naïve occupancy model'), col=c("red", "blue"), lty = c(1, 1), lwd = 2, cex = 0.8) # Turnover rate plot(1:(data$K-1), TO, type = "l", xlab = "Yearly Interval", ylab = "Probability", col = "red", xlim = c(0, data$K), ylim = c(0, 1), lwd = 2, lty = 1, frame.plot = FALSE, las = 1, main = "Turnover rate") lines(1:(data$K-1), out$mean$turnover, type = "l", col = "blue", lwd = 2, lty = 1) segments(1:(data$K-1), out$summary[38:46,3], 1:(data$K-1), out$summary[38:46,7], lwd = 2, col = "blue") # legend(0, 1, c('Truth', 'Estimate naïve occupancy model'), col=c("red", "blue"), lty = c(1, 1), lwd = 2, cex = 0.8) ############## # Exercise 3 ############## # Specify model in BUGS language sink("ImplicitDynocc.txt") cat(" model { # Specify priors for (k in 1:nyear){ psi[k] ~ dunif(0, 1) p[k] ~ dunif(0, 1) } # Ecological submodel: Define state conditional on parameters for (i in 1:nsite){ for (k in 1:nyear){ z[i,k] ~ dbern(psi[k]) } #k } #i # Observation model: Define observation conditional on state for (i in 1:nsite){ for (j in 1:nrep){ for (k in 1:nyear){ muy[i,j,k] <- z[i,k]*p[k] y[i,j,k] ~ dbern(muy[i,j,k]) } #k } #j } #i # Derived parameters: Sample occupancy and growth rate n.occ[1] <- sum(z[1:nsite,1]) for (k in 2:nyear){ n.occ[k] <- sum(z[1:nsite,k]) growthr[k] <- psi[k]/psi[k-1] } } ",fill = TRUE) sink() # Bundle data win.data <- list(y = y, nsite = dim(y)[1], nrep = dim(y)[2], nyear = dim(y)[3]) # Initial values zst <- apply(y, c(1, 3), max) # Observed occurrence as inits for z inits <- function(){ list(z = zst)} # Parameters monitored params <- c("psi", "p", "n.occ", "growthr") # MCMC settings ni <- 2500 nt <- 4 nb <- 500 nc <- 3 # Call WinBUGS from R out <- bugs(win.data, inits, params, "ImplicitDynocc.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out, dig = 2) ############## # Exercise 4 ############## data.fn <- function(R = 250, T = 3, psi = 0.2, p = 0.1){ y <- matrix(NA, nrow = R, ncol = T) z <- rbinom(n = R, size = 1, prob = psi) # Latent occurrence state for (j in 1:T){ y[,j] <- rbinom(n = R, size = 1, prob = z * p) } n.occ <- sum(z) n.occ.obs <- sum(apply(y, 1, max)) return(list(R=R, T=T, psi=psi, p=p, z=z, y=y, n.occ=n.occ, n.occ.obs=n.occ.obs)) } str(data <- data.fn()) # Define a function to do data augmentation, run analysis using Mh and return results all at once model.fn <- function(ni = 1200, nt = 2, nb = 200, nc = 3, data.file = data, debg = FALSE){ # Function arguments: # ni/nt/nb/nc -- MCMC settings # data.file -- name of the object with the simulated data # debg -- setting of DEBUG argument in bugs() # Specify model in BUGS language sink("model.txt") cat(" model { psi ~ dunif(0, 1) p ~ dunif(0, 1) for (i in 1:R){ z[i] ~ dbern(psi) p.eff[i] <- z[i] * p for (j in 1:T){ y[i,j] ~ dbern(p.eff[i]) } #j } #i occ.fs <- sum(z[]) } ",fill = TRUE) sink() # Bundle data win.data <- list(y = data$y, R = data$R, T = data$T) # Initial values zst <- apply(data$y, 1, max) # Observed occurrence as starting values for z inits <- function() list(z = zst) # Parameters monitored params <- c("psi", "p", "occ.fs") # Call WinBUGS from R out <- bugs(win.data, inits, params, "model.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug = debg, bugs.directory = bugs.dir, working.directory = getwd()) # Return stuff return(post.estimates = out$summary) } data <- data.fn(R = 250, T = 3, psi = 0.2, p = 0.1) estimates <- model.fn(ni = 2500, nt = 2, nb = 500, nc = 3, data.file = data, debg = TRUE) # Set up data structures to hold the results simreps <- 100 data.sets <- array(NA, dim = c(250, 3, simreps)) solutions <- array(NA, dim = c(4, 9, simreps)) rownames(solutions) <- rownames(estimates) colnames(solutions) <- colnames(estimates) # Run data generation/data analysis cycle simrep times for (i in 1:simreps){ cat(paste("\n\n*** SimRep", i, "***\n")) data <- data.fn(R = 250, T = 3, psi = 0.2, p = 0.1) data.sets[,,i] <- data$y estimates <- model.fn(ni = 2500, nt = 2, nb = 500, nc = 3, data.file = data, debg = FALSE) solutions[,,i] <- estimates } # Summarize simulation results par(mfrow = c(2, 1)) hist(solutions[1,1,], breaks = 25, col = "grey", xlab = "Estimates of psi", main = "psi: posterior mean", xlim = c(0, 1), las = 1) abline(v = 0.2, col = "red", lwd = 3) abline(v = mean(solutions[1,1,]), col = "blue", lwd = 3) hist(solutions[2,1,], breaks = 25, col = "grey", xlab = "Estimates of p", main = "p: posterior mean", xlim = c(0, 1), las = 1) abline(v = 0.1, col = "red", lwd = 3) abline(v = mean(solutions[2,1,]), col = "blue", lwd = 3) ############## # Exercise 5 ############## # Bundle data y <- as.matrix(owls[, 2:6]) y <- y + 1 DATE <- as.matrix(owls[, 7:11]) mn.date <- mean(DATE, na.rm = TRUE) sd.date <- sd(c(DATE), na.rm = TRUE) DATE <- (DATE - mn.date)/ sd.date DATE[is.na(DATE)]<- 0 win.data <- list(y = y, DATE = DATE, R = dim(y)[1], T = dim(y)[2]) # Specifiy model in BUGS language sink("model3.bug") cat(" model { # Priors psi ~ dunif(0, 1) r ~ dunif(0,1 ) for (t in 1:T){ for (s in 1:R){ # linear models on logit and multinomial logit scale logit(p2[s,t]) <- int.p2 + beta.p2 * DATE[s,t] lp3[2,s,t] <- int.p32 + beta.p32 * DATE[s,t] lp3[3,s,t] <- int.p33 + beta.p33 * DATE[s,t] p3[1,s,t] <- 1-p3[2,s,t]-p3[3,s,t] # calculate last p3 p3[2,s,t] <- exp(lp3[2,s,t]) / (1 + exp(lp3[2,s,t]) + exp(lp3[3,s,t])) # backtransformation to {0,1} scale p3[3,s,t] <- exp(lp3[3,s,t]) / (1 + exp(lp3[2,s,t]) + exp(lp3[3,s,t])) # backtransformation to {0,1} scale } # s } #t int.p2 ~ dnorm(0, 0.01) beta.p2 ~ dnorm(0, 0.01) int.p32 ~ dnorm(0, 0.01) beta.p32 ~ dnorm(0, 0.01) int.p33 ~ dnorm(0, 0.01) beta.p33 ~ dnorm(0, 0.01) # Define state vector for (s in 1:R){ phi[s,1] <- 1 - psi # Prob. of non-occupation phi[s,2] <- psi * (1-r) # Prob. of occupancy without repro. phi[s,3] <- psi * r # Prob. of occupancy and repro } # Define observation matrix # Order of indices: true state, time, observed state for (s in 1:R){ for (t in 1:T){ p[1,s,t,1] <- 1 p[1,s,t,2] <- 0 p[1,s,t,3] <- 0 p[2,s,t,1] <- 1-p2[s,t] p[2,s,t,2] <- p2[s,t] p[2,s,t,3] <- 0 p[3,s,t,1] <- p3[1,s,t] p[3,s,t,2] <- p3[2,s,t] p[3,s,t,3] <- p3[3,s,t] } #t } #s # State-space likelihood # State equation: model of true states (z) for (s in 1:R){ z[s] ~ dcat(phi[s,]) } # Observation equation for (s in 1:R){ for (t in 1:T){ y[s,t] ~ dcat(p[z[s],s,t,]) } #t } #s # Derived quantities for (s in 1:R){ occ1[s] <- equals(z[s], 1) occ2[s] <- equals(z[s], 2) occ3[s] <- equals(z[s], 3) } n.occ[1] <- sum(occ1[]) # Sites in state 1 n.occ[2] <- sum(occ2[]) # Sites in state 2 n.occ[3] <- sum(occ3[]) # Sites in state 3 } ",fill=TRUE) sink() # Initial values zst <- apply(y, 1, max, na.rm = TRUE) zst[zst == "-Inf"] <- 1 inits <- function(){list(z = zst, beta.p32 = 0, beta.p33 = 0)} # Parameters monitored params <- c("r", "psi", "n.occ", "int.p2", "beta.p2", "int.p32", "beta.p32", "int.p33", "beta.p33") # MCMC settings ni <- 5000 nt <- 1 nb <- 2000 nc <- 3 # Call WinBUGS from R (BRT < 1 min) out3 <- bugs(win.data, inits, params, "model3.bug", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, debug =TRUE, bugs.directory = bugs.dir, working.directory = getwd()) # Summarize posteriors print(out3, dig = 3)