diversity_analysis <- function(){ # code modified by Greg Jacobs, 8/13/2020 # Original code for WinBUGS downloaded from url: # https://www.mbr-pwrc.usgs.gov/site/communitymodeling/software-code/ ################################################################################# #Read in the occurence data library(tidyr) library(dplyr) #data1. <- read.csv("dat7.csv") %>% data1. <- dat7 %>% arrange(td) %>% mutate(Species = sub("female", "", spp), Point = site, Rep = as.numeric(as.factor(td)), Occ = ct>0) data1.$Occ <- 1*(data1.$ct>0) #How many sightings for each species total.count. <- tapply(data1.$Occ, data1.$Species, sum) data1. <- subset(data1., Species %in% names(total.count.)[total.count.>0]) %>% arrange(Point, Species) #Find the number of unique species (n) uspecies <- as.character(unique(data1.$Species)) n=length(uspecies) #Find the number of unique sampling locations (J) upoints = as.character(unique(data1.$Point)) J=length(upoints) #Reshape the data using the R package "reshape" library(reshape) junk.melt=melt(data1.,id.var=c("Species", "Point", "Rep"), measure.var="Occ") X=cast(junk.melt, Point ~ Rep ~ Species) #Collapse occurrence histories so that NAs are at the end (lose leverage on rep dimension). for(i in 1:n){ for(j in 1:J){ X[j,,i]<-X[j,c(which(!is.na(X[j,,i])),which(is.na(X[j,,i]))),i] } } #Create all zero encounter histories to add to the detection array X #as part of the data augmentation to account for additional #species (beyond the n observed species). #naug is the number of all zero encounter histories to be added naug = 10 #X.zero is a matrix of zeroes, including the NAs for when a point has not been sampled X.zero <- X[,,1]*0 #Xaug is the augmented version of X. The first n species were actually observed #and the n+1 through naug species are all zero encounter histories Xaug <- array(0, dim=c(dim(X)[1],dim(X)[2],dim(X)[3]+naug)) Xaug[,,(dim(X)[3]+1):dim(Xaug)[3]] = rep(X.zero, naug) dimnames(X)=NULL Xaug[,,1:dim(X)[3]] <- X #K is a vector of length J indicating the number of reps at each point j KK <- X.zero a=which(KK==0); KK[a] <- 1 K=apply(KK,1,sum, na.rm=TRUE) K=as.vector(K) #Assemble data list sp.data = list(n=n, naug=naug, J=J, K=K, X=Xaug) ############################################################ ############################################################ # Basic Model - Dorazio and Royle (2005; DOI: 10.1198/016214505000000015), Dorazio et al. (2006) ############################################################ ############################################################ #Write the basic model code to a text file (used to run WinBUGS) cat(" model{ #Define prior distributions for community-level model parameters omega ~ dunif(0,1) mu.u ~ dnorm(0, 0.6666667) mu.v ~ dnorm(0, 0.6666667) tau.u ~ dgamma(0.1,0.1) tau.v ~ dgamma(0.1,0.1) sigma.u <- sqrt(1/tau.u) sigma.v <- sqrt(1/tau.v) for (i in 1:(n+naug)) { #Create priors for species i from the community level prior distributions w[i] ~ dbern(omega) u[i] ~ dnorm(mu.u, tau.u) v[i] ~ dnorm(mu.v, tau.v) #Create a loop to estimate the Z matrix (true occurrence for species i at point j. for (j in 1:J) { logit(psi[j,i]) <- u[i] mu.psi[j,i] <- psi[j,i]*w[i] Z[j,i] ~ dbern(mu.psi[j,i]) #Create a loop to estimate detection for species i at point k during #sampling period k. for (k in 1:K[j]) { logit(p[j,k,i]) <- v[i] mu.p[j,k,i] <- p[j,k,i]*Z[j,i] X[j,k,i] ~ dbern(mu.p[j,k,i]) } } } # Derived diversity measures # Estimated total richness (R) n0 <- sum(w[(n+1):(n+naug)]) R <- n + n0 R_exp <- omega*(n+naug) # site-level richness (r) for(j in 1:J){ r[j]<- sum(Z[j, ]) } #j # species-level site-occupancy (S) for(i in 1:(n+naug)){ Nsites[i] <- sum(Z[ ,i]) } #i # Diversity metrics (alhpa, beta, zeta) alpha <- mean(r[ ]) beta <- R/alpha for (i in 1:(n+naug)){ Zall[i] <- prod(Z[,i]) } zeta <- sum(Zall[]) } ",file="basicmodel.txt") ###################### # Fit basic model ###################### #Specify the parameters to be monitored sp.params.basic = list('u', 'v', 'mu.u', 'mu.v', 'tau.u', 'tau.v', 'omega', 'R', 'r', 'Nsites', 'alpha', 'beta', 'zeta', 'R_exp', 'sigma.u', 'sigma.v') # initialize the z matrix get.z <- function(Xaug){ a <- apply(Xaug, c(1,3), function(x)sum(x,na.rm=T)) a[a>0] <- 1 a[a==0] <- NA return(a) } # function to initialize parameters get.inits <- function(n.chains){ inits <- list() for (i in 1:n.chains){ omegaGuess = runif(1, n/(n+naug), 1) psi.meanGuess = runif(1, .25, 1) initsi <- list() initsi$omega=omegaGuess initsi$w=c(rep(1, n), rbinom(naug, size=1, prob=omegaGuess)) initsi$u=rnorm(n+naug) initsi$v=rnorm(n+naug) initsi$Z = get.z(Xaug) initsi$.RNG.name = "base::Wichmann-Hill" initsi$.RNG.seed = i inits[[i]] <- initsi } return(inits) } #Run the model and call the results “fit” library(runjags) fit <- run.jags(model = "basicmodel.txt", monitor = unlist(sp.params.basic), data = sp.data, n.chains = 3, inits = get.inits(3), burnin = 10000, sample = 20000, thin = 1, method = "parallel", keep.jags.files = TRUE) return(fit) }