added new files, deleted some

COMING_UP.md

@@ -0,0 +1,11 @@
+
+The code is not complete yet. I expect the following updates to come up soon!
+
+- code for simulations (Appendix A)
+- code for JAGS
+- habitat specific model and code for Figure 6 
+- model optimizations for larger amount of sites
+- implementation of Maximum Likelihood model fitting
+
+Stay tuned!
+

code/example1-basic_model-dipper.R

@@ -0,0 +1,174 @@
+#######
+#
+# Extension of the Pradel (1996) model for transients and multiple sites with different temporal coverage
+# 
+# A simple example of the simplest - basic version of the model on the classic and well-known Dipper data set
+# (Warning: the Dipper data set might not meet the assumptions of Pradel model, notably the constant
+# sampled area requirement. It is only used as a demonstration of the model usage.)
+#
+# Telenský, T., Storch, D., Klvaňa, P., Reif, J. Extension of Pradel capture-recapture survival-recruitment model accounting for transients. 
+#	Methods in Ecology and Evolution. https://doi.org/10.1111/2041-210X.14262
+#
+# THIS IS THE MAIN SCRIPT TO RUN THE MODEL, START HERE
+
+library(coda)
+library(RMark) # dipper data set
+
+source("functions/data-capture_history.R")
+source("functions/marray.R")
+
+fast_model_run <- FALSE # very brief model run (few iterations) just for testing
+run_parallel <- TRUE # disable only for debugging
+
+### prepare the data
+
+data(dipper)
+
+# test the assumption that the capture history only contains 0s and 1s:
+stopifnot(length(grep("^[01]+$", dipper$ch, invert = TRUE)) == 0)
+
+# convert the string capture histories to matrix CH
+CH <- char_to_CH(dipper$ch)
+
+n.occasions <- ncol(CH)
+
+# create m-array
+data <- marray_multistage(CH)
+data$k <- n.occasions
+
+### run the model
+
+source("../model/model-pradel_with_transients-basic.nimble.R") # run the extended model, using Nimble
+
+save(out, file = "model-dipper.Rdata")
+
+### post-hoc calculation of the derived parameters. Performed on the MCMC samples
+### to ensure proper propagation of uncertainty of the primary model parameters to these derived quantities
+### (for more explanation see e.g. Kéry, M., & Schaub, M. (2012). Bayesian Population Analysis using WinBUGS. Elsevier.)
+
+source("functions/mcmc.R")
+
+mcmcsum(out$mcmc) # summary of the MCMC chains of the primary parameters
+
+mc <- as.mcmc.list(out$mcmc)
+mt <- as.matrix(mc, chains = TRUE, iters = TRUE)
+
+# extract the vital rates calculated by the model
+
+surv <- mt[,paste0('surv[',1:(n.occasions-1),']')]
+sen <- mt[,paste0('sen[',1:(n.occasions-1),']')]
+
+# calculate lambda = surv/sen
+
+lambda <- matrix(nrow = nrow(mt), ncol = n.occasions - 1)
+lambda <- surv/sen
+colnames(lambda) <- paste0("lambda[", 1:ncol(lambda), "]")
+mt <- cbind(mt, lambda)
+
+# calculate recruitment
+
+recr <- matrix(nrow = nrow(mt), ncol = n.occasions - 1)
+recr <- lambda - surv
+colnames(recr) <- paste0("recr", "[", 1:ncol(recr), "]")
+mt <- cbind(mt, recr)
+
+# calculate adult relative population index (equals 1 on the first occasion)
+
+pop.idx <- matrix(nrow = nrow(mt), ncol = n.occasions)
+pop.idx[,1] <- 1
+for (i in 2:n.occasions) {
+	pop.idx[,i] <- pop.idx[,i-1] * lambda[,i-1]
+}
+pop.idx.sum <- as.data.frame(cbind(time = 1:n.occasions, t(apply(pop.idx, 2, function (x) { c(Mean = mean(x), quantile(x, c(0.025, 0.975))) }))))
+
+# now create a summary of MCMC chains of all parameters, including the derived ones!
+mc <- matrix2mcmc(mt)
+mcmcsum(mc)
+
+### As an example of a simple follow-up analysis with proper propagation of uncertainty, we calculate correlations
+### of the demographic rates and density (population index). The follow-up analysis (here, correlation) is performed
+### for each MCMC sample from the posterior distributions of the parameter estimates. This ensures proper propagation
+### of uncertainty of the primary model parameters to the uncertainty of the resultant statistics (here, pearson correlation coefficient)
+### (for more explanation see e.g. Kéry, M., & Schaub, M. (2012). Bayesian Population Analysis using WinBUGS. Elsevier.)
+
+n.iter <- nrow(surv)
+r <- c()
+for (i in 1:n.iter) {
+	r[i] <- cor.test(surv[i,], recr[i,])$estimate
+}
+r.surv_recr <- c("r.surv_recr" = mean(r), quantile(r, c(0.025, 0.975)))
+
+n.iter <- nrow(lambda)
+r <- c()
+for (i in 1:n.iter) {
+	r[i] <- cor.test(lambda[i,], surv[i,])$estimate
+}
+r.lambda_surv <- c("r.lambda_surv" = mean(r), quantile(r, c(0.025, 0.975)))
+
+n.iter <- nrow(lambda)
+r <- c()
+for (i in 1:n.iter) {
+	r[i] <- cor.test(lambda[i,], recr[i,])$estimate
+}
+r.lambda_recr <- c("r.lambda_recr" = mean(r), quantile(r, c(0.025, 0.975)))
+
+n.iter <- nrow(lambda)
+r <- c()
+for (i in 1:n.iter) {
+	r[i] <- cor.test(lambda[i,], pop.idx[i,1:(n.occasions-1)])$estimate
+}
+r.density_lambda <- c("r.density_lambda" = mean(r), quantile(r, c(0.025, 0.975)))
+
+n.iter <- nrow(lambda)
+r <- c()
+for (i in 1:n.iter) {
+	r[i] <- cor.test(surv[i,], pop.idx[i,1:(n.occasions-1)])$estimate
+}
+r.density_surv <- c("r.density_surv" = mean(r), quantile(r, c(0.025, 0.975)))
+
+n.iter <- nrow(lambda)
+r <- c()
+for (i in 1:n.iter) {
+	r[i] <- cor.test(recr[i,], pop.idx[i,1:(n.occasions-1)])$estimate
+}
+r.density_recr <- c("r.density_recr" = mean(r), quantile(r, c(0.025, 0.975)))
+
+
+# print the results
+res <- rbind(r.surv_recr, r.lambda_surv, r.lambda_recr, r.density_lambda, r.density_surv, r.density_recr)
+colnames(res)[1] <- 'Mean'
+print(res)
+
+##### Talon plot
+
+library(talonplot)
+
+big_font <- TRUE
+
+talon_plot(
+	recr_pts = colMeans(recr), surv_pts = colMeans(surv),
+	recr_samples = recr, surv_samples = surv,
+	recr_CI_low = apply(recr, 2, quantile, 0.025),
+	recr_CI_high = apply(recr, 2, quantile, 0.975),
+	surv_CI_low = apply(surv, 2, quantile, 0.025),
+	surv_CI_high = apply(surv, 2, quantile, 0.975),
+	CI_type = "ellipse",
+	CI_transform = FALSE,
+	main = "Parus major",
+	big_font = big_font,
+	plot.samples = TRUE
+)
+
+#### Now, calculate the original Pradel model with MARK, to compare the values
+# Keep in mind that the results will differ because our extended model contains the transience parameter
+
+mod <- mark(dipper,
+	model.parameters = list(Phi = list(formula = ~0+time), p = list(formula = ~1), f = list(formula = ~0+time)), 
+	model="Pradrec")
+
+
+
+
+
+
+

code/fig2.R

@@ -0,0 +1,42 @@
+# Figure 2
+#
+# This script generates a part of Figure 2 - a graph showing the increasing amount of recaptured individuals
+# in a constant population.
+#
+# 
+
+all_cap <- 1
+rCt <- 0.7 # residency probability (within all captured animals)
+p <- 0.4 # capture probability
+sen <- 0.7 # seniority
+T <- 5 # number of temporal occasions
+
+# Calculate Equation 3 from the paper
+xi <- c()  
+xi[1] <- 1
+for (i in 2:T) {
+	xi[i] <- (1-sen) + sen*(1 - p)*xi[i-1]
+}
+
+# AC: allready captured before ("marked")
+AC <- rCt * (1 - xi) 
+
+
+par(cex = 1.3)
+plot(1:T, AC, type = "l", ylim = c(0,1), xlab = "", ylab = "", yaxt = "n", xaxt = "n", bty = "n")
+title(xlab = "time (occasions)", line = 2)
+axis(1, pos = 0)
+axis(2, at = c(0, rCt, 1), labels = c(0, "", 1), las = 1, pos = 1)
+lines(c(1,T), c(rCt, rCt))
+lines(c(1,T), c(1, 1))
+polygon(c(1:T,T,1), c(AC,1,1), density = 10, angle = 45)
+polygon(c(1:T,T,1), c(AC,0,0), density = 10, angle = -45)
+
+savePlot("res_tr_graf", "wmf")
+savePlot("res_tr_graf", "emf")
+
+#abline(h = rCt) # this plots beyond the vertical axes, unfortunatelly
+#abline(h = 1)
+
+
+

code/fig4.R

@@ -0,0 +1,3 @@
+
+# See script example2-extended_model-Czech_CES_bird_ringing.R to generate this Figure.
+