-
Notifications
You must be signed in to change notification settings - Fork 0
/
btgo.jags.txt
121 lines (108 loc) · 4.42 KB
/
btgo.jags.txt
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
######################################################################
# #
# Supplement to: #
# Murray et al (2017) Bar-tailed Godwit Analysis #
# #
# Created on 4/4/2017 by Nicholas Murray & Colin Studds #
# ------------------------------------------------------------------ #
# #
# Description: #
# JAGS code for Bar-tailed Godwit population analysis. #
# #
# Note: we use nested indexing [sub] to simultaneously fit the model #
# to the two subspecies considered in our analysis #
# #
# ------------------------------------------------------------------ #
# #
# The analysis was performed with: #
# JAGS 3.4.0 (64-bit) #
# R 3.2.2 (Fire Safety) #
# R library R2jags #
# #
######################################################################
model{
## PRIORS
for (s in 1:nsites){
ste[s] ~ dnorm(0, tau.site)
}
for (b in 1:nsub){
alpha[b] ~ dnorm(mu.alpha,tau.alpha)
trend[b] ~ dnorm(mu.trend,tau.trend)
for (c in 1:ncovs){
theta[b,c] ~ dnorm(mu.theta, tau.theta)
}
}
for (b in 1:nsub){
for (c in 1:ncovs){
ind[b,c] ~ dbern(0.5) # not included in final model
beta[b,c] <- theta[b,c] * ind[b,c] # modified in final model
}
}
for (k in 1:nyears){
beta.d[k] ~ dnorm(0,0.001)
}
## LIKELIHOOD
### Model for abundance
for (i in 1:nsites){
for (k in 1:nyears){
eps[i,k] ~ dnorm(0, tau.eps)
N[i,k] ~ dpois(lambda[i,k])
log(lambda[i,k]) <- ste[site[i]] +
alpha[sub[i]] +
trend[sub[i]] * year[k] +
beta[sub[i],1] * st.chl.std[i,k] + # not included in final model
beta[sub[i],2] * br.chl.std[i,k] + # not included in final model
beta[sub[i],3] * nb.chl.std[i,k] +
beta[sub[i],4] * st.rai.std[i,k] + # not included in final model
beta[sub[i],5] * nb.rai.std[i,k] + # not included in final model
beta[sub[i],6] * st.tem.std[i,k] +
beta[sub[i],7] * br.tem.std[i,k] +
eps[i,k]
## Observation model
for (j in 1:ncounts){
y[i,j,k] ~ dbin(p[i,j,k], N[i,k])
p[i,j,k] <- dunif(0,1)
expect[i,j,k] <- p[i,j,k] * N[i,k]
E[i,j,k] <- pow((y[i,j,k] - expect[i,j,k]),2) / (expect[i,j,k]+0.5)
y.new[i,j,k] ~ dbin(p[i,j,k], N[i,k])
E.new[i,j,k] <- pow((y.new[i,j,k] - expect[i,j,k]),2) / (expect[i,j,k]+0.5)
}
site.year.p[i,k] <- mean(p[i,,k])
}
}
for (k in 1:nyears){
total.N[k] <- sum(N[,k])
men.N[k] <- sum(N[1:5,k])
bau.N[k] <- sum(N[6:21,k])
detection[k] <- mean(site.year.p[,k])
men.detection[k] <- mean(site.year.p[1:5,k])
bau.detection[k] <- mean(site.year.p[6:21,k])
}
## DERIVED PARAMETERS
### Goodness of fit overall
fit <- sum(E[,,])
fit.new <- sum(E.new[,,])
### Goodness of fit
men.fit <- sum(E[1:5,,])
men.fit.new <- sum(E.new[1:5,,])
bau.fit <- sum(E[6:21,,])
bau.fit.new <- sum(E.new[6:21,,])
### Mean detection probability
mean.detection <- mean(detection[]) #overall
men.mean.detection <- mean(men.detection[]) #menzbieri
bau.mean.detection <- mean(bau.detection[]) #baueri
## HYPERPRIORS
mu.theta ~ dnorm(0,0.001)
tau.theta ~ dgamma(0.001,0.001)
sd.theta <- 1/sqrt(tau.theta)
tau.site ~ dgamma(0.001,0.001)
tau.alpha ~ dgamma(0.001,0.001)
tau.trend ~ dgamma(0.001,0.001)
tau.eps ~ dgamma(0.001,0.001)
sd.site <- 1 / sqrt(tau.site)
sd.alpha <- 1/sqrt(tau.alpha)
sd.trend <- 1/sqrt(tau.trend)
sd.eps <- 1/sqrt(tau.eps)
mu.alpha ~ dnorm(0,0.001)
mu.trend ~ dnorm(0,0.001)
}