diff --git a/.gitignore b/.gitignore index 7f9c92f..1bcccfa 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,6 @@ $RECYCLE.BIN/ # Operating System Files # ========================= .Rproj.user + +*OldFiles* +SimulationOutput/ diff --git a/DataAnalysis/CaseStudy/Model.R b/DataAnalysis/CaseStudy/Model.R new file mode 100644 index 0000000..d3a9195 --- /dev/null +++ b/DataAnalysis/CaseStudy/Model.R @@ -0,0 +1,343 @@ +#----------------# +#-Load Libraries-# +#----------------# + +library(nimble) +library(coda) + +#-----------# +#-Load Data-# +#-----------# + +load(file = "~/HMSNO/DataFormat/HMSNO.data.Rdata") +load(file = "~/HMSNO/DataFormat/HMSNO.con.Rdata") + +nspecs <- HMSNO.con$nspecs +parkS <- HMSNO.con$parkS +parkE <- HMSNO.con$parkE +nsite <- HMSNO.con$nsite +nstart <- HMSNO.con$nstart +nend <- HMSNO.con$nend +nparks <- max(parkE) +nyrs <- max(nend) + +#--------------# +#-NIMBLE model-# +#--------------# + +model.code <- nimbleCode({ + + #-Priors-# + + # Detection hyperparameters + mu.a0 ~ dunif(0, 1) # Community-level intercept on probability scale + mu.a0L <- logit(mu.a0) # Community-level intercept on logit scale + tau.a0 ~ dgamma(0.1, 0.1) # Community-level precision + + # Effect of effort (days sampled) + alpha1 ~ dnorm(0, 0.1) + + # Initial abundance hyperparameters + mu.b0 ~ dnorm(0, 0.1) # Community-level intercept on log scale + tau.b0 ~ dgamma(0.1, 0.1) # Community-level precision + + # Apparent survival hyperparameters + mu.o0 ~ dunif(0, 1) # Community-level intercept on probability scale + mu.o0L <- logit(mu.o0) # Community-level intercept on logit scale + tau.o0 ~ dgamma(0.1, 0.1) # Community-level precision + + # Gains hyperparameters + mu.g0 ~ dnorm(0, 0.1) # Community-level intercept on log scale + tau.g0 ~ dgamma(0.1, 0.1) # Community-level precision + + # Precison of park-level random effect on initial abundance + tau.eps.l ~ dgamma(0.1, 0.1) + + # Precison of park-level random effect on apparent survival + tau.eps.o ~ dgamma(0.1, 0.1) + + for(r in 1:nparks){ + + # Park-level random effect on initial abundance + eps.l[r] ~ dnorm(0, tau.eps.l) + + # Park-level random effect on apparent survival + eps.o[r] ~ dnorm(0, tau.eps.o) + + # Predicted community-level apparent survival for each park + logit(park.surv[r]) <- mu.o0L + eps.o[r] + + } # End r + + for(i in 1:nspecs){ + + # Species-specific intercept on detection + alpha0[i] ~ dnorm(mu.a0L, tau.a0) + + # Species-specific intercept on initial abundance + beta0[i] ~ dnorm(mu.b0, tau.b0) + + # Species-specific intercept on apparent survival + omega0[i] ~ dnorm(mu.o0L, tau.o0) + + # Species-specific intercept on gains + gamma0[i] ~ dnorm(mu.g0, tau.g0) + + for(r in parkS[i]:parkE[i]){ + + # Predicted population-level apparent survival + logit(pop.surv[r,i]) <- omega0[i] + eps.o[r] + + for(j in 1:nsite[r]){ + + # Linear predictor for initial abundance + log(lambda[j,r,i]) <- beta0[i] + eps.l[r] + + # Initial abundance + N[nstart[r],j,r,i] ~ dpois(lambda[j,r,i]) + + for(k in 1:nreps){ + + # Linear predictor for detection probability (in year 1) + logit(r[k,nstart[r],j,r,i]) <- alpha0[i] + alpha1 * days[k,nstart[r],j,i] + + # Observation process (in year 1) + y[k,nstart[r],j,r,i] ~ dbern(p[k,nstart[r],j,r,i]) + + # Site level detection (N-occupancy parameterization) (in year 1) + p[k,nstart[r],j,r,i] <- 1 - pow((1 - r[k,nstart[r],j,r,i]), N[nstart[r],j,r,i]) + + } # End k + + for(t in (nstart[r]+1):nend[r]){ + + for(k in 1:nreps){ + + # Linear predictor for detection probability (in year t+1) + logit(r[k,t,j,r,i]) <- alpha0[i] + alpha1 * days[k,t,j,i] + + # Observation process (in year t+1) + y[k,t,j,r,i] ~ dbern(p[k,t,j,r,i]) + + # Observation process (in year t+1) + p[k,t,j,r,i] <- 1 - pow((1 - r[k,t,j,r,i]), N[t,j,r,i]) + + } # End k + + # Linear predictor for apparent survival + logit(omega[t-1,j,r,i]) <- omega0[i] + eps.o[r] + + # Biological process + N[t,j,r,i] <- S[t-1,j,r,i] + G[t-1,j,r,i] + + # Apparent survival + S[t-1,j,r,i] ~ dbin(omega[t-1,j,r,i], N[t-1,j,r,i]) + + # Gains + G[t-1,j,r,i] ~ dpois(gamma[t-1,j,r,i]) + + # Linear predictor for gains + log(gamma[t-1,j,r,i]) <- gamma0[i] + + } # End t + + } # End j + + for(t in nstart[r]:nend[r]){ + + # Population (species-park) abundance per year + Nhat[t,r,i] <- sum(N[t,1:nsite[r],r,i]) + + } # End t + + } # End r + + } # End s + +}) + +#--------------# +#-Compile data-# +#--------------# + +Data <- list(y = HMSNO.data$y, days = HMSNO.data$days.scaled) + +#----------------# +#-Initial values-# +#----------------# + +dim2 <- dim1 <- dim(HMSNO.data$y)[-1] +dim2[1] <- dim2[1] - 1 + +Nst <- array(NA, dim = dim1) +Sst <- array(NA, dim = dim2) +Gst <- array(NA, dim = dim2) + +for(i in 1:nspecs){ + for(r in parkS[i]:parkE[i]){ + for(j in 1:nsite[r]){ + Nst[nstart[r],j,r,i] <- 10 + Sst[nstart[r]:(nend[r]-1),j,r,i] <- 5 + Gst[nstart[r]:(nend[r]-1),j,r,i] <- 2 + } + } +} + +alpha0.fun <- function(){ + alpha0 <- NULL + alpha0[1] <- runif(1, -1, 0) + alpha0[2] <- runif(1, -1.5, -1) + alpha0[3] <- runif(1, -1.5, -1) + alpha0[4] <- runif(1, -2.5, -1.5) + alpha0[5] <- runif(1, -1, -0.5) + alpha0[6] <- runif(1, -1.5, -0.5) + alpha0[7] <- runif(1, -0.7, -0.5) #Issues with this node + alpha0[8] <- runif(1, -2, -1.5) + alpha0[9] <- runif(1, -1, 0) + alpha0[10] <- runif(1, -3, -2) + alpha0[11] <- runif(1, -3, -2) + alpha0[12] <- runif(1, -1.1, -0.9) #Issues with this node + return(alpha0) +} + +beta0.fun <- function(){ + beta0 <- NULL + beta0[1] <- runif(1, 0, 1) + beta0[2] <- runif(1, -1, 1) + beta0[3] <- runif(1, 1, 3) + beta0[4] <- runif(1, -8, -4) + beta0[5] <- runif(1, 0, 2) + beta0[6] <- runif(1, 0, 2) + beta0[7] <- runif(1, -3, -1) #Issues with this node + beta0[8] <- runif(1, 1, 2) + beta0[9] <- runif(1, -3, -1) + beta0[10] <- runif(1, 0, 2) + beta0[11] <- runif(1, -5, -2) + beta0[12] <- runif(1, -2, 1) #Issues with this node + return(beta0) +} + +omega0.fun <- function(){ + omega0 <- NULL + omega0[1] <- runif(1, -0.5, 0.5) + omega0[2] <- runif(1, -2, 0) + omega0[3] <- runif(1, 0, 2) + omega0[4] <- runif(1, -1, 1) + omega0[5] <- runif(1, -1, 1) + omega0[6] <- runif(1, 0, 2) + omega0[7] <- runif(1, 1, 2) #Issues with this node + omega0[8] <- runif(1, 0, 2) + omega0[9] <- runif(1, -2, -1) + omega0[10] <- runif(1, -1, 1) + omega0[11] <- runif(1, -2, 2) + omega0[12] <- runif(1, 1, 3) + return(omega0) +} + +gamma0.fun <- function(){ + gamma0 <- NULL + gamma0[1] <- runif(1, 0, 0.5) + gamma0[2] <- runif(1, -0.2, 0.2) + gamma0[3] <- runif(1, -1.5, -0.5) + gamma0[4] <- runif(1, -2.5, -2) + gamma0[5] <- runif(1, -0.2, 0.2) + gamma0[6] <- runif(1, -2.5, -1.5) + gamma0[7] <- runif(1, -1.8, -1.4) + gamma0[8] <- runif(1, -1, -0.5) + gamma0[9] <- runif(1, -3, 2.5) + gamma0[10] <- runif(1, -3, -1) + gamma0[11] <- runif(1, -5, -3) + gamma0[12] <- runif(1, -2.5, -1.5) + return(gamma0) +} + +eps.l.fun <- function(){ + eps.l <- NULL + eps.l[1] <- runif(1, -2, 0) + eps.l[2] <- runif(1, 2, 4) + eps.l[3] <- runif(1, -2, 0) + eps.l[4] <- runif(1, 0, 2) + eps.l[5] <- runif(1, 0, 2) + eps.l[6] <- runif(1, -2, 0) + return(eps.l) +} + +eps.o.fun <- function(){ + eps.o <- NULL + eps.o[1] <- runif(1, -2, 0) + eps.o[2] <- runif(1, 0, 2) + eps.o[3] <- runif(1, -2, 0) + eps.o[4] <- runif(1, -1, 1) + eps.o[5] <- runif(1, -1, 1) + eps.o[6] <- runif(1, -2, 0) + return(eps.o) +} + +inits <- function()list(N=Nst, S=Sst, G=Gst, + mu.a0 = runif(1, 0.15, 0.3), tau.a0 = runif(1, 1, 3), + mu.b0 = runif(1, -2, 1), tau.b0 = runif(1, 0.2, 0.4), + mu.o0 = runif(1, 0.4, 0.7), tau.o0 = runif(1, 0, 1.5), + mu.g0 = runif(1, -2, -1), tau.g0 = runif(1, 0.25, 1), + gamma0 = gamma0.fun(), + alpha0 = alpha0.fun(), alpha1 = runif(1, 0.2, 0.4), + beta0 = beta0.fun(), + omega0 = omega0.fun(), + eps.l = eps.l.fun(), + eps.o = eps.o.fun(), + tau.eps.l = runif(1, 0, 1), + tau.eps.o = runif(1, 0, 1) +) + +#--------------------# +#-Parameters to save-# +#--------------------# + +params <- c("mu.a0", "tau.a0", "alpha0", "alpha1", + "mu.b0", "tau.b0", "beta0", "eps.l", "tau.eps.l", + "mu.o0", "tau.o0", "omega0", "eps.o", "tau.eps.o", + "mu.g0", "tau.g0", "gamma0", + "park.surv", "pop.surv", "Npark") + +#---------------# +#-MCMC settings-# +#---------------# + +model <- nimbleModel(model.code, constants = HMSNO.con, data = Data, inits = inits()) + +MCMCconf <- configureMCMC(model, monitors = params) + + +MCMCconf$addSampler(target = c("beta0[3]", "beta0[6]", "beta0[9]", "beta0[10]", "eps.l[1]"), + type = "AF_slice") + +MCMCconf$addSampler(target = c("beta0[1]", "beta0[2]", "beta0[4]", "beta0[5]", "beta0[11]", "eps.l[5]"), + type = "AF_slice") + +MCMCconf$addSampler(target = c("beta0[5]", "beta0[8]", "beta0[12]", "eps.l[6]"), + type = "AF_slice") + +MCMCconf$addSampler(target = c("omega0[3]", "omega0[6]", "omega0[9]", "omega0[10]", "eps.o[1]"), + type = "AF_slice") + +MCMCconf$addSampler(target = c("omega0[1]", "omega0[2]", "omega0[4]", "omega0[5]", "omega0[11]", "eps.o[5]"), + type = "AF_slice") + +MCMCconf$addSampler(target = c("omega0[5]", "omega0[8]", "omega0[12]", "eps.o[6]"), + type = "AF_slice") + +MCMC <- buildMCMC(MCMCconf) + +model.comp <- compileNimble(model, MCMC) + +ni <- 100000 +nb <- 75000 +nc <- 1 +nt <- 25 + +#Run NIMBLE model +out <- runMCMC(model.comp$MCMC, niter = ni, nburnin = nb, nchains = nc, thin = nt, samplesAsCodaMCMC = TRUE) + +#-Save output-# +ID <- paste("chain", length(list.files(path = "~/HMSNO/DataAnalysis/CaseStudy", pattern = "chain", full.names = FALSE)) + 1, sep="") +assign(ID, out) +save(list = ID, file = paste0(ID, ".Rds")) \ No newline at end of file diff --git a/DataAnalysis/CaseStudy/chain1.Rds b/DataAnalysis/CaseStudy/chain1.Rds new file mode 100644 index 0000000..7c54946 Binary files /dev/null and b/DataAnalysis/CaseStudy/chain1.Rds differ diff --git a/DataAnalysis/CaseStudy/chain2.Rds b/DataAnalysis/CaseStudy/chain2.Rds new file mode 100644 index 0000000..ebd7179 Binary files /dev/null and b/DataAnalysis/CaseStudy/chain2.Rds differ diff --git a/DataAnalysis/CaseStudy/chain3.Rds b/DataAnalysis/CaseStudy/chain3.Rds new file mode 100644 index 0000000..727207a Binary files /dev/null and b/DataAnalysis/CaseStudy/chain3.Rds differ diff --git a/DataAnalysis/CaseStudy/chain4.Rds b/DataAnalysis/CaseStudy/chain4.Rds new file mode 100644 index 0000000..d5962b1 Binary files /dev/null and b/DataAnalysis/CaseStudy/chain4.Rds differ diff --git a/DataAnalysis/CaseStudy/chain5.Rds b/DataAnalysis/CaseStudy/chain5.Rds new file mode 100644 index 0000000..47b49ae Binary files /dev/null and b/DataAnalysis/CaseStudy/chain5.Rds differ diff --git a/DataAnalysis/HMSNO.R b/DataAnalysis/HMSNO.R deleted file mode 100644 index f1dc63f..0000000 --- a/DataAnalysis/HMSNO.R +++ /dev/null @@ -1,231 +0,0 @@ -#-------------------# -#-Working Directory-# -#-------------------# - -setwd("./DataFormat") - -#----------------# -#-Load Libraries-# -#----------------# - -library(readxl) -library(dplyr) -library(tidyr) -library(reshape2) -library(jagsUI) - -#-----------# -#-Load Data-# -#-----------# - -Data <- NULL -files <- list.files(path = "./EditedData", pattern = "Data", full.names = TRUE) -for(i in 1:length(files)){ - sheets <- excel_sheets(files[i]) - for(j in 1:length(sheets)){ - Temp <- read_excel(files[i], sheet = sheets[j]) - Temp <- Temp %>% mutate(species = as.character(sheets[j]), park = as.factor(i)) %>% - select(park, `deployment ID`, Year, species, ls(Temp, pattern = "R"), density) - Data <- bind_rows(Data, Temp) - } -} - -covfile <- "./RawData/camera trap metadata.xlsx" -covsheets <- excel_sheets(covfile) -covsheets <- covsheets[-3] #Remove Nyungwe - -CovData <- NULL -for(i in 1:5){ - Temp <- read_excel(covfile, sheet = covsheets[i]) - Temp <- Temp %>% select(`deployment ID`, days, elevation, edge, Year) - CovData <- bind_rows(CovData, Temp) -} - -CovData <- bind_rows(CovData, read_excel("./EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 1) %>% - select(`deployment ID`, days, elevation, edge, Year)) %>% - bind_rows(., read_excel("./EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 2) %>% - select(`deployment ID`, days, elevation, edge, Year)) %>% - distinct() - -Data <- full_join(Data, CovData, by = c("deployment ID", "Year")) %>% drop_na() - -#-------------# -#-Format Data-# -#-------------# - -#Convert to NAs -Data[,5:10] <- apply(Data[,5:10], 2, function(y) as.numeric(gsub("-", NA, y))) - -#Remove cameras that were sampled in only 1 year -Data <- Data %>% group_by(`deployment ID`) %>% filter(n() > 2) %>% ungroup() - -#Convert to single vector format -Data2 <- melt(Data, id = c("park", "deployment ID", "Year", "species", "days", "elevation", "edge")) - -#Occupancy data -y <- Data2$value - -#Indices -nyrs <- Data %>% group_by(`deployment ID`) %>% summarize(min(Year), max(Year)) -nstart <- as.numeric(nyrs$`min(Year)`) - as.numeric(min(nyrs$`min(Year)`)) + 1 -nend <- 1 + as.numeric(max(nyrs$`max(Year)`)) - as.numeric(min(nyrs$`min(Year)`)) - - (as.numeric(max(nyrs$`max(Year)`)) - as.numeric(nyrs$`max(Year)`)) - -nyrs <- as.numeric(max(nyrs$`max(Year)`)) - as.numeric(min(nyrs$`min(Year)`)) + 1 - -nsites <- as.numeric(Data %>% summarize(n_distinct(`deployment ID`))) -nspec <- as.numeric(Data %>% summarize (n_distinct(species))) -nobs <- length(y) - -#Nested indices -yr <- as.numeric(as.factor(Data2$Year)) -site <- as.numeric(as.factor(Data2$`deployment ID`)) -spec <- as.numeric(as.factor(Data2$species)) - -#Covariates -Cov <- Data %>% group_by(`deployment ID`, Year) %>% - distinct(`deployment ID`, Year, .keep_all = TRUE) -Cov$elevation <- scale(as.numeric(Cov$elevation)) -Cov$edge <- scale(as.numeric(Cov$edge)) -Cov$days <- scale(as.numeric(Cov$days)) - -elev <- dcast(Cov %>% select(`deployment ID`, Year, elevation), `deployment ID` ~ Year) -elev <- t(as.matrix(elev[,-1])) -edge <- dcast(Cov %>% select(`deployment ID`, Year, edge), `deployment ID` ~ Year) -edge <- t(as.matrix(edge[,-1])) -days <- dcast(Cov %>% select(`deployment ID`, Year, days), `deployment ID` ~ Year) -days <- t(as.matrix(days[,-1])) -density <- Data %>% distinct(`deployment ID`, density) -density <- as.numeric(density$density) - -elev[is.na(elev)] <- 0 -edge[is.na(edge)] <- 0 -days[is.na(days)] <- 0 - -#--------------# -#-Compile Data-# -#--------------# - -#Data -jags.data <- list(y = y, nobs = nobs, nstart = nstart, nend = nend, nsites = nsites, nspec = nspec, - yr = yr, site = site, spec = spec, elev = elev, edge = edge, density = density, days = days) - -#Initial values -Nst <- array(NA, dim = c(nyrs, nsites, nspec)) -Sst <- array(NA, dim = c(nyrs-1, nsites, nspec)) -Gst <- array(NA, dim = c(nyrs-1, nsites, nspec)) -for(j in 1:nsites){ - Nst[nstart[j],j,] <- 10 - Sst[nstart[j]:(nend[j]-1),j,] <- 5 - Gst[nstart[j]:(nend[j]-1),j,] <- 2 -} - -alpha0.fun <- function(){ - alpha0 <- NULL - alpha0[1] <- runif(1, -0.2, 0.2) - alpha0[2] <- runif(1, -4, -2.75) - alpha0[3] <- runif(1, -2, -1.5) - alpha0[4] <- runif(1, -5, -3) - alpha0[5] <- runif(1, -1, -0.6) - alpha0[6] <- runif(1, -1.5, -1) - alpha0[7] <- runif(1, -2, -1.5) - alpha0[8] <- runif(1, -4, -2) - alpha0[9] <- runif(1, -2, -1.5) - return(alpha0) -} - -beta0.fun <- function(){ - beta0 <- NULL - beta0[1] <- runif(1, 0.4, 0.8) - beta0[2] <- runif(1, 1.5, 3.5) - beta0[3] <- runif(1, 0, 1) - beta0[4] <- runif(1, -7, -4) - beta0[5] <- runif(1, 0.5, 1) - beta0[6] <- runif(1, -0.5, 0) - beta0[7] <- runif(1, 1, 3) - beta0[8] <- runif(1, 0, 3) - beta0[9] <- runif(1, 0.5, 1) - return(beta0) -} - -beta4.fun <- function(){ - beta4 <- NULL - beta4[1] <- runif(1, -0.5, 0.5) - beta4[2] <- runif(1, 0, 0.5) - beta4[3] <- runif(1, -0.5, 0) - beta4[4] <- runif(1, -2, 2) - beta4[5] <- runif(1, 0, 0.6) - beta4[6] <- runif(1, 1, 1.5) - beta4[7] <- runif(1, 0, 2) - beta4[8] <- runif(1, 1, 2) - beta4[9] <- runif(1, 0.4, 0.6) - return(beta4) -} - -omega0.fun <- function(){ - omega0 <- NULL - omega0[1] <- runif(1, 1.4, 2.2) - omega0[2] <- runif(1, -5, -3) - omega0[3] <- runif(1, 0.5, 1) - omega0[4] <- runif(1, -2, -0.4) - omega0[5] <- runif(1, 0.6, 1) - omega0[6] <- runif(1, 2, 4) - omega0[7] <- runif(1, 0, 1) - omega0[8] <- runif(1, 0, 3) - omega0[9] <- runif(1, 2, 3) - return(omega0) -} - -# gamma0.fun <- function(){ -# gamma0 <- NULL -# gamma0[1] <- runif(1, 0, 1) -# gamma0[2] <- runif(1, -0.2, 0.2) -# gamma0[3] <- runif(1, -1.5, -0.5) -# gamma0[4] <- runif(1, -2.5, -1.5) -# gamma0[5] <- runif(1, -0.6, -0.2) -# gamma0[6] <- runif(1, -3, -2) -# gamma0[7] <- runif(1, -1.5, -0.5) -# gamma0[8] <- runif(1, -3, -1) -# gamma0[9] <- runif(1, -3, -2) -# return(gamma0) -# } - -inits <- function()list(N=Nst, S=Sst, G=Gst, - mu.a0 = runif(1, 0.1, 0.25), tau.a0 = runif(1, 0, 5), alpha0 = alpha0.fun(), alpha1 = runif(1, -0.25, 0.25), - mu.b0 = runif(1, 0, 3), tau.b0 = runif(1, 0, 1), beta0 = beta0.fun(), - mu.b1 = runif(1, -1, 1), tau.b1 = runif(1, 0, 10), beta1 = runif(nspec, -1, 1), - mu.b2 = runif(1, -1, 1), tau.b2 = runif(1, 0, 10), beta2 = runif(nspec, -0.5, 0.5), - mu.b3 = runif(1, -1, 1), tau.b3 = runif(1, 0, 10), beta3 = runif(nspec, -1, 1), - mu.b4 = runif(1, 0, 2), tau.b4 = runif(1, 0, 2), beta4 = beta4.fun(), - mu.o0 = runif(1, 0.5, 0.9), tau.o0 = runif(1, 0, 1), omega0 = omega0.fun(), - mu.o1 = runif(1, 0, 4), tau.o1 = runif(1, 0, 1), omega1 = runif(nspec, -5, 5), - mu.o2 = runif(1, -1, 1), tau.o2 = runif(1, 0, 10), omega2 = runif(nspec, -1, 1), - mu.o3 = runif(1, -1, 1), tau.o3 = runif(1, 0, 10), omega3 = runif(nspec, -1, 1), - gamma0 = runif(1, -4, 0) -) - -#Parameters to save -params <- c("mu.a0", "tau.a0", - "mu.b0", "tau.b0", "mu.b1", "tau.b1", "mu.b2", "tau.b2", "mu.b3", "tau.b3", "mu.b4", "tau.b4", - "mu.o0", "tau.o0", "mu.o1", "tau.o1", "mu.o2", "tau.o2", "mu.o3", "tau.o3", - "alpha0", "alpha1", "beta0", "beta1", "beta2", "beta3", "beta4", - "omega0", "omega1", "omega2", "omega3", "gamma0") - - -#MCMC settings -# ni <- 115000 -# na <- 1000 -# nb <- 100000 -# nc <- 3 -# nt <- 5 - -na <- 1000 -ni <- 100000 -nb <- 20000 -nc <- 1 -nt <- 40 - -#Run JAGS model -out <- jagsUI(jags.data, inits, params, model.file = "../MSdyn_cov.txt", - n.adapt = na, n.chains = nc, n.iter = ni, n.burnin = nb, n.thin = nt, - parallel = FALSE) diff --git a/DataAnalysis/HMSNO_nimble.R b/DataAnalysis/HMSNO_nimble.R deleted file mode 100644 index 7047319..0000000 --- a/DataAnalysis/HMSNO_nimble.R +++ /dev/null @@ -1,383 +0,0 @@ -#----------------# -#-Load Libraries-# -#----------------# - -library(nimble) -library(coda) - -#-----------# -#-Load Data-# -#-----------# - -load(file = "~/HMSNO/DataFormat/HMSNO.data.Rdata") -load(file = "~/HMSNO/DataFormat/HMSNO.con.Rdata") - -#---------# -#-Indices-# -#---------# - -nspec <- HMSNO.con$nspec -nparks <- max(HMSNO.data$park) -npark <- HMSNO.con$npark -nsites <- max(HMSNO.con$site) -nsite <- HMSNO.con$nsite -nyrs <- max(HMSNO.con$yr) -nstart <- HMSNO.data$nstart -nend <- HMSNO.data$nend -parkspec <- HMSNO.con$parkspec -sitespec <- HMSNO.con$sitespec - -HMSNO.con$nyrs <- nyrs - -#--------------# -#-NIMBLE model-# -#--------------# - -MSdyn.c <- nimbleCode({ - - mu.b0 ~ dnorm(0, 0.1) - tau.b0 ~ dgamma(0.1, 0.1) - mu.o0L <- logit(mu.o0) - mu.o0 ~ dunif(0, 1) - tau.o0 ~ dgamma(0.1, 0.1) - mu.o1 ~ dnorm(0, 0.1) - tau.o1 ~ dgamma(0.1, 0.1) - mu.o2 ~ dnorm(0, 0.1) - tau.o2 ~ dgamma(0.1, 0.1) - # mu.o3 ~ dnorm(0, 0.1) - # tau.o3 ~ dgamma(0.1, 0.1) - mu.g0 ~ dnorm(0, 0.1) - tau.g0 ~ dgamma(0.1, 0.1) - tau.phi ~ dgamma(0.1, 0.1) - # log(gamma) <- gamma0 - # gamma0 ~ dnorm(0, 0.1) - - alpha0 ~ dunif(0, 1) - alpha1 ~ dnorm(0, 0.1) - - for(t in 1:nyrs){ - eps.o[t] ~ T(dnorm(0, tau.phi), -5, 5) - }#end t - - for(s in 1:nspec){ - - omega0[s] ~ dnorm(mu.o0L, tau.o0) - omega1[s] ~ dnorm(mu.o1, tau.o1) - omega2[s] ~ dnorm(mu.o2, tau.o2) - # omega3[s] ~ dnorm(mu.o3, tau.o3) - gamma0[s] ~ dnorm(mu.g0, tau.g0) - log(gamma[s]) <- gamma0[s] - - for(j in sitespec[1,s]:sitespec[nsite[s],s]){ - logit(r[ns[j],j,s]) <- logit(alpha0) + alpha1 * days[ns[j],j] - N[ns[j],j,s] ~ dpois(lambda[park[j],s]) - for(t in (ns[j]+1):ne[j]){ - - logit(r[t,j,s]) <- alpha0 + alpha1 * days[t,j] - - logit(omega[t-1,j,s]) <- omega0[s] + omega1[s] * density[j] + omega2[s] * edge[t,j] - - N[t,j,s] <- S[t-1,j,s] + G[t-1,j,s] - S[t-1,j,s] ~ dbin(omega[t-1,j,s], N[t-1,j,s]) - G[t-1,j,s] ~ dpois(gamma[s]) - # G[t-1,j,s] ~ dpois(gamma) - - }#end t - }#end j - - for(i in 1:npark[s]){ - beta0[parkspec[i,s],s] ~ dnorm(mu.b0, tau.b0) - log(lambda[parkspec[i,s],s]) <- beta0[parkspec[i,s],s] - for(t in nyrstr[parkspec[i,s]]:nyrend[parkspec[i,s]]){ - Npark[t,parkspec[i,s],s] <- sum(N[t,nsitestr[parkspec[i,s]]:nsiteend[t,parkspec[i,s]],s]) - }#end t - }#end i - - }#end s - - - for(k in 1:nobs){ - y[k] ~ dbern(p[k]) - p[k] <- 1 - pow((1 - r[yr[k], site[k], spec[k]]), N[yr[k], site[k], spec[k]]) - }#end k - -}) - -#---------------# -#-Inital values-# -#---------------# - -Nst <- array(NA, dim = c(nyrs, nsites, nspec)) -Sst <- array(NA, dim = c(nyrs-1, nsites, nspec)) -Gst <- array(NA, dim = c(nyrs-1, nsites, nspec)) -for(s in 1:nspec){ - for(j in 1:nsite[s]){ - jj <- sitespec[j,s] - Nst[nstart[jj],jj,s] <- 10 - Sst[nstart[jj]:(nend[jj]-1),jj,s] <- 5 - Gst[nstart[jj]:(nend[jj]-1),jj,s] <- 2 - } -} -# for(j in 1:nsites){ -# Nst[nstart[j],j,] <- 10 -# Sst[nstart[j]:(nend[j]-1),j,] <- 5 -# Gst[nstart[j]:(nend[j]-1),j,] <- 2 -# } - -# alpha0.fun <- function(){ -# alpha0 <- NULL -# alpha0[1] <- runif(1, -1.5, -1) -# alpha0[2] <- runif(1, -1.5, -1) -# alpha0[3] <- runif(1, -2, -1.5) -# alpha0[4] <- runif(1, -2.5, -1.5) -# alpha0[5] <- runif(1, -1.5, -1) -# alpha0[6] <- runif(1, -1.5, -1) -# alpha0[7] <- runif(1, -2, -1.5) -# alpha0[8] <- runif(1, -4, -2) -# alpha0[9] <- runif(1, -2, -1.5) -# return(alpha0) -# } - -# beta0.fun <- function(){ -# beta0 <- NULL -# beta0[1] <- runif(1, 1.5, 2.5) -# beta0[2] <- runif(1, 0.5, 1.25) -# beta0[3] <- runif(1, 0, 0.75) -# beta0[4] <- runif(1, -5, -3) -# beta0[5] <- runif(1, 1.5, 2.25) -# beta0[6] <- runif(1, -0.25, 0.5) -# beta0[7] <- runif(1, 0.75, 1.5) -# beta0[8] <- runif(1, -0.75, 0) -# beta0[9] <- runif(1, -0.25, 0.5) -# return(beta0) -# } - -beta0.fun <- function(){ - beta0 <- matrix(NA, nrow = nparks, ncol = nspec) - # for(s in 1:nspec){ - # for(i in specPark[1:specNpark[s],s]){ - # beta0[i,s] <- runif(1, -1, 1) - # } - # } - for(s in 1:nspec){ - for(i in 1:npark[s]){ - ii <- parkspec[i,s] - beta0[ii,s] <- runif(1, -1, 2) - } - } - return(beta0) -} - -# beta4.fun <- function(){ -# beta4 <- NULL -# beta4[1] <- runif(1, 0, 1) -# beta4[2] <- runif(1, 0, 2) -# beta4[3] <- runif(1, -1, 0) -# beta4[4] <- runif(1, 0, 4) -# beta4[5] <- runif(1, 2, 3) -# beta4[6] <- runif(1, 1, 1.5) -# beta4[7] <- runif(1, 0, 2) -# beta4[8] <- runif(1, 1, 2) -# beta4[9] <- runif(1, 0.4, 0.6) -# return(beta4) -# } - -omega0.fun <- function(){ - omega0 <- NULL - omega0[1] <- runif(1, 0.2, 0.4) - omega0[2] <- runif(1, 0, 0.25) - omega0[3] <- runif(1, 1, 1.5) - omega0[4] <- runif(1, -1.5, -0.5) - omega0[5] <- runif(1, 0.4, 0.6) - omega0[6] <- runif(1, 1.9, 2.2) - omega0[7] <- runif(1, -0.25, 0) - omega0[8] <- runif(1, -0.5, 0) - omega0[9] <- runif(1, 2, 2.25) - return(omega0) -} - -gamma0.fun <- function(){ - gamma0 <- NULL - gamma0[1] <- runif(1, 0, 1) - gamma0[2] <- runif(1, -0.2, 0.2) - gamma0[3] <- runif(1, -1.5, -0.5) - gamma0[4] <- runif(1, -2.5, -1.5) - gamma0[5] <- runif(1, -0.6, -0.2) - gamma0[6] <- runif(1, -3, -2) - gamma0[7] <- runif(1, -1.5, -0.5) - gamma0[8] <- runif(1, -3, -1) - gamma0[9] <- runif(1, -3, -2) - return(gamma0) -} - -# inits <- function()list(N=Nst, S=Sst, G=Gst, -# mu.a0 = runif(1, 0.1, 0.25), tau.a0 = runif(1, 0, 5), alpha0 = alpha0.fun(), alpha1 = runif(1, -0.25, 0), -# mu.b0 = runif(1, 0, 3), tau.b0 = runif(1, 0, 1), beta0 = beta0.fun(), -# mu.b1 = runif(1, -1, 1), tau.b1 = runif(1, 0, 10), beta1 = runif(nspec, -1, 1), -# mu.b2 = runif(1, -1, 1), tau.b2 = runif(1, 0, 10), beta2 = runif(nspec, -0.5, 0.5), -# mu.b3 = runif(1, -1, 1), tau.b3 = runif(1, 0, 10), beta3 = runif(nspec, -1, 1), -# mu.b4 = runif(1, 0, 2), tau.b4 = runif(1, 0, 2), beta4 = beta4.fun(), -# mu.o0 = runif(1, 0.5, 0.9), tau.o0 = runif(1, 0, 1), omega0 = omega0.fun(), -# mu.o1 = runif(1, 0, 4), tau.o1 = runif(1, 0, 1), omega1 = runif(nspec, -5, 5), -# mu.o2 = runif(1, -1, 1), tau.o2 = runif(1, 0, 10), omega2 = runif(nspec, -1, 1), -# mu.o3 = runif(1, -1, 1), tau.o3 = runif(1, 0, 10), omega3 = runif(nspec, -1, 1), -# mu.g0 = runif(1, -2, 0), tau.g0 = runif(1, 0, 1), gamma0 = gamma0.fun() -# ) - -inits <- function()list(N=Nst, S=Sst, G=Gst, - alpha0 = runif(1, 0.15, 0.2), alpha1 = runif(1, -0.5, -0.45), - mu.b0 = runif(1, -0.3, 0.75), #sig.b0 = runif(1, 1, 2), beta0 = beta0.fun(), - beta0 = beta0.fun(), tau.b0 = runif(1, 0, 1), tau.phi = runif(1, 0, 1), - mu.o0 = runif(1, 0.5, 0.75), tau.o0 = runif(1, 0.5, 1.25), omega0 = omega0.fun(), - mu.o1 = runif(1, 0, 4), tau.o1 = runif(1, 0, 1), omega1 = runif(nspec, -5, 5), - mu.o2 = runif(1, -1, 1), tau.o2 = runif(1, 0, 10), omega2 = runif(nspec, -1, 1), - # mu.o3 = runif(1, -1, 1), tau.o3 = runif(1, 0, 10), omega3 = runif(nspec, -1, 1), - mu.g0 = runif(1, -2, 0), tau.g0 = runif(1, 0, 1), gamma0 = gamma0.fun() -) - -#Parameters to save -params <- c("mu.b0", "tau.b0", "mu.o0", "tau.o0", - "mu.o1", "tau.o1", "mu.o2", "tau.o2", - # "mu.o3", "tau.o3", - "mu.g0", "tau.g0", - "tau.phi", - "alpha0", "alpha1", "beta0", - "omega0", "omega1", "omega2", - # "omega3", - "gamma0", "Npark") - -#MCMC settings -MSdyn.m <- nimbleModel(MSdyn.c, constants = HMSNO.con, data = HMSNO.data, inits = inits()) - -MCMCconf <- configureMCMC(MSdyn.m, monitors = params) -#MCMCconf$printSamplers(1:122) - -# MCMCconf$removeSampler(c('alpha0', 'alpha1', 'gamma0', -# 'beta0', 'beta1', 'beta2', 'beta3', 'beta4', -# 'omega0', 'omega1', 'omega2', 'omega3')) -# -# MCMCconf$addSampler(target = c('mu.a0', 'alpha0', 'alpha1'), -# type = "AF_slice") -# -# MCMCconf$addSampler(target = c('mu.b0', 'beta0'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.b1', 'beta1'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.b2', 'beta2'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.b3', 'beta3'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.b4', 'beta4'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.o0', 'omega0'), -# type = "AF_slice") -# -# MCMCconf$addSampler(target = c('mu.o1', 'omega1'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.o2', 'omega2'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.o3', 'omega3'), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c('mu.g0', 'gamma0'), -# type = "RW_block") -# -# nFun <- function(node, i){ -# paste(node, "[", i, "]", sep = "") -# } -# -# for(i in 1:nspec){ -# MCMCconf$addSampler(target = c(nFun("beta0", i), nFun("beta1", i), -# nFun("beta2", i), nFun("beta3", i), -# nFun("beta4", i)), -# type = "RW_block") -# -# MCMCconf$addSampler(target = c(nFun("omega0", i), nFun("omega1", i), -# nFun("omega2", i), nFun("omega3", i)), -# type = "AF_slice") -# } -# -# MCMCconf$samplerExecutionOrder <- c(seq(1,22,1), seq(52, 41641, 1), seq(23, 51, 1)) - -MCMCconf$addSampler(target = c('alpha0', 'alpha1'), type = "AF_slice") - -species <- seq(1, nspec, 1) - -nFun <- function(node, i){ - paste(node, "[", i, "]", sep = "") -} - -omit <- NULL - -# MCMCconf$addSampler(target = c("beta0[1]", "beta0[2]", "beta0[3]", "beta0[4]", "beta0[5]", -# "beta0[6]", "beta0[7]", "beta0[8]", "beta0[9]"), -# type = "RW_block") - -MCMCconf$addSampler(target = c("omega0[1]", "omega0[2]", "omega0[3]", "omega0[4]", "omega0[5]", - "omega0[6]", "omega0[7]", "omega0[8]", "omega0[9]"), - type = "RW_block") - -for(i in 1:nspec){ - omit <- c(omit, i) - if(i != nspec){ - for(j in species[-omit]){ - # MCMCconf$addSampler(target = c(nFun("beta0", i), nFun("beta0", j)), - # type = "RW_block") - MCMCconf$addSampler(target = c(nFun("omega0", i), nFun("omega0", j)), - type = "RW_block") - MCMCconf$addSampler(target = c(nFun("gamma0", i), nFun("gamma0", j)), - type = "RW_block") - } - } - # MCMCconf$addSampler(target = c("mu.b0", nFun("beta0", i)), - # type = "RW_block") - MCMCconf$addSampler(target = c("mu.o0", nFun("omega0", i)), - type = "AF_slice") - MCMCconf$addSampler(target = c("mu.g0", nFun("gamma0", i)), - type = "AF_slice") - MCMCconf$addSampler(target = c(nFun("omega0", i), nFun("omega1", i), - nFun("omega2", i)), - type = "AF_slice") - # MCMCconf$addSampler(target = c(nFun("omega0", i), nFun("omega1", i), - # nFun("omega2", i), nFun("omega3", i)), - # type = "AF_slice") -} - -#MCMCconf$addSampler(target = c('mu.b0', 'beta0'), type = "RW_block") - -#MCMCconf$addSampler(target = c('mu.o0', 'omega0'), type = "AF_slice") - -MCMC <- buildMCMC(MCMCconf) - -MSdyn.comp <- compileNimble(MSdyn.m, MCMC) - -ni <- 10000 -nb <- 5000 -nc <- 3 -nt <- 5 - -#Run NIMBLE model -MSdyn.o <- runMCMC(MSdyn.comp$MCMC, niter = ni, nburnin = nb, nchains = nc, thin = nt, samplesAsCodaMCMC = TRUE) - -#-Save output-# -ID <- gsub(" ","_",Sys.time()) -ID <- gsub(":", "-", ID) -save(MSdyn.o, file = paste("output", ID, ".Rdata", sep="")) - -#traceplot(out[c(1:3)][,"Npark[1, 1, 1]"]) -#plot(unlist(output[c(1:5)][,"beta0[1]"]), unlist(output[c(1:5)][,"beta0[4]"]) - -for(i in 1:npark){ - plot(nyrstr[i], mean(unlist(output[c(1:5)][,paste("Npark[", nyrstr[i], ", ", i, ", 8]", sep = "")])), - xlim = c(1,9), ylim = c(0,500), xlab = "Abundance", ylab = "Year") - for(t in (nyrstr[i] + 1):nyrend[i]){ - points(t, mean(unlist(output[c(1:5)][,paste("Npark[", t, ", ", i, ", 8]", sep = "")]))) - } -} diff --git a/DataAnalysis/SSdyn_nimble.R b/DataAnalysis/SSdyn_nimble.R deleted file mode 100644 index ae7cd6a..0000000 --- a/DataAnalysis/SSdyn_nimble.R +++ /dev/null @@ -1,221 +0,0 @@ -#-------------------# -#-Working Directory-# -#-------------------# - -setwd("./DataFormat") - -#----------------# -#-Load Libraries-# -#----------------# - -library(readxl) -library(dplyr) -library(tidyr) -library(reshape2) -library(nimble) -library(compareMCMCs) -library(coda) - -#-----------# -#-Load Data-# -#-----------# - -Data <- NULL -files <- list.files(path = "./EditedData", pattern = "Data", full.names = TRUE) -for(i in 1:length(files)){ - sheets <- excel_sheets(files[i]) - for(j in 1:length(sheets)){ - Temp <- read_excel(files[i], sheet = sheets[j]) - Temp <- Temp %>% mutate(species = as.character(sheets[j]), park = as.factor(i)) %>% - select(park, `deployment ID`, Year, species, ls(Temp, pattern = "R"), density) - Data <- bind_rows(Data, Temp) - } -} - -covfile <- "./RawData/camera trap metadata.xlsx" -covsheets <- excel_sheets(covfile) -covsheets <- covsheets[-3] #Remove Nyungwe - -CovData <- NULL -for(i in 1:5){ - Temp <- read_excel(covfile, sheet = covsheets[i]) - Temp <- Temp %>% select(`deployment ID`, days, elevation, edge, Year) - CovData <- bind_rows(CovData, Temp) -} - -CovData <- bind_rows(CovData, read_excel("./EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 1) %>% - select(`deployment ID`, days, elevation, edge, Year)) %>% - bind_rows(., read_excel("./EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 2) %>% - select(`deployment ID`, days, elevation, edge, Year)) %>% - distinct() - -Data <- full_join(Data, CovData, by = c("deployment ID", "Year")) %>% drop_na() - -#-------------# -#-Format Data-# -#-------------# - -#Convert to NAs -Data[,5:10] <- apply(Data[,5:10], 2, function(y) as.numeric(gsub("-", NA, y))) - -#Remove cameras that were sampled in only 1 year -Data <- Data %>% group_by(`deployment ID`) %>% filter(n() > 2) %>% ungroup() - -#Filter for single-species -Data <- Data %>% filter(species == "nigrifrons") - -#Convert to single vector format -Data2 <- melt(Data, id = c("park", "deployment ID", "Year", "species", "days", "elevation", "edge")) - -#Occupancy data -y <- Data2$value - -#Indices -nyrs <- Data %>% group_by(`deployment ID`) %>% summarize(min(Year), max(Year)) -nstart <- as.numeric(nyrs$`min(Year)`) - as.numeric(min(nyrs$`min(Year)`)) + 1 -nend <- 1 + as.numeric(max(nyrs$`max(Year)`)) - as.numeric(min(nyrs$`min(Year)`)) - - (as.numeric(max(nyrs$`max(Year)`)) - as.numeric(nyrs$`max(Year)`)) - -nyrs <- as.numeric(max(nyrs$`max(Year)`)) - as.numeric(min(nyrs$`min(Year)`)) + 1 - -nsites <- as.numeric(Data %>% summarize(n_distinct(`deployment ID`))) -nobs <- length(y) - -#Nested indices -yr <- as.numeric(as.factor(Data2$Year)) -site <- as.numeric(as.factor(Data2$`deployment ID`)) - -#Covariates -Cov <- Data %>% group_by(`deployment ID`, Year) %>% - distinct(`deployment ID`, Year, .keep_all = TRUE) -Cov$elevation <- scale(as.numeric(Cov$elevation)) -Cov$edge <- scale(as.numeric(Cov$edge)) -Cov$days <- scale(as.numeric(Cov$days)) - -elev <- dcast(Cov %>% select(`deployment ID`, Year, elevation), `deployment ID` ~ Year) -elev <- t(as.matrix(elev[,-1])) -edge <- dcast(Cov %>% select(`deployment ID`, Year, edge), `deployment ID` ~ Year) -edge <- t(as.matrix(edge[,-1])) -days <- dcast(Cov %>% select(`deployment ID`, Year, days), `deployment ID` ~ Year) -days <- t(as.matrix(days[,-1])) -density <- Data %>% distinct(`deployment ID`, density) -density <- as.numeric(density$density) - -#--------------# -#-NIMBLE model-# -#--------------# - -SSdyn.c <- nimbleCode({ - - beta0 ~ dnorm(0, 0.1) - beta1 ~ dnorm(0, 0.1) - beta2 ~ dnorm(0, 0.1) - beta3 ~ dnorm(0, 0.1) - beta4 ~ dnorm(0, 0.1) - - alpha0 ~ dunif(0, 1) - alpha1 ~ dnorm(0, 0.1) - - omega0 ~ dunif(0, 1) - omega1 ~ dnorm(0, 0.1) - omega2 ~ dnorm(0, 0.1) - omega3 ~ dnorm(0, 0.1) - - gamma0 ~ dnorm(0, 0.1) - - for(j in 1:nsites){ - - log(lambda[j]) <- beta0 + beta1 * density[j] + beta2 * edge[nstart[j],j] + - beta3 * density[j] * edge[nstart[j],j] + beta4 * elev[nstart[j],j] - - logit(r[nstart[j],j]) <- logit(alpha0) + alpha1 * days[nstart[j],j] - N[nstart[j],j] ~ dpois(lambda[j]) - - for(t in (nstart[j]+1):nend[j]){ - - logit(r[t,j]) <- alpha0 + alpha1 * days[t,j] - - logit(omega[t-1,j]) <- logit(omega0) + omega1 * density[j] + omega2 * edge[t,j] + - omega3 * density[j] * edge[t,j] - - S[t-1,j] ~ dbin(omega[t-1,j], N[t-1,j]) - G[t-1,j] ~ dpois(gamma[t-1]) - N[t,j] <- S[t-1,j] + G[t-1,j] - - }#end t - }#end j - - for(t in 2:nyrs){ #check 9 - log(gamma[t-1]) <- gamma0 - }#end t - - for(k in 1:nobs){ - y[k] ~ dbern(p[k]) - p[k] <- 1 - pow((1 - r[yr[k], site[k]]), N[yr[k], site[k]]) - }#end k - -}) - -#--------------# -#-Compile Data-# -#--------------# - -#Data -SSdyn.d <- list(y = y, nobs = nobs, nstart = nstart, nend = nend, nsites = nsites, nyrs = nyrs, - yr = yr, site = site, elev = elev, edge = edge, density = density, days = days) - -#Initial values -Nst <- array(NA, dim = c(nyrs, nsites)) -Sst <- array(NA, dim = c(nyrs-1, nsites)) -Gst <- array(NA, dim = c(nyrs-1, nsites)) -for(j in 1:nsites){ - Nst[nstart[j],j] <- 10 - Sst[nstart[j]:(nend[j]-1),j] <- 5 - Gst[nstart[j]:(nend[j]-1),j] <- 2 -} - -inits <- function()list(N=Nst, S=Sst, G=Gst) - -#Parameters to save -params <- c("alpha0", "alpha1", "beta0", "beta1", "beta2", "beta3", "beta4", - "omega0", "omega1", "omega2", "omega3", "gamma0") - -#MCMC settings -SSdyn.m <- nimbleModel(SSdyn.c, constants = SSdyn.d, inits = inits()) - -MCMCconf <- configureMCMC(SSdyn.m) -MCMCconf$printSamplers(1:12) - -MCMCconf$removeSampler(c('beta0', 'beta1', 'beta2', 'beta3', 'beta4', - 'omega0', 'omega1', 'omega2', 'omega3')) - -MCMCconf$addSampler(target = c('beta0', 'beta1', 'beta2', 'beta3', 'beta4'), - type = "AF_slice") - -MCMCconf$addSampler(target = c('omega0', 'omega1', 'omega2', 'omega3'), - type = "AF_slice") - -MCMCconf$samplerExecutionOrder <- c(seq(3,4166,1),1,2) - -MCMC <- buildMCMC(MCMCconf) - -SSdyn.comp <- compileNimble(SSdyn.m, MCMC) - -ni <- 25000 -nb <- 20000 -nc <- 3 - -#Run NIMBLE model -SSdyn.o <- runMCMC(SSdyn.comp$MCMC, niter = ni, nburnin = nb, nchains = nc) -SSdyn.o2 <- nimbleMCMC(SSdyn.c, inits = inits, constants = SSdyn.d, - niter = ni, nburnin = nb) - -#------------# -#-Post check-# -#------------# - -samples <- as.mcmc.list(lapply(SSdyn.o, mcmc)) - -gelman.diag(samples) -effectiveSize(samples) -plot(samples) diff --git a/DataAnalysis/SimulationStudy/Simulation.R b/DataAnalysis/SimulationStudy/Simulation.R new file mode 100644 index 0000000..47778ce --- /dev/null +++ b/DataAnalysis/SimulationStudy/Simulation.R @@ -0,0 +1,406 @@ +#-----------# +#-Libraries-# +#-----------# + +library(nimble) +library(coda) + +#-----------# +#-Functions-# +#-----------# + +#Logit function +logit <- function(pp) +{ + log(pp) - log(1-pp) +} + +#Inverse logit +expit <- function(eta) +{ + 1/(1+exp(-eta)) +} + +#Species ID +spp.fun <- function(x,y){ + out <- c( + which.max(apply(y, MARGIN = 1, sum)), + which.min(apply(y, MARGIN = 1, sum)), + which.min(apply(x, MARGIN = 1, sum))) + return(out) +} + + +# "True" values +nYears <- 10 +nReps <- 3 +nSites <- 75 +nSpecies <- 30 + +mu.lambda <- runif(1,0.1,1.5) +mu.lambdaL <- log(mu.lambda) +mu.omega <- runif(1,0,1) +mu.omegaL <- logit(mu.omega) +mu.gamma <- runif(1,0,1) +mu.gammaL <- log(mu.gamma) +mu.r <- runif(1,0,1) +mu.rL <- logit(mu.r) + +sd.lambda <- 0.5 +sd.omega <- 0.5 +sd.gamma <- 0.5 +sd.r <- 0.5 + +lambda0 <- rnorm(nSpecies, mu.lambdaL, sd.lambda) +omega0 <- rnorm(nSpecies, mu.omegaL, sd.omega) +gamma0 <- rnorm(nSpecies, mu.gammaL, sd.gamma) +r0 <- rnorm(nSpecies, mu.rL, sd.r) + +lambda <- exp(lambda0) +# lambda <- lambda[order(lambda)] +omega <- expit(omega0) +# omega <- omega[order(omega)] +gamma <- exp(gamma0) +# gamma <- gamma[order(gamma)] +r <- expit(r0) +# r <- r[order(r)] + +#Simulate true abundances, N, for each location +N <- array(NA, dim = c(nSpecies, nSites, nYears)) +S <- G <- array(NA, dim = c(nSpecies, nSites, nYears-1)) + +#subsequent years follow the birth-death-immigration process +for(i in 1:nSpecies) { + N[i,,1] <- rpois(nSites, lambda[i]) + for(t in 2:nYears) { + S[i,,t-1] <- rbinom(nSites, N[i,,t-1], omega[i]) + G[i,,t-1] <- rpois(nSites, gamma[i]) + N[i,,t] <- S[i,,t-1] + G[i,,t-1] + } +} + +# Generate data vector y for the counts +y <- array(NA, c(nSpecies, nSites, nYears, nReps)) +for(i in 1:nSpecies){ + for(t in 1:nYears) { + for(k in 1:nReps) { + y[i,,t,k] <- rbinom(nSites, N[i,,t], r[i]) + } + } +} + +# Convert to occupancy data +y[which(y != 0)] <- 1 + + +# Common, rare, elusive +spp.ID <- spp.fun(N,y) + +#-------------------------# +#-Multispecies model code-# +#-------------------------# + +code <- nimbleCode({ + + mu.lambdaL ~ dnorm(0, 0.1) + log(mu.lambda) <- mu.lambdaL + tau.lambda ~ dgamma(0.1, 0.1) + sd.lambda <- 1/sqrt(tau.lambda) + mu.omegaL <- logit(mu.omega) + mu.omega ~ dunif(0, 1) + tau.omega ~ dgamma(0.1, 0.1) + sd.omega <- 1/sqrt(tau.omega) + mu.gammaL ~ dnorm(0, 0.1) + log(mu.gamma) <- mu.gammaL + tau.gamma ~ dgamma(0.1, 0.1) + sd.gamma <- 1/sqrt(tau.gamma) + mu.rL <- logit(mu.r) + mu.r ~ dunif(0, 1) + tau.r ~ dgamma(0.1, 0.1) + sd.r <- 1/sqrt(tau.r) + + for(i in 1:nSpecies){ + lambda0[i] ~ dnorm(mu.lambdaL, tau.lambda) + omega0[i] ~ dnorm(mu.omegaL, tau.omega) + gamma0[i] ~ dnorm(mu.gammaL, tau.gamma) + r0[i] ~ dnorm(mu.r, tau.r) + + log(lambda[i]) <- lambda0[i] + logit(omega[i]) <- omega0[i] + log(gamma[i]) <- gamma0[i] + logit(r[i]) <- r0[i] + + for(j in 1:nSites){ + + N[i,j,1] ~ dpois(lambda[i]) + + for(t in 2:nYears){ + + S[i,j,t-1] ~ dbin(omega[i], N[i,j,t-1]) + G[i,j,t-1] ~ dpois(gamma[i]) + N[i,j,t] <- S[i,j,t-1] + G[i,j,t-1] + + }#end t + + for(t in 1:nYears){ + + p[i,j,t] <- 1 - pow((1 - r[i]), N[i,j,t]) + + for(k in 1:nReps){ + + y[i,j,t,k] ~ dbern(p[i,j,t]) + + }#end k + }#end t + }#end j + }#end i +}) + +#--------------# +#-Compile data-# +#--------------# + +data <- list(y = y) + +constants <- list(nSpecies = nSpecies, nSites = nSites, nYears = nYears, nReps = nReps) + +#----------------# +#-Initial values-# +#----------------# + +inits <- function()list(N = N, S = S, G = G, + mu.lambdaL = mu.lambdaL, + # tau.lambda, + mu.omega = mu.omega, + # tau.omega, + mu.gammaL = mu.gammaL, + # tau.gamma, + mu.r = mu.r, + # tau.r, + lambda0 = lambda0, + omega0 = omega0, + gamma0 = gamma0, + r0 = r0 + ) + +#--------------------# +#-Parameters to save-# +#--------------------# + +params <- c("mu.lambda", "sd.lambda", + "mu.omega", "sd.omega", + "mu.gamma", "sd.gamma", + "mu.r", "sd.r", + "lambda", "omega", "gamma", "r") + +#---------------# +#-MCMC settings-# +#---------------# + +model <- nimbleModel(code = code, + constants = constants, + data = data, + inits = inits()) + +MCMCconf <- configureMCMC(model, monitors = params) + +MCMC <- buildMCMC(MCMCconf) + +compiled.model <- compileNimble(model, MCMC) + +ni <- 35000 +nb <- 10000 +nc <- 3 +nt <- 25 + +# ni <- 350 +# nb <- 100 +# nc <- 3 +# nt <- 1 + +#-----------# +#-Run model-# +#-----------# + +out <- runMCMC(compiled.model$MCMC, + niter = ni, nburnin = nb, + nchains = nc, thin = nt, + samplesAsCodaMCMC = TRUE) + +#----------------# +#-Compile output-# +#----------------# + +coefs <- c(outer(paste0(c("gamma", "lambda", "omega", "r"), "["), paste0(spp.ID, "]"), paste0), + "mu.gamma", "mu.lambda", "mu.omega", "mu.r", "sd.gamma", "sd.lambda", "sd.omega", "sd.r") + +coefs.names <- c(outer(c("gamma.", "lambda.", "omega.", "r."), c("common", "elusive", "rare"), paste0), + "mu.gamma", "mu.lambda", "mu.omega", "mu.r", "sd.gamma", "sd.lambda", "sd.omega", "sd.r") + +output <- cbind(sapply(coefs, FUN = function(x){eval(parse(text = x))}), + summary(out)$statistics[coefs,"Mean"], + summary(out)$statistics[coefs,"SD"], + gelman.diag(out)$psrf[coefs,1]) + +# output <- cbind(c(gamma, lambda, mu.gamma, mu.lambda, +# mu.omega, mu.r, omega, r, +# sd.gamma, sd.lambda, sd.omega, sd.r), +# summary(out)$statistics[,"Mean"], +# gelman.diag(out)$psrf[,1]) + +output <- as.data.frame(output) +# output$Params <- rownames(output) +output$Params <- coefs.names +rownames(output) <- NULL +colnames(output)[1:4] <- c("Truth", "Mean_MS", "SD_MS", "Rhat_MS") +output <- output[,c(5,1,2,3,4)] +output$Rhat_SS <- output$SD_SS <- output$Mean_SS <- NA + +#---------------------------# +#-Single-species model code-# +#---------------------------# + +code <- nimbleCode({ + + lambda0 ~ dnorm(0, 0.01) + omega ~ dunif(0, 1) + gamma0 ~ dnorm(0, 0.01) + r ~ dunif(0, 1) + + log(lambda) <- lambda0 + log(gamma) <- gamma0 + + for(j in 1:nSites){ + + N[j,1] ~ dpois(lambda) + + for(t in 2:nYears){ + + S[j,t-1] ~ dbin(omega, N[j,t-1]) + G[j,t-1] ~ dpois(gamma) + N[j,t] <- S[j,t-1] + G[j,t-1] + + }#end t + + for(t in 1:nYears){ + + p[j,t] <- 1 - pow((1 - r), N[j,t]) + + for(k in 1:nReps){ + + y[j,t,k] ~ dbern(p[j,t]) + + }#end k + }#end t + }#end j +}) + +#--------------# +#-Compile data-# +#--------------# + +iter <- 1 + +# for(i in 1:nSpecies){ +for(i in unique(spp.ID)){ + +data <- list(y = y[i,,,]) + +constants <- list(nSites = nSites, nYears = nYears, nReps = nReps) + +#----------------# +#-Initial values-# +#----------------# + +inits <- function()list(N = N[i,,], S = S[i,,], G = G[i,,], + lambda0 = lambda0[i], + omega = omega[i], + gamma0 = gamma0[i], + r = r[i] +) + +#--------------------# +#-Parameters to save-# +#--------------------# + +params <- c("gamma", "lambda", "omega", "r") + +#---------------# +#-MCMC settings-# +#---------------# + +model <- nimbleModel(code = code, + constants = constants, + data = data, + inits = inits()) + +MCMCconf <- configureMCMC(model, monitors = params) + +MCMC <- buildMCMC(MCMCconf) + +compiled.model <- compileNimble(model, MCMC) + +#-----------# +#-Run model-# +#-----------# + +out <- runMCMC(compiled.model$MCMC, + niter = ni, nburnin = nb, + nchains = nc, thin = nt, + samplesAsCodaMCMC = TRUE) + +#----------------# +#-Compile output-# +#----------------# + + +if(iter == 1){ + + output[1:4,6] <- summary(out)$statistics[,"Mean"] + output[1:4,7] <- summary(out)$statistics[,"SD"] + output[1:4,8] <- gelman.diag(out)$psrf[,1] + +}else{ + + if(length(unique(spp.ID) == 2)){ + + output[9:12,6] <- output[5:8,6] <- summary(out)$statistics[,"Mean"] + output[9:12,7] <- output[5:8,7] <- summary(out)$statistics[,"SD"] + output[9:12,8] <- output[5:8,8] <- gelman.diag(out)$psrf[,1] + + }else{ + + if(iter == 2){ + + output[5:8,6] <- summary(out)$statistics[,"Mean"] + output[5:8,7] <- summary(out)$statistics[,"SD"] + output[5:8,8] <- gelman.diag(out)$psrf[,1] + + }else{ + + output[9:12,6] <- summary(out)$statistics[,"Mean"] + output[9:12,7] <- summary(out)$statistics[,"SD"] + output[9:12,8] <- gelman.diag(out)$psrf[,1] + + } + } +} + +# output[c(i,i+nSpecies,i+4+nSpecies*2,i+4+nSpecies*3),5] <- summary(out)$statistics[,"Mean"] +# +# output[c(i,i+nSpecies,i+4+nSpecies*2,i+4+nSpecies*3),6] <- gelman.diag(out)$psrf[,1] +# +# print(i) + +iter <- iter + 1 + +print(iter) + +}#end i + +#-------------# +#-Save output-# +#-------------# + +ID <- length(list.files("./output_75_final/")) + 1 +save(output, file = paste("./output_75_final/output", ID, ".Rds", sep="")) diff --git a/DataFormat/DataFormat.R b/DataFormat/DataFormat.R index c67d7ff..5ba2771 100644 --- a/DataFormat/DataFormat.R +++ b/DataFormat/DataFormat.R @@ -1,8 +1,3 @@ -#TODO: - -#remove NAs by filtering for day > 1 -#nstart:nend would need to become ragged array skipping years when sampling doesn't occur at a site - #----------------# #-Load Libraries-# #----------------# @@ -10,13 +5,14 @@ library(readxl) library(tidyverse) library(reshape2) +library(sf) #-----------# #-Load Data-# #-----------# Data <- NULL -files <- list.files(path = "~/HMSNO/DataFormat/EditedData", pattern = "Data", full.names = TRUE) +files <- list.files(path = "~/HMSNO/DataFormat/RawData", pattern = "Data", full.names = TRUE) for(i in 1:length(files)){ sheets <- excel_sheets(files[i]) for(j in 1:length(sheets)){ @@ -27,6 +23,53 @@ for(i in 1:length(files)){ } } +#UDZ ungulates 2018-2019 +sheets <- excel_sheets("~/HMSNO/DataFormat/RawData/UDZ_ungulates_2018_2019.xlsx") +UDZdf <- NULL +for(j in 3:6){ + Temp <- read_excel("~/HMSNO/DataFormat/RawData/UDZ_ungulates_2018_2019.xlsx", + sheet = sheets[j]) + UDZdf <- bind_rows(UDZdf, Temp) +} + +Data <- bind_rows(Data, UDZdf %>% select(Trap...2, Year, Species, ls(Temp, pattern = "R")) %>% select(-R7) %>% + rename(`deployment ID` = Trap...2, species = Species) %>% + mutate(species = recode(species, "Cephalophus harveyi" = "harveyi", "Cephalophus spadix" = "spadix", + "Nesotragus moschatus" = "moschatus", "Tragelaphus scriptus" = "scriptus"), + park = as.factor(5), density = 0) %>% + relocate(park, `deployment ID`)) + +#Remaining ungulates +sheets <- excel_sheets("~/HMSNO/DataFormat/RawData/Remaining_ungulates.xlsx") +RUdf <- NULL +for(k in 1:length(sheets)){ + Temp <- read_excel("~/HMSNO/DataFormat/RawData/Remaining_ungulates.xlsx", + sheet = sheets[k]) + RUdf <- bind_rows(RUdf, Temp) +} + +Data <- bind_rows(Data, RUdf %>% select(park, `deployment ID`, Year, Species, ls(Temp, pattern = "R"), density) %>% + rename(species = Species) %>% + mutate(species = recode(species, "Tragelaphus scriptus" = "scriptus", + "Tragelaphus spekii" = "spekii"), + park = as.factor(park))) + +#Add missing UDZ data +sheets <- excel_sheets("~/HMSNO/DataFormat/RawData/Bushbuck_suni_UDZ.xlsx") +UDZdf <- NULL +for(k in 3:4){ + Temp <- read_excel("~/HMSNO/DataFormat/RawData/Bushbuck_suni_UDZ.xlsx", + sheet = sheets[k]) + Temp$park <- 5 + UDZdf <- bind_rows(UDZdf, Temp) +} + +Data <- bind_rows(Data, UDZdf %>% select(park, `deployment ID`, Year, species, ls(Temp, pattern = "R"), density) %>% + mutate(park = as.factor(park))) + + + +#Site data SiteData <- NULL file <- list.files(path = "~/HMSNO/DataFormat/RawData", pattern = "metadata", full.names = TRUE) sheets <- excel_sheets(file) @@ -37,189 +80,274 @@ for(j in 1:length(sheets)){ SiteData <- bind_rows(SiteData, Temp) } -SiteData <- bind_rows(SiteData, read_excel("~/HMSNO/DataFormat/EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 1) %>% +SiteData <- bind_rows(SiteData, read_excel("~/HMSNO/DataFormat/RawData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 1) %>% select(`deployment ID`, days, elevation, edge, Year)) %>% - bind_rows(., read_excel("~/HMSNO/DataFormat/EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 2) %>% + bind_rows(., read_excel("~/HMSNO/DataFormat/RawData/Data_Nyungwe_nig_sylv_April_1.xlsx", sheet = 2) %>% select(`deployment ID`, days, elevation, edge, Year)) %>% distinct() %>% select(`deployment ID`, days, elevation, edge, Year) +SiteData <- bind_rows(SiteData, SiteData %>% filter(grepl("CT-UDZ", `deployment ID`) & Year >= 2016) %>% + mutate(Year = ifelse(Year == 2016, 2018, 2019), + days = read_excel("~/HMSNO/DataFormat/RawData/UDZ_ungulates_2018_2019.xlsx", sheet = 4) %>% + select(trapdays) %>% .$trapdays)) + +#Add GPS data +sites <- NULL +sheets <- excel_sheets("~/HMSNO/DataFormat/RawData/trap locations and dates.xlsx") +for(i in 1:length(sheets)){ + Temp <- read_excel(path = "~/HMSNO/DataFormat/RawData/trap locations and dates.xlsx", sheet = sheets[i]) + Temp <- Temp %>% mutate(latitude = as.numeric(as.character(latitude)), + longitude = as.numeric(as.character(longitude)), + start_date = as.Date(start_date), + end_date = as.Date(end_date)) + sites <- bind_rows(sites, Temp) +} + +sites$deployment <- str_remove(sites$deployment, "CT-NNP-1-") +sites$deployment <- str_remove(sites$deployment, "NNP-2015-") +tmp <- sites$deployment[grep("TR", sites$deployment)] +tmp[-grep("-", tmp)] <- gsub("(....)(.*)","\\1-\\2",tmp[-grep("-", tmp)]) +sites$deployment[grep("TR", sites$deployment)] <- tmp +tmp <- sites$deployment[grep("TR-", sites$deployment)] +tmp <- stringr::str_replace(tmp, "-", "") +sites$deployment[grep("TR-", sites$deployment)] <- tmp +sites <- sites %>% rename(`deployment ID` = deployment, Year = year) +SiteData <- left_join(SiteData, sites, by = c("deployment ID", "Year")) + +#Set coordinate system +SiteData <- SiteData %>% group_by(`deployment ID`) %>% + fill(longitude, latitude, start_date, end_date, .direction = "downup") + +SiteData <- SiteData %>% + drop_na(longitude|latitude) %>% + st_as_sf(coords = c("longitude", "latitude")) %>% + st_set_crs(4326) + #-------------# #-Format Data-# #-------------# +#Turn dashes to NAs Data[,5:10] <- apply(Data[,5:10], 2, function(y) as.numeric(gsub("-", NA, y))) -park.levels <- c("1" = "2", "2" = "3", "3" = "4", "4" = "1", "5" = "6", "6" = "5") +#Convert CT-UDZ-1-14.1 to CT-UDZ-1-14 +Data$`deployment ID`[Data$`deployment ID` == "CT-UDZ-1-14.1"] <- "CT-UDZ-1-14" + +#Remove CT-UDZ-3-21 +Data <- Data %>% filter(`deployment ID` != "CT-UDZ-3-21") + +#Remove 2017 from VIR +Data <- Data %>% filter(!(park == 6 & Year == 2017)) + +#Reorder park levels +park.levels <- c("1" = "5", "2" = "6", "3" = "4", "4" = "1", "5" = "3", "6" = "2") Data <- Data %>% mutate(`deployment ID` = ifelse(species == "sylvilocutor" & Year == "2010" & park == 1, gsub("TEAM-001", "CT", `deployment ID`), `deployment ID`), `deployment ID` = ifelse(species == "sylvilocutor" & Year == "2010" & park == 1, gsub("-2010", "", `deployment ID`), `deployment ID`), - park = fct_recode(factor(park, levels = c("2","3","4","1","6","5")), !!!park.levels)) %>% + park = fct_recode(factor(park, levels = c("5","6","4","1","3","2")), !!!park.levels)) %>% arrange(park, Year, `deployment ID`) Data <- left_join(Data, SiteData, by = c("deployment ID", "Year")) +#Control for days sampled +Data$days[Data$days < 0] <- 0 #Set negative days to 0 +Data$days[Data$days > 30] <- 30 #Truncate days at max 30 days +Data$days[Data$days == 0] <- NA #Set 0 days to NA (ie, no sampling) +Data[is.na(Data$days),5:10] <- NA #Set all occ data to NA for sites with NA days (ie, no sampling) +Data <- Data %>% mutate(R6 = replace(R6, which(days>=25 & is.na(R6)), 0), #Replace values with NAs or 0s + R6 = replace(R6, which(days<25), NA), + R5 = replace(R5, which(days>=20 & is.na(R5)), 0), + R5 = replace(R5, which(days<20), NA), + R4 = replace(R4, which(days>=15 & is.na(R4)), 0), + R4 = replace(R4, which(days<15), NA), + R3 = replace(R3, which(days>=10 & is.na(R3)), 0), + R3 = replace(R3, which(days<10), NA), + R2 = replace(R2, which(days>=5 & is.na(R2)), 0), + R2 = replace(R2, which(days<5), NA), + R1 = replace(R1, which(days>0 & is.na(R1)), 0)) + +Data$days[is.na(Data$days)] <- 0 #Reset no sampling to zero days + Data <- Data %>% drop_na(days) %>% - mutate(siteID = as.numeric(factor(`deployment ID`, levels = unique(`deployment ID`))), - parkID = as.numeric(park), + mutate(parkID = as.numeric(park), yearID = as.numeric(as.factor(Year)), specID = as.numeric(as.factor(species))) %>% - arrange(parkID, yearID, siteID) - -Data.vec <- Data %>% - select(specID, parkID, siteID, yearID, ls(Data, pattern = "R")) %>% - melt(id = c("specID", "parkID", "siteID", "yearID")) - -#--------------# -#-Extract data-# -#--------------# - -# spec <- Data$specID -# site <- Data$siteID -# yr <- Data$yearID + group_by(parkID) %>% + mutate(siteID = as.numeric(factor(`deployment ID`, levels = unique(`deployment ID`)))) %>% + ungroup(parkID) %>% + arrange(parkID, specID, siteID, yearID) + +#---------# +#-Indices-# +#---------# + +#Number of species +nspecs <- max(Data$specID) + +#Number of parks +nparks <- max(Data$parkID) + +#First park for each park +parkS <- as.numeric(Data %>% group_by(specID) %>% + summarize(parkS = min(parkID) - min(Data$parkID) + 1) %>% + select(parkS) %>% .$parkS) + +#Last park for each park +parkE <- as.numeric(Data %>% group_by(specID) %>% + summarize(parkE = max(parkID) - min(Data$parkID) + 1) %>% + select(parkE) %>% .$parkE) +#Max number of sites +nsites <- max(Data$siteID) -spec <- Data.vec$specID -site <- Data.vec$siteID -yr <- Data.vec$yearID -park <- as.numeric(Data %>% group_by(siteID) %>% - distinct(parkID) %>% select(parkID) %>% .$parkID) +#Number of sites/park +nsite <- as.numeric(Data %>% group_by(parkID) %>% + summarize(nsite = n_distinct(siteID)) %>% + select(nsite) %>% .$nsite) -nspec <- max(Data$specID) -nparks <- max(park) -nsites <- max(Data$siteID) +#Number of years nyrs <- max(Data$yearID) -nreps <- 6 -nobs <- dim(Data.vec)[1] - -npark <- as.numeric(Data %>% group_by(specID) %>% - summarize(npark = n_distinct(parkID)) %>% - select(npark) %>% .$npark) -nsite <- as.numeric(Data %>% group_by(specID) %>% - summarize(nsite = n_distinct(siteID)) %>% - select(nsite) %>% .$nsite) - -parkspec <- matrix(NA, nrow = nparks, ncol = nspec) -sitespec <- matrix(NA, nrow = nsites, ncol = nspec) - -tmp1 <- as.matrix(Data %>% group_by(specID) %>% - distinct(parkID) %>% arrange(specID, parkID)) -tmp2 <- as.matrix(Data %>% group_by(specID) %>% - distinct(siteID) %>% arrange(specID, siteID)) -k <- 1 -v <- 1 -for(s in 1:nspec){ - for(i in 1:npark[s]){ - parkspec[i,s] <- tmp1[k,1] - k <- k + 1 - } - for(j in 1:nsite[s]){ - sitespec[j,s] <- tmp2[v,1] - v <- v + 1 - } -} -rm(tmp1, tmp2, k, v) - -nstart <- as.numeric(Data %>% group_by(parkID, siteID) %>% - summarize(nstart = min(yearID) - min(Data$yearID) + 1) %>% - select(nstart) %>% .$nstart) -nend <- as.numeric(Data %>% group_by(parkID, siteID) %>% - summarize(nend = max(yearID) - min(Data$yearID) + 1) %>% - select(nend) %>% .$nend) -nyrstr <- as.numeric(Data %>% group_by(parkID) %>% - summarize(nyrstr = min(yearID) - min(Data$yearID) + 1) %>% - select(nyrstr) %>% .$nyrstr) -nyrend <- as.numeric(Data %>% group_by(parkID) %>% - summarize(nyrend = max(yearID) - min(Data$yearID) + 1) %>% - select(nyrend) %>% .$nyrend) -nsitestr <- as.numeric(Data %>% group_by(parkID) %>% - summarize(nsitestr = min(siteID)) %>% - select(nsitestr) %>% .$nsitestr) -nsiteend <- matrix(data = NA, nrow = nyrs, ncol = nparks) -for(t in 1:nyrs){ - nsiteend[t,] <- as.numeric(Data %>% group_by(parkID) %>% - summarize(nsiteend = max(siteID)) %>% - select(nsiteend) %>% .$nsiteend) -} +#First year for each park +nstart <- as.numeric(Data %>% group_by(parkID) %>% + summarize(nstart = min(yearID) - min(Data$yearID) + 1) %>% + select(nstart) %>% .$nstart) -for(t in nyrstr[3]:nyrend[3]){ - nsiteend[t,3] <- as.numeric(Data %>% filter(parkID == 3 & yearID == t)%>% - group_by(parkID) %>% - summarize(nsiteend = max(siteID)) %>% - select(nsiteend) %>% .$nsiteend) -} +#Last year for each park +nend <- as.numeric(Data %>% group_by(parkID) %>% + summarize(nend = max(yearID) - min(Data$yearID) + 1) %>% + select(nend) %>% .$nend) -#SKIP specific parks for a given species... - -#SKIP specific years or reps for a given site... -# as.numeric(Data %>% filter(days > 1) %>% -# group_by(park, `deployment ID`) %>% -# summarize(nstart = min(Year) - min(Data$Year) + 1) %>% -# select(nstart) %>% .$nstart) +#Number of replicates +nreps <- 6 -occ <- Data.vec$value +#Nested indices +yr <- Data$yearID +site <- Data$siteID +park <- Data$parkID +spec <- Data$specID -# occ <- array(NA, dim = c(nreps, nyrs, nsites, nspec)) -# -# for(s in 1:nspec){ -# for(j in 1:nsites){ -# #for(t in ifelse(nstart[j] == nend[j], nstart[j], nstart[j]:nend[j])){ -# for(t in nstart[j]:nend[j]){ -# occ[1:nreps,t,j,s] <- rep(0,nreps) -# } -# } -# } -# for(i in 1:dim(Data)[1]){ -# for(k in 5:(nreps+5)){ -# if(is.na(Data[i,k])){ -# occ[k-4,yr[i],site[i],1:nspec] <- rep(NA, nspec) -# } -# } -# occ[1:6,yr[i],site[i],spec[i]] <- as.numeric(Data[i,5:10]) -# } +#--------------# +#-Extract data-# +#--------------# +#Occupancy +y <- array(NA, dim = c(nreps, nyrs, nsites, nparks, nspecs)) +for(i in 1:dim(Data)[1]){ + y[1:6,yr[i],site[i],park[i],spec[i]] <- as.numeric(Data[i,5:10]) +} -Covariatedf <- Data %>% group_by(siteID, yearID) %>% - select(days, edge, elevation) %>% distinct(siteID, yearID, .keep_all = TRUE) +#Covariates +Cov <- Data %>% group_by(parkID, siteID, yearID) %>% + select(days, edge, elevation, geometry, start_date, end_date) %>% + distinct(parkID, siteID, yearID, .keep_all = TRUE) -days <- edge <- elevation <- array(NA, dim = c(nyrs, nsites)) -for(i in 1:dim(Covariatedf)[1]){ - days[Covariatedf$yearID[i], Covariatedf$siteID[i]] <- Covariatedf$days[i] - edge[Covariatedf$yearID[i], Covariatedf$siteID[i]] <- Covariatedf$edge[i] - elevation[Covariatedf$yearID[i], Covariatedf$siteID[i]] <- Covariatedf$elevation[i] +#Fix start_date years that do not match year ID +for(i in 1:dim(Cov)[1]){ + if(as.numeric(as.factor(format(Cov$start_date[i], "%Y"))) != Cov$yearID[i]){ + lubridate::year(Cov$start_date[i]) <- Cov$yearID[i] + 2008 + } } -density <- Data %>% group_by(siteID) %>% - select(density) %>% distinct() %>% .$density +Cov <- Cov %>% mutate(end_date_old = end_date, + end_date = start_date + 30) + + +# edge <- elevation <- array(NA, dim = c(nyrs, nsites, nparks)) +# days <- array(NA, dim = c(nreps, nyrs, nsites, nparks)) + +for(i in 1:dim(Cov)[1]){ + + if(Cov$days[i] > 25){ + tmp <- c(rep(5,5), (Cov$days[i] - 25)) + }else{ + if(Cov$days[i] <= 25 & Cov$days[i] > 20){ + tmp <- c(rep(5,4), (Cov$days[i] - 20), 0) + }else{ + if(Cov$days[i] <= 20 & Cov$days[i] > 15){ + tmp <- c(rep(5,3), (Cov$days[i] - 15), rep(0, 2)) + }else{ + if(Cov$days[i] <= 15 & Cov$days[i] > 10){ + tmp <- c(rep(5,2), (Cov$days[i] - 10), rep(0, 3)) + }else{ + if(Cov$days[i] <= 10 & Cov$days[i] > 5){ + tmp <- c(1, (Cov$days[i] - 5), rep(0, 4)) + }else{ + if(Cov$days[i] <= 5 & Cov$days[i] > 0){ + tmp <- c(Cov$days[i], rep(0, 5)) + }else{ + if(Cov$days[i] == 0 | is.na(Cov$days[i])){ + tmp <- rep(0, 6) + } + } + } + } + } + } + } + + days[ , Cov$yearID[i], Cov$siteID[i], Cov$parkID[i]] <- tmp + # edge[Cov$yearID[i], Cov$siteID[i], Cov$parkID[i]] <- Cov$edge[i] + # elevation[Cov$yearID[i], Cov$siteID[i], Cov$parkID[i]] <- Cov$elevation[i] +} -days <- (days - mean(days, na.rm = TRUE))/sd(days, na.rm = TRUE) -edge <- (edge - mean(edge, na.rm = TRUE))/sd(edge, na.rm = TRUE) -elevation <- (elevation - mean(elevation, na.rm = TRUE))/sd(elevation, na.rm = TRUE) +days.scaled <- (days - mean(days, na.rm = TRUE))/sd(days, na.rm = TRUE) +# edge.scaled <- (edge - mean(edge, na.rm = TRUE))/sd(edge, na.rm = TRUE) +# elevation.scaled <- (elevation - mean(elevation, na.rm = TRUE))/sd(elevation, na.rm = TRUE) +# density <- as.numeric(Data %>% group_by(park) %>% distinct(density) %>% select(density) %>% .$density) -Effort <- Data %>% - group_by(parkID, yearID) %>% - summarize(nsite = n_distinct(siteID), days = sum(days)) +#----------# +#-Rainfall-# +#----------# +# Geomdf <- st_as_sf(Cov) +# +# Geomdf <- Geomdf %>% mutate( +# annual_precip = 0, +# ) +# +# year <- as.character(rep(2009:2020, each = 12)) +# month <- as.character(rep(sprintf('%0.2d', 1:12), length(2009:2020))) +# dates <- paste0(year, "-", month) +# year_ID <- as.numeric(as.factor(year)) +# +# +# #Code to download CHIRPS data +# for(i in 1:length(year)){ +# url <- paste0("https://data.chc.ucsb.edu/products/CHIRPS-2.0/africa_monthly/tifs/chirps-v2.0.", year[i], ".", month[i], ".tif.gz") +# httr::GET(url = url, +# httr::write_disk(path = paste0(getwd(), "/DataFormat/CHIRPS/", year[i], month[i], ".tif.gz"), overwrite = TRUE)) +# } +# +# filenames <- list.files("~/HMSNO/DataFormat/CHIRPS/", pattern = "chirps", full.names = TRUE) +# for(i in 1:length(filenames)){ +# rainfall <- raster::raster(filenames[i]) +# raster::values(rainfall)[raster::values(rainfall) < 0] <- NA +# Geomdf$cell <- tabularaster::cellnumbers(rainfall, Geomdf)$cell_ +# Geomdf <- Geomdf %>% mutate(annual_precip = ifelse(yearID == year_ID[i], raster::extract(rainfall, cell) + annual_precip, annual_precip)) +# } +# +# precip <- array(NA, dim = c(22, nsites, nparks)) +# +# for(i in 1:dim(Geomdf)[1]){ +# precip[Geomdf$yearID[i], Geomdf$siteID[i], Geomdf$parkID[i]] <- Geomdf$annual_precip[i] +# } #--------------# #-Compile data-# #--------------# -HMSNO.data <- list(y = occ, nstart = nstart, nend = nend, - park = park, elevation = elevation, edge = edge, density = density, days = days) +HMSNO.data <- list(y = y, + days.scaled = days.scaled) -HMSNO.con <- list(nspec = nspec, nobs = nobs, nsitestr = nsitestr, nsiteend = nsiteend, - nyrstr = nyrstr, nyrend = nyrend, ns = nstart, ne = nend, nsite = nsite, npark = npark, - parkspec = parkspec, sitespec = sitespec, spec = spec, site = site, yr = yr) +HMSNO.con <- list(nspecs = nspecs, parkS = parkS, parkE = parkE, + nsite = nsite, nreps = nreps, nstart = nstart, nend = nend) save(HMSNO.data, file = "~/HMSNO/DataFormat/HMSNO.data.Rdata") -save(HMSNO.con, file = "~/HMSNO/DataFormat/HMSNO.con.Rdata") - -save(Effort, file = "~/HMSNO/DataFormat/Effort.Rdata") - +save(HMSNO.con, file = "~/HMSNO/DataFormat/HMSNO.con.Rdata") \ No newline at end of file diff --git a/DataFormat/EditedData/Data_Bwindi_nig_sylv_Nov18.xlsx b/DataFormat/EditedData/Data_Bwindi_nig_sylv_Nov18.xlsx deleted file mode 100644 index a9e34d9..0000000 Binary files a/DataFormat/EditedData/Data_Bwindi_nig_sylv_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/EditedData/Data_Korup_olg_syl_mon_Nov18.xlsx b/DataFormat/EditedData/Data_Korup_olg_syl_mon_Nov18.xlsx deleted file mode 100644 index 8fcf66b..0000000 Binary files a/DataFormat/EditedData/Data_Korup_olg_syl_mon_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/EditedData/Data_Nouabale_Ndoki_all_Nov18.xlsx b/DataFormat/EditedData/Data_Nouabale_Ndoki_all_Nov18.xlsx deleted file mode 100644 index 417e9a0..0000000 Binary files a/DataFormat/EditedData/Data_Nouabale_Ndoki_all_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx b/DataFormat/EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx deleted file mode 100644 index 9ff08ee..0000000 Binary files a/DataFormat/EditedData/Data_Nyungwe_nig_sylv_April_1.xlsx and /dev/null differ diff --git a/DataFormat/EditedData/Data_UDZ_harveyi_spadix_Nov18.xlsx b/DataFormat/EditedData/Data_UDZ_harveyi_spadix_Nov18.xlsx deleted file mode 100644 index 87dab86..0000000 Binary files a/DataFormat/EditedData/Data_UDZ_harveyi_spadix_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/EditedData/Data_Virunga_nigrifrons_Nov18.xlsx b/DataFormat/EditedData/Data_Virunga_nigrifrons_Nov18.xlsx deleted file mode 100644 index 51df0f8..0000000 Binary files a/DataFormat/EditedData/Data_Virunga_nigrifrons_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/EditedData/Old_Nyungwe_nig_sylv_Nov18.xlsx b/DataFormat/EditedData/Old_Nyungwe_nig_sylv_Nov18.xlsx deleted file mode 100644 index d78cd5f..0000000 Binary files a/DataFormat/EditedData/Old_Nyungwe_nig_sylv_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/Effort.Rdata b/DataFormat/Effort.Rdata deleted file mode 100644 index e250c0f..0000000 Binary files a/DataFormat/Effort.Rdata and /dev/null differ diff --git a/DataFormat/HMSNO.con.Rdata b/DataFormat/HMSNO.con.Rdata index 9f3c36f..bbb7bcd 100644 Binary files a/DataFormat/HMSNO.con.Rdata and b/DataFormat/HMSNO.con.Rdata differ diff --git a/DataFormat/HMSNO.data.Rdata b/DataFormat/HMSNO.data.Rdata index 3b2bd42..a670871 100644 Binary files a/DataFormat/HMSNO.data.Rdata and b/DataFormat/HMSNO.data.Rdata differ diff --git a/DataFormat/MSdny.Rdata b/DataFormat/MSdny.Rdata deleted file mode 100644 index 9d164c9..0000000 Binary files a/DataFormat/MSdny.Rdata and /dev/null differ diff --git a/DataFormat/RawData/Bushbuck_suni_UDZ.xlsx b/DataFormat/RawData/Bushbuck_suni_UDZ.xlsx new file mode 100644 index 0000000..f64df24 Binary files /dev/null and b/DataFormat/RawData/Bushbuck_suni_UDZ.xlsx differ diff --git a/DataFormat/RawData/Data_Bwindi_nig_sylv_Nov18.xlsx b/DataFormat/RawData/Data_Bwindi_nig_sylv_Nov18.xlsx index 8de29e6..a0bedb9 100644 Binary files a/DataFormat/RawData/Data_Bwindi_nig_sylv_Nov18.xlsx and b/DataFormat/RawData/Data_Bwindi_nig_sylv_Nov18.xlsx differ diff --git a/DataFormat/RawData/Data_Korup_olg_syl_mon_Nov18.xlsx b/DataFormat/RawData/Data_Korup_olg_syl_mon_Nov18.xlsx index 28be223..8fcf66b 100644 Binary files a/DataFormat/RawData/Data_Korup_olg_syl_mon_Nov18.xlsx and b/DataFormat/RawData/Data_Korup_olg_syl_mon_Nov18.xlsx differ diff --git a/DataFormat/RawData/Data_Nouabale_Ndoki_all_Nov18.xlsx b/DataFormat/RawData/Data_Nouabale_Ndoki_all_Nov18.xlsx deleted file mode 100644 index 45adc3a..0000000 Binary files a/DataFormat/RawData/Data_Nouabale_Ndoki_all_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/RawData/Data_Nouabale_Ndoki_all_Nov18_Aug10.xlsx b/DataFormat/RawData/Data_Nouabale_Ndoki_all_Nov18_Aug10.xlsx new file mode 100644 index 0000000..0585548 Binary files /dev/null and b/DataFormat/RawData/Data_Nouabale_Ndoki_all_Nov18_Aug10.xlsx differ diff --git a/DataFormat/RawData/Data_Nyungwe_nig_sylv_April_1.xlsx b/DataFormat/RawData/Data_Nyungwe_nig_sylv_April_1.xlsx index 9003e33..b33ad96 100644 Binary files a/DataFormat/RawData/Data_Nyungwe_nig_sylv_April_1.xlsx and b/DataFormat/RawData/Data_Nyungwe_nig_sylv_April_1.xlsx differ diff --git a/DataFormat/RawData/Data_Nyungwe_nig_sylv_Nov18.xlsx b/DataFormat/RawData/Data_Nyungwe_nig_sylv_Nov18.xlsx deleted file mode 100644 index 30d60a3..0000000 Binary files a/DataFormat/RawData/Data_Nyungwe_nig_sylv_Nov18.xlsx and /dev/null differ diff --git a/DataFormat/RawData/Data_UDZ_harveyi_spadix_Nov18.xlsx b/DataFormat/RawData/Data_UDZ_harveyi_spadix_Nov18.xlsx index 53b95be..35d9bb8 100644 Binary files a/DataFormat/RawData/Data_UDZ_harveyi_spadix_Nov18.xlsx and b/DataFormat/RawData/Data_UDZ_harveyi_spadix_Nov18.xlsx differ diff --git a/DataFormat/RawData/Data_Virunga_nigrifrons_Nov18.xlsx b/DataFormat/RawData/Data_Virunga_nigrifrons_Nov18.xlsx index 501987f..51df0f8 100644 Binary files a/DataFormat/RawData/Data_Virunga_nigrifrons_Nov18.xlsx and b/DataFormat/RawData/Data_Virunga_nigrifrons_Nov18.xlsx differ diff --git a/DataFormat/RawData/Duiker background.xlsx b/DataFormat/RawData/Duiker background.xlsx deleted file mode 100644 index 2c588af..0000000 Binary files a/DataFormat/RawData/Duiker background.xlsx and /dev/null differ diff --git a/DataFormat/RawData/Remaining_ungulates.xlsx b/DataFormat/RawData/Remaining_ungulates.xlsx new file mode 100644 index 0000000..5c85b22 Binary files /dev/null and b/DataFormat/RawData/Remaining_ungulates.xlsx differ diff --git a/DataFormat/RawData/UDZ_ungulates_2018_2019.xlsx b/DataFormat/RawData/UDZ_ungulates_2018_2019.xlsx new file mode 100644 index 0000000..d4f8060 Binary files /dev/null and b/DataFormat/RawData/UDZ_ungulates_2018_2019.xlsx differ diff --git a/DataFormat/RawData/trap locations and dates.xlsx b/DataFormat/RawData/trap locations and dates.xlsx new file mode 100644 index 0000000..254c1a4 Binary files /dev/null and b/DataFormat/RawData/trap locations and dates.xlsx differ diff --git a/MSdyn.txt b/MSdyn.txt deleted file mode 100644 index 107a228..0000000 --- a/MSdyn.txt +++ /dev/null @@ -1,47 +0,0 @@ -model{ - -mu.b0 ~ dnorm(0, 0.1) -tau.b0 ~ dgamma(0.1, 0.1) -mu.a0 ~ dnorm(0, 0.1) -tau.a0 ~ dgamma(0.1, 0.1) -mu.o0 ~ dnorm(0, 0.1) -tau.o0 ~ dgamma(0.1, 0.1) -mu.g0 ~ dnorm(0, 0.1) -tau.g0 ~ dgamma(0.1, 0.1) - -for(s in 1:nspec){ -log(lambda[s]) <- beta0[s] -beta0[s] ~ dnorm(mu.b0, tau.b0) - -alpha0[s] ~ dnorm(mu.a0, tau.a0) -omega0[s] ~ dnorm(mu.o0, tau.o0) -gamma0[s] ~ dnorm(mu.g0, tau.g0) - -for(j in 1:nsites){ - -logit(r[1,j,s]) <- alpha0[s] -N[1,j,s] ~ dpois(lambda[s]) - -for(t in 2:nyrs){ - -logit(r[t,j,s]) <- alpha0[s] -logit(omega[t-1,j,s]) <- omega0[s] -S[t-1,j,s] ~ dbin(omega[t-1,j,s], N[t-1,j,s]) -G[t-1,j,s] ~ dpois(gamma[t-1,s]) -N[t,j,s] <- S[t-1,j,s] + G[t-1,j,s] - -}#end t -}#end j - -for(t in 2:nyrs){ -log(gamma[t-1,s]) <- gamma0[s] -}#end t - -}#end s - -for(k in 1:nobs){ -y[k] ~ dbern(p[k]) -p[k] <- 1 - pow((1 - r[yr[k], site[k], spec[k]]), N[yr[k], site[k], spec[k]]) -}#end k - -} \ No newline at end of file diff --git a/MSdyn_cov.txt b/MSdyn_cov.txt deleted file mode 100644 index b0d8e4f..0000000 --- a/MSdyn_cov.txt +++ /dev/null @@ -1,79 +0,0 @@ -model{ - -mu.b0 ~ dnorm(0, 0.1) -tau.b0 <- 1/(sig.b0 * sig.b0) -sig.b0 ~ dt(0, pow(2.5,-2), 1) T(0,) -mu.b1 ~ dnorm(0, 0.1) -tau.b1 ~ dgamma(0.1, 0.1) -mu.b2 ~ dnorm(0, 0.1) -tau.b2 ~ dgamma(0.1, 0.1) -mu.b3 ~ dnorm(0, 0.1) -tau.b3 ~ dgamma(0.1, 0.1) -mu.b4 ~ dnorm(0, 0.1) -tau.b4 ~ dgamma(0.1, 0.1) -mu.a0 ~ dnorm(0, 0.368) -tau.a0 ~ dgamma(0.1, 0.1) -mu.o0 ~ dnorm(0, 0.368) -tau.o0 <- 1/(sig.o0 * sig.o0) -sig.o0 ~ dt(0, pow(2.5,-2), 1) T(0,) -mu.o1 ~ dnorm(0, 0.368) -tau.o1 <- 1/(sig.o1 * sig.o1) -sig.o1 ~ dt(0, pow(2.5,-2), 1) T(0,) -mu.o2 ~ dnorm(0, 0.368) -tau.o2 <- 1/(sig.o2 * sig.o2) -sig.o2 ~ dt(0, pow(2.5,-2), 1) T(0,) -mu.o3 ~ dnorm(0, 0.368) -tau.o3 <- 1/(sig.o3 * sig.o3) -sig.o3 ~ dt(0, pow(2.5,-2), 1) T(0,) -mu.g0 ~ dnorm(0, 0.1) -tau.g0 <- 1/(sig.g0 * sig.g0) -sig.g0 ~ dt(0, pow(2.5,-2), 1) T(0,) - -alpha1 ~ dnorm(0, 0.368) - -for(s in 1:nspec){ - -beta0[s] ~ dnorm(mu.b0, tau.b0) -beta1[s] ~ dnorm(mu.b1, tau.b1) -beta2[s] ~ dnorm(mu.b2, tau.b2) -beta3[s] ~ dnorm(mu.b3, tau.b3) -beta4[s] ~ dnorm(mu.b4, tau.b4) -alpha0[s] ~ dnorm(mu.a0, tau.a0) -omega0[s] ~ dnorm(mu.o0, tau.o0) -omega1[s] ~ dnorm(mu.o1, tau.o1) -omega2[s] ~ dnorm(mu.o2, tau.o2) -omega3[s] ~ dnorm(mu.o3, tau.o3) -gamma0[s] ~ dnorm(mu.g0, tau.g0) - -for(j in 1:nsites){ - -logit(r[nstart[j],j,s]) <- alpha0[s] + alpha1 * days[nstart[j],j] - -N[nstart[j],j,s] ~ dpois(lambda[j,s]) -log(lambda[j,s]) <- beta0[s] + beta1[s] * density[j] + beta2[s] * edge[nstart[j],j] + - beta3[s] * density[j] * edge[nstart[j],j] + beta4[s] * elev[nstart[j],j] - -for(t in (nstart[j]+1):nend[j]){ - -logit(r[t,j,s]) <- alpha0[s] + alpha1 * days[t,j] -logit(omega[t-1,j,s]) <- omega0[s] + omega1[s] * density[j] + omega2[s] * edge[t,j] + - omega3[s] * density[j] * edge[t,j] -S[t-1,j,s] ~ dbin(omega[t-1,j,s], N[t-1,j,s]) -G[t-1,j,s] ~ dpois(gamma[t-1,s]) -N[t,j,s] <- S[t-1,j,s] + G[t-1,j,s] - -}#end t -}#end j - -for(t in 2:9){ -log(gamma[t-1,s]) <- gamma0[s] -}#end t - -}#end s - -for(k in 1:nobs){ -y[k] ~ dbern(p[k]) -p[k] <- 1 - pow((1 - r[yr[k], site[k], spec[k]]), N[yr[k], site[k], spec[k]]) -}#end k - -} \ No newline at end of file diff --git a/MSstatic.txt b/MSstatic.txt deleted file mode 100644 index 36fc9d5..0000000 --- a/MSstatic.txt +++ /dev/null @@ -1,31 +0,0 @@ -model{ - -mu.a0 ~ dnorm(0, 0.1) -tau.a0 ~ dgamma(0.1, 0.1) - -mu.b0 ~ dnorm(0, 0.1) -tau.b0 ~ dgamma(0.1, 0.1) - -for(j in 1:nsite){ -gamma[j] ~ dnorm(0, tau.g) -} - -tau.g ~ dgamma(0.1, 0.1) - -for(i in 1:nspec){ - -alpha0[i] ~ dnorm(mu.a0, tau.a0) -beta0[i] ~ dnorm(mu.b0, tau.b0) - -for(t in 1:nyrs){ -for(j in 1:nsite){ -y[i,t,j] ~ dbin(v[t,j], p[i,t,j]) -p[i,t,j] <- 1 - pow((1 - r[i,t,j]), N[i,t,j]) -logit(r[i,t,j]) <- alpha0[i] -N[i,t,j] <- dpois(lambda[i,t,j]) -log(lambda[i,t,j]) <- beta0[i] + gamma[j] -}#end j -}#end t -}#end i - -} \ No newline at end of file diff --git a/PostAnalysis/Figure1.tiff b/PostAnalysis/Figure1.tiff new file mode 100644 index 0000000..351976f Binary files /dev/null and b/PostAnalysis/Figure1.tiff differ diff --git a/PostAnalysis/Figure2.tiff b/PostAnalysis/Figure2.tiff new file mode 100644 index 0000000..d12faea Binary files /dev/null and b/PostAnalysis/Figure2.tiff differ diff --git a/PostAnalysis/Figure3.tiff b/PostAnalysis/Figure3.tiff new file mode 100644 index 0000000..7e81e04 Binary files /dev/null and b/PostAnalysis/Figure3.tiff differ diff --git a/PostAnalysis/PostAnalysis.R b/PostAnalysis/PostAnalysis.R index 2d8af76..76e76f5 100644 --- a/PostAnalysis/PostAnalysis.R +++ b/PostAnalysis/PostAnalysis.R @@ -1,286 +1,577 @@ +#-----------# #-Libraries-# +#-----------# library(coda) -library(ggplot2) +library(tidyverse) +library(sf) library(ggthemes) -library(dplyr) library(reshape2) +library(grid) +library(gridExtra) +library(readxl) +library(nimble) +library(abind) +#-----------# #-Load data-# +#-----------# -load(file = "~/HMSNO/MSdynout.Rdata") load(file = "~/HMSNO/DataFormat/HMSNO.data.Rdata") load(file = "~/HMSNO/DataFormat/HMSNO.con.Rdata") -load(file = "~/HMSNO/DataFormat/Effort.Rdata") - - -# Effort <- Data %>% -# mutate(siteID = as.numeric(as.factor(`deployment ID`)), -# year = as.integer(as.factor(Year))) %>% -# group_by(park, year) %>% -# summarize(nsite = n_distinct(siteID), days = sum(days)) -# Effort$park <- as.factor(Effort$park) - -# Effort <- data.frame(park = HMSNO.data$park, t(HMSNO.data$days)) %>% -# melt(id = "park") %>% -# mutate(year = as.numeric(variable)) %>% -# group_by(park, year) %>% -# summarize(nsite = n(), days = sum(value, na.rm = TRUE)) -# -# group_by(park) - -nspec <- HMSNO.con$nspec -nparks <- max(HMSNO.data$park) -nyrs <- max(HMSNO.con$yr) -nyrstr <- HMSNO.con$nyrstr -nyrend <- HMSNO.con$nyrend -parkspec <- HMSNO.con$parkspec -#-Summary-# - -names <- factor(c("Peters's", "Bay", "Harvey's", "White-bellied", "Blue", - "Black-fronted", "Ogilby's", "Abbott's", "Yellow-backed", "Community"), - levels = c("Abbott's", "Bay", "Black-fronted", "Blue", - "Harvey's", "Ogilby's", "Peters's","White-bellied", - "Yellow-backed", "Community")) - -output <- MSdyn.o -nchains <- length(output) -out <- summary(output[c(1:nchains)][,543:578]) -ni <- length(unlist(output[c(1)][,1])) * nchains - -N <- array(NA, dim = c(nspec,nparks,nyrs,ni)) - -for(s in 1:nspec){ - for(i in 1:nparks){ - for(t in nyrstr[i]:nyrend[i]){ - N[s,i,t,] <- unlist(output[c(1:nchains)][,paste("Npark[", t, ", ", i, ", ", s, "]", sep = "")]) + +attach(HMSNO.data) +attach(HMSNO.con) + +#-------------------# +#-Load model output-# +#-------------------# + +pattern <- "chain" + +#Retrospective analysis +files <- list.files(path = "~/HMSNO/DataAnalysis/CaseStudy", pattern = pattern, full.names = TRUE) + +nc <- 5 + +for(i in 1:nc){ + load(files[i]) +} + +out <- mcmc.list(mget(ls()[grep(pattern, ls())])) + +rm(list = ls()[grep(pattern, ls())]) + +#------------# +#-Parameters-# +#------------# + +params <- attr(out[[1]], "dimnames")[[2]][!grepl("Npark|park.surv|pop.surv|park.gain|pop.gain", attr(out[[1]], "dimnames")[[2]])] + +#-------------# +#-Convergence-# +#-------------# + +Rhat <- gelman.diag(out[c(1:nc)][,params]) +if(all(Rhat[[1]][,1] < 1.1)){ + print("Converged") +}else{ + tmp <- as.numeric(which(Rhat[[1]][,1] > 1.1)) + print("Not converged") + print(params[tmp]) + traceplot(out[c(1:nc)][,params[tmp]]) +} + +#------------# +#-Extraction-# +#------------# + +for(i in 1:6){ + params <- c(params, paste0("park.surv[", i, "]")) +} + +for(s in 1:nspecs){ + for(i in parkS[s]:parkE[s]){ + for(t in nstart[i]:nend[i]){ + params <- c(params, paste0("Npark[", t, ", ", i, ", ", s, "]")) } + params <- c(params, paste0("pop.surv[", i, ", ", s, "]")) } } + +spec.names <- c("C. callipygus", "C. dorsalis", "C. harveyi", "C. leucogaster", + "P. monticola", "N. moschatus", "C. nigrifrons", "C. ogilbyi", + "T. scriptus", "C. spadix", "T. spekii", "C. silvicultor") + +park.names <- c("UDZ", "VNP", "NFNP", "BIF", "NNNP", "KRP") + +year.names <- function(x){ + return(format(as.Date("2008", format = "%Y") + lubridate::years(x), "%Y")) +} + +ni <- nc * length(out[[1]][,1]) + +#---------# +#-Results-# +#---------# +N <- array(NA, dim = c(nspecs,max(parkE),11,ni)) +park.surv <- array(NA, dim = c(max(parkE),ni)) +pop.surv <- array(NA, dim = c(nspecs,max(parkE),ni)) + +for(i in 1:6){ + park.surv[i,] <- unlist(out[c(1:nc)][,paste("park.surv[", i, "]", sep = "")]) +} + +for(s in 1:nspecs){ + for(i in parkS[s]:parkE[s]){ + pop.surv[s,i,] <- unlist(out[c(1:nc)][,paste("pop.surv[", i, ", ", s, "]", sep = "")]) + for(t in nstart[i]:nend[i]){ + N[s,i,t,] <- unlist(out[c(1:nc)][,paste("Npark[", t, ", ", i, ", ", s, "]", sep = "")]) + } + } +} meanN <- apply(N, MARGIN = c(1,2,3), mean, na.rm = TRUE) N97.5 <- apply(N, MARGIN = c(1,2,3), quantile, probs = 0.975, na.rm = TRUE) N2.5 <- apply(N, MARGIN = c(1,2,3), quantile, probs = 0.025, na.rm = TRUE) +meanpark.surv <- apply(park.surv, MARGIN = 1, mean, na.rm = TRUE) +park.surv97.5 <- apply(park.surv, MARGIN = 1, quantile, probs = 0.975, na.rm = TRUE) +park.surv75 <- apply(park.surv, MARGIN = 1, quantile, probs = 0.75, na.rm = TRUE) +park.surv25 <- apply(park.surv, MARGIN = 1, quantile, probs = 0.25, na.rm = TRUE) +park.surv2.5 <- apply(park.surv, MARGIN = 1, quantile, probs = 0.025, na.rm = TRUE) + +meanpop.surv <- apply(pop.surv, MARGIN = c(1,2), mean, na.rm = TRUE) +pop.surv97.5 <- apply(pop.surv, MARGIN = c(1,2), quantile, probs = 0.975, na.rm = TRUE) +pop.surv75 <- apply(pop.surv, MARGIN = c(1,2), quantile, probs = 0.75, na.rm = TRUE) +pop.surv25 <- apply(pop.surv, MARGIN = c(1,2), quantile, probs = 0.25, na.rm = TRUE) +pop.surv2.5 <- apply(pop.surv, MARGIN = c(1,2), quantile, probs = 0.025, na.rm = TRUE) + data <- reshape2::melt(meanN, varnames = c("species", "parkID", "yearID"), value.name = "abundance") -data$species <- factor(data$species, - labels = c("Peters's", "Bay", "Harvey's", "White-bellied", "Blue", - "Black-fronted", "Ogilby's", "Abbott's", "Yellow-backed")) -data$park <- as.factor(data$park) - -high <- reshape2::melt(N97.5, varnames = c("species", "parkID", "yearID"), value.name = "97.5%")$`97.5%` -low <- reshape2::melt(N2.5, varnames = c("species", "parkID", "yearID"), value.name = "2.5%")$`2.5%` -data$`97.5%` <- high -data$`2.5%` <- low - -data <- data %>% - left_join(., Effort, by = c("parkID", "yearID")) %>% - mutate(abundance_site = abundance/nsite, - `97.5%_site` = `97.5%`/nsite, - `2.5%_site` = `2.5%`/nsite, - abundance_days = abundance/days, - `97.5%_days` = `97.5%`/days, - `2.5%_days` = `2.5%`/days, - parkID = as.factor(parkID)) - -# data2 <- data.frame(t(rep(NA, dim(data)[2]))) -# -# for(s in 1:nspec){ -# for(i in parkspec[,s]){ -# data2 <- rbind(data2, data %>% filter(species == names[s] & parkID == i)) -# } -# } - -data %>% - #filter(park != "4" & park != "6") %>% - ggplot(., aes(x = yearID, y = abundance)) + - geom_line(aes(col = parkID)) + - geom_ribbon(aes(x = yearID, ymin = `2.5%`, ymax = `97.5%`, fill = parkID), alpha = 0.5) + - facet_wrap(. ~ species) + - theme_few() - -data %>% - #filter(park != "4" & park != "6") %>% - ggplot(., aes(x = yearID, y = abundance_site)) + - geom_line(aes(col = parkID)) + - geom_ribbon(aes(x = yearID, ymin = `2.5%_site`, ymax = `97.5%_site`, fill = parkID), alpha = 0.5) + - facet_wrap(. ~ species) + - theme_few() + - labs(y = "Abundance") +data$species <- factor(data$species, labels = spec.names) +data$park <- factor(data$park, labels = park.names) +data$year <- as.numeric(year.names(data$yearID)) +data$`97.5%` <- reshape2::melt(N97.5, varnames = c("species", "parkID", "yearID"), value.name = "97.5%")$`97.5%` +data$`2.5%` <- reshape2::melt(N2.5, varnames = c("species", "parkID", "yearID"), value.name = "2.5%")$`2.5%` +data$effort <- nsite[data$parkID] + +data <- data %>% + mutate(abundance_site = abundance/effort, + `97.5%_site` = `97.5%`/effort, + `2.5%_site` = `2.5%`/effort) + +data$park <- factor(data$park, levels = c("KRP", "NNNP", "UDZ", "BIF", "NFNP", "VNP")) + +parkdata <- data.frame(reshape2::melt(meanpark.surv, value.name = "mean"), + reshape2::melt(park.surv2.5, value.name = "2.5"), + reshape2::melt(park.surv25, value.name = "25"), + reshape2::melt(park.surv75, value.name = "75"), + reshape2::melt(park.surv97.5, value.name = "97.5")) +parkdata$park <- factor(park.names) +colnames(parkdata)[2:5] <- c("l2.5", "l25", "u75", "u97.5") + +popdata <- reshape2::melt(meanpop.surv, varnames = c("species", "parkID"), value.name = "mean") +popdata$`l2.5` <- reshape2::melt(pop.surv2.5, varnames = c("species", "parkID"), value.name = "2.5%")$`2.5` +popdata$`l25` <- reshape2::melt(pop.surv25, varnames = c("species", "parkID"), value.name = "25%")$`25` +popdata$`u75` <- reshape2::melt(pop.surv75, varnames = c("species", "parkID"), value.name = "75%")$`75` +popdata$`u97.5` <- reshape2::melt(pop.surv97.5, varnames = c("species", "parkID"), value.name = "97.5%")$`97.5` + +popdata$species <- factor(popdata$species, labels = spec.names) +popdata$park <- factor(popdata$parkID, labels = park.names) +popdata <- popdata[complete.cases(popdata),] +popdata <- popdata[,-2] + +#Summary statistics +measurements <- c("Initial", "Final", "PopGrowth") + +#Extract summary statistics +Data.pop <- array(NA, dim = c(nspecs,max(parkE),length(measurements),ni)) +for(s in 1:nspecs){ + for(i in parkS[s]:parkE[s]){ + Data.pop[s,i,1,] <- N[s,i,nstart[i],]/nsite[i] + Data.pop[s,i,2,] <- N[s,i,nend[i],]/nsite[i] + tmp2 <- 1 + for(t in (nstart[i]+1):nend[i]){ + tmp1 <- N[s,i,t,]/N[s,i,t-1,] + tmp2 <- tmp1 * tmp2 + } + Data.pop[s,i,3,] <- tmp2^(1/(nend[i]-nstart[i])) + } +} -data %>% - #filter(park != "4" & park != "6") %>% - ggplot(., aes(x = yearID, y = abundance)) + - geom_line(aes(col = species)) + - geom_ribbon(aes(x = yearID, ymin = `2.5%`, ymax = `97.5%`, fill = species), alpha = 0.5) + - facet_wrap(. ~ parkID) + - theme_few() - # theme(panel.background = element_rect(fill = "transparent", color = NA), - # axis.text.x = element_text(angle = 50, hjust = 0.5, vjust = 0.5), - # legend.position = "none") + - # labs(y ="Detection probability", x = expression()) - -Fig1 <- data %>% - mutate(year = factor(as.character(yearID))) %>% - #filter(park != "4" & park != "6") %>% - ggplot(., aes(x = year, y = abundance_site)) + - geom_line(aes(col = species)) + - geom_line(aes(x = yearID, col = species)) + - geom_ribbon(aes(x = yearID, ymin = `2.5%_site`, ymax = `97.5%_site`, fill = species), alpha = 0.5) + - facet_wrap(. ~ parkID, scales='free_x', labeller = labeller(parkID = c("1" = "Korup", - "2" = "Nouabale-Ndoki", - "3" = "Nyungwe", - "4" = "Bwindi", - "5" = "Virunga", - "6" = "Udzungwa"))) + - labs(y = "Abundance per site", x = "Year") + - theme_few() - -Fig2 <- data %>% - mutate(year = factor(as.character(yearID + 2008))) %>% - #filter(park != "4" & park != "6") %>% - ggplot(., aes(x = year, y = abundance_site)) + +#Calculate quantiles from MCMC +Data.mean <- apply(Data.pop, MARGIN = c(1,2,3), mean, na.rm = TRUE) +Data.97.5 <- apply(Data.pop, MARGIN = c(1,2,3), quantile, probs = 0.975, na.rm = TRUE) +Data.2.5 <- apply(Data.pop, MARGIN = c(1,2,3), quantile, probs = 0.025, na.rm = TRUE) + +#Format summary statistics into dataframe +sum.df <- full_join(full_join(dcast(melt(Data.mean, varnames = c("species", "parkID", "measurements")), formula = species + parkID ~ measurements), + dcast(melt(Data.2.5, varnames = c("species", "parkID", "measurements")), formula = species + parkID ~ measurements), + by = c("species", "parkID")), + dcast(melt(Data.97.5, varnames = c("species", "parkID", "measurements")), formula = species + parkID ~ measurements), + by = c("species", "parkID")) +colnames(sum.df) <- c("species", "park", paste0(measurements, rep(c("_mean", "_lower", "_upper"), each = length(measurements)))) +sum.df$species <- factor(sum.df$species, labels = spec.names) +sum.df$park <- factor(sum.df$park, labels = park.names) +sum.df <- sum.df[complete.cases(sum.df),] +sum.df <- sum.df[,c("species", "park", paste0(rep(measurements, each = 3), c("_mean", "_lower", "_upper")))] + +#Fill in missing values +for(k in 1:dim(sum.df)[1]){ + if(is.infinite(sum.df[k,"PopGrowth_mean"])){ + i <- as.numeric(sum.df[k,"species"]) + r <- as.numeric(sum.df[k,"park"]) + sum.df[k,"PopGrowth_mean"] <- mean(Data.pop[i,r,3,!is.infinite(Data.pop[i,r,3,])], na.rm = TRUE) + } + if(is.infinite(sum.df[k,"PopGrowth_lower"])){ + i <- as.numeric(sum.df[k,"species"]) + r <- as.numeric(sum.df[k,"park"]) + sum.df[k,"PopGrowth_lower"] <- quantile(Data.pop[i,r,3,!is.infinite(Data.pop[i,r,3,])], probs = 0.025, na.rm = TRUE) + } + if(is.infinite(sum.df[k,"PopGrowth_upper"])){ + i <- as.numeric(sum.df[k,"species"]) + r <- as.numeric(sum.df[k,"park"]) + sum.df[k,"PopGrowth_upper"] <- quantile(Data.pop[i,r,3,!is.infinite(Data.pop[i,r,3,])], probs = 0.975, na.rm = TRUE) + } +} + +#Initial abundances minimum +sum.df[which.min(sum.df$Initial_mean), c(1:5)] + +#Initial abundances maximum +sum.df[which.max(sum.df$Initial_mean), c(1:5)] + +#Final abundances minimum +sum.df[which.min(sum.df$Final_mean), c(1:2, 6:8)] + +#Final abundances maximum +sum.df[which.max(sum.df$Final_mean), c(1:2, 6:8)] + +#Population growth rate increasing +sum(sum.df$PopGrowth_mean > 1) + +sum(sum.df$PopGrowth_mean < 0.8) + +#Population growth rate minimum +sum.df[which.min(sum.df$PopGrowth_mean), c(1:2, 9:11)] + +#Population growth rate maximum +sum.df[which.max(sum.df$PopGrowth_mean), c(1:2, 9:11)] + +#Population growth rate decreasing +sum.df[which(sum.df$PopGrowth_upper<1), c(1:2, 9:11)] + +#NNNP +sum.df %>% filter(park == "NNNP") + +sum.df %>% filter(park == "VNP") + +sum.df %>% filter(park == "UDZ") + +sum.df %>% filter(park == "BIF") + +sum.df %>% filter(park == "NFNP") + +sum.df %>% filter(park == "KRP") + +#Save table +sum.df %>% arrange(park) %>% select(park, species, PopGrowth_mean, PopGrowth_lower, PopGrowth_upper) %>% + mutate(PopGrowth_mean = round(PopGrowth_mean, digits = 2), + PopGrowth_lower = round(PopGrowth_lower, digits = 2), + PopGrowth_upper = round(PopGrowth_upper, digits = 2)) %>% + write.csv(., file = "./SupportingInformation/AppS3TableS1.csv") + +#-------------------# +#-Covariate effects-# +#-------------------# + +output <- summary(out[c(1:nc)][,params]) + +omega0 <- data.frame(species = factor(c(spec.names, "community"), levels = c(spec.names, "community")), mean = c(plogis(output$statistics[grep("omega0", params),"Mean"]), output$statistics["mu.o0","Mean"]), + l2.5 = c(plogis(output$quantiles[grep("omega0", params),"2.5%"]), output$quantiles["mu.o0","2.5%"]), + l25 = c(plogis(output$quantiles[grep("omega0", params),"25%"]), output$quantiles["mu.o0","25%"]), + u75 = c(plogis(output$quantiles[grep("omega0", params),"75%"]), output$quantiles["mu.o0","75%"]), + u97.5 = c(plogis(output$quantiles[grep("omega0", params),"97.5%"]), output$quantiles["mu.o0","97.5%"])) + +popdata <- full_join(popdata, omega0 %>% mutate(park = "TEAM")) +popdata <- full_join(popdata, parkdata %>% mutate(species = "community")) +popdata$park <- factor(popdata$park, levels = c("KRP", "NNNP", "UDZ", "BIF", "NFNP", "VNP", "TEAM")) +popdata$species <- factor(popdata$species, levels = c(spec.names, "community")) + +#-Gains-# + +gamma0 <- data.frame(species = factor(c(spec.names, "community"), levels = c(spec.names, "community")), mean = c(exp(output$statistics[grep("gamma0", params),"Mean"]), exp(output$statistics["mu.g0","Mean"])), + l2.5 = c(exp(output$quantiles[grep("gamma0", params),"2.5%"]), exp(output$quantiles["mu.g0","2.5%"])), + l25 = c(exp(output$quantiles[grep("gamma0", params),"25%"]), exp(output$quantiles["mu.g0","25%"])), + u75 = c(exp(output$quantiles[grep("gamma0", params),"75%"]), exp(output$quantiles["mu.g0","75%"])), + u97.5 = c(exp(output$quantiles[grep("gamma0", params),"97.5%"]), exp(output$quantiles["mu.g0","97.5%"]))) + +data %>% drop_na(.) %>% select(species, park, year, abundance_site, `2.5%_site`, `97.5%_site`) %>% + arrange(park, species, year) %>% + mutate(abundance_site = round(abundance_site, digits = 2), + `2.5%_site` = round(`2.5%_site`, digits = 2), + `97.5%_site` = round(`97.5%_site`, digits = 2), + abundance = paste0(abundance_site, " (CI ", `2.5%_site`, ", ", `97.5%_site`, ")" )) %>% + select(park, species, year, abundance) %>% + write.csv(., file = "./SupportingInformation/AppS3TableS2.csv") + +#Table S3 +full_join(omega0 %>% mutate(mean = format(round(mean, digits=2), nsmall = 2), + l2.5 = format(round(l2.5, digits=2), nsmall = 2), + u97.5 = format(round(u97.5, digits=2), nsmall = 2), + `Annual apparent survival` = paste0(mean, " (CI ", l2.5, ", ", u97.5, ")" )) %>% + select(species, `Annual apparent survival`), + gamma0 %>% mutate(mean = format(round(mean, digits=2), nsmall = 2), + l2.5 = format(round(l2.5, digits=2), nsmall = 2), + u97.5 = format(round(u97.5, digits=2), nsmall = 2), + Gains = paste0(mean, " (CI ", l2.5, ", ", u97.5, ")" )) %>% + select(species, Gains), + by = "species") %>% + write.csv(., file = "./SupportingInformation/AppS3TableS3.csv") + +#Table S4 +popdata %>% filter(!(park == "TEAM") & species == "community") %>% + mutate(mean = format(round(mean, digits=2), nsmall = 2), + l2.5 = format(round(l2.5, digits=2), nsmall = 2), + u97.5 = format(round(u97.5, digits=2), nsmall = 2), + `Annual apparent survival` = paste0(mean, " (CI ", l2.5, ", ", u97.5, ")" )) %>% + select(park, `Annual apparent survival`) %>% + write.csv(., file = "./SupportingInformation/AppS3TableS4.csv") + +#TEAM community mean +popdata[popdata$species == "community" & popdata$park == "TEAM",] + +#Range of species survival +sd.o0 <- mean(1/sqrt(unlist(out[c(1:nc)][,"tau.o0"]))) +sd.o0; quantile(1/sqrt(unlist(out[c(1:nc)][,"tau.o0"])), c(0.025, 0.975)) + +plogis(logit(output$statistics["mu.o0","Mean"]) + sd.o0) - +plogis(logit(output$statistics["mu.o0","Mean"]) - sd.o0) + +#Range of park survival +sd.eps.o <- mean(1/sqrt(unlist(out[c(1:nc)][,"tau.eps.o"]))) +sd.eps.o; quantile(1/sqrt(unlist(out[c(1:nc)][,"tau.eps.o"])), c(0.025, 0.975)) + +plogis(logit(output$statistics["mu.o0","Mean"]) + sd.eps.o) - + plogis(logit(output$statistics["mu.o0","Mean"]) - sd.eps.o) + + +#Species vs park magnitude +data.frame(species = factor(spec.names), + mean = output$statistics[grep("omega0", params),"Mean"] - logit(output$statistics["mu.o0","Mean"]), + l2.5 = output$quantiles[grep("omega0", params),"2.5%"] - logit(output$quantiles["mu.o0","2.5%"]), + l25 = output$quantiles[grep("omega0", params),"25%"] - logit(output$quantiles["mu.o0","25%"]), + u75 = output$quantiles[grep("omega0", params),"75%"] - logit(output$quantiles["mu.o0","75%"]), + u97.5 = output$quantiles[grep("omega0", params),"97.5%"] - logit(output$quantiles["mu.o0","97.5%"])) + +data.frame(park = factor(park.names), + mean = output$statistics[grep("eps.o", params),"Mean"][-7], + l2.5 = output$quantiles[grep("eps.o", params),"2.5%"][-7], + l25 = output$quantiles[grep("eps.o", params),"25%"][-7], + u75 = output$quantiles[grep("eps.o", params),"75%"][-7], + u97.5 = output$quantiles[grep("eps.o", params),"97.5%"][-7]) + +#Maximum species survival +max(popdata[popdata$species != "community",'mean']) + +popdata[popdata$species != "community",] + +popdata[popdata$species == "community" & popdata$park != "TEAM",] + +#TEAM community gains +gamma0[gamma0$species == 'community',] + +#----------# +#-Figure 2-# +#----------# + +pos <- position_dodge(0.75, preserve = "total") + +data$park <- factor(data$park, labels = c("KRP" = "Korup National Park", + "NNNP" = "Nouabale-Ndoki National Park", + "UDZ" = "Udzungwa National Park", + "BIF" = "Bwindi National Park", + "NFNP" = "Nyungwe Forest National Park", + "VNP" = "Volcanoes National Park")) + +Figure2 <- ggplot(data %>% drop_na(.), aes(x = year, y = abundance_site)) + + geom_point(col = "grey") + geom_line(aes(col = species)) + - geom_line(aes(x = yearID, col = species)) + - geom_ribbon(aes(x = yearID, ymin = `2.5%_site`, ymax = `97.5%_site`, fill = species), alpha = 0.5) + - facet_wrap(. ~ parkID, scales='free', labeller = labeller(parkID = c("1" = "Korup", - "2" = "Nouabale-Ndoki", - "3" = "Nyungwe", - "4" = "Bwindi", - "5" = "Virunga", - "6" = "Udzungwa"))) + - labs(y = "Abundance per site", x = "Year") + - theme_few() - -# alpha0 <- data.frame(species = names, mean = c(plogis(out$statistics[1:9,"Mean"]), out$statistics[13,"Mean"]), -# l2.5 = c(plogis(out$quantiles[1:9,"2.5%"]), out$quantiles[13,"2.5%"]), -# l25 = c(plogis(out$quantiles[1:9,"25%"]), out$quantiles[13,"25%"]), -# u75 = c(plogis(out$quantiles[1:9,"75%"]), out$quantiles[13,"75%"]), -# u97.5 = c(plogis(out$quantiles[1:9,"97.5%"]), out$quantiles[13,"97.5%"])) - -omega0 <- data.frame(species = names, mean = c(plogis(out$statistics[6:14,"Mean"]), out$statistics[3,"Mean"]), - l2.5 = c(plogis(out$quantiles[6:14,"2.5%"]), out$quantiles[3,"2.5%"]), - l25 = c(plogis(out$quantiles[6:14,"25%"]), out$quantiles[3,"25%"]), - u75 = c(plogis(out$quantiles[6:14,"75%"]), out$quantiles[3,"75%"]), - u97.5 = c(plogis(out$quantiles[6:14,"97.5%"]), out$quantiles[3,"97.5%"])) - -omega1 <- data.frame(species = names, mean = c(out$statistics[15:23,"Mean"], out$statistics[4,"Mean"]), - l2.5 = c(out$quantiles[15:23,"2.5%"], out$quantiles[4,"2.5%"]), - l25 = c(out$quantiles[15:23,"25%"], out$quantiles[4,"25%"]), - u75 = c(out$quantiles[15:23,"75%"], out$quantiles[4,"75%"]), - u97.5 = c(out$quantiles[15:23,"97.5%"], out$quantiles[4,"97.5%"])) - -omega2 <- data.frame(species = names, mean = c(out$statistics[24:32,"Mean"], out$statistics[5,"Mean"]), - l2.5 = c(out$quantiles[24:32,"2.5%"], out$quantiles[5,"2.5%"]), - l25 = c(out$quantiles[24:32,"25%"], out$quantiles[5,"25%"]), - u75 = c(out$quantiles[24:32,"75%"], out$quantiles[5,"75%"]), - u97.5 = c(out$quantiles[24:32,"97.5%"], out$quantiles[5,"97.5%"])) - -# FigA0 <- ggplot(alpha0) + -# geom_errorbar(aes(x = species, ymin = mean, ymax = mean), -# width = 0.275) + -# geom_errorbar(aes(x = species, ymin = l2.5, ymax = u97.5), -# width = 0, size = 1.25) + -# geom_errorbar(aes(x = species, ymin = l25, ymax = u75), -# width = 0, size = 3) + -# coord_cartesian(ylim = c(0, 0.65)) + -# scale_y_continuous(expand = c(0, 0)) + -# geom_hline(yintercept = 0, alpha = 0.75) + -# theme_few() + -# theme(panel.background = element_rect(fill = "transparent", color = NA), -# axis.text.x = element_text(angle = 50, hjust = 0.5, vjust = 0.5), -# legend.position = "none") + -# labs(y ="Detection probability", x = expression()) - -FigO0 <- ggplot(omega0) + - geom_errorbar(aes(x = species, ymin = mean, ymax = mean), - width = 0.275) + - geom_errorbar(aes(x = species, ymin = l2.5, ymax = u97.5), - width = 0, size = 1.25) + - geom_errorbar(aes(x = species, ymin = l25, ymax = u75), - width = 0, size = 3) + - coord_cartesian(ylim = c(0, 1)) + + geom_ribbon(aes(x = year, ymin = `2.5%_site`, ymax = `97.5%_site`, fill = species), alpha = 0.5) + + facet_wrap(. ~ park, scales = 'free', ncol = 2, dir = "v") + + scale_x_continuous(breaks = c(2010, 2012, 2014, 2016, 2018, 2020)) + + theme_few() + + theme( + legend.title = element_blank(), + legend.position = "bottom", + legend.text = element_text(face = "italic")) + + labs(y = "Relative abundance/site\n", x = "Year") + + +tiff(file = "~/HMSNO/PostAnalysis/Figure2.tiff", res = 600, width = 6, height = 6, units = "in") +Figure2 +dev.off() + +#----------# +#-Figure 3-# +#----------# + +Fig3A <- popdata %>% filter(park == "TEAM" | species == "community") %>% + mutate(xname = ifelse(park == "TEAM", as.character(species), as.character(park))) %>% + mutate(xname = factor(xname, levels = c(spec.names, "community", park.names))) %>% + ggplot(., aes(x = xname, color = park)) + + geom_boxplot(aes(ymin = l2.5, lower = l25, middle = mean, upper = u75, ymax = u97.5), + stat = "identity", size = 0.75, width = 0.5) + + scale_color_manual(values = c("#66c2a5", "#fc8d62", "#8da0cb", "#e78ac3", "#a6d854", "#ffd92f", "black")) + theme_few() + theme(panel.background = element_rect(fill = "transparent", color = NA), - axis.text.x = element_text(angle = 50, hjust = 0.5, vjust = 0.5), - legend.position = "none") + - labs(y ="Apparent survival probability", x = expression()) - - -FigO1 <- ggplot(omega1) + - geom_errorbar(aes(x = species, ymin = mean, ymax = mean), - width = 0.275) + - geom_errorbar(aes(x = species, ymin = l2.5, ymax = u97.5), - width = 0, size = 1.25) + - geom_errorbar(aes(x = species, ymin = l25, ymax = u75), - width = 0, size = 3) + - coord_cartesian(ylim = c(-7, 7)) + - geom_hline(yintercept = 0, alpha = 0.75) + + axis.text.x = element_blank(), + axis.ticks.x = element_blank(), + legend.title = element_blank(), + strip.placement = "outside", + strip.background = element_blank(), + strip.text = element_text(size = 12)) + + labs(y = "Apparent survival", x = expression()) + + geom_vline(xintercept=(13 - 0.5), linetype = "dotdash") + +Fig3B <- gamma0 %>% + ggplot(., aes(x = species)) + + geom_boxplot(aes(ymin = l2.5, lower = l25, middle = mean, upper = u75, ymax = u97.5), + stat = "identity", size = 0.75, width = 0.5, col = "black") + theme_few() + theme(panel.background = element_rect(fill = "transparent", color = NA), - axis.text.x = element_text(angle = 50, hjust = 0.5, vjust = 0.5), - legend.position = "none") + - labs(y ="Effect of density", x = expression()) - -FigO2 <- ggplot(omega2) + - geom_errorbar(aes(x = species, ymin = mean, ymax = mean), - width = 0.275) + - geom_errorbar(aes(x = species, ymin = l2.5, ymax = u97.5), - width = 0, size = 1.25) + - geom_errorbar(aes(x = species, ymin = l25, ymax = u75), - width = 0, size = 3) + - coord_cartesian(ylim = c(-1.5, 1.5)) + - geom_hline(yintercept = 0, alpha = 0.75) + + axis.text.x = element_text(angle = 50, hjust = 0.5, vjust = 0.5, face = "italic"), + legend.title = element_blank(), + strip.placement = "outside", + strip.background = element_blank(), + strip.text = element_text(size = 12)) + + labs(y = "Individuals gained \nper site", x = expression()) + + geom_vline(xintercept=(13 - 0.5), linetype = "dotdash") + +Fig3A <- ggplotGrob(Fig3A) +Fig3B <- ggplotGrob(Fig3B) + +Fig3B$heights <- Fig3A$heights + +Legend <- Fig3A$grobs[[15]] + +Fig3A$grobs[[15]] <- nullGrob() + +Fig3B$widths <- Fig3A$widths + +Fig3A$widths[9] <- unit(0, "cm") + +Fig3B$widths[9] <- sum(unit(5.4, "cm"), unit(2.5, "point")) + +tiff(file = "~/HMSNO/PostAnalysis/Figure3.tiff", res = 600, width = 8, height = 5, units = "in") +grid.arrange(arrangeGrob(Fig3A, Fig3B), bottom = "\n \n", left = "") +grid.draw(textGrob("A", x = unit(0.5, "cm"), y = unit(12.25, "cm"), + gp=grid::gpar(fontsize=14, fontface = 2))) +grid.draw(textGrob("B", x = unit(0.5, "cm"), y = unit(6.75, "cm"), + gp=grid::gpar(fontsize=14, fontface = 2))) +pushViewport(viewport(x = unit(17, "cm"), y = unit(4.75, "cm"))) +grid.draw(Legend) +dev.off() + +#------------------# +#-Simulation study-# +#------------------# + +#List of all filenames of 75 site output +filenames_75 <- list.files(path = "~/HMSNO/DataAnalysis/SimulationStudy/SimulationOutput/output_75", full.names = TRUE) +filenames_25 <- list.files(path = "~/HMSNO/DataAnalysis/SimulationStudy/SimulationOutput/output_25", full.names = TRUE) + +filenames <- c(filenames_75, filenames_25) + +#File length of each scenario +n75 <- length(filenames_75) +n25 <- length(filenames_25) + +#Load first file +load(filenames[1]) + +#Initialize vector for all output +out <- output[,2:7] + +#Unconverged index +ii <- NULL + +#Harvout parameters from files +for(i in 2:(n75 + n25)){ + load(filenames[i]) + if(max(output[,c(5,8)], na.rm = TRUE) < 1.1){ + out <- abind(out, output[,2:7], along = 3) + }else{ + ii <- c(ii, i) + } +} + +jj <- seq(1,c(n75+n25))[-ii] +nn75 <- max(which(jj<500)) +nn25 <- length(jj) + +#-Summary dataframe-# + +ms75_75 <- apply((out[1:12,2,1:nn75] - out[1:12,1,1:nn75])/out[1:12,1,1:nn75], MARGIN = 1, FUN = quantile, probs = 0.75, na.rm = TRUE) +ms50_75 <- apply((out[1:12,2,1:nn75] - out[1:12,1,1:nn75])/out[1:12,1,1:nn75], MARGIN = 1, FUN = quantile, probs = 0.5, na.rm = TRUE) +ms25_75 <- apply((out[1:12,2,1:nn75] - out[1:12,1,1:nn75])/out[1:12,1,1:nn75], MARGIN = 1, FUN = quantile, probs = 0.25, na.rm = TRUE) +msmax_75 <- ((ms75_75 - ms25_75) * 1.5) + ms75_75 +msmin_75 <- ms25_75 - ((ms75_75 - ms25_75) * 1.5) + +ss75_75 <- apply((out[1:12,5,1:nn75] - out[1:12,1,1:nn75])/out[1:12,1,1:nn75], MARGIN = 1, FUN = quantile, probs = 0.75, na.rm = TRUE) +ss50_75 <- apply((out[1:12,5,1:nn75] - out[1:12,1,1:nn75])/out[1:12,1,1:nn75], MARGIN = 1, FUN = quantile, probs = 0.5, na.rm = TRUE) +ss25_75 <- apply((out[1:12,5,1:nn75] - out[1:12,1,1:nn75])/out[1:12,1,1:nn75], MARGIN = 1, FUN = quantile, probs = 0.25, na.rm = TRUE) +ssmax_75 <- ((ss75_75 - ss25_75) * 1.5) + ss75_75 +ssmin_75 <- ss25_75 - ((ss75_75 - ss25_75) * 1.5) + +ms75_25 <- apply((out[1:12,2,(nn75+1):nn25] - out[1:12,1,(nn75+1):nn25])/out[1:12,1,(nn75+1):nn25], MARGIN = 1, FUN = quantile, probs = 0.75, na.rm = TRUE) +ms50_25 <- apply((out[1:12,2,(nn75+1):nn25] - out[1:12,1,(nn75+1):nn25])/out[1:12,1,(nn75+1):nn25], MARGIN = 1, FUN = quantile, probs = 0.5, na.rm = TRUE) +ms25_25 <- apply((out[1:12,2,(nn75+1):nn25] - out[1:12,1,(nn75+1):nn25])/out[1:12,1,(nn75+1):nn25], MARGIN = 1, FUN = quantile, probs = 0.25, na.rm = TRUE) +msmax_25 <- ((ms75_25 - ms25_25) * 1.5) + ms75_25 +msmin_25 <- ms25_25 - ((ms75_25 - ms25_25) * 1.5) + +ss75_25 <- apply((out[1:12,5,(nn75+1):nn25] - out[1:12,1,(nn75+1):nn25])/out[1:12,1,(nn75+1):nn25], MARGIN = 1, FUN = quantile, probs = 0.75, na.rm = TRUE) +ss50_25 <- apply((out[1:12,5,(nn75+1):nn25] - out[1:12,1,(nn75+1):nn25])/out[1:12,1,(nn75+1):nn25], MARGIN = 1, FUN = quantile, probs = 0.5, na.rm = TRUE) +ss25_25 <- apply((out[1:12,5,(nn75+1):nn25] - out[1:12,1,(nn75+1):nn25])/out[1:12,1,(nn75+1):nn25], MARGIN = 1, FUN = quantile, probs = 0.25, na.rm = TRUE) +ssmax_25 <- ((ss75_25 - ss25_25) * 1.5) + ss75_25 +ssmin_25 <- ss25_25 - ((ss75_25 - ss25_25) * 1.5) + +data <- data.frame(rbind(cbind(msmin_75, ms25_75, ms50_75, ms75_75, msmax_75), + cbind(ssmin_75, ss25_75, ss50_75, ss75_75, ssmax_75), + cbind(msmin_25, ms25_25, ms50_25, ms75_25, msmax_25), + cbind(ssmin_25, ss25_25, ss50_25, ss75_25, ssmax_25))) + +params <- as.factor(output[,1]) +data$params <- rep(params[1:12], 4) +data$framework <- rep(rep(c("multi-species", "single-species"), each = 12), 2) +data$sites <- rep(c("75", "25"), each = 24) +colnames(data)[1:5] <- c("qmin","q25","q50","q75","qmax") +data <- data %>% separate(params, c("params", "species.type")) +data$framework <- factor(data$framework, labels = c("multi-species" = "Multi-species", + "single-species" = "Single-species")) +data$species.type <- factor(data$species.type, labels = c("common" = "Common", + "elusive" = "Elusive", + "rare" = "Rare")) +data$params <- factor(data$params, + labels = c("gamma" = "Gains", + "lambda" = "Initial Abundance", + "omega" = "Apparent Survival", + "r" = "Detection Probability")) +data$params <- factor(data$params, levels = c("Detection Probability", "Initial Abundance", "Apparent Survival", "Gains")) +data$sites <- factor(data$sites, labels = c("25" = "25 sites", + "75" = "75 sites")) +data$species.type <- factor(data$species.type, levels = c("Common", "Rare", "Elusive")) + +#----------# +#-Figure 1-# +#----------# + +Figure1 <- ggplotGrob(ggplot(data, aes(x = sites, fill = framework)) + + geom_boxplot(aes(ymin = qmin, lower = q25, middle = q50, upper = q75, ymax = qmax), + stat = "identity", size = 0.75) + + scale_fill_manual(values = c("white", "grey")) + + facet_grid(params ~ species.type, scale = "free", + labeller = label_wrap_gen(width = 2, multi_line = TRUE)) + + geom_hline(yintercept = 0, col = "black", size = 1) + theme_few() + theme(panel.background = element_rect(fill = "transparent", color = NA), - axis.text.x = element_text(angle = 50, hjust = 0.5, vjust = 0.5), - legend.position = "none") + - labs(y ="Effect of distance", x = expression()) - -#-Single-species-# - -# species <- c("callipygus", "dorsalis", "harveyi", "leucogaster", "monticola", "nigrifrons") -# names <- c("Peters's", "Bay", "Harvey's", "White-bellied", "Blue", "Black-fronted") -# -# alpha0 <- beta0 <- beta1 <- beta2 <- beta3 <- -# beta4 <- gamma0 <- omega0 <- omega1 <- omega2 <- omega3 <- rep(NA, 6) -# -# for(i in 1:length(species)){ -# load(file = paste("Z:/HMSNO/", species[i], ".Rdata", sep="")) -# output <- summary(as.mcmc.list(out)) -# alpha0 <- rbind(alpha0, c("mean" = output$statistics[1,1], output$quantiles[1,])) -# beta0 <- rbind(beta0, c("mean" = output$statistics[3,1], output$quantiles[3,])) -# beta1 <- rbind(beta1, c("mean" = output$statistics[4,1], output$quantiles[4,])) -# beta2 <- rbind(beta2, c("mean" = output$statistics[5,1], output$quantiles[5,])) -# beta3 <- rbind(beta3, c("mean" = output$statistics[6,1], output$quantiles[6,])) -# beta4 <- rbind(beta4, c("mean" = output$statistics[7,1], output$quantiles[7,])) -# gamma0 <- rbind(gamma0, c("mean" = output$statistics[8,1], output$quantiles[8,])) -# omega0 <- rbind(omega0, c("mean" = output$statistics[9,1], output$quantiles[9,])) -# omega1 <- rbind(omega1, c("mean" = output$statistics[10,1], output$quantiles[10,])) -# omega2 <- rbind(omega2, c("mean" = output$statistics[11,1], output$quantiles[11,])) -# omega3 <- rbind(omega3, c("mean" = output$statistics[12,1], output$quantiles[12,])) -# } -# -# alpha0 <- alpha0[-1,]; beta0 <- beta0[-1,]; beta1 <- beta1[-1,]; beta2 <- beta2[-1,]; beta3 <- beta3[-1,]; beta4 <- beta4[-1,]; -# gamma0 <- gamma0[-1,]; omega0 <- omega0[-1,]; omega1 <- omega1[-1,]; omega2 <- omega2[-1,]; omega3 <- omega3[-1,] -# -# out <- list(alpha0, beta0, beta1, beta2, beta3, beta4, gamma0, omega0, omega1, omega2, omega3) -# -# for(i in 1:11){ -# data.frame(species = names, out[[i]]) %>% -# ggplot(.) + -# geom_errorbar(aes(x = species, ymin = mean, ymax = mean), -# width = 0.275) + -# geom_errorbar(aes(x = species, ymin = X2.5., ymax = X97.5.), -# width = 0, size = 1.25) + -# geom_errorbar(aes(x = species, ymin = X25., ymax = X75.), -# width = 0, size = 3) + -# geom_hline(yintercept = 0, alpha = 0.75) + -# theme_few() + -# theme(panel.background = element_rect(fill = "transparent", color = NA), -# axis.text.x = element_text(angle = 50, hjust = 0.5, vjust = 0.5), -# legend.position = "none") + -# labs(y = expression(), x = expression()) -# -# } + legend.title = element_blank(), + legend.position = "bottom", + strip.text = element_text(size = 12)) + + labs(y = "Relative bias", x = expression())) + +loc <- Figure1$layout[grep("panel", Figure1$layout$name),1:4] +loc <- loc %>% arrange(t,l) + +for(i in 1:dim(loc)[1]){ + Figure1 <- gtable::gtable_add_grob(Figure1, grid::textGrob(LETTERS[i], x = unit(0.1, "in"), y = unit(1.5, "in"), + gp = grid::gpar(fontsize = 14, + fontface = 2)), + t = loc$t[i], l = loc$l[i]) +} + + +tiff(file = "~/HMSNO/PostAnalysis/Figure1.tiff", res = 600, width = 6, height = 8, units = "in") +grid::grid.draw(Figure1) +dev.off() \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..4505b64 --- /dev/null +++ b/README.md @@ -0,0 +1,29 @@ +# [Quantifying the conservation status and abundance trends of wildlife communities with detection-nondetection data]() + +### [Matthew T. Farr](https://farrmt.github.io/), Timothy O'Brien, Charles B. Yackulic, & [Elise F. Zipkin](https://zipkinlab.org/) + +### Conservation Biology + +### Code/Data DOI: [![DOI]()]() + +### Please contact the first author for questions about the code or data: Matthew T. Farr (matthewtfarr@gmail.com) +__________________________________________________________________________________________________________________________________________ + +## Abstract: +Effective conservation requires understanding species’ abundance patterns and demographic rates across space and time. Ideally, such knowledge should be available for whole communities, as variation in species’ dynamics can elucidate factors leading to biodiversity losses. However, collecting data to simultaneously estimate abundance and demographic rates is often prohibitively time-intensive and expensive for communities of species. We developed a “multi-species dynamic N-occupancy model” to estimate unbiased, community-wide relative abundance and demographic rates. Our model uses detection-nondetection data (e.g., repeated presence-absence surveys) to estimate both species- and community-level parameters as well as the effects of environmental factors. We conducted a simulation study that validated our modeling framework, demonstrating how and when such an approach can be valuable. Using data from a network of camera traps across tropical equatorial Africa, we then used our model to evaluate the statuses and trends of a forest-dwelling antelope community. We estimated relative abundance, rates of recruitment (i.e., reproduction and immigration), and apparent survival probabilities for each species’ local population. Our analysis indicated that the antelope community was fairly stable in this region (although 17% of populations [species-park combinations] declined over the study period), with variation in apparent survival linked more closely to differences among national parks rather than individual species’ life histories. The multi-species dynamic N-occupancy model requires only detection-nondetection data to evaluate the population dynamics of multiple sympatric species and can thus be a valuable tool for conservation efforts seeking to understand the reasons behind recent biodiversity loss. + +## Repository Directory + +### [DataAnalysis](./DataAnalysis): Contains code for modeling, analysis, and results for both the simulation and case studies +### [DataFormat](./DataFormat): Contains raw data, code to format raw data for analysis, and formatted data +### [PostAnalysis](./PostAnalysis): Contains code to create tables & figures +### [SupportingInformation](./SupportingInformation): Contains supporting information and code to generate supporting information +### [PublishedPDF](): PDF of published paper + +## Data +See the following subdirectories for data and metadata: [DataFromatting/RawData](./DataFormatting/RawData) + +## Code +See the following subdirectories for code and metadata: [DataAnalysis](./DataAnalysis), [DataFormatting](./DataFormatting), [PostAnalysis](./PostAnalysis), [SupportingInformation](./SupportingInformation) + + diff --git a/SS.txt b/SS.txt deleted file mode 100644 index 15d440c..0000000 --- a/SS.txt +++ /dev/null @@ -1,25 +0,0 @@ -model{ - -for(j in 1:nsites){ -gamma[j] ~ dnorm(0, tau.g) -} - -tau.g ~ dgamma(0.1, 0.1) - -alpha0 ~ dnorm(0, 0.1) -beta0 ~ dnorm(0, 0.1) - -for(t in 1:nyrs){ -for(j in 1:nsites){ -logit(r[t,j]) <- alpha0 -N[t,j] ~ dpois(lambda[t,j]) -log(lambda[t,j]) <- beta0 + gamma[j] -}#end j -}#end t - -for(k in 1:nobs){ -y[k] ~ dbern(p[k]) -p[k] <- 1 - pow((1 - r[yr[k], site[k]]), N[yr[k], site[k]]) -} - -} \ No newline at end of file diff --git a/SSdyn.txt b/SSdyn.txt deleted file mode 100644 index f0713ca..0000000 --- a/SSdyn.txt +++ /dev/null @@ -1,35 +0,0 @@ -model{ - -log(lambda) <- beta0 -beta0 ~ dnorm(0, 0.1) - -alpha0 ~ dnorm(0, 0.1) -omega0 ~ dnorm(0, 0.1) -gamma0 ~ dnorm(0, 0.1) - -for(j in 1:nsites){ - -logit(r[1,j]) <- alpha0 -N[1,j] ~ dpois(lambda) - -for(t in 2:nyrs){ - -logit(r[t,j]) <- alpha0 -logit(omega[t-1,j]) <- omega0 -S[t-1,j] ~ dbin(omega[t-1,j], N[t-1,j]) -G[t-1,j] ~ dpois(gamma[t-1]) -N[t,j] <- S[t-1,j] + G[t-1,j] - -}#end t -}#end j - -for(t in 2:nyrs){ -log(gamma[t-1]) <- gamma0 -} - -for(k in 1:nobs){ -y[k] ~ dbern(p[k]) -p[k] <- 1 - pow((1 - r[yr[k], site[k]]), N[yr[k], site[k]]) -} - -} \ No newline at end of file diff --git a/SSstatic.txt b/SSstatic.txt deleted file mode 100644 index e110286..0000000 --- a/SSstatic.txt +++ /dev/null @@ -1,22 +0,0 @@ -model{ - -for(j in 1:nsite){ -gamma[j] ~ dnorm(0, tau.g) -} - -tau.g ~ dgamma(0.1, 0.1) - -alpha0 ~ dnorm(0, 0.1) -beta0 ~ dnorm(0, 0.1) - -for(t in 1:nyrs){ -for(j in 1:nsite){ -y[t,j] ~ dbin(v[t,j], p[t,j]) -p[t,j] <- 1 - pow((1 - r[t,j]), N[t,j]) -logit(r[t,j]) <- alpha0 -N[t,j] <- dpois(lambda[t,j]) -log(lambda[t,j]) <- beta0 + gamma[j] -}#end j -}#end t - -} \ No newline at end of file