#################################################################### # Code for workshop at BES 2017, Ghent, Belgium: An accessible introduction to Integrated Population Models for ecological and evolutionary research # Organiser: Mitch D. Weegman # Facilitators: David J. Hodgson & David N. Koons # Note: Attendees, please download Program R, the jagsUI package in R, and JAGS # Code based on Kery & Schaub (2012), revised by Mitch D. Weegman & David N. Koons (2017) #################################################################### #################################################################### # State-Space example: for IPM workshop at BES 2017, Ghent, Belgium, # 11-14 December 2017 #################################################################### library(jagsUI) # Population count data popcount <- c(32, 42, 64, 85, 82, 78, 73, 69, 79) # Specify model in BUGS language sink("ssm.jags") cat(" model { # Priors and constraints logN[1] ~ dnorm(3, 0.01) # Prior for initial population size meanr ~ dnorm(0.1, 25) # Prior for mean growth rate sigmaproc ~ dunif(0, 1) # Prior for sd of state process sigma2proc <- pow(sigmaproc, 2) tauproc <- pow(sigmaproc, -2) sigmaobs ~ dunif(0, 1) # Prior for sd of observation process sigma2obs <- pow(sigmaobs, 2) tauobs <- pow(sigmaobs, -2) # Likelihood # State process for (t in 1:(T-1)){ r[t] ~ dnorm(meanr, tauproc) logN[t+1] <- logN[t] + r[t] } # Observation process for (t in 1:T) { y[t] ~ dnorm(logN[t], tauobs) } # Derived population abundances on real scale for (t in 1:T) { N[t] <- exp(logN[t]) } } ",fill = TRUE) sink() # Bundle data jags.data <- list(y = log(popcount), T = length(popcount)) # Initial values inits <- function(){list(sigmaproc = runif(1, 0, 1), meanr = rnorm(1,0.1,30), logN = c(rnorm(1, 3, 0.5), rep(NA, (length(popcount)-1))), sigmaobs = runif(1, 0, 1))} # Parameters monitored parameters <- c("r", "meanr", "sigma2obs", "sigma2proc", "N") # MCMC settings ni <- 20000 nt <- 10 nb <- 10000 nc <- 3 # Call JAGS from R (JRT <1 min) ssm <- jags(jags.data, inits, parameters, "ssm.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, parallel = FALSE) # Summarize posteriors print(ssm, digits = 3) # Draw figure fitted <- lower <- upper <- numeric() year <- 1:length(popcount) n.years <- length(popcount) for (i in 1:n.years){ fitted[i] <- mean(ssm$sims.list$N[,i]) lower[i] <- quantile(ssm$sims.list$N[,i], 0.025) upper[i] <- quantile(ssm$sims.list$N[,i], 0.975)} m1 <- min(c(fitted, popcount, lower), na.rm = TRUE) m2 <- max(c(fitted, popcount, upper), na.rm = TRUE) par(mar = c(4.5, 4, 1, 1)) plot(0, 0, ylim = c(m1, m2), xlim = c(1, n.years), ylab = "Abundance", xlab = "Year", col = "black", type = "l", lwd = 2, axes = FALSE, frame = FALSE) axis(2, las = 1) axis(1, at = 1:n.years, labels = year) polygon(x = c(1:n.years, n.years:1), y = c(lower, upper[n.years:1]), col = "gray90", border = "gray90") points(popcount, type = "l", col = "black", lwd = 2) points(fitted, type = "l", col = "blue", lwd = 2) legend(x = 1, y = 100, legend = c("Counts", "Estimates"), lty = c(1, 1), lwd = c(2, 2), col = c("black", "blue"), bty = "n", cex = 1) ################################################################## # Fecundity example: for IPM workshop at BES 2017, Ghent, Belgium, # 11-14 December 2017 ################################################################## # Fecundity data J <- c(189, 274, 398, 538, 520, 476, 463, 438) # Number of offspring R <- c(28, 36, 57, 77, 81, 83, 77, 72) # Number of surveyed broods # Specify model in BUGS language sink("fec.jags") cat(" model { # Define the priors for the parameters B0 ~ dnorm(0,0.0001) # prior for intercept B1 ~ dnorm(0,0.0001) # prior for covariate slope # Constrain parameters using regression on the link scale for (t in 1:nyears){ log(f[t]) <- B0 + B1*X[t] # Fecundity } # The likelihood for (t in 1:nyears){ J[t] ~ dpois(rho[t]) rho[t] <- R[t] * f[t] } } ",fill = TRUE) sink() # Bundle data jags.data <- list(J = J, R = R, X = 1:length(J), nyears = length(J)) # Initial values inits <- function(){list(B0 = rnorm(1,0,0.5), B1 = rnorm(1,0,0.5))} # Parameters monitored parameters <- c("B0", "B1", "f") # MCMC settings ni <- 20000 nt <- 10 nb <- 10000 nc <- 3 # Call JAGS from R (JRT <1 min) fec <- jags(jags.data, inits, parameters, "fec.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, parallel = FALSE) # Summarize posteriors print(fec, digits = 3) # Draw figure fitted <- lower <- upper <- numeric() year <- 1:length(J) n.years <- length(J) for (i in 1:n.years){ fitted[i] <- mean(fec$sims.list$f[,i]) lower[i] <- quantile(fec$sims.list$f[,i], 0.025) upper[i] <- quantile(fec$sims.list$f[,i], 0.975)} m1 <- min(c(fitted, J/R, lower), na.rm = TRUE) m2 <- max(c(fitted, J/R, upper), na.rm = TRUE) par(mar = c(4.5, 4, 1, 1)) plot(0, 0, ylim = c(m1, m2), xlim = c(1, n.years), ylab = "Fecundity", xlab = "Year", col = "black", type = "l", lwd = 2, axes = FALSE, frame = FALSE) axis(2, las = 1) axis(1, at = 1:n.years, labels = year) polygon(x = c(1:n.years, n.years:1), y = c(lower, upper[n.years:1]), col = "gray90", border = "gray90") points(J/R, type = "l", col = "black", lwd = 2) points(fitted, type = "l", col = "blue", lwd = 2) legend(x = 5, y = 7, legend = c("Data", "Estimates"), lty = c(1, 1), lwd = c(2, 2), col = c("black", "blue"), bty = "n", cex = 1) ################################################################## # Survival example: for IPM workshop at BES 2017, Ghent, Belgium, # 11-14 December 2017 ################################################################## marray <- 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) # Specify model in BUGS language sink("cjs.jags") cat(" model { # Priors and constraints for (t in 1:(n.occasions-1)){ logit(phi[t]) <- mu.phi + epsilon.phi[t] epsilon.phi[t] ~ dnorm(0, tau.phi) logit(p[t]) <- mu.p + epsilon.p[t] epsilon.p[t] ~ dnorm(0, tau.p) } mean.phi ~ dunif(0, 1) # Prior for mean survival mu.phi <- log(mean.phi / (1-mean.phi)) # Logit transformation sigma.phi ~ dunif(0, 5) # Prior for standard deviation tau.phi <- pow(sigma.phi, -2) mean.p ~ dunif(0, 1) # Prior for mean recapture mu.p <- log(mean.p / (1-mean.p)) # Logit transformation sigma.p ~ dunif(0, 5) # Prior for standard deviation tau.p <- pow(sigma.p, -2) # Define the multinomial likelihood for (t in 1:(n.occasions-1)){ marr[t,1:n.occasions] ~ dmulti(pr[t,], rel[t]) } # Define the cell probabilities of the m-array: # Main diagonal for (t in 1:(n.occasions-1)){ q[t] <- 1-p[t] 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 jags.data <- list(marr = marray, n.occasions = dim(marray)[2], rel = rowSums(marray)) # Initial values inits <- function(){list(mean.phi = runif(1, 0, 1), sigma.phi = runif(1, 0, 5), mean.p = runif(1, 0, 1), sigma.p = runif(1, 0, 5))} # Parameters monitored parameters <- c("phi", "p", "mean.phi", "mean.p", "sigma.phi", "sigma.p") # MCMC settings ni <- 20000 nt <- 10 nb <- 10000 nc <- 3 # Call JAGS from R (JRT <1 min) cjs <- jags(jags.data, inits, parameters, "cjs.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb) # Summarize posteriors print(cjs, digits = 3) # Draw figure par(mar = c(4.5, 4, 1, 1)) lower <- upper <- numeric() T <- dim(marray)[2]-1 for (t in 1:T){ lower[t] <- quantile(cjs$sims.list$phi[,t], 0.025) upper[t] <- quantile(cjs$sims.list$phi[,t], 0.975) } plot(y = cjs$mean$phi, x = (1:T)+0.5, type = "b", pch = 16, ylim = c(0.1, 0.45), ylab = "Annual survival probability", xlab = "", axes = F) axis(1, at = seq(1,(T+1),2), labels = seq(1990,1998,2)) axis(1, at = 1:(T+1), labels = rep("", T+1), tcl = -0.25) axis(2, las = 1) mtext("Year", 1, line = 2.25) segments((1:T)+0.5, lower, (1:T)+0.5, upper) segments(1, cjs$mean$mean.phi, T+1, cjs$mean$mean.phi, lty = 2, col = "red", lwd = 2) segments(1, quantile(cjs$sims.list$mean.phi,0.025), T+1, quantile(cjs$sims.list$mean.phi, 0.025), lty = 2, col = "red") segments(1, quantile(cjs$sims.list$mean.phi, 0.975), T+1, quantile(cjs$sims.list$mean.phi, 0.975), lty = 2, col = "red") #################################################################### ##Hoopoe code for an IPM; IPM workshop at BES 2017, Ghent, Belgium, # 11-14 December 2017 #################################################################### nyears <- 9 # Number of years # Capture-recapture data: m-arrays for juveniles and adults (these are males and females together; apparent surival is equivalent between sexes) 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 count data popcount <- c(32, 42, 64, 85, 82, 78, 73, 69, 79) # Fecundity data J <- c(189, 274, 398, 538, 520, 476, 463, 438) # Number of offspring R <- c(28, 36, 57, 77, 81, 83, 77, 72) # Number of surveyed broods # Specify model in BUGS language sink("ipm.test.jags") 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 = 1 year # - Female-based model with prebreeding census # - All vital rates are assumed to fluctuate 'randomly' over time (environmental stochasticity) # - Explicit estimation of immigration #------------------------------------------------------------- #---------------------------------------- # 1. Define the priors for the parameters #---------------------------------------- # Initial population sizes n1 ~ dnorm(100, 0.0001)I(0,) # 1-year old individuals nadlocal ~ dnorm(100, 0.0001)I(0,) # Local adults >= 2 years nadimm ~ dnorm(100, 0.0001)I(0,) # Immigrants N1[1] <- round(n1) # Rounding to a whole number is needed Nadlocal[1] <- round(nadlocal) # when implementing branching processes Nadimm[1] <- round(nadimm) # for demographic stochasticity (below) # Mean demographic parameters (on appropriate link-function scale) l.mphij ~ dnorm(0, 0.0001) # apparent survival of juveniles l.mphia ~ dnorm(0, 0.0001) # apparent survival of 1 yr olds and adults l.mfec ~ dnorm(0, 0.0001) # fecundity (number of offspring produced per female) l.mim ~ dnorm(0, 0.0001) # immigrants l.p ~ dnorm(0, 0.0001) # recapture probability (assumed to be constant) # Precision of standard deviations for temporal process variability in each # demographic parameter (environmental stochasticity) 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 temporal process error terms (environmental stochasticity) # (The bounding of the s.d. above helps with convergence, and we can check # posterior results if the bounding is too restrictive) for (t in 1:(nyears-1)){ epsilon.phij[t] ~ dnorm(0, tau.phij) epsilon.phia[t] ~ dnorm(0, tau.phia) epsilon.fec[t] ~ dnorm(0, tau.fec) epsilon.im[t] ~ dnorm(0, tau.im) } #------------------------------------------------------------- # 2. Constrain parameters using mixed models on the link scale #------------------------------------------------------------- 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(f[t]) <- l.mfec + epsilon.fec[t] # Fecundity log(omega[t]) <- l.mim + epsilon.im[t] # Immigration rate logit(p[t]) <- l.p # Recapture probability } #----------------------- # 3. Derived parameters #----------------------- mphij <- exp(l.mphij)/(1+exp(l.mphij)) # Mean juvenile apparent survival probability mphia <- exp(l.mphia)/(1+exp(l.mphia)) # Mean adult apparent survival probability mfec <- exp(l.mfec) # Mean fecundity (productivity) mim <- exp(l.mim) # 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 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]) # branching process for demographic stochasticity Nadlocal[t] ~ dbin(phia[t-1], Ntot[t-1]) # b.p. for demographic stochasticity mimm[t] <- Ntot[t-1] * omega[t-1] Nadimm[t] ~ dpois(mimm[t]) # b.p. for demographic stochasticity } # 4.1.2 Observation process for (t in 1:nyears){ Ntot[t] <- Nadlocal[t] + Nadimm[t] + N1[t] y[t] ~ dpois(Ntot[t]) # models observation (sampling) error } # 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]) } # 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 fecundity data: Poisson regression for (t in 1:(nyears-1)){ J[t] ~ dpois(rho[t]) rho[t] <- R[t] * f[t] } } ",fill = TRUE) sink() # Bundle data jags.data <- list(nyears = nyears, marray.j = marray.j, marray.a = marray.a, y = popcount, J = J, R = R, r.j = rowSums(marray.j), r.a = rowSums(marray.a)) # Initial values inits <- function(){list(l.mphij = rnorm(1, 0.2, 0.5), l.mphia = rnorm(1, 0.2, 0.5), l.mfec = rnorm(1, 0.2, 0.5), l.mim = rnorm(1, 0.2, 0.5), l.p = rnorm(1, 0.2, 1), sig.phij = runif(1, 0.1, 10), sig.phia = runif(1, 0.1, 10), sig.fec = runif(1, 0.1, 10), sig.im = runif(1, 0.1, 10), n1 = round(runif(1, 1, 50), 0), nadlocal = round(runif(1, 5, 50), 0), nadimm = round(runif(1, 1, 50), 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", "Nadlocal", "Nadimm", "Ntot") # MCMC settings ni <- 200000 nt <- 10 nb <- 100000 nc <- 3 # Call JAGS from R (BRT 5 min) ipm.hoopoe <- jags(jags.data, inits, parameters, "ipm.test.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, parallel = TRUE) # Summarize posteriors print(ipm.hoopoe, digits = 3) # Check traceplots for convergence (chains should be mixed) par(mfrow = c(2, 2)) traceplot(ipm.hoopoe) # Draw figures par(mfrow = c(2, 2), cex.axis = 1.2, cex.lab = 1.2, mar = c(3, 6, 1.5, 2), las = 1) lower <- upper <- numeric() year <- 2002:2010 for (i in 1:nyears){ lower[i] <- quantile(ipm.hoopoe$sims.list$Ntot[,i], 0.025) upper[i] <- quantile(ipm.hoopoe$sims.list$Ntot[,i], 0.975)} m1 <- min(c(ipm.hoopoe$mean$Ntot, popcount, lower), na.rm = T) m2 <- max(c(ipm.hoopoe$mean$Ntot, popcount, 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(popcount, type = "l", col = "grey30", lwd = 2) points(ipm.hoopoe$mean$Ntot, type = "l", col = "blue", lwd = 2) legend(x = 2, y = 25, 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-1 for (t in 1:T){ lower[t] <- quantile(ipm.hoopoe$sims.list$phij[,t], 0.025) upper[t] <- quantile(ipm.hoopoe$sims.list$phij[,t], 0.975)} par(mgp=c(3.8,1,0)) plot(y = ipm.hoopoe$mean$phij, x = (1:T)+0.5, xlim= c(1, 9), type = "b", pch = 16, ylim = c(0, 0.6), ylab = "Annual survival probability", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2) axis(2) axis(1, at = 1:(T+1), labels = 2002:2010) segments((1:T)+0.5, lower, (1:T)+0.5, upper, lwd = 2) segments(1, ipm.hoopoe$mean$mphij, T+1, ipm.hoopoe$mean$mphij, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mphij, 0.025), T+1, quantile(ipm.hoopoe$sims.list$mphij, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mphij, 0.975), T+1, quantile(ipm.hoopoe$sims.list$mphij, 0.975), lty = 2, col = "red") for (t in 1:T){ lower[t] <- quantile(ipm.hoopoe$sims.list$phia[,t], 0.025) upper[t] <- quantile(ipm.hoopoe$sims.list$phia[,t], 0.975)} points(y=ipm.hoopoe$mean$phia, 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.hoopoe$mean$mphia, T+1, ipm.hoopoe$mean$mphia, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mphia, 0.025), T+1, quantile(ipm.hoopoe$sims.list$mphia, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mphia, 0.975), T+1, quantile(ipm.hoopoe$sims.list$mphia, 0.975), lty = 2, col = "red") legend(x = 4.5, y = 0.61, legend = c("Adults", "Juveniles"), pch = c(1, 16), bty = "n") lower <- upper <- numeric() T <- nyears-1 for (t in 1:T){ lower[t] <- quantile(ipm.hoopoe$sims.list$f[,t], 0.025) upper[t] <- quantile(ipm.hoopoe$sims.list$f[,t], 0.975)} plot(y=ipm.hoopoe$mean$f, x = (1:T), type = "b", pch = 16, ylim = c(5, 9), ylab = "Fecundity (fledgling / female)", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2) axis(2) axis(1, at = 1:T, labels = 2003:2010) segments((1:T), lower, (1:T), upper) segments(1, ipm.hoopoe$mean$mfec, T, ipm.hoopoe$mean$mfec, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mfec, 0.025), T, quantile(ipm.hoopoe$sims.list$mfec, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mfec, 0.975), T, quantile(ipm.hoopoe$sims.list$mfec, 0.975), lty = 2, col = "red") lower <- upper <- numeric() T <- nyears-1 for (t in 1:T){ lower[t] <- quantile(ipm.hoopoe$sims.list$omega[,t], 0.025) upper[t] <- quantile(ipm.hoopoe$sims.list$omega[,t], 0.975)} plot(y = ipm.hoopoe$mean$omega, x = (1:T)+0.5, xlim = c(1, 9), type = "b", pch = 16, ylim = c(0, 1.1), ylab = "Immigration rate", xlab = "", axes = F, cex = 1.5, frame = F, lwd = 2) axis(2) axis(1, at = 1:(T+1), labels = 2002:2010) segments((1:T)+0.5, lower, (1:T)+0.5, upper) segments(1, ipm.hoopoe$mean$mim, T+1, ipm.hoopoe$mean$mim, lty = 2, lwd = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mim, 0.025), T+1, quantile(ipm.hoopoe$sims.list$mim, 0.025), lty = 2, col = "red") segments(1, quantile(ipm.hoopoe$sims.list$mim, 0.975), T+1, quantile(ipm.hoopoe$sims.list$mim, 0.975), lty = 2, col = "red") # Calculate means and 95% credible intervals from posteriors nyears <- 9 lambda.h <- lam.lower.h <- lam.upper.h <- numeric() Fitted.h <- lower.h <- upper.h <- matrix(NA, nrow = nyears-1, ncol = 4) for (i in 1:(nyears-1)){ lambda.h[i] <- mean(ipm.hoopoe$sims.list$lambda[,i]) lam.lower.h[i] <- quantile(ipm.hoopoe$sims.list$lambda[,i], 0.025) lam.upper.h[i] <- quantile(ipm.hoopoe$sims.list$lambda[,i], 0.975) } for (i in 1:(nyears-1)){ Fitted.h[i,1] <- mean(ipm.hoopoe$sims.list$phij[,i]) lower.h[i,1] <- quantile(ipm.hoopoe$sims.list$phij[,i], 0.025) upper.h[i,1] <- quantile(ipm.hoopoe$sims.list$phij[,i], 0.975) } for (i in 1:(nyears-1)){ Fitted.h[i,2] <- mean(ipm.hoopoe$sims.list$phia[,i]) lower.h[i,2] <- quantile(ipm.hoopoe$sims.list$phia[,i], 0.025) upper.h[i,2] <- quantile(ipm.hoopoe$sims.list$phia[,i], 0.975) } for (i in 1:(nyears-1)){ Fitted.h[i,3] <- mean(ipm.hoopoe$sims.list$f[,i]) lower.h[i,3] <- quantile(ipm.hoopoe$sims.list$f[,i], 0.025) upper.h[i,3] <- quantile(ipm.hoopoe$sims.list$f[,i], 0.975) } for (i in 1:(nyears-1)){ Fitted.h[i,4] <- mean(ipm.hoopoe$sims.list$omega[,i]) lower.h[i,4] <- quantile(ipm.hoopoe$sims.list$omega[,i], 0.025) upper.h[i,4] <- quantile(ipm.hoopoe$sims.list$omega[,i], 0.975) } # Calculate correlation coefficients correl.h <- matrix(NA, ncol = 4, nrow = 30000) for (i in 1:30000){ correl.h[i,1] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$phij[i,]) correl.h[i,2] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$phia[i,]) correl.h[i,3] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$f[i,]) correl.h[i,4] <- cor(ipm.hoopoe$sims.list$lambda[i,], ipm.hoopoe$sims.list$omega[i,]) } # Credible intervals of correlation coefficients quantile(correl.h[,1], c(0.05, 0.5, 0.95), na.rm = TRUE) quantile(correl.h[,2], c(0.05, 0.5, 0.95), na.rm = TRUE) quantile(correl.h[,3], c(0.05, 0.5, 0.95), na.rm = TRUE) quantile(correl.h[,4], c(0.05, 0.5, 0.95), na.rm = TRUE) # Compute the posterior modes of correlation coefficients m <- density(correl.h[,1], na.rm = TRUE) m$x[which(m$y==max(m$y))] m <- density(correl.h[,2], na.rm = TRUE) m$x[which(m$y==max(m$y))] m <- density(correl.h[,3], na.rm = TRUE) m$x[which(m$y==max(m$y))] m <- density(correl.h[,4], na.rm = TRUE) m$x[which(m$y==max(m$y))] # Probability that correlation coefficients (r) > 0 sum(correl.h[!is.na(correl.h[,1]),1]>0)/30000 sum(correl.h[!is.na(correl.h[,2]),2]>0)/30000 sum(correl.h[!is.na(correl.h[,3]),3]>0)/30000 sum(correl.h[!is.na(correl.h[,4]),4]>0)/30000 # Draw figures par(mfrow = c(2, 2), mar = c(4.5, 4, 1.5, 1), mgp=c(3, 1, 0), las = 1, cex = 1) linecol <- c("grey70") plot(y = lambda.h, Fitted.h[,1], type = "n", xlim = c(0.05, 0.25), ylim = c(0.6, 1.8), ylab = "Population growth rate", xlab = "Juvenile survival", frame = FALSE, pch = 19) segments(Fitted.h[,1], lam.lower.h, Fitted.h[,1], lam.upper.h, col = linecol) segments(lower.h[,1], lambda.h, upper.h[,1], lambda.h, col = linecol) points(y = lambda.h, Fitted.h[,1], pch = 19, col = "black") text(x = 0.13, y = 0.75, "r = 0.78 (0.14, 0.90)", pos = 4, font = 3, cex = 0.8) text(x = 0.13, y = 0.65, "P(r>0) = 0.98", pos = 4, font = 3, cex = 0.8) par(mar = c(4.5, 4, 1.5, 1)) plot(y = lambda.h, Fitted.h[,2], type = "n", xlim = c(0.3, 0.51), ylim = c(0.6, 1.8), ylab = "", xlab = "Adult survival", frame.plot = FALSE, pch = 19) segments(Fitted.h[,2], lam.lower.h, Fitted.h[,2], lam.upper.h, col = linecol) segments(lower.h[,2], lambda.h, upper.h[,2], lambda.h, col = linecol) points(y = lambda.h, Fitted.h[,2], pch = 19, col = "black") text(x = 0.295, y = 1.75, "r = 0.25 (-0.46, 0.73)", pos = 4, font = 3, cex = 0.8) text(x = 0.295, y = 1.65, "P(r>0) = 0.70", pos = 4, font = 3, cex = 0.8) par(mar = c(4.5, 4, 2, 1)) plot(y = lambda.h, Fitted.h[,3], type = "n", xlim = c(5, 8.5), ylim = c(0.6, 1.8), ylab = "Population growth rate", xlab = "Fecundity", frame.plot = FALSE, pch = 19) segments(Fitted.h[,3], lam.lower.h, Fitted.h[,3], lam.upper.h, col = linecol) segments(lower.h[,3], lambda.h, upper.h[,3], lambda.h, col = linecol) points(y=lambda.h, Fitted.h[,3], pch = 19, col = "black") text(x = 5, y = 1.75, "r = 0.70 (0.19, 0.87)", pos = 4, font = 3, cex = 0.8) text(x = 5, y = 1.65, "P(r>0) = 0.99", pos = 4, font = 3, cex = 0.8) par(mar = c(4.5, 4, 2, 1)) plot(y = lambda.h, Fitted.h[,4], type = "n", xlim = c(0, 0.8), ylim = c(0.6, 1.8), ylab = "", xlab = "Immigration rate", frame.plot = FALSE, pch = 19) segments(Fitted.h[,4], lam.lower.h, Fitted.h[,4], lam.upper.h, col = linecol) segments(lower.h[,4], lambda.h, upper.h[,4], lambda.h, col = linecol) points(y=lambda.h, Fitted.h[,4], pch = 19, col = "black") text(x = 0.35, y = 0.8, "r = 0.84 (-0.30, 0.93)", pos = 4, font = 3, cex = 0.8) text(x = 0.35, y = 0.7, "P(r>0) = 0.86", pos = 4, font = 3, cex = 0.8)