# Supplemental Material

# Designing a long-term occupancy monitoring plan for a cryptic 
# reptile

# This is the R Code for simulating occupancy survey data given underlying 
# changes in true occupancy, while in the presence of process (desert 
# tortoise availability for detection) and sampling noise (detection 
# probability < 1.0).

library(unmarked)
library(truncnorm)


# Set your working directory.  This is where ‘Full_simulations_parms.txt’
# is stored and where all simulated datasets will be exported to. See below
# for example format for Full_simulations_parms.txt.


setwd("C:/Users/…")

# Begin simulation code.  Creates blank datasets to be filled in by the 
# for-loops, generates simulation parameters for all combinations in
# ‘Full_simulation_parms.txt’, simulates the datasets, then exports to 
# the working directory.

 combinations.dp=read.delim("Full_simulation_parms.txt")
 
 n.row=nrow(combinations.dp)
 names(combinations.dp)
 niter=150    # The number of desired datasets to create for each combination
 
 psi.10=rep(NA,niter)
 smoothed.occu=rep(NA,niter)
 smoothed.se=rep(NA,niter)
 n.years=10
 psi0=0.5   #Sets original psi in year 1 as 0.5 for all datasets
 psi <- rep(NA,n.years)
 psi.true <- rep(NA,n.years)
 psi.logit <- rep(NA,n.years)
 psi.logit.noise <- rep(NA,n.years)
 psi[1] <- psi0
 psi.true[1] <- psi0
ann.noise <- rep(NA,n.years-1)
year.app <- rep(NA,n.years)

 for(m in 1:n.row){
     n.sites=combinations.dp$n.sites[m]
     n.visits=combinations.dp$n.visits[m]
     r.change=combinations.dp$r.change[m]
     name=paste(as.character(n.sites),"Nsites_", as.character(n.visits),"Nvisits_",as.character(r.change),"Rchange",sep="")
 
 site <- 1:n.sites
 year <- 1:n.years
 change <-1:r.change

 muZ <- z <- array(dim=c(n.sites, n.years))
 y <- array(NA, dim=c(n.sites, n.visits, n.years))
 prob <- array(NA, dim=c(n.sites,n.years))


# Generate detection probability p values from a uniform distribution,
# with the bounds determined by the observed range of p from 
# Harju & Cambrin (2019)
for(t in 1:niter){
 p <- runif(n = n.years, min = 0.185, max = 0.734)   
 

 #Generate latent occupancy states for all sites
 #First year
 z[,1] <- rbinom(n.sites,1,psi[1])
 
 #Later years
 for(i in 1:n.sites){
     for(k in 2:n.years){
	# Use truncated normal dist based on observed deviation in
	# occupancy from Harju & Cambrin (2019)
         ann.noise[k] <- rtruncnorm(n=1,a=-1.53,b=1.647, mean=0, sd=1.063) 
         psi.true[k] <-psi.true[k-1]*r.change
         psi.logit[k] <- log(psi.true[k]/(1-psi.true[k]))
         psi.logit.noise[k] <- psi.logit[k]+ann.noise[k]
         psi[k] <- exp(psi.logit.noise[k])/(1+exp(psi.logit.noise[k]))
         z[i,k] <- rbinom(1,1,psi[k])
     }
 }
 
 
 #Generate detection/non-detection data
 for(i in 1:n.sites){
     for(k in 1:n.years){
         # p[k] <- runif(1,0.185,0.734)
         prob[i,k] <- z[i,k] * p[k]
         for(j in 1:n.visits){
             y[i,j,k] <- rbinom(1, 1, prob[i,k])
         }
     }
 }
 
 # Calculate annual apparent occupancy (including both process noise and 
 # sampling noise from rbinom function.  psi[k] has the latent occupancy
 # rate with process noise only)
 for(i in 1:n.years){
     year.app[i] <- mean(z[,i])
 }
 
 # Run dynamic occupancy model to analyze occupancy trend 
 # using package ‘unmarked’
 yy <- matrix (y, n.sites, n.visits*n.years)
 year <- matrix(c('01','02','03','04','05','06','07','08','09','10'),nrow(yy), n.years, byrow=TRUE)
 
 simUMF <- unmarkedMultFrame(y = yy,numPrimary=n.years)
 m0 <- colext(psiformula= ~1, gammaformula = ~ 1, epsilonformula = ~ 1,pformula = ~ 1, data = simUMF, method="BFGS")
 m0 <- nonparboot(m0, B = 100)  # 100 nonparametric bootstrap estimates of se
 #print(name)
 #print(paste("Iteration Number:",i))
 #print(paste("Row number:",m))
 
 # Extract estimated proportion of sites occupied and boostrapped se
 m0.out <- data.frame(cbind(n.sites=rep(n.sites,n.years),n.visits=rep(n.visits,n.years),r.change=rep(r.change,n.years),psi=psi, smoothed=m0@smoothed.mean[2,], SE=m0@smoothed.mean.bsse[2,],ann.noise=ann.noise))

write.table(m0.out, file=paste(name, "occu",t, "txt", sep=".")) 

# Extract estimated colonization and extinction rates
coloniz <- summary(m0@estimates@estimates$col)
extinct <- summary(m0@estimates@estimates$ext)

m0.colext <- data.frame(cbind(n.sites=n.sites,n.visits=n.visits,r.change=r.change,col.est=coloniz[,1],col.se=coloniz[,2], col.pval=coloniz[,4],ext.est=extinct[,1],ext.se=extinct[,2],ext.pval=extinct[,4]))

write.table(m0.colext, file=paste(name, "colext",t, "txt", sep="."))
}}


# Code for importing all .txt files in a directory, combining into a single
# file, then writing out the single file back to the directory

out.file.occu <- ""
 file.names.occu <- Sys.glob("*occu*.txt")
 
 for(i in 1:length(file.names.occu)){
     file.occu <- read.table(file.names.occu[i], header=TRUE, sep=" ", stringsAsFactors=FALSE)
     out.file.occu <- rbind(out.file.occu, file.occu)
 }

# out.file may have a blank first row, need to remove
sims.out.occu <- out.file.occu #[-c(1),]

# Make sure all columns are numeric
sims.out.occu[] <- lapply(sims.out.occu,function(x) as.numeric(x))
# Data check
str(sims.out.occu)

# Calculate confidence intervals
sims.out.occu$model.psi <- as.numeric(sims.out.occu$smoothed)
sims.out.occu$SE <- as.numeric(sims.out.occu$SE)
sims.out.occu$l95cl <- sims.out.occu$model.psi - (sims.out.occu$SE*1.96)
sims.out.occu$u95cl <- sims.out.occu$model.psi + (sims.out.occu$SE*1.96)
sims.out.occu$l90cl <- sims.out.occu$model.psi - (sims.out.occu$SE*1.645)
sims.out.occu$u90cl <- sims.out.occu$model.psi + (sims.out.occu$SE*1.645)
sims.out.occu$l85cl <- sims.out.occu$model.psi - (sims.out.occu$SE*1.44)
sims.out.occu$u85cl <- sims.out.occu$model.psi + (sims.out.occu$SE*1.44)


# Now load and combine the colonization and extinction estimates
out.file.colext <- ""
 file.names.colext <- Sys.glob("*colext*.txt")


for(i in 1:length(file.names.colext)){
     file.colext <- read.table(file.names.colext[i],header=TRUE, sep=" ", stringsAsFactors=FALSE)
     out.file.colext <- rbind(out.file.colext, file.colext)
 }

# out.file may have a blank first row, need to remove
sims.out.colext <- out.file.colext #[-c(1),]

sims.out.colext[] <- lapply(sims.out.colext,function(x) as.numeric(x))

# Data check
str(sims.out.colext)

# Calculate confidence intervals
sims.out.colext$col.l95cl <- sims.out.colext$col.est - (sims.out.colext$col.se*1.96)
sims.out.colext$col.u95cl <- sims.out.colext$col.est + (sims.out.colext$col.se*1.96)
sims.out.colext$col.l90cl <- sims.out.colext$col.est - (sims.out.colext$col.se *1.645)
sims.out.colext$col.u90cl <- sims.out.colext$col.est + (sims.out.colext$col.se *1.645)
sims.out.colext$col.l85cl <- sims.out.colext$col.est - (sims.out.colext$col.se *1.44)
sims.out.colext$col.u85cl <- sims.out.colext$col.est + (sims.out.colext$col.se *1.44)
 
sims.out.colext$ext.l95cl <- sims.out.colext$ext.est - (sims.out.colext$ext.se*1.96)
sims.out.colext$ext.u95cl <- sims.out.colext$ext.est + (sims.out.colext$ext.se*1.96)
sims.out.colext$ext.l90cl <- sims.out.colext$ext.est - (sims.out.colext$ext.se *1.645)
sims.out.colext$ext.u90cl <- sims.out.colext$ext.est + (sims.out.colext$ext.se *1.645)
sims.out.colext$ext.l85cl <- sims.out.colext$ext.est - (sims.out.colext$ext.se *1.44)
sims.out.colext$ext.u85cl <- sims.out.colext$ext.est + (sims.out.colext$ext.se *1.44)

# Derived patterns

# Occupancy estimates
# Create a new variable using mutate that follows logical ifelse arguments   # for the occu estimates
library(plyr)
sims.out.occu <- mutate(sims.out.occu,r.change.type=ifelse(r.change < 1,"decreasing", ifelse(r.change==1,"stable","increasing")))

# Calculate 0/1 to see if CI’s exclude original psi
sims.out.occu <- mutate(sims.out.occu,ci.95.exclude.origpsi=ifelse(r.change.type=="decreasing", ifelse(u95cl>0.5,0,1), ifelse(r.change.type=="increasing",ifelse(l95cl<0.5,0,1),ifelse(l95cl<0.5&u95cl>0.5,0,1))))  

sims.out.occu <- mutate(sims.out.occu,ci.90.exclude.origpsi=ifelse(r.change.type=="decreasing", ifelse(u90cl>0.5,0,1), ifelse(r.change.type=="increasing",ifelse(l90cl<0.5,0,1),ifelse(l90cl<0.5&u90cl>0.5,0,1))))  

sims.out.occu <- mutate(sims.out.occu,ci.85.exclude.origpsi=ifelse(r.change.type=="decreasing", ifelse(u85cl>0.5,0,1), ifelse(r.change.type=="increasing",ifelse(l85cl<0.5,0,1),ifelse(l85cl<0.5&u85cl>0.5,0,1))))  

# Calculate true psi as function of modeled annual rate of change
sims.out.occu <- mutate(sims.out.occu,true.psi=ifelse(r.change==0.96,0.5*0.96^9,ifelse(r.change==0.97,0.5*0.97^9,ifelse(r.change==0.98,0.5*0.98^9,ifelse(r.change==0.99,0.5*0.99^9,ifelse(r.change==1,0.5,ifelse(r.change==1.01,0.5*1.01^9,ifelse(r.change==1.02,0.5*1.02^9,ifelse(r.change==1.03,0.5*1.03^9,ifelse(r.change==1.04,0.5*1.04^9,NA))))))))))


# Colonization and extinction estimates

# Add column labeling R.change as decreasing, stable, or increasing
sims.out.colext <- mutate(sims.out.colext,r.change.type=ifelse(r.change < 1,"decreasing", ifelse(r.change==1,"stable","increasing")))

# Calculate 0/1 to see if CI’s exclude col or ext = 0.0
# col
sims.out.colext <- mutate(sims.out.colext,ci.95.col.exclude0=ifelse(r.change.type=="decreasing", ifelse(col.u95cl>0,0,1), ifelse(r.change.type=="increasing",ifelse(col.l95cl<0,0,1),ifelse(col.l95cl<0&col.u95cl>0,0,1))))  

sims.out.colext <- mutate(sims.out.colext,ci.90.col.exclude0=ifelse(r.change.type=="decreasing", ifelse(col.u90cl>0,0,1), ifelse(r.change.type=="increasing",ifelse(col.l90cl<0,0,1),ifelse(col.l90cl<0&col.u90cl>0,0,1))))  

sims.out.colext <- mutate(sims.out.colext,ci.85.col.exclude0=ifelse(r.change.type=="decreasing", ifelse(col.u85cl>0,0,1), ifelse(r.change.type=="increasing",ifelse(col.l85cl<0,0,1),ifelse(col.l85cl<0&col.u85cl>0,0,1))))  

# ext
sims.out.colext <- mutate(sims.out.colext,ci.95.ext.exclude0=ifelse(r.change.type=="decreasing", ifelse(ext.u95cl>0,0,1), ifelse(r.change.type=="increasing",ifelse(ext.l95cl<0,0,1),ifelse(ext.l95cl<0&ext.u95cl>0,0,1))))  

sims.out.colext <- mutate(sims.out.colext,ci.90.ext.exclude0=ifelse(r.change.type=="decreasing", ifelse(ext.u90cl>0,0,1), ifelse(r.change.type=="increasing",ifelse(ext.l90cl<0,0,1),ifelse(ext.l90cl<0&ext.u90cl>0,0,1))))  

sims.out.colext <- mutate(sims.out.colext,ci.85.ext.exclude0=ifelse(r.change.type=="decreasing", ifelse(ext.u85cl>0,0,1), ifelse(r.change.type=="increasing",ifelse(ext.l85cl<0,0,1),ifelse(ext.l85cl<0&ext.u85cl>0,0,1))))  




# Input parameters
# This is the sample input code for the file ‘Full_simulation_parms.txt’ that
# tells the simulation loop what combination of variables to assess.  This 
# needs to be a .txt file with all of these description lines removed (i.e., # the very first line of the .txt file is the column headers below). Be sure # to remove the ellipse from the last row (included here to truncate the
# parameter list for display purposes).

id	n.sites	n.visits	r.change
1	10	3	1.03
2	10	5	1.03
3	10	7	1.03
4	20	3	1.03
5	20	5	1.03
6	20	7	1.03
7	30	3	1.03
8	30	5	1.03
9	30	7	1.03
10	40	3	1.03
11	40	5	1.03
12	40	7	1.03
13	50	3	1.03
14	50	5	1.03
15	50	7	1.03
16	60	3	1.03
17	60	5	1.03
18	60	7	1.03
19	70	3	1.03
20	70	5	1.03
21	70	7	1.03
22	80	3	1.03
23	80	5	1.03
24	80	7	1.03
25	90	3	1.03
26	90	5	1.03
27	90	7	1.03
28	100	3	1.03
29	100	5	1.03
30	100	7	1.03
31	10	3	1.02
32	10	5	1.02
33	10	7	1.02
34	20	3	1.02
35	20	5	1.02
36	20	7	1.02
37	30	3	1.02
38	30	5	1.02
39	30	7	1.02
40	40	3	1.02
41	40	5	1.02
42	40	7	1.02
…