#R code for (1) simulating 400 datasets of bird occurrence across three temporal replicates within a closed period at 50 survey sites for a community of 40 species of birds when species detection probabilities declined with the site covariate and averaged ~0.65 (scenario 3 in text) and (2) analyzing these datasets to evaluate four estimators of local species richness: raw species counts (naïve richness), second order jackknife, Chao, and multi-species occupancy model.
#It's not elegant but it will work.
#Packages needed
library(R2jags)
library(SPECIES)

############ LOOP FOR REPLICATING DATASETS AND ASSESSING BIAS ##########
#file full analysis p~site covariate neg
#abundance parameters per original code set to range from 0.14-7.4 individuals per sp per site increasing or decreasing with siteCov
#alpha.lambda = -2 beta.lambda = 4 (increasing abundance)
#alpha.lambda = 2 beta.lambda = -4 (decreasing abundance)
#detection parameters set so p decreases from 0.88-0.12 with siteCov
#alpha.p = 2 beta.p = -4

nsim <- 2 # how many times to simulate dataset and run model

# Create a bunch of empty vectors to store results from replicated datasets

m.bias.naive <- vector("list",nsim)
sd.bias.naive <- vector("list",nsim)
baye.p.naive <- vector("list",nsim)
m.bias.jack <- vector("list",nsim)
sd.bias.jack <- vector("list",nsim)
baye.p.jack <- vector("list",nsim)
m.bias.chao <- vector("list",nsim)
sd.bias.chao <- vector("list",nsim)
baye.p.chao <- vector("list",nsim)
m.bias.msocc <- vector("list",nsim)
sd.bias.msocc <- vector("list",nsim)
baye.p.msocc <- vector("list",nsim)

# Means of Model intercepts and effect sizes (beta) of each estimator vs. siteCov (to assess bias)
m.alpha.t <- vector("list",nsim)
m.alpha.n <- vector("list",nsim)
m.alpha.j <- vector("list",nsim)
m.alpha.c <- vector("list",nsim)
m.alpha.msocc <- vector("list",nsim)
m.beta.t <- vector("list",nsim)
m.beta.n <- vector("list",nsim)
m.beta.j <- vector("list",nsim)
m.beta.c <- vector("list",nsim)
m.beta.msocc <- vector("list",nsim)

p.bias.msocc <- vector("list",nsim)

# Variability of Model intercepts and effect sizes (beta) of each estimator vs. siteCov
sd.alpha.t <- vector("list",nsim)
sd.alpha.n <- vector("list",nsim)
sd.alpha.j <- vector("list",nsim)
sd.alpha.c <- vector("list",nsim)
sd.alpha.msocc <- vector("list",nsim)
sd.beta.t <- vector("list",nsim)
sd.beta.n <- vector("list",nsim)
sd.beta.j <- vector("list",nsim)
sd.beta.c <- vector("list",nsim)
sd.beta.msocc <- vector("list",nsim)

#Deviance of summary models for each estimator (spread of residuals for glm(richness~siteCov))
#for each simulation collect the deviance value for each of the 4 post-processing regressions
dev.t <- vector("list", nsim)
dev.n <- vector("list", nsim)
dev.j <- vector("list", nsim)
dev.c <- vector("list", nsim)
dev.msocc <- vector("list", nsim)

system.time(for (k in 1:nsim){      # Stick the simulations in a for loop and replicate them nsim times
  
  #data are arranged in a 3-d matrix 
  n.site = 50
  T = 3
  n.species = 40
  siteCov = (sort(runif(n.site,0,1))) #Create a site covariate  
  
  #It would shorten the code to design a loop to simulate the species observations; but here I simulated each species'
  #true site-specific abundance, occurance, and build individual observation matrices (Site x T),
  #then combine them all into a n.site x 3 x n.species data array for analysis with JAGS
  
  #70% of the species (n=28) will have a positive relationship with the the covariate, the rest negative
  
  #Species 1. Positive relationship of abundance with SiteCov and negative relationship with detection
  
  #Abundance
  alpha.lambda1 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda1 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam1 <- exp(alpha.lambda1+beta.lambda1*siteCov)
  
  N1 <- rpois(n = n.site, lambda = lam1) # give it some noise
  
  table(N1)        # Distribution of abundances across sites
  tru.occ1<-sum(N1 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p1 <- rnorm(1,2,0.1)     	# Intercept
  beta.p1 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob1 <- exp(alpha.p1 + beta.p1 * siteCov) / (1 + exp(alpha.p1 + beta.p1 * siteCov))
  
  #Create observations
  y1 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y1[,j]<- rbinom(n=n.site, size=N1, prob=det.prob1)
  }
  #naive occupancy
  m<- apply(y1,1,max)
  m[m>1]<-1
  naive.occ1<-mean(m)
  
  occ1 <- rep(0,n.site) #create a vector 40 units long and fill it with 0's
  occ1[N1>=1]<-1
  
  #----------------------species 1 end-----------------------------
  
  #Species 2.
  
  #Abundance
  alpha.lambda2 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda2 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam2 <- exp(alpha.lambda2+beta.lambda2*siteCov)
  
  N2 <- rpois(n = n.site, lambda = lam2) # give it some noise
  
  table(N2)      	# Distribution of abundances across sites
  (tru.occ2<-sum(N2 > 0) / n.site	)		# Empirical occupancy
  
  #Detection
  alpha.p2 <- rnorm(1,1,0.1) 			# Intercept
  beta.p2 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob2 <- exp(alpha.p2 + beta.p2 * siteCov) / (1 + exp(alpha.p2 + beta.p2 * siteCov))
  
  #Create observations
  y2 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y2[,j]<- rbinom(n=n.site, size=N2, prob=det.prob2)
  }
  #naive occupancy
  m<- apply(y2,1,max)
  m[m>1]<-1
  naive.occ2<-mean(m)
  
  occ2 <- rep(0,n.site)
  occ2[N2>=1]<-1
  #----------------------species 2 end-----------------------------
  
  #Species 3. Negative relationship of abundance with SiteCov and negative relationship with detection
  
  #Abundance
  alpha.lambda3 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda3 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam3 <- exp(alpha.lambda3+beta.lambda3*siteCov)
  
  N3 <- rpois(n = n.site, lambda = lam3) # give it some noise
  
  table(N3)        # Distribution of abundances across sites
  (tru.occ3<-sum(N3 > 0) / n.site	)		# Empirical occupancy
  
  #Detection
  alpha.p3 <- rnorm(1,2,0.1) 			# Intercept
  beta.p3 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob3 <- exp(alpha.p3 + beta.p3 * siteCov) / (1 + exp(alpha.p3 + beta.p3 * siteCov))
  
  #Create observations
  y3 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y3[,j]<- rbinom(n=n.site, size=N3, prob=det.prob3)
  }
  #naive occupancy
  m<- apply(y3,1,max)
  m[m>1]<-1
  naive.occ3<-mean(m)
  
  occ3<- rep(0,n.site)
  occ3[N3>=1]<-1
  
  #------------End species 3---------------------------------------
  
  #Species 4
  #Abundance
  alpha.lambda4 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda4 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam4 <- exp(alpha.lambda4+beta.lambda4*siteCov)
  
  N4 <- rpois(n = n.site, lambda = lam4) # give it some noise
  
  table(N4)      	# Distribution of abundances across sites
  tru.occ4<-sum(N4 > 0) / n.site			# Empirical occupancy
  
  #Detection
  alpha.p4 <- rnorm(1,2,0.1) 			# Intercept
  beta.p4 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob4 <- exp(alpha.p4 + beta.p4 * siteCov) / (1 + exp(alpha.p4 + beta.p4 * siteCov))
  
  #Create observations
  y4 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y4[,j]<- rbinom(n=n.site, size=N4, prob=det.prob4)
  }
  #naive occupancy
  m<- apply(y4,1,max)
  m[m>1]<-1
  naive.occ4<-mean(m)
  
  occ4<- rep(0,n.site)
  occ4[N4>=1]<-1
  #--------------------End species 4-------------------------------------------
  
  #Species 5
  #Abundance
  alpha.lambda5 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda5 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam5 <- exp(alpha.lambda5+beta.lambda5*siteCov)
  
  N5 <- rpois(n = n.site, lambda = lam5) # give it some noise
  
  table(N5)        # Distribution of abundances across sites
  tru.occ5<-sum(N5 > 0) / n.site			# Empirical occupancy
  
  #Detection
  alpha.p5 <- rnorm(1,2,0.1) 			# Intercept
  beta.p5 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob5 <- exp(alpha.p5 + beta.p5 * siteCov) / (1 + exp(alpha.p5 + beta.p5 * siteCov))
  
  #Create observations
  y5 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y5[,j]<- rbinom(n=n.site, size=N5, prob=det.prob5)
  }
  #naive occupancy
  m<- apply(y5,1,max)
  m[m>1]<-1
  naive.occ5<-mean(m)
  
  occ5<- rep(0,n.site)
  occ5[N5>=1]<-1
  #----------------End species 5---------------------------------------
  
  #Species 6
  #Abundance
  alpha.lambda6 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda6 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam6 <- exp(alpha.lambda6+beta.lambda6*siteCov)
  
  N6 <- rpois(n = n.site, lambda = lam6) # give it some noise
  
  table(N6)        # Distribution of abundances across sites
  tru.occ6<-sum(N6 > 0) / n.site  		# Empirical occupancy
  
  #Detection
  alpha.p6 <- rnorm(1,2,0.1) 			# Intercept
  beta.p6 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob6 <- exp(alpha.p6 + beta.p6 * siteCov) / (1 + exp(alpha.p6 + beta.p6 * siteCov))
  
  #Create observations
  y6 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y6[,j]<- rbinom(n=n.site, size=N6, prob=det.prob6)
  }
  #naive occupancy
  m<- apply(y6,1,max)
  m[m>1]<-1
  naive.occ6<-mean(m)
  
  #true occurance and each site
  occ6<- rep(0,n.site)
  occ6[N6>=1]<-1
  #----------------------------------End species 6-------------------------------------
  
  #Species 7
  #Abundance
  alpha.lambda7 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda7 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam7 <- exp(alpha.lambda7+beta.lambda7*siteCov)
  
  N7 <- rpois(n = n.site, lambda = lam7) # give it some noise
  
  table(N7)        # Distribution of abundances across sites
  tru.occ7<-sum(N7 > 0) / n.site    	# Empirical occupancy
  
  #Detection
  alpha.p7 <- rnorm(1,2,0.1) 			# Intercept
  beta.p7 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob7 <- exp(alpha.p7 + beta.p7 * siteCov) / (1 + exp(alpha.p7 + beta.p7 * siteCov))
  
  #Create observations
  y7 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y7[,j]<- rbinom(n=n.site, size=N7, prob=det.prob7)
  }
  #naive occupancy
  m<- apply(y7,1,max)
  m[m>1]<-1
  naive.occ7<-mean(m)
  
  #true occurance and each site
  occ7<- rep(0,n.site)
  occ7[N7>=1]<-1
  #----------------------------------End species 7-------------------------------------
  #Species 8
  #Abundance
  alpha.lambda8 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda8 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam8 <- exp(alpha.lambda8+beta.lambda8*siteCov)
  
  N8 <- rpois(n = n.site, lambda = lam8) # give it some noise
  
  table(N8)        # Distribution of abundances across sites
  tru.occ8<-sum(N8 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p8 <- rnorm(1,2,0.1) 			# Intercept
  beta.p8 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob8 <- exp(alpha.p8 + beta.p8 * siteCov) / (1 + exp(alpha.p8 + beta.p8 * siteCov))
  
  #Create observations
  y8 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y8[,j]<- rbinom(n=n.site, size=N8, prob=det.prob8)
  }
  #naive occupancy
  m<- apply(y8,1,max)
  m[m>1]<-1
  naive.occ8<-mean(m)
  
  #true occurance and each site
  occ8<- rep(0,n.site)
  occ8[N8>=1]<-1
  #----------------------------------End species 8-------------------------------------
  
  #Species 9
  #Abundance
  alpha.lambda9 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda9 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam9 <- exp(alpha.lambda9+beta.lambda9*siteCov)
  
  N9 <- rpois(n = n.site, lambda = lam9) # give it some noise
  
  table(N9)        # Distribution of abundances across sites
  tru.occ9<-sum(N9 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p9 <- rnorm(1,2,0.1) 			# Intercept
  beta.p9 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob9 <- exp(alpha.p9 + beta.p9 * siteCov) / (1 + exp(alpha.p9 + beta.p9 * siteCov))
  
  #Create observations
  y9 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y9[,j]<- rbinom(n=n.site, size=N9, prob=det.prob9)
  }
  #naive occupancy
  m<- apply(y9,1,max)
  m[m>1]<-1
  naive.occ9<-mean(m)
  
  #true occurance and each site
  occ9<- rep(0,n.site)
  occ9[N9>=1]<-1
  #----------------------------------End species 9-------------------------------------
  
  #Species 10
  #Abundance
  alpha.lambda10 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda10 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam10 <- exp(alpha.lambda10+beta.lambda10*siteCov)
  
  N10 <- rpois(n = n.site, lambda = lam10) # give it some noise
  
  table(N10)        # Distribution of abundances across sites
  tru.occ10<-sum(N10 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p10 <- rnorm(1,2,0.1) 			# Intercept
  beta.p10 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob10 <- exp(alpha.p10 + beta.p10 * siteCov) / (1 + exp(alpha.p10 + beta.p10 * siteCov))
  
  #Create observations
  y10 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y10[,j]<- rbinom(n=n.site, size=N10, prob=det.prob10)
  }
  #naive occupancy
  m<- apply(y10,1,max)
  m[m>1]<-1
  naive.occ10<-mean(m)
  
  #true occurance and each site
  occ10<- rep(0,n.site)
  occ10[N10>=1]<-1
  #----------------------------------End species 10-------------------------------------
  #Species 11
  #Abundance
  alpha.lambda11 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda11 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam11 <- exp(alpha.lambda11+beta.lambda11*siteCov)
  
  N11 <- rpois(n = n.site, lambda = lam11) # give it some noise
  
  table(N11)        # Distribution of abundances across sites
  tru.occ11<-sum(N11 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p11 <- rnorm(1,2,0.1)   		# Intercept
  beta.p11 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob11 <- exp(alpha.p11 + beta.p11 * siteCov) / (1 + exp(alpha.p11 + beta.p11 * siteCov))
  
  #Create observations
  y11 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y11[,j]<- rbinom(n=n.site, size=N11, prob=det.prob11)
  }
  #naive occupancy
  m<- apply(y11,1,max)
  m[m>1]<-1
  naive.occ11<-mean(m)
  
  #true occurance and each site
  occ11<- rep(0,n.site)
  occ11[N11>=1]<-1
  #----------------------------------End species 11-----------
  #Species 12
  #Abundance
  alpha.lambda12 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda12 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam12 <- exp(alpha.lambda12+beta.lambda12*siteCov)
  
  N12 <- rpois(n = n.site, lambda = lam12) # give it some noise
  
  table(N12)        # Distribution of abundances across sites
  tru.occ12<-sum(N12 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p12 <- rnorm(1,2,0.1)   		# Intercept
  beta.p12 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob12 <- exp(alpha.p12 + beta.p12 * siteCov) / (1 + exp(alpha.p12 + beta.p12 * siteCov))
  
  #Create observations
  y12 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y12[,j]<- rbinom(n=n.site, size=N12, prob=det.prob12)
  }
  #naive occupancy
  m<- apply(y12,1,max)
  m[m>1]<-1
  naive.occ12<-mean(m)
  
  #true occurance and each site
  occ12<- rep(0,n.site)
  occ12[N12>=1]<-1
  #----------------------------------End species 12-----------
  
  #Species 13
  #Abundance
  alpha.lambda13 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda13 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam13 <- exp(alpha.lambda13+beta.lambda13*siteCov)
  
  N13 <- rpois(n = n.site, lambda = lam13) # give it some noise
  
  table(N13)        # Distribution of abundances across sites
  tru.occ13<-sum(N13 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p13 <- rnorm(1,2,0.1)   		# Intercept
  beta.p13 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob13 <- exp(alpha.p13 + beta.p13 * siteCov) / (1 + exp(alpha.p13 + beta.p13 * siteCov))
  
  #Create observations
  y13 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y13[,j]<- rbinom(n=n.site, size=N13, prob=det.prob13)
  }
  #naive occupancy
  m<- apply(y13,1,max)
  m[m>1]<-1
  naive.occ13<-mean(m)
  
  #true occurance and each site
  occ13<- rep(0,n.site)
  occ13[N13>=1]<-1
  #----------------------------------End species 13-----------
  
  #Species 14
  #Abundance
  alpha.lambda14 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda14 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam14 <- exp(alpha.lambda14+beta.lambda14*siteCov)
  
  N14 <- rpois(n = n.site, lambda = lam14) # give it some noise
  
  table(N14)        # Distribution of abundances across sites
  tru.occ14<-sum(N14 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p14 <- rnorm(1,2,0.1)   		# Intercept
  beta.p14 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob14 <- exp(alpha.p14 + beta.p14 * siteCov) / (1 + exp(alpha.p14 + beta.p14 * siteCov))
  
  #Create observations
  y14 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y14[,j]<- rbinom(n=n.site, size=N14, prob=det.prob14)
  }
  #naive occupancy
  m<- apply(y14,1,max)
  m[m>1]<-1
  naive.occ14<-mean(m)
  
  #true occurance and each site
  occ14<- rep(0,n.site)
  occ14[N14>=1]<-1
  #----------------------------------End species 14-----------
  
  #Species 15
  #Abundance
  alpha.lambda15 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda15 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam15 <- exp(alpha.lambda15+beta.lambda15*siteCov)
  
  N15 <- rpois(n = n.site, lambda = lam15) # give it some noise
  
  table(N15)        # Distribution of abundances across sites
  tru.occ15<-sum(N15 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p15 <- rnorm(1,2,0.1)   		# Intercept
  beta.p15 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob15 <- exp(alpha.p15 + beta.p15 * siteCov) / (1 + exp(alpha.p15 + beta.p15 * siteCov))
  
  #Create observations
  y15 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y15[,j]<- rbinom(n=n.site, size=N15, prob=det.prob15)
  }
  #naive occupancy
  m<- apply(y15,1,max)
  m[m>1]<-1
  naive.occ15<-mean(m)
  
  #true occurance and each site
  occ15<- rep(0,n.site)
  occ15[N15>=1]<-1
  #----------------------------------End species 15-----------
  
  #Species 16
  #Abundance
  alpha.lambda16 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda16 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam16 <- exp(alpha.lambda16+beta.lambda16*siteCov)
  
  N16 <- rpois(n = n.site, lambda = lam16) # give it some noise
  
  table(N16)        # Distribution of abundances across sites
  tru.occ16<-sum(N16 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p16 <- rnorm(1,2,0.1)   		# Intercept
  beta.p16 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob16 <- exp(alpha.p16 + beta.p16 * siteCov) / (1 + exp(alpha.p16 + beta.p16 * siteCov))
  
  #Create observations
  y16 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y16[,j]<- rbinom(n=n.site, size=N16, prob=det.prob16)
  }
  #naive occupancy
  m<- apply(y16,1,max)
  m[m>1]<-1
  naive.occ16<-mean(m)
  
  #true occurance and each site
  occ16<- rep(0,n.site)
  occ16[N16>=1]<-1
  #----------------------------------End species 16-----------
  
  #Species 17
  #Abundance
  alpha.lambda17 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda17 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam17 <- exp(alpha.lambda17+beta.lambda17*siteCov)
  
  N17 <- rpois(n = n.site, lambda = lam17) # give it some noise
  
  table(N17)        # Distribution of abundances across sites
  tru.occ17<-sum(N17 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p17 <- rnorm(1,2,0.1)   		# Intercept
  beta.p17 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob17 <- exp(alpha.p17 + beta.p17 * siteCov) / (1 + exp(alpha.p17 + beta.p17 * siteCov))
  
  #Create observations
  y17 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y17[,j]<- rbinom(n=n.site, size=N17, prob=det.prob17)
  }
  #naive occupancy
  m<- apply(y17,1,max)
  m[m>1]<-1
  naive.occ17<-mean(m)
  
  #true occurance and each site
  occ17<- rep(0,n.site)
  occ17[N17>=1]<-1
  #----------------------------------End species 17-----------
  
  #Species 18
  #Abundance
  alpha.lambda18 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda18 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam18 <- exp(alpha.lambda18+beta.lambda18*siteCov)
  
  N18 <- rpois(n = n.site, lambda = lam18) # give it some noise
  
  table(N18)        # Distribution of abundances across sites
  tru.occ18<-sum(N18 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p18 <- rnorm(1,2,0.1)   		# Intercept
  beta.p18 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob18 <- exp(alpha.p18 + beta.p18 * siteCov) / (1 + exp(alpha.p18 + beta.p18 * siteCov))
  
  #Create observations
  y18 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y18[,j]<- rbinom(n=n.site, size=N18, prob=det.prob18)
  }
  #naive occupancy
  m<- apply(y18,1,max)
  m[m>1]<-1
  naive.occ18<-mean(m)
  
  #true occurance and each site
  occ18<- rep(0,n.site)
  occ18[N18>=1]<-1
  #----------------------------------End species 18-----------
  
  #Species 19
  #Abundance
  alpha.lambda19 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda19 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam19 <- exp(alpha.lambda19+beta.lambda19*siteCov)
  
  N19 <- rpois(n = n.site, lambda = lam19) # give it some noise
  
  table(N19)        # Distribution of abundances across sites
  tru.occ19<-sum(N19 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p19 <- rnorm(1,2,0.1)   		# Intercept
  beta.p19 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob19 <- exp(alpha.p19 + beta.p19 * siteCov) / (1 + exp(alpha.p19 + beta.p19 * siteCov))
  
  #Create observations
  y19 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y19[,j]<- rbinom(n=n.site, size=N19, prob=det.prob19)
  }
  #naive occupancy
  m<- apply(y19,1,max)
  m[m>1]<-1
  naive.occ19<-mean(m)
  
  #true occurance and each site
  occ19<- rep(0,n.site)
  occ19[N19>=1]<-1
  #----------------------------------End species 19-----------
  
  #Species 20
  #Abundance
  alpha.lambda20 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda20 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam20 <- exp(alpha.lambda20+beta.lambda20*siteCov)
  
  N20 <- rpois(n = n.site, lambda = lam20) # give it some noise
  
  table(N20)        # Distribution of abundances across sites
  tru.occ20<-sum(N20 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p20 <- rnorm(1,2,0.1)   		# Intercept
  beta.p20 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob20 <- exp(alpha.p20 + beta.p20 * siteCov) / (1 + exp(alpha.p20 + beta.p20 * siteCov))
  plot(siteCov, det.prob1, ylim = c(0,1), main = "", xlab = "", ylab = "Detection probability")
  
  #Create observations
  y20 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y20[,j]<- rbinom(n=n.site, size=N20, prob=det.prob20)
  }
  #naive occupancy
  m<- apply(y20,1,max)
  m[m>1]<-1
  naive.occ20<-mean(m)
  
  #true occurance and each site
  occ20<- rep(0,n.site)
  occ20[N20>=1]<-1
  #----------------------------------End species 20-----------
  
  #Species 21
  #Abundance
  alpha.lambda21 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda21 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam21 <- exp(alpha.lambda21+beta.lambda21*siteCov)
  
  N21 <- rpois(n = n.site, lambda = lam21) # give it some noise
  
  table(N21)        # Distribution of abundances across sites
  tru.occ21<-sum(N21 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p21 <- rnorm(1,2,0.1)     	# Intercept
  beta.p21 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob21 <- exp(alpha.p21 + beta.p21 * siteCov) / (1 + exp(alpha.p21 + beta.p21 * siteCov))
  
  #Create observations
  y21 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y21[,j]<- rbinom(n=n.site, size=N21, prob=det.prob21)
  }
  #naive occupancy
  m<- apply(y21,1,max)
  m[m>1]<-1
  naive.occ21<-mean(m)
  
  #true occurance and each site
  occ21<- rep(0,n.site)
  occ21[N21>=1]<-1
  #----------------------------------End species 21-----------
  
  #Species 22
  #Abundance
  alpha.lambda22 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda22 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam22 <- exp(alpha.lambda22+beta.lambda22*siteCov)
  
  N22 <- rpois(n = n.site, lambda = lam22) # give it some noise
  
  table(N22)        # Distribution of abundances across sites
  tru.occ22<-sum(N22 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p22 <- rnorm(1,2,0.1)     	# Intercept
  beta.p22 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob22 <- exp(alpha.p22 + beta.p22 * siteCov) / (1 + exp(alpha.p22 + beta.p22 * siteCov))
  
  #Create observations
  y22 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y22[,j]<- rbinom(n=n.site, size=N22, prob=det.prob22)
  }
  #naive occupancy
  m<- apply(y22,1,max)
  m[m>1]<-1
  naive.occ22<-mean(m)
  
  #true occurance and each site
  occ22<- rep(0,n.site)
  occ22[N22>=1]<-1
  #----------------------------------End species 22-----------
  
  #Species 23
  #Abundance
  alpha.lambda23 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda23 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam23 <- exp(alpha.lambda23+beta.lambda23*siteCov)
  
  N23 <- rpois(n = n.site, lambda = lam23) # give it some noise
  
  table(N23)        # Distribution of abundances across sites
  tru.occ23<-sum(N23 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p23 <- rnorm(1,2,0.1)     	# Intercept
  beta.p23 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob23 <- exp(alpha.p23 + beta.p23 * siteCov) / (1 + exp(alpha.p23 + beta.p23 * siteCov))
  
  #Create observations
  y23 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y23[,j]<- rbinom(n=n.site, size=N23, prob=det.prob23)
  }
  #naive occupancy
  m<- apply(y23,1,max)
  m[m>1]<-1
  naive.occ23<-mean(m)
  
  #true occurance and each site
  occ23<- rep(0,n.site)
  occ23[N23>=1]<-1
  #----------------------------------End species 23-----------
  
  #Species 24
  #Abundance
  alpha.lambda24 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda24 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam24 <- exp(alpha.lambda24+beta.lambda24*siteCov)
  
  N24 <- rpois(n = n.site, lambda = lam24) # give it some noise
  
  table(N24)        # Distribution of abundances across sites
  tru.occ24<-sum(N24 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p24 <- rnorm(1,2,0.1)     	# Intercept
  beta.p24 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob24 <- exp(alpha.p24 + beta.p24 * siteCov) / (1 + exp(alpha.p24 + beta.p24 * siteCov))
  
  #Create observations
  y24 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y24[,j]<- rbinom(n=n.site, size=N24, prob=det.prob24)
  }
  #naive occupancy
  m<- apply(y24,1,max)
  m[m>1]<-1
  naive.occ24<-mean(m)
  
  #true occurance and each site
  occ24<- rep(0,n.site)
  occ24[N24>=1]<-1
  #----------------------------------End species 24-----------
  
  #Species 25
  #Abundance
  alpha.lambda25 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda25 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam25 <- exp(alpha.lambda25+beta.lambda25*siteCov)
  
  N25 <- rpois(n = n.site, lambda = lam25) # give it some noise
  
  table(N25)        # Distribution of abundances across sites
  tru.occ25<-sum(N25 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p25 <- rnorm(1,2,0.1)     	# Intercept
  beta.p25 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob25 <- exp(alpha.p25 + beta.p25 * siteCov) / (1 + exp(alpha.p25 + beta.p25 * siteCov))
  
  #Create observations
  y25 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y25[,j]<- rbinom(n=n.site, size=N25, prob=det.prob25)
  }
  #naive occupancy
  m<- apply(y25,1,max)
  m[m>1]<-1
  naive.occ25<-mean(m)
  
  #true occurance and each site
  occ25<- rep(0,n.site)
  occ25[N25>=1]<-1
  #----------------------------------End species 25-----------
  
  #Species 26
  #Abundance
  alpha.lambda26 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda26 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam26 <- exp(alpha.lambda26+beta.lambda26*siteCov)
  
  N26 <- rpois(n = n.site, lambda = lam26) # give it some noise
  
  table(N26)        # Distribution of abundances across sites
  tru.occ26<-sum(N26 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p26 <- rnorm(1,2,0.1)     	# Intercept
  beta.p26 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob26 <- exp(alpha.p26 + beta.p26 * siteCov) / (1 + exp(alpha.p26 + beta.p26 * siteCov))
  
  #Create observations
  y26 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y26[,j]<- rbinom(n=n.site, size=N26, prob=det.prob26)
  }
  #naive occupancy
  m<- apply(y26,1,max)
  m[m>1]<-1
  naive.occ26<-mean(m)
  
  #true occurance and each site
  occ26<- rep(0,n.site)
  occ26[N26>=1]<-1
  #----------------------------------End species 26-----------
  
  #Species 27
  #Abundance
  alpha.lambda27 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda27 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam27 <- exp(alpha.lambda27+beta.lambda27*siteCov)
  
  N27 <- rpois(n = n.site, lambda = lam27) # give it some noise
  
  table(N27)        # Distribution of abundances across sites
  tru.occ27<-sum(N27 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p27 <- rnorm(1,2,0.1)     	# Intercept
  beta.p27 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob27 <- exp(alpha.p27 + beta.p27 * siteCov) / (1 + exp(alpha.p27 + beta.p27 * siteCov))
  
  #Create observations
  y27 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y27[,j]<- rbinom(n=n.site, size=N27, prob=det.prob27)
  }
  #naive occupancy
  m<- apply(y27,1,max)
  m[m>1]<-1
  naive.occ27<-mean(m)
  
  #true occurance and each site
  occ27<- rep(0,n.site)
  occ27[N27>=1]<-1
  #----------------------------------End species 27-----------
  
  #Species 28
  #Abundance
  alpha.lambda28 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda28 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam28 <- exp(alpha.lambda28+beta.lambda28*siteCov)
  
  N28 <- rpois(n = n.site, lambda = lam28) # give it some noise
  
  table(N28)        # Distribution of abundances across sites
  tru.occ28<-sum(N28 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p28 <- rnorm(1,2,0.1)     	# Intercept
  beta.p28 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob28 <- exp(alpha.p28 + beta.p28 * siteCov) / (1 + exp(alpha.p28 + beta.p28 * siteCov))
  
  #Create observations
  y28 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y28[,j]<- rbinom(n=n.site, size=N28, prob=det.prob28)
  }
  #naive occupancy
  m<- apply(y28,1,max)
  m[m>1]<-1
  naive.occ28<-mean(m)
  
  #true occurance and each site
  occ28<- rep(0,n.site)
  occ28[N28>=1]<-1
  #----------------------------------End species 28-----------
  
  #Species 29
  #Abundance
  alpha.lambda29 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda29 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam29 <- exp(alpha.lambda29+beta.lambda29*siteCov)
  
  N29 <- rpois(n = n.site, lambda = lam29) # give it some noise
  
  table(N29)        # Distribution of abundances across sites
  tru.occ29<-sum(N29 > 0) / n.site      # Empirical occupancy
  
  #Detection
  alpha.p29 <- rnorm(1,2,0.1)     	# Intercept
  beta.p29 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob29 <- exp(alpha.p29 + beta.p29 * siteCov) / (1 + exp(alpha.p29 + beta.p29 * siteCov))
  
  #Create observations
  y29 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y29[,j]<- rbinom(n=n.site, size=N29, prob=det.prob29)
  }
  #naive occupancy
  m<- apply(y29,1,max)
  m[m>1]<-1
  naive.occ29<-mean(m)
  
  #true occurance and each site
  occ29<- rep(0,n.site)
  occ29[N29>=1]<-1
  #----------------------------------End species 29-----------
  
  #Species 30
  #Abundance
  alpha.lambda30 <- rnorm(1,-2,0.3)  #Logit linear intercept for occurrence
  beta.lambda30 <- rnorm(1,4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam30 <- exp(alpha.lambda30+beta.lambda30*siteCov)
  
  N30 <- rpois(n = n.site, lambda = lam30) # give it some noise
  
  table(N30)        # Distribution of abundances across sites
  (tru.occ30<-sum(N30 > 0) / n.site	)		# Empirical occupancy
  
  #Detection
  alpha.p30 <- rnorm(1,2,0.1) 			# Intercept
  beta.p30 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob30 <- exp(alpha.p30 + beta.p30 * siteCov) / (1 + exp(alpha.p30 + beta.p29 * siteCov))
  
  #Create observations
  y30 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y30[,j]<- rbinom(n=n.site, size=N30, prob=det.prob30)
  }
  #naive occupancy
  m<- apply(y30,1,max)
  m[m>1]<-1
  naive.occ30<-mean(m)
  
  occ30 <- rep(0,n.site)
  occ30[N30>=1]<-1
  #----------------------species 30 end----
  #Species 31
  #Abundance
  alpha.lambda31 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda31 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam31 <- exp(alpha.lambda31+beta.lambda31*siteCov)
  
  N31 <- rpois(n = n.site, lambda = lam31) # give it some noise
  
  table(N31)        # Distribution of abundances across sites
  (tru.occ31<-sum(N31 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p31 <- rnorm(1,2,0.1) 			# Intercept
  beta.p31 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob31 <- exp(alpha.p31 + beta.p31 * siteCov) / (1 + exp(alpha.p31 + beta.p31 * siteCov))
  
  #Create observations
  y31 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y31[,j]<- rbinom(n=n.site, size=N31, prob=det.prob31)
  }
  #naive occupancy
  m<- apply(y31,1,max)
  m[m>1]<-1
  naive.occ31<-mean(m)
  
  occ31 <- rep(0,n.site)
  occ31[N31>=1]<-1
  #----------------------species 31 end------------------------------------------
  
  #Species 32
  #Abundance
  alpha.lambda32 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda32 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam32 <- exp(alpha.lambda32+beta.lambda32*siteCov)
  
  N32 <- rpois(n = n.site, lambda = lam32) # give it some noise
  
  table(N32)        # Distribution of abundances across sites
  (tru.occ32<-sum(N32 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p32 <- rnorm(1,2,0.1) 			# Intercept
  beta.p32 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob32 <- exp(alpha.p32 + beta.p32 * siteCov) / (1 + exp(alpha.p32 + beta.p32 * siteCov))
  
  #Create observations
  y32 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y32[,j]<- rbinom(n=n.site, size=N32, prob=det.prob32)
  }
  #naive occupancy
  m<- apply(y32,1,max)
  m[m>1]<-1
  naive.occ32<-mean(m)
  
  occ32 <- rep(0,n.site)
  occ32[N32>=1]<-1
  #----------------------species 32 end----
  
  #Species 33
  #Abundance
  alpha.lambda33 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda33 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam33 <- exp(alpha.lambda33+beta.lambda33*siteCov)
  
  N33 <- rpois(n = n.site, lambda = lam33) # give it some noise
  
  table(N33)        # Distribution of abundances across sites
  (tru.occ33<-sum(N33 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p33 <- rnorm(1,2,0.1) 			# Intercept
  beta.p33 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob33 <- exp(alpha.p33 + beta.p33 * siteCov) / (1 + exp(alpha.p33 + beta.p33 * siteCov))
  
  #Create observations
  y33 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y33[,j]<- rbinom(n=n.site, size=N33, prob=det.prob33)
  }
  #naive occupancy
  m<- apply(y33,1,max)
  m[m>1]<-1
  naive.occ33<-mean(m)
  
  occ33 <- rep(0,n.site)
  occ33[N33>=1]<-1
  #----------------------species 33 end----
  
  #Species 34
  #Abundance
  alpha.lambda34 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda34 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam34 <- exp(alpha.lambda34+beta.lambda34*siteCov)
  
  N34 <- rpois(n = n.site, lambda = lam34) # give it some noise
  
  table(N34)        # Distribution of abundances across sites
  (tru.occ34<-sum(N34 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p34 <- rnorm(1,2,0.1) 			# Intercept
  beta.p34 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob34 <- exp(alpha.p34 + beta.p34 * siteCov) / (1 + exp(alpha.p34 + beta.p34 * siteCov))
  
  #Create observations
  y34 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y34[,j]<- rbinom(n=n.site, size=N34, prob=det.prob34)
  }
  #naive occupancy
  m<- apply(y34,1,max)
  m[m>1]<-1
  naive.occ34<-mean(m)
  
  occ34 <- rep(0,n.site)
  occ34[N34>=1]<-1
  #----------------------species 34 end----
  
  #Species 35
  #Abundance
  alpha.lambda35 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda35 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam35 <- exp(alpha.lambda35+beta.lambda35*siteCov)
  
  N35 <- rpois(n = n.site, lambda = lam35) # give it some noise
  
  table(N35)        # Distribution of abundances across sites
  (tru.occ35<-sum(N35 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p35 <- rnorm(1,2,0.1) 			# Intercept
  beta.p35 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob35 <- exp(alpha.p35 + beta.p35 * siteCov) / (1 + exp(alpha.p35 + beta.p35 * siteCov))
  
  #Create observations
  y35 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y35[,j]<- rbinom(n=n.site, size=N35, prob=det.prob35)
  }
  #naive occupancy
  m<- apply(y35,1,max)
  m[m>1]<-1
  naive.occ35<-mean(m)
  
  occ35 <- rep(0,n.site)
  occ35[N35>=1]<-1
  #----------------------species 35 end----
  
  #Species 36
  #Abundance
  alpha.lambda36 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda36 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam36 <- exp(alpha.lambda36+beta.lambda36*siteCov)
  
  N36 <- rpois(n = n.site, lambda = lam36) # give it some noise
  
  table(N36)        # Distribution of abundances across sites
  (tru.occ36<-sum(N36 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p36 <- rnorm(1,2,0.1) 			# Intercept
  beta.p36 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob36 <- exp(alpha.p36 + beta.p36 * siteCov) / (1 + exp(alpha.p36 + beta.p36 * siteCov))
  
  #Create observations
  y36 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y36[,j]<- rbinom(n=n.site, size=N36, prob=det.prob36)
  }
  #naive occupancy
  m<- apply(y36,1,max)
  m[m>1]<-1
  naive.occ36<-mean(m)
  
  occ36 <- rep(0,n.site)
  occ36[N36>=1]<-1
  #----------------------species 36 end----
  
  #Species 37
  #Abundance
  alpha.lambda37 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda37 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam37 <- exp(alpha.lambda37+beta.lambda37*siteCov)
  
  N37 <- rpois(n = n.site, lambda = lam37) # give it some noise
  
  table(N37)        # Distribution of abundances across sites
  (tru.occ37<-sum(N37 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p37 <- rnorm(1,2,0.1) 			# Intercept
  beta.p37 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob37 <- exp(alpha.p37 + beta.p37 * siteCov) / (1 + exp(alpha.p37 + beta.p37 * siteCov))
  
  #Create observations
  y37 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y37[,j]<- rbinom(n=n.site, size=N37, prob=det.prob37)
  }
  #naive occupancy
  m<- apply(y37,1,max)
  m[m>1]<-1
  naive.occ37<-mean(m)
  
  occ37 <- rep(0,n.site)
  occ37[N37>=1]<-1
  #----------------------species 37 end----
  
  #Species 38
  #Abundance
  alpha.lambda38 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda38 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam38 <- exp(alpha.lambda38+beta.lambda38*siteCov)
  
  N38 <- rpois(n = n.site, lambda = lam38) # give it some noise
  
  table(N38)        # Distribution of abundances across sites
  (tru.occ38<-sum(N38 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p38 <- rnorm(1,2,0.1) 			# Intercept
  beta.p38 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob38 <- exp(alpha.p38 + beta.p38 * siteCov) / (1 + exp(alpha.p38 + beta.p38 * siteCov))
  
  #Create observations
  y38 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y38[,j]<- rbinom(n=n.site, size=N38, prob=det.prob38)
  }
  #naive occupancy
  m<- apply(y38,1,max)
  m[m>1]<-1
  naive.occ38<-mean(m)
  
  occ38 <- rep(0,n.site)
  occ38[N38>=1]<-1
  #----------------------species 38 end----
  
  #Species 39
  #Abundance
  alpha.lambda39 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda39 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam39 <- exp(alpha.lambda39+beta.lambda39*siteCov)
  
  N39 <- rpois(n = n.site, lambda = lam39) # give it some noise
  
  table(N39)        # Distribution of abundances across sites
  (tru.occ39<-sum(N39 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p39 <- rnorm(1,2,0.1) 			# Intercept
  beta.p39 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob39 <- exp(alpha.p39 + beta.p39 * siteCov) / (1 + exp(alpha.p39 + beta.p39 * siteCov))
  
  #Create observations
  y39 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y39[,j]<- rbinom(n=n.site, size=N39, prob=det.prob39)
  }
  #naive occupancy
  m<- apply(y39,1,max)
  m[m>1]<-1
  naive.occ39<-mean(m)
  
  occ39 <- rep(0,n.site)
  occ39[N39>=1]<-1
  #----------------------species 39 end----
  
  #Species 40
  #Abundance
  alpha.lambda40 <- rnorm(1,2,0.3)  #Logit linear intercept for occurrence
  beta.lambda40 <- rnorm(1,-4,0.4)  #Logit linear slope for siteCov on occurrence
  
  lam40 <- exp(alpha.lambda40+beta.lambda40*siteCov)
  
  N40 <- rpois(n = n.site, lambda = lam40) # give it some noise
  
  table(N40)        # Distribution of abundances across sites
  (tru.occ40<-sum(N40 > 0) / n.site  )		# Empirical occupancy
  
  #Detection
  alpha.p40 <- rnorm(1,2,0.1) 			# Intercept
  beta.p40 <- rnorm(1,-4,0.3)				# Linear effect of vegetation
  det.prob40 <- exp(alpha.p40 + beta.p40 * siteCov) / (1 + exp(alpha.p40 + beta.p40 * siteCov))
  
  #Create observations
  y40 <- array(dim=c(n.site,T))
  for(j in 1:T){
    y40[,j]<- rbinom(n=n.site, size=N40, prob=det.prob40)
  }
  #naive occupancy
  m<- apply(y40,1,max)
  m[m>1]<-1
  naive.occ40<-mean(m)
  
  occ40 <- rep(0,n.site)
  occ40[N40>=1]<-1
  #----------------------species 40 end----
  
  #Data arrays of observations to feed into Bugs 
  #Put all species in a multi-dimensional array
  #array grouped by species - counts of birds
  data.array<-array(dim=c(n.site,T,40))  #This is the format to use.
  data.array[,,1]<-y1
  data.array[,,2]<-y2
  data.array[,,3]<-y3
  data.array[,,4]<-y4
  data.array[,,5]<-y5
  data.array[,,6]<-y6
  data.array[,,7]<-y7
  data.array[,,8]<-y8
  data.array[,,9]<-y9
  data.array[,,10]<-y10
  data.array[,,11]<-y11
  data.array[,,12]<-y12
  data.array[,,13]<-y13
  data.array[,,14]<-y14
  data.array[,,15]<-y15
  data.array[,,16]<-y16
  data.array[,,17]<-y17
  data.array[,,18]<-y18
  data.array[,,19]<-y19
  data.array[,,20]<-y20
  data.array[,,21]<-y21
  data.array[,,22]<-y22
  data.array[,,23]<-y23
  data.array[,,24]<-y24
  data.array[,,25]<-y25
  data.array[,,26]<-y26
  data.array[,,27]<-y27
  data.array[,,28]<-y28
  data.array[,,29]<-y29
  data.array[,,30]<-y30
  data.array[,,31]<-y31
  data.array[,,32]<-y32
  data.array[,,33]<-y33
  data.array[,,34]<-y34
  data.array[,,35]<-y35
  data.array[,,36]<-y36
  data.array[,,37]<-y37
  data.array[,,38]<-y38
  data.array[,,39]<-y39
  data.array[,,40]<-y40
  
  #Degrade Counts to Detected/Not detected for input into MS-OCC model
  data.array[data.array>=1]<-1
  
  #Here are the true values of all parameters
  
  #True species richness
  site.spp.occ <- cbind(occ1,occ2,occ3,occ4,occ5,occ6,occ7,occ8,occ9,occ10,
                        occ11,occ12,occ13,occ14,occ15,occ16,occ17,occ18,occ19,occ20,
                        occ21,occ22,occ23,occ24,occ25,occ26,occ27,occ28,occ29,occ30,
                        occ31,occ32,occ33,occ34,occ35,occ36,occ37,occ38,occ39,occ40)
  tru.rich <- rowSums(site.spp.occ)  #true species richness at each site
  
  
  ##Naive richness - based on whether each species was observed during the first survey (which was created randomly)
  #first visit observation
  occ.1 <- rep(0,n.site)
  occ.1[y1[,1]>=1]<-1
  occ.2 <- rep(0,n.site)
  occ.2[y2[,1]>=1]<-1
  occ.3 <- rep(0,n.site)
  occ.3[y3[,1]>=1]<-1
  occ.4 <- rep(0,n.site)
  occ.4[y4[,1]>=1]<-1
  occ.5 <- rep(0,n.site)
  occ.5[y5[,1]>=1]<-1
  occ.6 <- rep(0,n.site)
  occ.6[y6[,1]>=1]<-1
  occ.7 <- rep(0,n.site)
  occ.7[y7[,1]>=1]<-1
  occ.8 <- rep(0,n.site)
  occ.8[y8[,1]>=1]<-1
  occ.9 <- rep(0,n.site)
  occ.9[y9[,1]>=1]<-1
  occ.10 <- rep(0,n.site)
  occ.10[y10[,1]>=1]<-1
  occ.11 <- rep(0,n.site)
  occ.11[y11[,1]>=1]<-1
  occ.12 <- rep(0,n.site)
  occ.12[y12[,1]>=1]<-1
  occ.13 <- rep(0,n.site)
  occ.13[y13[,1]>=1]<-1
  occ.14 <- rep(0,n.site)
  occ.14[y14[,1]>=1]<-1
  occ.15 <- rep(0,n.site)
  occ.15[y15[,1]>=1]<-1
  occ.16 <- rep(0,n.site)
  occ.16[y16[,1]>=1]<-1
  occ.17 <- rep(0,n.site)
  occ.17[y17[,1]>=1]<-1
  occ.18 <- rep(0,n.site)
  occ.18[y18[,1]>=1]<-1
  occ.19 <- rep(0,n.site)
  occ.19[y19[,1]>=1]<-1
  occ.20 <- rep(0,n.site)
  occ.20[y20[,1]>=1]<-1
  occ.21 <- rep(0,n.site)
  occ.21[y21[,1]>=1]<-1
  occ.22 <- rep(0,n.site)
  occ.22[y22[,1]>=1]<-1
  occ.23 <- rep(0,n.site)
  occ.23[y23[,1]>=1]<-1
  occ.24 <- rep(0,n.site)
  occ.24[y24[,1]>=1]<-1
  occ.25 <- rep(0,n.site)
  occ.25[y25[,1]>=1]<-1
  occ.26 <- rep(0,n.site)
  occ.26[y26[,1]>=1]<-1
  occ.27 <- rep(0,n.site)
  occ.27[y27[,1]>=1]<-1
  occ.28 <- rep(0,n.site)
  occ.28[y28[,1]>=1]<-1
  occ.29 <- rep(0,n.site)
  occ.29[y29[,1]>=1]<-1
  occ.30 <- rep(0,n.site)
  occ.30[y30[,1]>=1]<-1
  occ.31 <- rep(0,n.site)
  occ.31[y31[,1]>=1]<-1
  occ.32 <- rep(0,n.site)
  occ.32[y32[,1]>=1]<-1
  occ.33 <- rep(0,n.site)
  occ.33[y33[,1]>=1]<-1
  occ.34 <- rep(0,n.site)
  occ.34[y34[,1]>=1]<-1
  occ.35 <- rep(0,n.site)
  occ.35[y35[,1]>=1]<-1
  occ.36 <- rep(0,n.site)
  occ.36[y36[,1]>=1]<-1
  occ.37 <- rep(0,n.site)
  occ.37[y37[,1]>=1]<-1
  occ.38 <- rep(0,n.site)
  occ.38[y38[,1]>=1]<-1
  occ.39 <- rep(0,n.site)
  occ.39[y39[,1]>=1]<-1
  occ.40 <- rep(0,n.site)
  occ.40[y40[,1]>=1]<-1
  
  #Observed species richness
  site.spp.occ.obs <- cbind(occ.1,occ.2,occ.3,occ.4,occ.5,occ.6,occ.7,occ.8,occ.9,occ.10,
                            occ.11,occ.12,occ.13,occ.14,occ.15,occ.16,occ.17,occ.18,occ.19,occ.20,
                            occ.21,occ.22,occ.23,occ.24,occ.25,occ.26,occ.27,occ.28,occ.29,occ.30,
                            occ.31,occ.32,occ.33,occ.34,occ.35,occ.36,occ.37,occ.38,occ.39,occ.40)
  
  naive.rich <- rowSums(site.spp.occ.obs)  #true species richness at each site
  tot.naive.rich <- colSums(site.spp.occ.obs)  #next 3 lines calculate observed total species richness at all sites
  tot.naive.rich[tot.naive.rich>=1]<-1
  sum(tot.naive.rich)
  
  #True intercepts and slopes of species occupancy vs. siteCov
  #I had to generate data as counts for the jacknife estiamtors, not as presence/absence,
  #so I fit a logistic regression on the true occurance data to yield unbiased estimate of the parameters of occupancy vs. siteCov
  
  #Assemble table of true alphas and betas for each species
  occ.table <- array(dim=c(n.species,2))
  colnames(occ.table)<-c("alpha", "beta")
  
  occ.table[1,1] <- as.numeric(coef(glm(occ1~siteCov, family=binomial))[1])  #stick the alpha for species1 in the table
  occ.table[1,2] <- as.numeric(coef(glm(occ1~siteCov, family=binomial))[2])
  occ.table[2,1] <- as.numeric(coef(glm(occ2~siteCov, family=binomial))[1])
  occ.table[2,2] <- as.numeric(coef(glm(occ2~siteCov, family=binomial))[2])
  occ.table[3,1] <- as.numeric(coef(glm(occ3~siteCov, family=binomial))[1])
  occ.table[3,2] <- as.numeric(coef(glm(occ3~siteCov, family=binomial))[2])
  occ.table[4,1] <- as.numeric(coef(glm(occ4~siteCov, family=binomial))[1])
  occ.table[4,2] <- as.numeric(coef(glm(occ4~siteCov, family=binomial))[2])
  occ.table[5,1] <- as.numeric(coef(glm(occ5~siteCov, family=binomial))[1])
  occ.table[5,2] <- as.numeric(coef(glm(occ5~siteCov, family=binomial))[2])
  occ.table[6,1] <- as.numeric(coef(glm(occ6~siteCov, family=binomial))[1])
  occ.table[6,2] <- as.numeric(coef(glm(occ6~siteCov, family=binomial))[2])
  occ.table[7,1] <- as.numeric(coef(glm(occ7~siteCov, family=binomial))[1])
  occ.table[7,2] <- as.numeric(coef(glm(occ7~siteCov, family=binomial))[2])
  occ.table[8,1] <- as.numeric(coef(glm(occ8~siteCov, family=binomial))[1])
  occ.table[8,2] <- as.numeric(coef(glm(occ8~siteCov, family=binomial))[2])
  occ.table[9,1] <- as.numeric(coef(glm(occ9~siteCov, family=binomial))[1])
  occ.table[9,2] <- as.numeric(coef(glm(occ9~siteCov, family=binomial))[2])
  occ.table[10,1] <- as.numeric(coef(glm(occ10~siteCov, family=binomial))[1])
  occ.table[10,2] <- as.numeric(coef(glm(occ10~siteCov, family=binomial))[2])
  occ.table[11,1] <- as.numeric(coef(glm(occ11~siteCov, family=binomial))[1])
  occ.table[11,2] <- as.numeric(coef(glm(occ11~siteCov, family=binomial))[2])
  occ.table[12,1] <- as.numeric(coef(glm(occ12~siteCov, family=binomial))[1])
  occ.table[12,2] <- as.numeric(coef(glm(occ12~siteCov, family=binomial))[2])
  occ.table[13,1] <- as.numeric(coef(glm(occ13~siteCov, family=binomial))[1])
  occ.table[13,2] <- as.numeric(coef(glm(occ13~siteCov, family=binomial))[2])
  occ.table[14,1] <- as.numeric(coef(glm(occ14~siteCov, family=binomial))[1])
  occ.table[14,2] <- as.numeric(coef(glm(occ14~siteCov, family=binomial))[2])
  occ.table[15,1] <- as.numeric(coef(glm(occ15~siteCov, family=binomial))[1])
  occ.table[15,2] <- as.numeric(coef(glm(occ15~siteCov, family=binomial))[2])
  occ.table[16,1] <- as.numeric(coef(glm(occ16~siteCov, family=binomial))[1])
  occ.table[16,2] <- as.numeric(coef(glm(occ16~siteCov, family=binomial))[2])
  occ.table[17,1] <- as.numeric(coef(glm(occ17~siteCov, family=binomial))[1])
  occ.table[17,2] <- as.numeric(coef(glm(occ17~siteCov, family=binomial))[2])
  occ.table[18,1] <- as.numeric(coef(glm(occ18~siteCov, family=binomial))[1])
  occ.table[18,2] <- as.numeric(coef(glm(occ18~siteCov, family=binomial))[2])
  occ.table[19,1] <- as.numeric(coef(glm(occ19~siteCov, family=binomial))[1])
  occ.table[19,2] <- as.numeric(coef(glm(occ19~siteCov, family=binomial))[2])
  occ.table[20,1] <- as.numeric(coef(glm(occ20~siteCov, family=binomial))[1])
  occ.table[20,2] <- as.numeric(coef(glm(occ20~siteCov, family=binomial))[2])
  occ.table[21,1] <- as.numeric(coef(glm(occ21~siteCov, family=binomial))[1])
  occ.table[21,2] <- as.numeric(coef(glm(occ21~siteCov, family=binomial))[2])
  occ.table[22,1] <- as.numeric(coef(glm(occ22~siteCov, family=binomial))[1])
  occ.table[22,2] <- as.numeric(coef(glm(occ22~siteCov, family=binomial))[2])
  occ.table[23,1] <- as.numeric(coef(glm(occ23~siteCov, family=binomial))[1])
  occ.table[23,2] <- as.numeric(coef(glm(occ23~siteCov, family=binomial))[2])
  occ.table[24,1] <- as.numeric(coef(glm(occ24~siteCov, family=binomial))[1])
  occ.table[24,2] <- as.numeric(coef(glm(occ24~siteCov, family=binomial))[2])
  occ.table[25,1] <- as.numeric(coef(glm(occ25~siteCov, family=binomial))[1])
  occ.table[25,2] <- as.numeric(coef(glm(occ25~siteCov, family=binomial))[2])
  occ.table[26,1] <- as.numeric(coef(glm(occ26~siteCov, family=binomial))[1])
  occ.table[26,2] <- as.numeric(coef(glm(occ26~siteCov, family=binomial))[2])
  occ.table[27,1] <- as.numeric(coef(glm(occ27~siteCov, family=binomial))[1])
  occ.table[27,2] <- as.numeric(coef(glm(occ27~siteCov, family=binomial))[2])
  occ.table[28,1] <- as.numeric(coef(glm(occ28~siteCov, family=binomial))[1])
  occ.table[28,2] <- as.numeric(coef(glm(occ28~siteCov, family=binomial))[2])
  occ.table[29,1] <- as.numeric(coef(glm(occ29~siteCov, family=binomial))[1])
  occ.table[29,2] <- as.numeric(coef(glm(occ29~siteCov, family=binomial))[2])
  occ.table[30,1] <- as.numeric(coef(glm(occ30~siteCov, family=binomial))[1])
  occ.table[30,2] <- as.numeric(coef(glm(occ30~siteCov, family=binomial))[2])
  occ.table[31,1] <- as.numeric(coef(glm(occ31~siteCov, family=binomial))[1])
  occ.table[31,2] <- as.numeric(coef(glm(occ31~siteCov, family=binomial))[2])
  occ.table[32,1] <- as.numeric(coef(glm(occ32~siteCov, family=binomial))[1])
  occ.table[32,2] <- as.numeric(coef(glm(occ32~siteCov, family=binomial))[2])
  occ.table[33,1] <- as.numeric(coef(glm(occ33~siteCov, family=binomial))[1])
  occ.table[33,2] <- as.numeric(coef(glm(occ33~siteCov, family=binomial))[2])
  occ.table[34,1] <- as.numeric(coef(glm(occ34~siteCov, family=binomial))[1])
  occ.table[34,2] <- as.numeric(coef(glm(occ34~siteCov, family=binomial))[2])
  occ.table[35,1] <- as.numeric(coef(glm(occ35~siteCov, family=binomial))[1])
  occ.table[35,2] <- as.numeric(coef(glm(occ35~siteCov, family=binomial))[2])
  occ.table[36,1] <- as.numeric(coef(glm(occ36~siteCov, family=binomial))[1])
  occ.table[36,2] <- as.numeric(coef(glm(occ36~siteCov, family=binomial))[2])
  occ.table[37,1] <- as.numeric(coef(glm(occ37~siteCov, family=binomial))[1])
  occ.table[37,2] <- as.numeric(coef(glm(occ37~siteCov, family=binomial))[2])
  occ.table[38,1] <- as.numeric(coef(glm(occ38~siteCov, family=binomial))[1])
  occ.table[38,2] <- as.numeric(coef(glm(occ38~siteCov, family=binomial))[2])
  occ.table[39,1] <- as.numeric(coef(glm(occ39~siteCov, family=binomial))[1])
  occ.table[39,2] <- as.numeric(coef(glm(occ39~siteCov, family=binomial))[2])
  occ.table[40,1] <- as.numeric(coef(glm(occ40~siteCov, family=binomial))[1])
  occ.table[40,2] <- as.numeric(coef(glm(occ40~siteCov, family=binomial))[2])
  
  
  #Combine all the Counts from the first surveys into a site x species data table
  obs.counts <- cbind(y1[,1], y2[,1], y3[,1], y4[,1],  y5[,1], y6[,1], y7[,1], y8[,1], y9[,1], y10[,1],
                      y11[,1], y12[,1], y13[,1], y14[,1], y15[,1], y16[,1], y17[,1], y18[,1], y19[,1], y20[,1],
                      y21[,1], y22[,1], y23[,1], y24[,1], y25[,1], y26[,1], y27[,1], y28[,1], y29[,1], y30[,1],
                      y31[,1], y32[,1], y33[,1], y34[,1], y35[,1], y36[,1], y37[,1], y38[,1], y39[,1], y40[,1])
  
  #Build a function to transform data to right format for SPECIES and get jackknife estimate for each site
  jack.func <- function(i){
    j <- as.data.frame(table(obs.counts[i,]))  #Manipulate data into the right format for package SPECIES
    j$Var1<-strtoi(j[,1])                      #converts a character string to integer to pass to jackknife
    ifelse(j[2,1]>1,NA, ifelse(nrow(j)<3,NA,jackknife(j[-1,],k=2)$Nhat))
    # 2nd order jackknife estimator (delete the first row of 0 observations of species), also the jackknife needs the frequecies in 1:n ordered array. cases where there are no singletons make it croak
  }       #A jackknife estimator for some sites can't be estimated due to low number of observations,
  #for these cases, assign "NA" to the data frame output
  
  #Estimate jackknife species richness at each site and stick in a data array
  jack.rich <- array(dim=c(50,1))
  for (i in 1:50){
    jack.rich[i,] = jack.func(i)
  }
  
  #Estimate asymptotic species richness according to Chao and stick in an array.
  chao.func <- function(i){
    chao <- as.data.frame(table(obs.counts[i,]))  #Manipulate data into the right format for package SPECIES
    chao$Var1<-strtoi(chao[,1])                      #converts a character string to integer to pass to chao1984
    ifelse(nrow(chao)<3,NA,chao1984(chao[-1,])$Nhat)
    #chao estimator (delete the first row of 0 observations of species)
  }
  
  chao.rich <- array(dim=c(50,1))
  for (i in 1:50){
    chao.rich[i,] = chao.func(i)
  }
  chao.rich<-ifelse(chao.rich>200,NA,chao.rich)  #Replace all Inf solutions with NA
  
  ###############################ANALYSIS IN JAGS######################################
  #Analysis with Doyle-Royle Multi-species occupancy model 
  nrep=3
  nspec=40
  
  Y <- data.array
  
  # Bundle data
  win.data <- list(Y = Y, nsite = dim(Y)[1], nrep = dim(Y)[2], nspec = dim(Y)[3], 
                   siteCov = siteCov)
  
  #########Multispecies Occupancy Model for Species Richness Comparison######################
  
  ####No data augmentation - Dorazio et al. 2005
  
  # Specify model in BUGS language
  sink("model11.txt")
  cat("
      model {
      # Priors
      for(k in 1:nspec){ # For each species
      alpha0[k] ~ dnorm(mu.alpha0, tau.alpha0)
      alpha1[k] ~ dnorm(mu.alpha1, tau.alpha1)
      beta0[k] ~ dnorm(mu.beta0, tau.beta0)
      beta1[k] ~ dnorm(mu.beta1, tau.beta1)
      }
      
      #Hyperpriors
      mu.alpha0 ~ dunif(-5,5)
      tau.alpha0 <- 1 / (sd.alpha0 * sd.alpha0)
      sd.alpha0 ~ dunif(0, 10)
      mu.alpha1 ~ dunif(-5, 5)
      tau.alpha1 <- 1 / (sd.alpha1 * sd.alpha1)
      sd.alpha1 ~ dunif(0, 5)
      mu.beta0 ~ dunif(-5,5)
      tau.beta0 <- 1 / (sd.beta0 * sd.beta0)
      sd.beta0 ~ dunif(0, 5)
      mu.beta1 ~ dunif(-5, 5)
      tau.beta1 <- 1 / (sd.beta1 * sd.beta1)
      sd.beta1 ~ dunif(0, 5)
      
      
      # Likelihood
      for (k in 1:nspec) {                    # Loop over sites
      for(i in 1:nsite){                       # Loop over species
      logit(psi[i,k]) <- alpha0[k] + alpha1[k] * siteCov[i]
      z[i,k] ~ dbern(psi[i,k])                 # True occupancy z at site i for species k   
      }
      }
      
      # Observation model for the actual observations
      for (k in 1:nspec) {                    # Loop over sites
      for(i in 1:nsite){                       # Loop over species
      for (j in 1:nrep){                     # Loop over surveys  
      logit(p[i,j,k]) <- beta0[k] + beta1[k]*siteCov[i]
      mu.p[i,j,k] <- z[i,k] * p[i,j,k]    
      Y[i,j,k] ~ dbern(mu.p[i,j,k])
      }
      }
      }
      
      
      # Derived quantities
      for(k in 1:(nspec)){
      occ.fs[k] <- sum(z[,k])        # Number of occupied sites among the 50
      }
      
      for (i in 1:nsite) {
      Nsite[i] <- sum(z[i,])         # Number of occurring species at each site
      }
      }
      ",fill = TRUE)
  sink()
  
  # Initial values
  zst <- apply(Y, c(1,3), max) # Observed occurrence as starting values for z
  zst[is.na(zst)] <- 1
  #wst <- c(rep(1, nspec), rep(0, nz))
  
  inits <- function() list(z = zst)  #, lpsi = rnorm(n = nspec), 
  #betalpsi1 = rnorm(n = nspec), lp = rnorm(n = nspec), betalp1 = rnorm(n = nspec))
  
  # Parameters monitored
  params <-
    c("alpha0", "sd.alpha0", "alpha1", "sd.alpha1" ,
      "beta0", "sd.beta0","beta1", "sd.beta1", "occ.fs", "Nsite")
  
  # MCMC settings
  ni <- 10000
  nt <- 4
  nb <- 2000
  nc <- 3
  
  
  # Call JAGS from R
  
  out <- jags(win.data, inits, params, "model11.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb) 
  
  
  #################################################################################
  #Evaualting bias in estiating local richness  - put this stuff into an overall loop to build distributions of   #bias
  Nsite <- out$BUGSoutput$mean$Nsite
  
  #Bias of MSOCC
  bias.msocc <- out$BUGSoutput$mean$Nsite - tru.rich
  m.bias.msocc[k] <- mean(bias.msocc)
  sd.bias.msocc[k] <- sd(bias.msocc)
  baye.p.msocc[k] <-mean(tru.rich > out$BUGSoutput$mean$Nsite)  #Baysian P-value (proportion of simulations where the true richness
  #was greater than the estimated richness- values close to 0 or 1
  #indicate signficant bias)
  
  #Bias of the Jackknife estimates
  bias.jack <- jack.rich - tru.rich
  m.bias.jack[k] <- mean(bias.jack, na.rm=TRUE)
  sd.bias.jack[k] <- sd(bias.jack, na.rm=TRUE)   #spread of the bias  
  baye.p.jack[k] <-mean(tru.rich > jack.rich, na.rm=T) 
  
  #Bias of the Chao estimates
  bias.chao <- chao.rich - tru.rich
  m.bias.chao[k] <- mean(bias.chao, na.rm=TRUE)
  sd.bias.chao[k] <- sd(bias.chao, na.rm=TRUE)
  baye.p.chao[k] <-mean(tru.rich > chao.rich, na.rm=T) 
  
  #Bias of the naive estiates
  bias.naive <- naive.rich - tru.rich
  m.bias.naive[k] <- mean(bias.naive, na.rm=TRUE)
  sd.bias.naive[k] <- sd(bias.naive, na.rm=TRUE)
  baye.p.naive[k] <-mean(tru.rich > naive.rich) 
  
  
  ###Bias in the slopes of effects  -  Beta of estimator-siteCov minus Betat of true relationship
  
  rich.cov = glm(tru.rich~siteCov)
  coeff.t = as.vector(coef(rich.cov))
  m.alpha.t[k] = coeff.t[1]                   # "m." = means , "sd." standard deviation
  sd.alpha.t[k]= summary(rich.cov)$coefficients[1,2]
  m.beta.t[k] = coeff.t[2]
  sd.beta.t[k]= summary(rich.cov)$coefficients[2,2]
  
  nrich.cov = glm(naive.rich~siteCov)
  coeff.n = as.vector(coef(nrich.cov))
  m.alpha.n[k] = coeff.n[1]
  sd.alpha.n[k]= summary(nrich.cov)$coefficients[1,2]
  m.beta.n[k] = coeff.n[2]
  sd.beta.n[k]= summary(nrich.cov)$coefficients[2,2]
  
  jrich.cov = glm(jack.rich~siteCov)
  coeff.j = as.vector(coef(jrich.cov))
  m.alpha.j[k] = coeff.j[1]
  sd.alpha.j[k]= summary(jrich.cov)$coefficients[1,2]
  m.beta.j[k] = coeff.j[2]
  sd.beta.j[k]= summary(jrich.cov)$coefficients[2,2]
  
  crich.cov = glm(chao.rich~siteCov)
  coeff.c = as.vector(coef(crich.cov))
  m.alpha.c[k] = coeff.c[1]
  sd.alpha.c[k]= summary(crich.cov)$coefficients[1,2]
  m.beta.c[k] = coeff.c[2]
  sd.beta.c[k]= summary(crich.cov)$coefficients[2,2]
  
  msoccrich.cov = glm(Nsite~siteCov)
  coeff.msocc = as.vector(coef(msoccrich.cov))
  m.alpha.msocc[k] = coeff.msocc[1]
  sd.alpha.msocc[k]= summary(msoccrich.cov)$coefficients[1,2]
  m.beta.msocc[k] = coeff.msocc[2]
  sd.beta.msocc[k]= summary(msoccrich.cov)$coefficients[2,2]
  
  
  #For each site determine whether the true richness fell outside the 95% crI of the estimated local richness
  q.Nsite <- array(dim=c(50,2))
  for(i in 1:length(out$BUGSoutput$mean$Nsite)){
    q.Nsite[i,] <- quantile(out$BUGSoutput$sims.list$Nsite[,i],probs=c(0.025,0.975))
  }
  c_Nsite<-ifelse(tru.rich>=q.Nsite[,1]&tru.rich<=q.Nsite[,2],0,1)
  
  #Probability of estimate being biased
  p.bias.msocc[k] <- mean(c_Nsite)  # Probability of bias = 0.02
  
  
  
}  ) #This will be the end of the simulations

########################################

print(out,dig=3)

#####################################################################################
#Some example plots to visualize output

#Plot true species richness, naive species richness, and jackknife richness against site Covariate

par(mfrow=c(1,1), cex.lab=2, las=1, cex.axis=2, mar=c(4,4.5,1,1))
plot(siteCov, tru.rich, ylim = c(0,60), main = "", xlab = "Site Covariate", ylab = "Species Richness",pch=16, col="red")
rich.cov=glm(tru.rich~siteCov)
lines(siteCov, predict(rich.cov), ylim=c(0,100), main="", las=1, lwd=2, col="red")

points(siteCov, naive.rich, ylim = c(0,40), main = "", pch=16)
nrich.cov=glm(naive.rich~siteCov)
lines(siteCov, predict(nrich.cov), ylim=c(0,40), main="", las=1, lwd=2)

points(siteCov, jack.rich, ylim = c(0,40), main = "", pch=16, col="blue")
jrich.cov = glm(jack.rich~siteCov)
coeff= as.vector(coef(jrich.cov[1]))
j.pred <- coeff[1] + coeff[2]*siteCov
lines(siteCov, j.pred, ylim=c(0,40), main="", las=1, lwd=2, col="blue")

points(siteCov, chao.rich, ylim = c(0,40), main = "", pch=16, col="orange")
crich.cov = glm(chao.rich~siteCov)    #there are some Inf values for chao.rich so an error message here
coeff= as.vector(coef(crich.cov[1]))
c.pred <- coeff[1] + coeff[2]*siteCov
lines(siteCov, c.pred, ylim=c(0,40), main="", las=1, lwd=2, col="orange")

points(siteCov, Nsite, col="green", pch=16)                              #MS OCC predictions
lines(siteCov, predict(lm(Nsite~siteCov)), col="green", lwd=2)

legend("topleft", c("True", "Naive", "Jackknife", "Chao", "MSOCC"), lty=c(1,1,1,1), lwd=2, 
       col=c("red", "black", "blue", "orange", "green"), cex=1.3)

#----------------------------Evaluate Bias of Estimators----------------------------------------------------------

par(mfrow = c(4,2), mai=c(0.2,0.2,0.2,0.2), mar=c(4,4.5,1,1), oma=c(1,1,0,0), las=1,cex.lab=2, cex.axis=2)

hist(unlist(m.bias.naive), xlim=c(-20,20), main="", xlab="", ylab="Naive")
abline(v=0, col="red", lwd=3)
hist(unlist(sd.bias.naive), xlim=c(0,30), main="", xlab="",ylab="")

hist(unlist(m.bias.jack), xlim=c(-20,20), main="", xlab="",ylab="Jackknife")
abline(v=0, col="red", lwd=3)
hist(unlist(sd.bias.jack), xlim=c(0,30), main="", xlab="",ylab="")

hist(unlist(m.bias.chao), xlim=c(-20,20), main="", xlab="",ylab="Chao")
abline(v=0, col="red", lwd=3)
hist(unlist(sd.bias.chao), xlim=c(0,30), main="", xlab="",ylab="")

hist(unlist(m.bias.msocc), xlim=c(-20,20), main="", xlab="No. of species", ylab="MSOCC")
abline(v=0, col="red", lwd=3)
hist(unlist(sd.bias.msocc), xlim=c(0,20), main="", xlab="Standard deviation", ylab="")

# Means of Model intercepts and effect sizes (beta) of each estimator vs. siteCov (to assess bias)

##Bias in the model intercept (alpha)
par(mfrow = c(4,2), mai=c(0.2,0.2,0.2,0.2), mar=c(4,4.5,1,1), oma=c(1,1,0,0), las=1,cex.lab=2, cex.axis=2)

hist(unlist(m.alpha.n), xlim=c(0,50), main="", xlab="", ylab="Naive")
abline(v=mean(unlist(m.alpha.t)), col="red", lwd=3)
hist(unlist(sd.alpha.n), xlim=c(0,10), main="", xlab="",ylab="")

hist(unlist(m.alpha.j), xlim=c(0,50), main="", xlab="",ylab="Jackknife")
abline(v=mean(unlist(m.alpha.t)), col="red", lwd=3)
hist(unlist(sd.alpha.j), xlim=c(0,10), main="", xlab="",ylab="")

hist(unlist(m.alpha.c), xlim=c(0,50), main="", xlab="",ylab="Chao")
abline(v=mean(unlist(m.alpha.t)), col="red", lwd=3)
hist(unlist(sd.alpha.c), xlim=c(0,10), main="", xlab="",ylab="")

hist(unlist(m.alpha.msocc), xlim=c(0,50), main="", xlab="",ylab="MSOCC")
abline(v=mean(unlist(m.alpha.t)), col="red", lwd=3)
hist(unlist(sd.alpha.msocc), xlim=c(0,10), main="", xlab="",ylab="")

#Bias in slope
par(mfrow = c(4,2), mai=c(0.2,0.2,0.2,0.2), mar=c(4,4.5,1,1), oma=c(1,1,0,0), las=1,cex.lab=2, cex.axis=2)

hist(unlist(m.beta.n), xlim=c(-20,50), main="", xlab="", ylab="Naive")
abline(v=mean(unlist(m.beta.t)), col="red", lwd=3)
hist(unlist(sd.beta.n), xlim=c(0,20), main="", xlab="", ylab="")

hist(unlist(m.beta.j), xlim=c(-20,50), main="", xlab="", ylab="Jackknife")
abline(v=mean(unlist(m.beta.t)), col="red", lwd=3)
hist(unlist(sd.beta.j), xlim=c(0,20), main="", xlab="", ylab="")

hist(unlist(m.beta.c), xlim=c(-20,50), main="", xlab="", ylab="Chao")
abline(v=mean(unlist(m.beta.t)), col="red", lwd=3)
hist(unlist(sd.beta.c), xlim=c(0,20), main="", xlab="", ylab="")

hist(unlist(m.beta.msocc), xlim=c(-20,50), main="", xlab="Slope coefficient", ylab="MSOCC")
abline(v=mean(unlist(m.beta.t)), col="red", lwd=3)
hist(unlist(sd.beta.msocc), xlim=c(0,20), main="", xlab="Standard deviation", ylab="")


#####Comparison of the Occupancy - site covariate relationships for each species
#Plot the relationships between occupancy and siteCovariate for model
#tiff("D:\\Users\\lmcnew\\Documents\\CAE\\Manuscripts\\Species Richness Comparison\\Figures\\Psi.siteCov, P.siteCov negative\\No augmentation\\No group effect\\occupancy 400.tif", width=12, height=8, units="in", res=400) 

par(mfrow = c(2,1),mar=c(4,5,1,1), oma=c(1,1,1,1), las=1, cex.lab=2, cex.axis=2)

pred.psi.m <- array(NA, dim = c(length(siteCov), 40))
for(i in 1:length(siteCov)){
  pred.psi.m[i,] <- plogis(out$BUGSoutput$summary[51:90,1] + 
                             out$BUGSoutput$summary[91:130,1]*siteCov[i])
}

plot(siteCov, plogis(out$BUGSoutput$summary[51,1] + 
                       out$BUGSoutput$summary[91,1] * siteCov), ylim = c(0,1),
     lwd = 1, col = "black", type = "l", xlab = "",  ylab = expression(psi), las=1,cex.lab=1.5, cex.axis=1.5)  

as.data.frame(pred.psi.m)

for(i in 1:40){
  lines(siteCov, pred.psi.m[,i], col = i, lwd = 1)
}
for(i in 31:40){
  lines(siteCov, pred.psi.m[,i], col = i, lwd = 3)
} #Thick line for species with declining detect
for(i in 2:3){
  lines(siteCov, pred.psi.m[,i], col = i,lwd = 3)
} #Thick line for species with declining detection

###and compare that with truth.
##True relationship between species-specific site occupancy and siteCov

pred.psi <- array(NA, dim = c(length(siteCov), 40))
for(i in 1:length(siteCov)){
  pred.psi[i,] <- plogis(occ.table[,1]+occ.table[,2]*siteCov[i])
}


plot(siteCov, pred.psi[,1], type="l", xlab="Site covariate", ylab="True occupancy",  ylim=c(0,1),las=1, cex.lab=1.5, cex.axis=1.5)
for(i in 1:40){
  lines(siteCov, pred.psi[,i], col =i, lwd=1)
}
for(i in 31:40){
  lines(siteCov, pred.psi[,i], col =i, lwd=3)
} #Thick line for species with declining detection

for(i in 2:3){
  lines(siteCov, pred.psi[,i], col =i, lwd=3)
}  #Thick line for species with declining detection

###Both intercepts and slopes are pulled toward the mean

###Summary stats for bias of species richness from estimators
bias.msocc.95CI<-quantile(unlist(m.bias.msocc),c(0.025,0.975))
bias.msocc.95CI
bias.chao.95CI<-quantile(unlist(m.bias.chao),c(0.025,0.975))
bias.chao.95CI
bias.jack.95CI<-quantile(unlist(m.bias.jack),c(0.025,0.975))
bias.jack.95CI
bias.naive.95CI<-quantile(unlist(m.bias.naive),c(0.025,0.975))
bias.naive.95CI
bias.msocc.mean<-mean(unlist(m.bias.msocc))
bias.msocc.mean
bias.chao.mean<-mean(unlist(m.bias.chao))
bias.chao.mean
bias.jack.mean<-mean(unlist(m.bias.jack))
bias.jack.mean
bias.naive.mean<-mean(unlist(m.bias.naive))
bias.naiv