-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathexample2-extended_model-Czech_CES_bird_ringing.R
169 lines (126 loc) · 5.49 KB
/
example2-extended_model-Czech_CES_bird_ringing.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
#######
#
# Extension of the Pradel (1996) model for transients and multiple sites with different temporal coverage
#
# Code for Chapter 6 | Case study: Constant Effort Sites mist-netting scheme in the Czech Republic
# of the manuscript:
#
# 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)
fast_model_run <- FALSE # very brief model run (few iterations) just for testing
run_parallel <- TRUE # disable only for debugging
source("../data/ces-data.R", chdir = TRUE) # load CES data
source("functions/data-prepare.R") # functions to prepare data for the model
### prepare the data
YRS <- min(ces$year):max(ces$year) # year span
#sp_id <- 12510 # Euring species id, 12510 = Acrocephalus scirpaceus
sp_id <- 14640 # Euring species id, 14640 = Parus major
ces$site <- ces$site_ces #
ces$time <- ces$year # use year as temporal occasion
# build global table of visits (long format)
visits <- visits(ces[ces$year %in% YRS,])
# process the individual captures and prepare the data in the required format:
filter <- ces$euring == sp_id & ces$age >= 4 & ces$year %in% YRS # filter adults of given species
data <- prepare_data(
cap = ces[filter,],
visits = visits,
known_resid_first_cap = function(cap_sel) recaptured_on_first_occasion(cap_sel, visit.col = "visit")
)
### run the model
source("../model/model-pradel_with_transients-extended.nimble.R") # run the extended model, using Nimble
save(out, file = paste0("model-", sp_id, ".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)
n.occasions <- length(YRS)
# 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(Year = YRS, 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
)