-
Notifications
You must be signed in to change notification settings - Fork 0
/
LeopardFrog_Occupancy.R
331 lines (257 loc) · 12.4 KB
/
LeopardFrog_Occupancy.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
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
######################################################################################
# R code used to model metapopulation dynamics of lowland leopard frogs
# in southern Arizona based on 22 years of detection-nondetection data.
# Code adapted from Chandler et al. 2015 (Journal of Applied Ecology)
#
# Reference:
# Zylstra, E.R., D.E. Swann, B.R. Hossack, E. Muths, and R.J. Steidl. Drought-mediated
# extinction of an arid-land amphibian: insights from a spatially-explicit dynamic
# occupancy model. Ecological Applications.
######################################################################################
#--- Load libraries -----------------------------------------------------------------#
#install.packages(c("reshape2","jagsUI"))
library(reshape2)
library(jagsUI)
#--- Load data ----------------------------------------------------------------------#
#set working directory
setwd("C:/...")
rm(list=ls(all=TRUE))
#Data stored in Open Science Framework: https://doi.org/10.17605/OSF.IO/43PVJ
surveys <- read.table("LeopardFrog_Occupancy_SurveyData.csv", header=TRUE, sep=",", na.strings=c(NA,""))
site.cov <- read.table("LeopardFrog_Occupancy_SiteXSeasonCovariates.csv", header=TRUE, sep=",", na.strings=c(NA,""))
seas.cov <- read.table("LeopardFrog_Occupancy_SeasonCovariates.csv", header=TRUE, sep=",", na.strings=c(NA,""))
neighbors <- read.table("LeopardFrog_Occupancy_Neighbors.csv", header=FALSE, sep=",", na.strings=c(NA,""))
distances <- read.table("LeopardFrog_Occupancy_Distances.csv", header=FALSE, sep=",", na.strings=c(NA,""))
#--- Format and standardize survey-level covariates ---------------------------------#
#Index of surface-water availability
water.orig <- surveys$water
water <- (water.orig - mean(water.orig))/sd(water.orig)
#Indicator for fall sampling periods
fall <- surveys$fall
#Indicator for surveys led by less-experienced observers
obs <- surveys$inexp.obs
#--- Format and standardize site-level covariates -----------------------------------#
#Index of surface-water reliability
rel.orig <- site.cov$reliability
rel <- (rel.orig - mean(rel.orig))/sd(rel.orig)
#Area of catchment basin
area.orig <- site.cov$basin.area
area <- (area.orig - mean(area.orig))/sd(area.orig)
#Number of pools
npools.orig <- site.cov$n.pools
npools <- (npools.orig - mean(npools.orig))/sd(npools.orig)
#Elevation
elev.orig <- site.cov$elevation
elev <- (elev.orig - mean(elev.orig))/sd(elev.orig)
#--- Format and standardize covariates that vary among sites and seasons ------------#
#Palmer Drought Severity Index (PDSI), averaged over six months prior to each winter/summer
pdsi.orig <- site.cov[,grep('pdsi6mo',names(site.cov))]
#standardize using the mean and SD among canyon-level observations:
pdsi6mo.mn <- -1.797
pdsi6mo.sd <- 2.298
pdsi <- (pdsi.orig - pdsi6mo.mn)/pdsi6mo.sd
#Palmer Drought Severity Index (PDSI), averaged over prior larval period
pdsiJA.orig <- site.cov[,grep('pdsiJA',names(site.cov))]
#standardize using the mean and SD among canyon-level observations:
pdsiJA.mn <- -1.865
pdsiJA.sd <- 2.327
pdsiJA <- (pdsiJA.orig - pdsiJA.mn)/pdsiJA.sd
#Percent of 30-year precipitation normals during current winter/summer season
precip.orig <- site.cov[,grep('precip',names(site.cov))]
#standardize using the mean and SD among canyon-level observations:
precip.mn <- 88.629
precip.sd <- 36.444
precip <- (precip.orig - precip.mn)/precip.sd
#Sediment levels during each winter/summer season
sed <- site.cov[,grep('sed',names(site.cov))]
#--- Neighbor index and distance matrices ------------------------------------------------#
#Indices for neighbors of each site
nsi <- as.matrix(neighbors)
#Squared distances between a site and each of its neighbors
distSq <- as.matrix(distances^2)
#--- User-defined functions ---------------------------------------------------------#
#Function to create a matrix with information about known latent states, z[i,t]
known.state.dynocc <- function(obsarray){
state <- apply(obsarray,c(1,3),function(x) ifelse(sum(is.na(x))==length(x),NA,max(x,na.rm=T)))
state[state==0] <- NA
return(state)
}
#Function to create a matrix with initial values for unknown latent states, z[i,t]
inits.state.dynocc <- function(obsarray){
state <- apply(obsarray,c(1,3),function(x) ifelse(sum(is.na(x))==length(x),NA,max(x,na.rm=T)))
#Initial value of 1 whenever occupancy state is unknown
state[state==1] <- 2
state[is.na(state) | state==0] <- 1
state[state==2] <- NA
return(state)
}
#Creating arrays of observed survey data:
y.array <- acast(surveys[,c('y','site','period','rep')], site ~ rep ~ period, value.var="y")
#--- Bundle data for JAGS -----------------------------------------------------------#
#Create matrix of detection covariates: season, water, observer
water2 <- water*water
cov.p <- cbind(fall, water, water2, obs)
ncov.p <- ncol(cov.p)
#create matrix of site covariates for initial occupancy: elev, npools, area
elev2 <- elev*elev
cov.psi <- cbind(elev, elev2, npools, area)
ncov.psi <- ncol(cov.psi)
#See LeopardFrog_Occupancy_JAGSdata.pdf for detailed decriptions of each object in jags.data
jags.data <- list(y=as.vector(surveys$y),
nSites=max(surveys$site),
nSites1=sum(site.cov$watershed=='South'),
nSampPeriod=max(surveys$period),
nObs=nrow(surveys),
ncov.p=ncov.p,
cov.p=cov.p,
ncov.psi=ncov.psi,
cov.psi=cov.psi,
site=as.vector(surveys$site),
period=as.vector(surveys$period),
distSq=distSq,
nNeighbors=as.vector(site.cov$n.Neighbors),
nsi=nsi,
winter=as.vector(seas.cov$winter),
pdsi=as.matrix(pdsi),
pdsiJA=as.matrix(pdsiJA),
precip=as.matrix(precip),
sed=as.matrix(sed),
rel=as.vector(rel),
area=as.vector(area),
z=known.state.dynocc(y.array))
#--- Parameters to monitor ----------------------------------------------------------#
params <- c('p0','psi0','eps0','rho0','theta',
'beta.p0','beta.psi0','beta.eps0','beta.rho0',
'beta.p','beta.psi',
'beta.eps1.win','beta.eps2.pdsi','beta.eps3.winpdsi',
'beta.eps4.sed','beta.eps5.area','beta.eps6.rel',
'beta.rho1.win','beta.rho2.pdsiJA','beta.rho3.winpdsiJA','beta.rho4.precip',
'beta.rho5.winprecip','beta.rho6.sed','beta.rho7.area','beta.rho8.rel',
'PAO')
#--- Initial values -----------------------------------------------------------------#
inits <- function(){list(beta.p0=runif(1,-2,2),beta.psi0=runif(1,-2,2),
beta.eps0=runif(1,-2,2),beta.rho0=runif(1,-2,2),
theta=runif(1,0,5),
beta.p=runif(ncov.p,-2,2),beta.psi=runif(ncov.psi,-2,2),
beta.eps1.win=runif(1,-2,2),
beta.eps2.pdsi=runif(1,-2,2),
beta.eps3.winpdsi=runif(1,-2,2),
beta.eps4.sed=runif(1,-2,2),
beta.eps5.area=runif(1,-2,2),
beta.eps6.rel=runif(1,-2,2),
beta.rho1.win=runif(1,-2,2),
beta.rho2.pdsiJA=runif(1,-2,2),
beta.rho3.winpdsiJA=runif(1,-2,2),
beta.rho4.precip=runif(1,-2,2),
beta.rho5.winprecip=runif(1,-2,2),
beta.rho6.sed=runif(1,-2,2),
beta.rho7.area=runif(1,-2,2),
beta.rho8.rel=runif(1,-2,2),
z=inits.state.dynocc(y.array))}
#--- Create JAGS model file ---------------------------------------------------------#
sink("LeopardFrog_Occupancy_JAGSmodel.txt")
cat("
model{
#Priors
# Detection
beta.p0 ~ dlogis(0,1)
for(i in 1:ncov.p){
beta.p[i] ~ dnorm(0,0.1)
} #i
# Initial occupancy
beta.psi0 ~ dlogis(0,1)
for(i in 1:ncov.psi){
beta.psi[i] ~ dnorm(0,0.1)
} #i
# Extinction
beta.eps0 ~ dlogis(0,1)
beta.eps1.win ~ dnorm(0,0.1)
beta.eps2.pdsi ~ dnorm(0,0.1)
beta.eps3.winpdsi ~ dnorm(0,0.1)
beta.eps4.sed ~ dnorm(0,0.1)
beta.eps5.area ~ dnorm(0,0.1)
beta.eps6.rel ~ dnorm(0,0.1)
# Baseline colonization
beta.rho0 ~ dlogis(0,1)
beta.rho1.win ~ dnorm(0,0.1)
beta.rho2.pdsiJA ~ dnorm(0,0.1)
beta.rho3.winpdsiJA ~ dnorm(0,0.1)
beta.rho4.precip ~ dnorm(0,0.1)
beta.rho5.winprecip ~ dnorm(0,0.1)
beta.rho6.sed ~ dnorm(0,0.1)
beta.rho7.area ~ dnorm(0,0.1)
beta.rho8.rel ~ dnorm(0,0.1)
theta ~ dunif(0,15)
#State model
for(i in 1:nSites){
logit(psi[i]) <- beta.psi0 + inprod(beta.psi[], cov.psi[i,])
z[i,1] ~ dbern(psi[i])
for(t in 2:nSampPeriod){
# Extinction probability at site i between sampling period t-1 and t
logit(eps[i,t-1]) <- beta.eps0 +
beta.eps1.win*winter[t-1] +
beta.eps2.pdsi*pdsi[i,t-1] +
beta.eps3.winpdsi*winter[t-1]*pdsi[i,t-1] +
beta.eps4.sed*sed[i,t-1] +
beta.eps5.area*area[i] +
beta.eps6.rel*rel[i]
# Baseline colonization probability (probability a site is colonized by a frog
# from a coincident site) between sampling period t-1 and t:
logit(rho[i,t-1]) <- beta.rho0 +
beta.rho1.win*winter[t-1] +
beta.rho2.pdsiJA*pdsiJA[i,t-1] +
beta.rho3.winpdsiJA*winter[t-1]*pdsiJA[i,t-1] +
beta.rho4.precip*precip[i,t-1] +
beta.rho5.winprecip*winter[t-1]*precip[i,t-1] +
beta.rho6.sed*sed[i,t-1] +
beta.rho7.area*area[i] +
beta.rho8.rel*rel[i]
for(n in 1:nNeighbors[i]){
# Pairwise colonization probability
gammapair[i,n,t-1] <- rho[i,t-1]*exp(-distSq[i,n]/(2*theta^2))*z[nsi[i,n],t-1]
# We used a Gaussian dispersal kernel, but a different kernel could be
# specified by replacing exp(-distSq[i,n]/(2*theta^2)) with another kernel
# (e.g., exponential: lambda*exp(-dist[i,n]*lambda)
gammapairsinv[i,n,t-1] <- 1 - gammapair[i,n,t-1]
} #n
# Colonization probability
gamma[i,t-1] <- 1 - prod(gammapairsinv[i,1:nNeighbors[i],t-1])
# Probability of being occupied at time t
# Use the following if you don't want to include a rescue effect:
# Ez[i,t-1] <- gamma[i,t-1]*(1-z[i,t-1])+(1-eps[i,t-1])*z[i,t-1]
# Use the following if you want to include a pseudo-rescue effect:
Ez[i,t-1] <- gamma[i,t-1]*(1-z[i,t-1])+(1-eps[i,t-1]*(1-gamma[i,t-1]))*z[i,t-1]
z[i,t] ~ dbern(Ez[i,t-1])
} #t
} #i
#Observation model (data in long form)
for(w in 1:nObs){
logit(p[w]) <- beta.p0 + inprod(beta.p[], cov.p[w,])
y[w] ~ dbern(z[site[w],period[w]]*p[w])
} #w
#Derived parameters
logit(p0) <- beta.p0
logit(psi0) <- beta.psi0
logit(eps0) <- beta.eps0
logit(rho0) <- beta.rho0
for(t in 1:nSampPeriod){
# Estimated occupancy in the southern watershed
PAO[1,t] <- sum(z[1:nSites1,t])/nSites1
# Estimated occupancy in the northern watershed
PAO[2,t] <- sum(z[(nSites1+1):nSites,t])/(nSites-nSites1)
} #t
} #model
",fill=TRUE)
sink()
#--- Run model in JAGS --------------------------------------------------------------#
ni <- 83000
nt <- 40
nc <- 3
nb <- 3000
na <- 2000
out <- jags(data=jags.data, inits=inits, parameters.to.save=params,
model.file="LeopardFrog_Occupancy_JAGSmodel.txt",
n.chains=nc, n.adapt=na, n.iter=ni, n.burnin=nb, n.thin=nt,
parallel=T)
print(out)
plot(out,ask=T)