forked from steffenoppel/TRAL_IPM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathTRAL_IPM_marray_age_recruit_immat_FINAL.jags
384 lines (281 loc) · 18.7 KB
/
TRAL_IPM_marray_age_recruit_immat_FINAL.jags
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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
model {
#-------------------------------------------------
# integrated population model for the Gough TRAL population
# - age structured model with 10 age classes
# - adult survival based on CMR ringing data
# - pre breeding census, female-based assuming equal sex ratio & survival
# - productivity based on all areas incu and chick counts
# - simplified population process with assumed age at recruiting = 10 (Wanless et al. 2009)
# - adult breeders skipping when unsuccessful at rate of 22-32%, all successful breeders will skip
# - linked population process with SUM OF count data
# - v4 includes 3 scenarios of future projection: no change, improved fecundity, reduced adult survival
# - marray_v1 uses marray for survival estimation to speed up computation time
# - marray_simplified adjusts population process to make number of bird returning completely random
# - marray_simplified_v2 further simplifies return process to tie to recapture probability
# - marray_simplified_v3 includes carrying capacity for future projections and re-inserts breeding success
# - marray_simplified_v4 calculates Ntot and lambda based on Ntot (rather than just breeding pop).
# - marray_simplified_v5 includes two values for p (high and low monitoring effort years) and sets p.juv to 0 for first year.
# - marray_simplified_v5 also calculates mean propensity and recruitment from annual values for good monitoring years.
# - marray_age_recruit allows for varying juvenile recapture probability with age.
# - marray_age_recruit_immat creates a loop over immature age classes rather than specify each age separately
# - marray_age_recruit_immat_2008 sets juv surv to mean and starts in 2008 rather than 2004
# -------------------------------------------------
#-------------------------------------------------
# 1. PRIORS FOR ALL DATA SETS
#-------------------------------------------------
# -------------------------------------------------
# 1.1. Priors and constraints FOR FECUNDITY
# -------------------------------------------------
for (t in 1:T){
ann.fec[t] ~ dbeta(32,68) ## dnorm(0.32,10) T(0.001,0.999) ## Informative Priors on fecundity based on Wanless et al 2009
} #t
# -------------------------------------------------
# 1.2. Priors and constraints FOR POPULATION COUNTS
# -------------------------------------------------
for (s in 1:n.sites){ ### start loop over every study area
for (t in 1:T){ ### start loop over every year
sigma.obs[s,t] ~ dexp(0.1) #Prior for SD of observation process (variation in detectability)
tau.obs[s,t]<-pow(sigma.obs[s,t],-2)
}
}
# -------------------------------------------------
# 1.3. Priors and constraints FOR SURVIVAL
# -------------------------------------------------
### RECAPTURE PROBABILITY
for (gy in 1:2){ ## for good and poor monitoring years
mean.p.juv[gy] ~ dunif(0, 1) # Prior for mean juvenile recapture - should be higher than 20% if they survive!
mean.p.ad[gy] ~ dunif(0, 1) # Prior for mean adult recapture - should be higher than 20%
mu.p.juv[gy] <- log(mean.p.juv[gy] / (1-mean.p.juv[gy])) # Logit transformation
mu.p.ad[gy] <- log(mean.p.ad[gy] / (1-mean.p.ad[gy])) # Logit transformation
}
agebeta ~ dunif(0,1) # Prior for shape of increase in juvenile recapture probability with age
## RANDOM TIME EFFECT ON RESIGHTING PROBABILITY OF JUVENILES
for (t in 1:(n.occasions-1)){
for (j in 1:t){ ## zero by definition (these are never actually used)
p.juv[t,j] <- 0
}
for (j in (t+1):(n.occasions-1)){
logit(p.juv[t,j]) <- mu.p.juv[goodyear[j]] + agebeta*(j - t) + eps.p[j]
}
}
## PRIORS FOR RANDOM EFFECTS
sigma.p ~ dexp(1) # Prior for standard deviation
tau.p <- pow(sigma.p, -2)
### SURVIVAL PROBABILITY
mean.phi.juv ~ dbeta(75.7,24.3) # Prior for mean juvenile survival first year 0.757, second year 0.973 in Laysan albatross
mean.phi.ad ~ dbeta(91,9) # Prior for mean adult survival - should be higher than 70%
mu.juv <- log(mean.phi.juv / (1-mean.phi.juv)) # Logit transformation
mu.ad <- log(mean.phi.ad / (1-mean.phi.ad)) # Logit transformation
## PRIORS FOR RANDOM EFFECTS
sigma.phi ~ dexp(1) # Prior for standard deviation
tau.phi <- pow(sigma.phi, -2)
## RANDOM TIME EFFECT ON SURVIVAL AND ADULT RECAPTURE
for (j in 1:(n.occasions-1)){
logit(phi.juv[j]) <- mu.juv + eps.phi[j]*juv.poss[j]
logit(phi.ad[j]) <- mu.ad + eps.phi[j]
eps.phi[j] ~ dnorm(0, tau.phi)
logit(p.ad[j]) <- mu.p.ad[goodyear[j]] + eps.p[j]
eps.p[j] ~ dnorm(0, tau.p)
}
#-------------------------------------------------
# 2. LIKELIHOODS AND ECOLOGICAL STATE MODEL
#-------------------------------------------------
# -------------------------------------------------
# 2.1. System process: female based matrix model
# -------------------------------------------------
### INITIAL VALUES FOR COMPONENTS FOR YEAR 1 - based on deterministic multiplications
## ADJUSTED BASED ON PAST POPULATION SIZES WIT CHICK COUNTS SINCE 1999
IM[1,1,1] ~ dnorm(324,20) T(0,) ### number of 1-year old survivors is low because few chicks hatched in 2003 - CAN BE MANIPULATED
IM[1,1,2] <- 0
IM[1,1,3] <- IM[1,1,1] - IM[1,1,2]
IM[1,2,1] ~ dnorm(257,20) T(0,) ### number of 2-year old survivors is very low because very few chicks hatched in 2002 - CAN BE MANIPULATED
IM[1,2,2] <- IM[1,2,1]*p.juv.recruit.f[2]
IM[1,2,3] <- IM[1,1,1] - IM[1,1,2]
IM[1,3,1] ~ dnorm(462,20) T(0,) ### number of 3-year old survivors is higher because many chicks hatched in 2001 - CAN BE MANIPULATED
IM[1,3,2] <- IM[1,3,1]*p.juv.recruit.f[3]
IM[1,3,3] <- IM[1,1,1] - IM[1,1,2]
IM[1,4,1] ~ dnorm(207,20) T(0,) ### number of 4-year old survivors is very low because few chicks hatched in 2000 - CAN BE MANIPULATED
IM[1,4,2] <- IM[1,4,1]*p.juv.recruit.f[4]
IM[1,4,3] <- IM[1,1,1] - IM[1,1,2]
IM[1,5,1] ~ dnorm(700,10) T(0,) ### number of 5-year old survivors is huge because a lot of chicks hatched in 1999 - CAN BE MANIPULATED
IM[1,5,2] <- IM[1,5,1]*p.juv.recruit.f[5]
IM[1,5,3] <- IM[1,1,1] - IM[1,1,2]
IM[1,6,1] ~ dnorm(225,20) T(0,) ### very uncertain number of of 6-year old survivors because no data from 1998 or previously - CAN BE MANIPULATED
IM[1,6,2] <- IM[1,6,1]*p.juv.recruit.f[6]
IM[1,6,3] <- IM[1,1,1] - IM[1,1,2]
for(age in 7:30) {
IM[1,age,1] ~ dbin(pow(mean.phi.ad,(age-1)), round(IM[1,age-1,3]))
IM[1,age,2] <- IM[1,age,1]*p.juv.recruit.f[age]
IM[1,age,3] <- IM[1,age,1] - IM[1,age,2]
}
N.recruits[1] <- sum(IM[1,,2]) ### number of this years recruiters - irrelevant in year 1 as already included in Ntot.breed prior
Ntot.breed[1] ~ dnorm(1869,100) T(0,) ### sum of counts is 1869
JUV[1] ~ dnorm(510,100) T(0,) ### sum of chicks is 510
N.atsea[1] ~ dnorm(530,20) T(0,) ### unknown number - CAN BE MANIPULATED
Ntot[1]<-sum(IM[1,,3]) + Ntot.breed[1]+N.atsea[1] ## total population size is all the immatures plus adult breeders and adults at sea - does not include recruits in Year 1
### FOR EVERY SUBSEQUENT YEAR POPULATION PROCESS
for (tt in 2:T){
## THE PRE-BREEDING YEARS ##
## define recruit probability for various ages ##
for (age in 1:30) {
logit(p.juv.recruit[age,tt])<-mu.p.juv[2] + eps.p[tt+24] + (agebeta * age)
}
## IMMATURE MATRIX WITH 3 columns:
# 1: survivors from previous year
# 2: recruits in current year
# 3: unrecruited in current year (available for recruitment next year)
nestlings[tt] <- ann.fec[tt] * 0.5 * Ntot.breed[tt] ### number of locally produced FEMALE chicks
JUV[tt] ~ dpois(nestlings[tt]) ### need a discrete number otherwise dbin will fail, dpois must be >0
IM[tt,1,1] ~ dbin(phi.juv[tt+24], max(1,round(JUV[tt-1]))) ### number of 1-year old survivors
IM[tt,1,2] <- 0
IM[tt,1,3] <- IM[tt,1,1] - IM[tt,1,2]
for(age in 2:30) {
IM[tt,age,1] ~ dbin(phi.ad[tt+24], max(1,round(IM[tt-1,age-1,3])))
IM[tt,age,2] <- min(round(IM[tt,age-1,3]),IM[tt,age,1])*p.juv.recruit[age,tt]
IM[tt,age,3] <- IM[tt,age,1] - IM[tt,age,2]
}
N.recruits[tt] <- sum(IM[tt,,2]) ### number of this years recruiters
## THE BREEDING POPULATION ##
# Ntot.breed comprised of first-time breeders, previous skippers, and previous unsuccessful breeders
# simplified in simplified_v2 to just adult survivors with p.ad as proportion returning
N.ad.surv[tt] ~ dbin(phi.ad[tt+24], round(Ntot.breed[tt-1]+N.atsea[tt-1])) ### previous year's adults that survive
N.breed.ready[tt] ~ dbin(p.ad[tt+24], N.ad.surv[tt]) ### number of available breeders is proportion of survivors that returns
Ntot.breed[tt]<- round(N.breed.ready[tt]+N.recruits[tt]) ### number of counted breeders is sum of old breeders returning and first recruits
N.atsea[tt] <- round(N.ad.surv[tt]-N.breed.ready[tt]) ### potential breeders that remain at sea
### THE TOTAL TRAL POPULATION ###
Ntot[tt]<-sum(IM[tt,,3]) + Ntot.breed[tt]+N.atsea[tt] ## total population size is all the immatures plus adult breeders and adults at sea
} # tt
# -------------------------------------------------
# 2.2. Observation process for population counts: state-space model of annual counts
# -------------------------------------------------
for (s in 1:n.sites){ ### start loop over every study area
## Observation process
for (t in 1:T){
y.count[t,s] ~ dnorm(Ntot.breed[t]*prop.sites[s], tau.obs[s,t]) # Distribution for random error in observed numbers (counts)
} # run this loop over t= nyears
} ## end site loop
# -------------------------------------------------
# 2.3. Likelihood for fecundity: Logistic regression from the number of surveyed broods
# -------------------------------------------------
for (t in 1:(T-1)){
J[t] ~ dbin(ann.fec[t], R[t])
} # close loop over every year in which we have fecundity data
# -------------------------------------------------
# 2.4. Likelihood for adult and juvenile survival from CMR
# -------------------------------------------------
# Define the multinomial likelihood
for (t in 1:(n.occasions-1)){
marr.j[t,1:n.occasions] ~ dmulti(pr.j[t,], r.j[t])
marr.a[t,1:n.occasions] ~ dmulti(pr.a[t,], r.a[t])
}
# Define the cell probabilities of the m-arrays
# Main diagonal
for (t in 1:(n.occasions-1)){
q.ad[t] <- 1-p.ad[t] # Probability of non-recapture
for(j in 1:(n.occasions-1)){
q.juv[t,j] <- 1 - p.juv[t,j]
}
pr.j[t,t] <- 0
pr.a[t,t] <- phi.ad[t]*p.ad[t]
# Above main diagonal
for (j in (t+1):(n.occasions-1)){
pr.j[t,j] <- phi.juv[t]*prod(phi.ad[(t+1):j])*prod(q.juv[t,t:(j-1)])*p.juv[t,j]
pr.a[t,j] <- prod(phi.ad[t:j])*prod(q.ad[t:(j-1)])*p.ad[j]
} #j
# Below main diagonal
for (j in 1:(t-1)){
pr.j[t,j] <- 0
pr.a[t,j] <- 0
} #j
} #t
# Last column: probability of non-recapture
for (t in 1:(n.occasions-1)){
pr.j[t,n.occasions] <- 1-sum(pr.j[t,1:(n.occasions-1)])
pr.a[t,n.occasions] <- 1-sum(pr.a[t,1:(n.occasions-1)])
} #t
#-------------------------------------------------
# 3. DERIVED PARAMETERS FOR OUTPUT REPORTING
#-------------------------------------------------
## DERIVED POPULATION GROWTH RATE PER YEAR
for (t in 1:(T-1)){
lambda[t]<-Ntot[t+1]/max(1,Ntot[t]) ## division by 0 creates invalid parent value
} ## end year loop
## DERIVED MEAN FECUNDITY
mean.fec <- mean(ann.fec)
pop.growth.rate <- exp((1/(T-1))*sum(log(lambda[1:(T-1)]))) # Geometric mean
#-------------------------------------------------
# 4. PROJECTION INTO FUTURE
#-------------------------------------------------
## includes 3 scenarios
## scenario 1: projection with no changes in demography
## scenario 2: successful mouse eradication in 2021 - fecundity doubles
## scenario 3: increasing mouse impacts on adult survival (adult survival decreases by 10%)
## recruit probability
for (age in 1:30) {
logit(p.juv.recruit.f[age])<-mu.p.juv[2] + (agebeta * age)
}
# -------------------------------------------------
# 4.1. System process for future
# -------------------------------------------------
## LOOP OVER EACH SCENARIO
for(scen in 1:n.scenarios){
### ~~~~~~~~~~ COPY POPULATIONS FROM LAST YEAR OF DATA SERIES FOR FIRST FUTURE YEAR ~~~~~~~~~###
## IMMATURE MATRIX WITH 3 columns:
# 1: survivors from previous year
# 2: recruits in current year
# 3: unrecruited in current year (available for recruitment next year)
nestlings.f[scen,1] ~ dbin(fut.fec.change[scen]*mean.fec*0.5,round(Ntot.breed.f[scen,1])) ### number of locally produced FEMALE chicks based on average fecundity
IM.f[scen,1,1,1] ~ dbin(mean.phi.juv, max(1,round(JUV[T]))) ### number of 1-year old survivors
IM.f[scen,1,1,2] <- 0
IM.f[scen,1,1,3] <- IM.f[scen,1,1,1] - IM.f[scen,1,1,2]
for(age in 2:30) {
IM.f[scen,1,age,1] ~ dbin(mean.phi.ad, max(1,round(IM[T,age-1,3])))
IM.f[scen,1,age,2] <- min(round(IM[T,age-1,3]),IM.f[scen,1,age,1])*p.juv.recruit.f[age]
IM.f[scen,1,age,3] <- IM.f[scen,1,age,1] - IM.f[scen,1,age,2]
}
N.recruits.f[scen,1] <- sum(IM.f[scen,1,,2]) ### number of this years recruiters
N.ad.surv.f[scen,1] ~ dbin(mean.phi.ad, round(Ntot.breed[T]+N.atsea[T])) ### previous year's adults that survive
N.breed.ready.f[scen,1] ~ dbin(mean.p.ad[2], round(N.ad.surv.f[scen,1])) ### number of available breeders is proportion of survivors that returns, with fecundity INCLUDED in return probability
Ntot.breed.f[scen,1]<- round(N.breed.ready.f[scen,1]+N.recruits.f[scen,1]) ### number of counted breeders is sum of old breeders returning and first recruits
N.atsea.f[scen,1] <- round(N.ad.surv.f[scen,1]-N.breed.ready.f[scen,1]) ### potential breeders that remain at sea
N.succ.breed.f[scen,1] ~ dbin(mean.fec, round(Ntot.breed.f[scen,1])) ### these birds will remain at sea because they bred successfully
### THE TOTAL TRAL POPULATION ###
Ntot.f[scen,1]<-sum(IM.f[scen,1,,3])+Ntot.breed.f[scen,1]+N.atsea.f[scen,1] ## total population size is all the immatures plus adult breeders and adults at sea
### ~~~~~~~~~~ LOOP OVER ALL SUBSEQUENT FUTURE YEARS ~~~~~~~~~###
for (tt in 2:FUT.YEAR){
## INCLUDE CARRYING CAPACITY OF 2500 breeding pairs (slightly more than maximum ever counted)
carr.capacity[scen,tt] ~ dnorm(2500,5) T(0,)
## THE PRE-BREEDING YEARS ##
## because it goes for 30 years, all pops must be safeguarded to not become 0 because that leads to invald parent error
## IMMATURE MATRIX WITH 3 columns:
# 1: survivors from previous year
# 2: recruits in current year
# 3: unrecruited in current year (available for recruitment next year)
nestlings.f[scen,tt] ~ dbin(fut.fec.change[scen]*mean.fec*0.5,round(Ntot.breed.f[scen,tt])) ### number of locally produced FEMALE chicks based on average fecundity
IM.f[scen,tt,1,1] ~ dbin(mean.phi.juv, max(1,round(nestlings.f[scen,tt-1]))) ### number of 1-year old survivors
IM.f[scen,tt,1,2] <- 0
IM.f[scen,tt,1,3] <- IM.f[scen,tt,1,1] - IM.f[scen,tt,1,2]
for(age in 2:30) {
IM.f[scen,tt,age,1] ~ dbin(mean.phi.ad, max(1,round(IM.f[scen,tt-1,age-1,3])))
IM.f[scen,tt,age,2] <- min(round(IM.f[scen,tt-1,age-1,3]),IM.f[scen,tt,age,1])*p.juv.recruit.f[age]
IM.f[scen,tt,age,3] <- IM.f[scen,tt,age,1] - IM.f[scen,tt,age,2]
}
N.recruits.f[scen,tt] <- sum(IM.f[scen,tt,,2]) ### number of this years recruiters
## THE BREEDING POPULATION ##
N.ad.surv.f[scen,tt] ~ dbin(fut.surv.change[tt,scen]*mean.phi.ad, round((Ntot.breed.f[scen,tt-1]-N.succ.breed.f[scen,tt-1])+N.atsea.f[scen,tt-1])) ### previous year's adults that survive
N.prev.succ.f[scen,tt] ~ dbin(fut.surv.change[tt,scen]*mean.phi.ad, round(N.succ.breed.f[scen,tt-1])) ### these birds will remain at sea because tey bred successfully
N.breed.ready.f[scen,tt] ~ dbin(min(0.99,(mean.p.ad[2]/(1-mean.fec))), max(1,round(N.ad.surv.f[scen,tt]))) ### number of available breeders is proportion of survivors that returns, with fecundity partialled out of return probability
Ntot.breed.f[scen,tt]<- min(carr.capacity[scen,tt],round(N.breed.ready.f[scen,tt]+N.recruits.f[scen,tt])) ### number of counted breeders is sum of old breeders returning and first recruits
N.succ.breed.f[scen,tt] ~ dbin(fut.fec.change[scen]*mean.fec, round(Ntot.breed.f[scen,tt])) ### these birds will remain at sea because tey bred successfully
N.atsea.f[scen,tt] <- round(N.ad.surv.f[scen,tt]-N.breed.ready.f[scen,tt]+N.prev.succ.f[scen,tt]) ### potential breeders that remain at sea
### THE TOTAL TRAL POPULATION ###
Ntot.f[scen,tt]<-sum(IM.f[scen,tt,,3])+Ntot.breed.f[scen,tt]+N.atsea.f[scen,tt] ## total population size is all the immatures plus adult breeders and adults at sea
} ### end future loop
## CALCULATE ANNUAL POP GROWTH RATE ##
for (fut2 in 1:(FUT.YEAR-1)){
fut.lambda[scen,fut2] <- Ntot.f[scen,fut2+1]/max(1,Ntot.f[scen,fut2]) ### inserted safety to prevent denominator being 0
} # fut2
## DERIVED MEAN FUTURE GROWTH RATE
fut.growth.rate[scen] <- exp((1/(FUT.YEAR-1))*sum(log(fut.lambda[scen,1:(FUT.YEAR-1)]))) # Geometric mean
} # end future projection scenarios
} ## end model loop