## Code to run the two species models ## Adapted from Brodie et al. 2018, Biotropica ## Written seperately for each species pair # The data for species 1 are in data0 # The data for species 2 are in data1 library(R2jags) library(readxl) # Set working directory to where your "spp_occ75" and "site_covs" Excel files are. #--------------------------------- BTM/GEN ------------------------------------------------------------------------------------ # rm(list=ls(all=TRUE)) # Note: have to delete NAs from excel sheet for the counts to be read as integers data0 <- read_xlsx("spp_occ75.xlsx", sheet = "BTM") data1 <- read_xlsx("spp_occ75.xlsx", sheet = "GEN") # readxl makes it into a tibble, which messes up code below data0 <- as.data.frame(data0) data1 <- as.data.frame(data1) rownames(data0) = data0[,1] # makes the station names the row names data0 = data0[,-1] # removes the "station" column data0 = as.matrix(data0) # creates a matrix rownames(data1) = data1[,1] data1 = data1[,-1] data1 = as.matrix(data1) ################## COVARIATES ####################################### PSIcovars <- read_xlsx("site_covs.xlsx", sheet = "PsiCov") PSIcovars <- as.data.frame(PSIcovars) DETcovars <- read_xlsx("site_covs.xlsx", sheet="DetCov") DETcovars <- as.data.frame(DETcovars) dates <- read_xlsx("site_covs.xlsx", sheet = "dates") dates <- as.data.frame(dates) # same thing, make station names the row names rownames(DETcovars)=DETcovars[,1] DETcovars=DETcovars[,-1] # and delete the station column rownames(PSIcovars)=PSIcovars[,1] PSIcovars=PSIcovars[,-1] rownames(dates)=dates[,1] dates=dates[,-1] dates=as.matrix(dates) #standardize continuous data columns PSIcovars<-as.data.frame(cbind(scale(PSIcovars[1:4], center = TRUE, scale = TRUE))) DETcovars<-cbind(DETcovars[1:2], cbind(scale(DETcovars[3], center = TRUE, scale = TRUE))) # ^check column names #standardize julian date matrix mdates=mean(as.vector(dates),na.rm=T) sddates = sd(as.vector(dates),na.rm=T) dates1 <- (dates-mdates)/sddates for (i in 1:dim(dates1)[1]) { a <- which(dates[i,]>0) dates1[i,-a]=0} #detection paired = DETcovars$paired trail = DETcovars$ontrail slope <- DETcovars$slope dates = dates1 #occupancy water = PSIcovars$dist_allriv human = PSIcovars$human_rate ndvi500 = PSIcovars$ndvi_500buff settle = PSIcovars$dist_settle #Number of sites nsite = dim(data0)[1] #Number of reps nrep = dim(data0)[2] ################################# BAYESIAN MODELING ############################# sp.data = list(y0=data0, y1=data1, nsite=nsite, nrep=dim(data0)[2], trail=trail,water=water,ndvi500=ndvi500,paired=paired,settle=settle, human=human,dates=dates,slope=slope) # Specify the parameters to be monitored sp.params = c('a0','a1','a2','a3','a4','a5','b0','b1','b2','b3','b4','N0','N1','lambda0','lambda1') # Specify the initial values sp.inits = function() { list(a0=rnorm(2), a1=rnorm(2), a2=rnorm(2), a3=rnorm(2), a4=rnorm(2), a5=rnorm(1), b0=rnorm(2), b1=rnorm(2), b2=rnorm(2), b3=rnorm(2), b4=rnorm(2), N0 = as.vector(apply(data0,1,max, na.rm=T) + 5), N1 = as.vector(apply(data1,1,max, na.rm=T) + 5) ) } # MCMC settings (choose one) ni <- 300; nt <- 2; nb <- 100; nc <- 3 # quick test ni <- 20000; nt <- 10; nb <- 10000; nc <- 3 # examining parameter values ni <- 100000; nt <- 20; nb <- 50000; nc <- 3 # thorough run # BTM/GEN Sys.time() out_btmgen <- jags(data=sp.data, inits=sp.inits, parameters.to.save = sp.params, model.file = "TE_2spp.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb) save.image(file="btmgen.RData") # To look at parameter values, mean, sd, and CIs # This output used for Table 2, Figures 5 & 6 out_btmgen_sum <- out_btmgen$BUGSoutput$summary Sys.time() #--------------------------------- BTM/CIV ------------------------------------------------------------------------------------ #--------------------------------------------------------------------------------------------------------------------- #rm(list=ls(all=TRUE)) #have to delete NAs from excel sheet for the counts to be read as integers data0 <- read_xlsx("spp_occ75.xlsx", sheet = "BTM") data1 <- read_xlsx("spp_occ75.xlsx", sheet = "CIV") data0 <- as.data.frame(data0) #readxl makes it into a tibble, which messes up code below data1 <- as.data.frame(data1) rownames(data0) = data0[,1] #makes the station names the row names data0 = data0[,-1] #removes the "station" column data0 = as.matrix(data0) #creates a matrix rownames(data1) = data1[,1] data1 = data1[,-1] data1 = as.matrix(data1) ################## COVARIATES ####################################### PSIcovars <- read_xlsx("site_covs.xlsx", sheet = "PsiCov") PSIcovars <- as.data.frame(PSIcovars) DETcovars <- read_xlsx("site_covs.xlsx", sheet="DetCov") DETcovars <- as.data.frame(DETcovars) dates <- read_xlsx("site_covs.xlsx", sheet = "dates") dates <- as.data.frame(dates) # same thing, make station names the row names rownames(DETcovars)=DETcovars[,1] DETcovars=DETcovars[,-1] # and delete the station column rownames(PSIcovars)=PSIcovars[,1] PSIcovars=PSIcovars[,-1] rownames(dates)=dates[,1] dates=dates[,-1] dates=as.matrix(dates) #standardize continuous data columns PSIcovars<-as.data.frame(cbind(scale(PSIcovars[1:4], center = TRUE, scale = TRUE))) DETcovars<-cbind(DETcovars[1:2], cbind(scale(DETcovars[3], center = TRUE, scale = TRUE))) # ^check column names #standardize julian date matrix mdates=mean(as.vector(dates),na.rm=T) sddates = sd(as.vector(dates),na.rm=T) dates1 <- (dates-mdates)/sddates for (i in 1:dim(dates1)[1]) { a <- which(dates[i,]>0) dates1[i,-a]=0} #detection paired = DETcovars$paired trail = DETcovars$ontrail slope <- DETcovars$slope dates = dates1 #occupancy water = PSIcovars$dist_allriv human = PSIcovars$human_rate ndvi500 = PSIcovars$ndvi_500buff settle = PSIcovars$dist_settle #Number of sites nsite = dim(data0)[1] #Number of reps nrep = dim(data0)[2] ################################# BAYESIAN MODELING ############################# sp.data = list(y0=data0, y1=data1, nsite=nsite, nrep=dim(data0)[2], trail=trail,water=water,ndvi500=ndvi500,paired=paired,settle=settle, human=human,dates=dates,slope=slope) # Specify the parameters to be monitored sp.params = c('a0','a1','a2','a3','a4','a5','b0','b1','b2','b3','b4','N0','N1','lambda0','lambda1') # Specify the initial values sp.inits = function() { list(a0=rnorm(2), a1=rnorm(2), a2=rnorm(2), a3=rnorm(2), a4=rnorm(2), a5=rnorm(1), b0=rnorm(2), b1=rnorm(2), b2=rnorm(2), b3=rnorm(2), b4=rnorm(2), N0 = as.vector(apply(data0,1,max, na.rm=T) + 5), N1 = as.vector(apply(data1,1,max, na.rm=T) + 5) ) } # MCMC settings (choose one) ni <- 300; nt <- 2; nb <- 100; nc <- 3 # quick test ni <- 20000; nt <- 10; nb <- 10000; nc <- 3 # examining parameter values ni <- 100000; nt <- 20; nb <- 50000; nc <- 3 # thorough run # BTM/CIV Sys.time() out_btmciv <- jags(data=sp.data, inits=sp.inits, parameters.to.save = sp.params, model.file = "TE_2spp.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb) save.image(file="btmciv.RData") out_btmciv_sum <- out_btmciv$BUGSoutput$summary Sys.time() #--------------------------------- GEN/CIV ------------------------------------------------------------------------------------ #--------------------------------------------------------------------------------------------------------------------- #rm(list=ls(all=TRUE)) #have to delete NAs from excel sheet for the counts to be read as integers data0 <- read_xlsx("spp_occ75.xlsx", sheet = "GEN") data1 <- read_xlsx("spp_occ75.xlsx", sheet = "CIV") data0 <- as.data.frame(data0) #readxl makes it into a tibble, which messes up code below data1 <- as.data.frame(data1) rownames(data0)=data0[,1] #makes the station names the row names data0=data0[,-1] #removes the "station" column data0=as.matrix(data0) #creates a matrix rownames(data1)=data1[,1] data1=data1[,-1] data1=as.matrix(data1) ################## COVARIATES ####################################### PSIcovars <- read_xlsx("site_covs.xlsx", sheet = "PsiCov") PSIcovars <- as.data.frame(PSIcovars) DETcovars <- read_xlsx("site_covs.xlsx", sheet="DetCov") DETcovars <- as.data.frame(DETcovars) dates <- read_xlsx("site_covs.xlsx", sheet = "dates") dates <- as.data.frame(dates) # same thing, make station names the row names rownames(DETcovars)=DETcovars[,1] DETcovars=DETcovars[,-1] # and delete the station column rownames(PSIcovars)=PSIcovars[,1] PSIcovars=PSIcovars[,-1] rownames(dates)=dates[,1] dates=dates[,-1] dates=as.matrix(dates) #standardize continuous data columns PSIcovars<-as.data.frame(cbind(scale(PSIcovars[1:4], center = TRUE, scale = TRUE))) DETcovars<-cbind(DETcovars[1:2], cbind(scale(DETcovars[3], center = TRUE, scale = TRUE))) # ^check column names #standardize julian date matrix mdates=mean(as.vector(dates),na.rm=T) sddates = sd(as.vector(dates),na.rm=T) dates1 <- (dates-mdates)/sddates for (i in 1:dim(dates1)[1]) { a <- which(dates[i,]>0) dates1[i,-a]=0} #detection paired = DETcovars$paired trail = DETcovars$ontrail slope <- DETcovars$slope dates = dates1 #occupancy water = PSIcovars$dist_allriv human = PSIcovars$human_rate ndvi500 = PSIcovars$ndvi_500buff settle = PSIcovars$dist_settle #Number of sites nsite = dim(data0)[1] #Number of reps nrep = dim(data0)[2] ################################# BAYESIAN MODELING ############################# sp.data = list(y0=data0, y1=data1, nsite=nsite, nrep=dim(data0)[2], trail=trail,water=water,ndvi500=ndvi500,paired=paired,settle=settle, human=human,dates=dates,slope=slope) # Specify the parameters to be monitored sp.params = c('a0','a1','a2','a3','a4','a5','b0','b1','b2','b3','b4','N0','N1','lambda0','lambda1') # Specify the initial values sp.inits = function() { list(a0=rnorm(2), a1=rnorm(2), a2=rnorm(2), a3=rnorm(2), a4=rnorm(2), a5=rnorm(1), b0=rnorm(2), b1=rnorm(2), b2=rnorm(2), b3=rnorm(2), b4=rnorm(2), N0 = as.vector(apply(data0,1,max, na.rm=T) + 5), N1 = as.vector(apply(data1,1,max, na.rm=T) + 5) ) } # MCMC settings (choose one) ni <- 300; nt <- 2; nb <- 100; nc <- 3 # quick test ni <- 20000; nt <- 10; nb <- 10000; nc <- 3 # examining parameter values ni <- 100000; nt <- 20; nb <- 50000; nc <- 3 # thorough run # GEN/CIV Sys.time() out_genciv <- jags(data=sp.data, inits=sp.inits, parameters.to.save = sp.params, model.file = "TE_2spp.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb) out_genciv_sum <- out_genciv$BUGSoutput$summary save.image(file="genciv.RData") Sys.time()