## Change to match your computer setwd("~/Desktop/R/RBirdCorrelation") ## Required package library("unmarked") library("MuMIn") ## Get survey data and covariates data<-read.csv("BrettFurnas9SurveysPerSite.csv",as.is=T) data$Occupancy[382]<-1# Was a 3, is this a typo? data$Canopy.Cover<- data$Canopy.Cover/100 survey<-get.survey(data) ## Create unmarked frame for a single species species: umf.comy<-get.umf(survey,"Common Myna") umf.cowa<-get.umf(survey,"Common Waxbill") umf.ggfd<-get.umf(survey,"Grey-green Fruit Dove") umf.moki<-get.umf(survey,"Mo'orean Kingfisher") umf.rbfi<-get.umf(survey,"Red-browed Firetail") umf.rvbu<-get.umf(survey,"Red-vented Bulbul") umf.rjfo<-get.umf(survey,"Red Jungle Fowl") umf.silv<-get.umf(survey,"Silvereye") umf.zedo<-get.umf(survey,"Zebra Dove") ##FOR MO'OREAN KINGFISHER model.moki.D1<-occu(~1 ~ 1, data=umf.moki) model.moki.D2<-occu(~before + after ~ 1, data=umf.moki) model.moki.D3<-occu(~ canopy ~ 1, data=umf.moki) model.moki.D4<-occu(~ before + after + canopy ~ 1, data=umf.moki) model.moki.D5<-occu(~ habitat ~ 1, data=umf.moki) model.moki.D6<-occu(~ before + after + habitat~ 1, data=umf.moki) model.moki.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.moki) models.set.mokiD<- model.sel('Intercept Only'= model.moki.D1, 'before + after'= model.moki.D2, 'canopy'= model.moki.D3, 'before + after + canopy'= model.moki.D4) models.set.mokiD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.moki.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Mo'orean Kingfisher", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ### For Common Waxbill model.cowa.D1<-occu(~1 ~ 1, data=umf.cowa) model.cowa.D2<-occu(~before + after ~ 1, data=umf.cowa) model.cowa.D3<-occu(~ canopy ~ 1, data=umf.cowa) model.cowa.D4<-occu(~ before + after + canopy ~ 1, data=umf.cowa) model.cowa.D5<-occu(~ habitat ~ 1, data=umf.cowa) model.cowa.D6<-occu(~ before + after + habitat~ 1, data=umf.cowa) model.cowa.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.cowa) models.set.cowaD<- model.sel('Intercept Only'= model.cowa.D1, 'before + after'= model.cowa.D2, 'canopy'= model.cowa.D3, 'before + after + canopy'= model.cowa.D4, 'habitat' = model.cowa.D5, 'before + after +habitat'=model.cowa.D6, 'before + after + habitat + canopy'=model.cowa.D7) models.set.cowaD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.cowa.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Common Waxbill", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ###FOR COMMON MYNA model.comy.D1<-occu(~1 ~ 1, data=umf.comy) model.comy.D2<-occu(~before + after ~ 1, data=umf.comy) model.comy.D3<-occu(~ canopy ~ 1, data=umf.comy) model.comy.D4<-occu(~ before + after + canopy ~ 1, data=umf.comy) model.comy.D5<-occu(~ habitat ~ 1, data=umf.comy) model.comy.D6<-occu(~ before + after + habitat~ 1, data=umf.comy) model.comy.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.comy) models.set.comyD<- model.sel('Intercept Only'= model.comy.D1, 'before + after'= model.comy.D2, 'canopy'= model.comy.D3, 'before + after + canopy'= model.comy.D4, 'habitat' = model.comy.D5, 'before + after +habitat'=model.comy.D6, 'before + after + habitat + canopy'=model.comy.D7) models.set.comyD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.comy.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Common Myna", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ### FOR Grey-green Fruit Dove model.ggfd.D1<-occu(~1 ~ 1, data=umf.ggfd) model.ggfd.D2<-occu(~before + after ~ 1, data=umf.ggfd) model.ggfd.D3<-occu(~ canopy ~ 1, data=umf.ggfd) model.ggfd.D4<-occu(~ before + after + canopy ~ 1, data=umf.ggfd) model.ggfd.D5<-occu(~ habitat ~ 1, data=umf.ggfd) model.ggfd.D6<-occu(~ before + after + habitat~ 1, data=umf.ggfd) model.ggfd.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.ggfd) models.set.ggfdD<- model.sel('Intercept Only'= model.ggfd.D1, 'before + after'= model.ggfd.D2, 'canopy'= model.ggfd.D3, 'before + after + canopy'= model.ggfd.D4, 'habitat' = model.ggfd.D5, 'before + after +habitat'=model.ggfd.D6, 'before + after + habitat + canopy'=model.ggfd.D7) models.set.ggfdD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.ggfd.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Grey-green Fruit Dove", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ### FOR Red-vented Bulbul model.rvbu.D1<-occu(~1 ~ 1, data=umf.rvbu) model.rvbu.D2<-occu(~before + after ~ 1, data=umf.rvbu) model.rvbu.D3<-occu(~ canopy ~ 1, data=umf.rvbu) model.rvbu.D4<-occu(~ before + after + canopy ~ 1, data=umf.rvbu) model.rvbu.D5<-occu(~ habitat ~ 1, data=umf.rvbu) model.rvbu.D6<-occu(~ before + after + habitat~ 1, data=umf.rvbu) model.rvbu.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.rvbu) models.set.rvbuD<- model.sel('Intercept Only'= model.rvbu.D1, 'before + after'= model.rvbu.D2, 'canopy'= model.rvbu.D3, 'before + after + canopy'= model.rvbu.D4, 'habitat' = model.rvbu.D5, 'before + after +habitat'=model.rvbu.D6) 'before + after + habitat + canopy'=model.rvbu.D7) models.set.rvbuD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.rvbu.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main= "Detectibility of Red-vented Bulbul", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ### FOR Red-browed Firetail model.rbfi.D1<-occu(~1 ~ 1, data=umf.rbfi) model.rbfi.D2<-occu(~before + after ~ 1, data=umf.rbfi) model.rbfi.D3<-occu(~ canopy ~ 1, data=umf.rbfi) model.rbfi.D4<-occu(~ before + after + canopy ~ 1, data=umf.rbfi) model.rbfi.D5<-occu(~ habitat ~ 1, data=umf.rbfi) model.rbfi.D6<-occu(~ before + after + habitat~ 1, data=umf.rbfi) model.rbfi.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.rbfi) models.set.rbfiD<- model.sel('Intercept Only'= model.rbfi.D1, 'before + after'= model.rbfi.D2, 'canopy'= model.rbfi.D3, 'before + after + canopy'= model.rbfi.D4) models.set.rbfiD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.rbfi.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Red-browed Firetail", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ## FOR Red Jungle Fowl model.rjfo.D1<-occu(~1 ~ 1, data=umf.rjfo) model.rjfo.D2<-occu(~before + after ~ 1, data=umf.rjfo) model.rjfo.D3<-occu(~ canopy ~ 1, data=umf.rjfo) model.rjfo.D4<-occu(~ before + after + canopy ~ 1, data=umf.rjfo) model.rjfo.D5<-occu(~ habitat ~ 1, data=umf.rjfo) model.rjfo.D6<-occu(~ before + after + habitat~ 1, data=umf.rjfo) model.rjfo.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.rjfo) models.set.rjfoD<- model.sel('Intercept Only'= model.rjfo.D1, 'before + after'= model.rjfo.D2, 'canopy'= model.rjfo.D3, 'before + after + canopy'= model.rjfo.D4) models.set.rjfoD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.rjfo.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Red Jungle Fowl", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ## FOR Silvereye model.silv.D1<-occu(~1 ~ 1, data=umf.silv) model.silv.D2<-occu(~before + after ~ 1, data=umf.silv) model.silv.D3<-occu(~ canopy ~ 1, data=umf.silv) model.silv.D4<-occu(~ before + after + canopy ~ 1, data=umf.silv) model.silv.D5<-occu(~ habitat ~ 1, data=umf.silv) model.silv.D6<-occu(~ before + after + habitat~ 1, data=umf.silv) model.silv.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.silv) models.set.silvD<- model.sel('Intercept Only'= model.silv.D1, 'before + after'= model.silv.D2, 'canopy'= model.silv.D3, 'before + after + canopy'= model.silv.D4) models.set.silvD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.silv.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Silvereye", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2, col = "blue") points(x.fit, y.fit$upper, type = "l", lty=2, col = "blue") ## FOR zebra dove model.zedo.D1<-occu(~1 ~ 1, data=umf.zedo) model.zedo.D2<-occu(~before + after ~ 1, data=umf.zedo) model.zedo.D3<-occu(~ canopy ~ 1, data=umf.zedo) model.zedo.D4<-occu(~ before + after + canopy ~ 1, data=umf.zedo) model.zedo.D5<-occu(~ habitat ~ 1, data=umf.zedo) model.zedo.D6<-occu(~ before + after + habitat~ 1, data=umf.zedo) model.zedo.D7<-occu(~ before + after + habitat + canopy~ 1, data=umf.zedo) models.set.zedoD<- model.sel('Intercept Only'= model.zedo.D1, 'before + after'= model.zedo.D2, 'canopy'= model.zedo.D3, 'before + after + canopy'= model.zedo.D4) models.set.zedoD # Here are the results x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<- predict(model.zedo.D3, newdata= data.frame(canopy=x.fit.c), type = 'det') plot(x.fit, y.fit$Predicted, type="l", ylim=c(0,1), ylab="Detection Probability", xlab="Canopy Cover", lwd = 2, main = "Detectibility of Zebra Dove", col = "blue") points(x.fit, y.fit$lower, type = "l", lty=2) points(x.fit, y.fit$upper, type = "l", lty=2) ######## ######## ######## RUN THESE FIRST #### FUNCTIONS get.survey<-function(data){ data$Survey2<-NA data$Survey2[data$Day==1 &data$Time=="Before"]<-1 data$Survey2[data$Day==1 &data$Time=="During"]<-2 data$Survey2[data$Day==1 &data$Time=="After"]<-3 data$Survey2[data$Day==2 &data$Time=="Before"]<-4 data$Survey2[data$Day==2 &data$Time=="During"]<-5 data$Survey2[data$Day==2 &data$Time=="After"]<-6 data$Survey2[data$Day==3 &data$Time=="Before"]<-7 data$Survey2[data$Day==3 &data$Time=="During"]<-8 data$Survey2[data$Day==3 &data$Time=="After"]<-9 species<-sort(unique(data$Species)) n.species<-length(species) sites<-sort(unique(data$Site)) n.sites<-length(sites) surveys<-1:9 n.surveys<-length(surveys) # detection histories y<-array(NA,dim=c(n.sites,n.surveys,n.species)) dimnames(y)<-list(sites,surveys,species) for(i in 1:n.species){ for(j in 1:n.sites){ for(k in 1:n.surveys){ temp.1<-subset(data,Species==species[i] & Site==sites[j] & Survey2==surveys[k]) y[j,k,i]<-temp.1$Occupancy } # k } # j } # i # site covariates habitat<-elevation<-canopy<-rep(NA,n.sites) for(j in 1:n.sites){ temp.2<-subset(data,Site==sites[j]) habitat[j]<-unique(temp.2$Habitat) elevation[j]<-unique(temp.2$Elevation) canopy[j]<-unique(temp.2$Canopy.Cover) # survey covariates (could use survey date or temperature) } # j # other survey covariates time<-array(c(1,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1),dim=c(9,3)) dimnames(time)<-list(1:9,c("B","D","A")) # standardize continuous covariates standards<-array(NA,dim=c(2,2)) dimnames(standards)<-list(c("elevation","canopy"),c("mn","sd")) standards["elevation","mn"]<-mean(elevation) standards["elevation","sd"]<-sd(elevation) standards["canopy","mn"]<-mean(canopy) standards["canopy","sd"]<-sd(canopy) elevation<-(elevation-standards["elevation","mn"])/standards["elevation","sd"] canopy<-(canopy-standards["canopy","mn"])/standards["canopy","sd"] # output out<-list(y,habitat,elevation,canopy,time,standards) names(out)<-c("y","habitat","elevation","canopy","time","standards") return(out) } #### Get unmarked frame for occupancy modeling get.umf<-function(survey,spp){ y<-survey[["y"]][,,spp] umf<- unmarkedFrameOccu(y=y, siteCovs = data.frame(habitat=factor(survey[["habitat"]]), elevation=survey[["elevation"]], canopy=survey[["canopy"]]), obsCovs = list( before=t(array(survey[["time"]][,1],dim=c(9,dim(y)[[1]]))), during=t(array(survey[["time"]][,2],dim=c(9,dim(y)[[1]]))), after=t(array(survey[["time"]][,3],dim=c(9,dim(y)[[1]]))) ) ) return(umf) } ## WHAT YOU USE FOR X.FIT AND Y.FIT x.fit<-seq(1,100,1) x.fit.c<-(x.fit-survey[["standards"]]["canopy","mn"])/survey[["standards"]]["canopy","sd"] y.fit<-data.frame(canopy=x.fit.c)