-
Notifications
You must be signed in to change notification settings - Fork 2
/
EGVU_PVA_2023_2scen.jags
136 lines (105 loc) · 7.42 KB
/
EGVU_PVA_2023_2scen.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
model {
#-------------------------------------------------
# 1. PRIORS FOR ALL DATA SETS
#-------------------------------------------------
# Priors and constraints
N.est.count[1,1] ~ dunif(20,40) ## draw random value from a uniform distribution between 50 and 200 for initial population size
N.est.count[1,2] ~ dunif(60,100) ## draw random value from a uniform distribution between 50 and 200 for initial population size
N.est.count[1,3] ~ dunif(30,60) ## draw random value from a uniform distribution between 50 and 200 for initial population size
N.est.count[1,4] ~ dunif(40,70) ## draw random value from a uniform distribution between 50 and 200 for initial population size
for (t in 1:2) { ###allow two separate time trends
mean.lambda.count[t] ~ dunif(0.1,2) #Prior for mean growth rate
sigma.proc.count[t] ~ dunif(0.1,10) #Prior for SD of state process (annual variation in pop size)
sigma2.proc.count[t] <-pow(sigma.proc.count[t] ,2)
tau.proc.count[t] <-pow(sigma.proc.count[t] ,-2)
}
# Priors and constraints FOR POPULATION COUNTS OBSERVATION
for (s in 1:countries){ ### start loop over every country
sigma.obs.count[s] ~ dunif(0.1,100) #Prior for SD of observation process (variation in detectability)
tau.obs.count[s]<-pow(sigma.obs.count[s],-2)
}
##### Likelihood function
## State process for entire time series
for (t in 1:(T.count-1)){
for (s in 1:countries){ ### start loop over every country
lambda.count[t,s] ~ dnorm(mean.lambda.count[period[t]], tau.proc.count[period[t]]) # Distribution for random error of growth rate
N.est.count[t+1,s]<-N.est.count[t,s]*lambda.count[t,s] # Linear predictor (population size based on past pop size and change rate)
}
}
## Observation process
for (t in 1:T.count){
for (s in 1:countries){ ### start loop over every country
y.count[t,s] ~ dnorm(N.est.count[t,s], tau.obs.count[s]) # Distribution for random error in observed numbers (counts)
}
}
##### Derived quantity
for (t in 1:T.count){
Ntot[t]<-sum(N.est.count[t,])
}
# Assessing the fit of the state-space model
# 1. Compute fit statistic for observed data
# Discrepancy measure: mean absolute error
for (t in 1:T.count){
C.exp[t] <- Ntot[t] # Expected counts
C.act[t] <- sum(y.count[t,])
Dssm.obs[t] <- abs((C.act[t] - C.exp[t]) / C.act[t]) # Discrepancy measure
}
Dmape.obs <- sum(Dssm.obs)
# 2. Compute fit statistic for replicate data
# Discrepancy measure: mean absolute error
for (t in 1:T.count){
for (s in 1:countries){ ### start loop over every country
C.rep.c[t,s] ~ dnorm(N.est.count[t,s], tau.obs.count[s]) # Distribution for random error in observed numbers (counts)
}
C.rep[t] <- sum(C.rep.c[t,]) # Generate replicate data
Dssm.rep[t] <- abs((C.rep[t] - C.exp[t]) / C.rep[t]) # Discrepancy measure
}
Dmape.rep <- sum(Dssm.rep)
# -------------------------------------------------
# 2. specify parameters for projection
# -------------------------------------------------
mean.fec ~ dunif(0.914,1.109)
for (scen in 1:2) { ### explore two scenarios with pre-LIFE estimates and post-LIFE estimates
phi.ad[scen] ~ dbeta(ad.phi.mean[scen],ad.phi.sd[scen])
phi.juv[scen] ~ dbeta(juv.phi.mean[scen],juv.phi.sd[scen])
phi.capt.juv[scen] ~ dbeta(sec.capt.phi.mean[scen],sec.capt.phi.sd[scen])
phi.sec[scen] ~ dbeta(sec.phi.mean[scen],sec.phi.sd[scen])
phi.third[scen] ~ dbeta(third.phi.mean[scen],third.phi.sd[scen])
# -------------------------------------------------
# 3. PREDICTION INTO THE FUTURE
# -------------------------------------------------
## POPULATION PROCESS
## need to set up population structure
nestlings.f[1,scen] <- round(Nterr.f[1,scen]*mean.fec * 0.5) ##JUV[T.count]
N1.f[1,scen] ~ dbin((mean.fec * 0.5*phi.juv[scen]),round(Ntot[T.count-1]))
N2.f[1,scen] ~ dbin((mean.fec * 0.5*phi.juv[scen]*phi.sec[scen]),round(Ntot[T.count-2]))
N3.f[1,scen] ~ dbin((mean.fec * 0.5*phi.juv[scen]*phi.sec[scen]*phi.third[scen]),round(Ntot[T.count-3]))
N4.f[1,scen] ~ dbin((mean.fec * 0.5*phi.juv[scen]*phi.sec[scen]*phi.third[scen]*phi.ad[scen]),round(Ntot[T.count-4]))
N5.f[1,scen] ~ dbin((mean.fec * 0.5*phi.juv[scen]*phi.sec[scen]*phi.third[scen]*phi.ad[scen]*phi.ad[scen]),max(1,round(Ntot[T.count-5])))
# N6.f[1] ~ dbin(phi.ad[scen],round(Ntot[T.count]*0.44+(mean.fec * 0.5*phi.juv[scen]*phi.sec[scen]*phi.third[scen]*phi.ad[scen]*phi.ad[scen]*phi.ad[scen]*Ntot[T.count-6]))) ### number of 6-year or older (adult) birds
# Nterr.f[1] <- round((N4.f[1,scen] * 0.024) + (N5.f[1] * 0.124) + (N6.f[1,scen]))
N6.f[1,scen] <- round(Nterr.f[1,scen]-(N4.f[1,scen] * 0.024) - (N5.f[1,scen] * 0.124))
Nterr.f[1,scen] <- Ntot[T.count]
N2wild.f[1,scen] <- 0
N2released.f[1,scen] <- 0
for (fut in 2:PROJECTION){
### probabilistic formulation
nestlings.f[fut,scen] <- (mean.fec * 0.5 * Nterr.f[fut,scen]) ### number of local recruits calculated as REDUCED fecundity times number of territorial pairs
N1.f[fut,scen] ~ dbin(phi.juv[scen],round(nestlings.f[fut-1,scen])) ### +rescued[fut] number of 1-year old survivors
N2wild.f[fut,scen] ~ dbin(phi.sec[scen],max(1,round(N1.f[fut-1,scen]))) ### number of 2-year old wild survivors
N2released.f[fut,scen] ~ dbin(phi.capt.juv[scen],round(capt.release[fut])) ### +rescued[fut] number of 1-year old survivors
N2.f[fut,scen] <- N2wild.f[fut,scen] + N2released.f[fut,scen] ### sum of the N1 cohort derived from wild and captive-bred juveniles
N3.f[fut,scen] ~ dbin(phi.third[scen],max(1,round(N2.f[fut-1,scen]))) ### number of 3-year old survivors
N4.f[fut,scen] ~ dbin(phi.ad[scen],max(1,round(N3.f[fut-1,scen]))) ### number of 4-year old survivors
N5.f[fut,scen] ~ dbin(phi.ad[scen],max(1,round(N4.f[fut-1,scen]))) ### number of 5-year old survivors
N6.f[fut,scen] ~ dbin(phi.ad[scen],max(1,round((N5.f[fut-1,scen]+N6.f[fut-1,scen])))) ### number of 6-year or older (adult) birds
Nterr.f[fut,scen] <- round((N4.f[fut,scen] * 0.024) + (N5.f[fut,scen] * 0.124) + (N6.f[fut,scen]))
} # fut
for (fut2 in 1:(PROJECTION-1)){
lambda.t.f[fut2,scen] <- Nterr.f[fut2+1,scen] / max(Nterr.f[fut2,scen],1)
loglambda.t.f[fut2,scen]<-log(lambda.t.f[fut2,scen])## for calculating geometric mean of overall population growth rate
} # fut2
#### FUTURE POPULATION GROWTH RATE #########
fut.lambda[scen]<-exp((1/(PROJECTION-1))*sum(loglambda.t.f[1:(PROJECTION-1),scen])) # Geometric mean
} # end of future projection scenarios
} # END of the JAGS model