library(RPresence) library(gplots) library(sp) library(rgdal) # R wrapper around GDAL/OGR library(ggplot2) # for general plotting library(ggmap) # for fortifying shapefiles getwd() # set working directory, where files are stored setwd("C:/Users/koust/Dropbox (Snow Leopard Trust)/Papers in Prep/Sanjay Gubbi") ########################################################### # Load saved models to avoid the need to run them again ########################################################### load("Sanjaymodels2.RData") #including new analyses using theta! # Load leopard data used to run the occupancy analysis load("leopData.RData") #Load saved data ########################################################### # load data that uses 1km transect segments ########################################################### leop.hist10<-read.csv("LeopardHistory_10.csv", header = TRUE) leop.site.cov<-read.csv("SiteCov_sanjay.csv", header = TRUE) #long leop.survey.cov<-read.csv("SamplingCovs_10.csv") # standardise site covariates to z-scores leop.site.cov$std.habitat<-(leop.site.cov$Natural.Habitat-mean(leop.site.cov$Natural.Habitat))/ sd(leop.site.cov$Natural.Habitat) leop.site.cov$std.irrigated<-(leop.site.cov$Irrigated-mean(leop.site.cov$Irrigated))/ sd(leop.site.cov$Irrigated) leop.site.cov$std.rainfed<-(leop.site.cov$Rainfed-mean(leop.site.cov$Rainfed))/ sd(leop.site.cov$Rainfed) # check headers and names of standardized site covariates names(leop.site.cov) summary(leop.site.cov) # Standardize survey covariates leop.survey.cov$Dogs.std<-(leop.survey.cov$Dogs_Prop-mean(leop.survey.cov$Dogs_Prop))/ sd(leop.survey.cov$Dogs_Prop) leop.survey.cov$Mining.std<-(leop.survey.cov$Mining_prop-mean(leop.survey.cov$Mining_prop))/ sd(leop.survey.cov$Mining_prop) leop.survey.cov$DryDec.std<-(leop.survey.cov$DryDeciduous_prop-mean(leop.survey.cov$DryDeciduous_prop))/ sd(leop.survey.cov$DryDeciduous_prop) leop.survey.cov$Eucaly.std<-(leop.survey.cov$Eucalyptus_prop-mean(leop.survey.cov$Eucalyptus_prop))/ sd(leop.survey.cov$Eucalyptus_prop) leop.survey.cov$Grassland.std<-(leop.survey.cov$Grassland_prop-mean(leop.survey.cov$Grassland_prop))/ sd(leop.survey.cov$Grassland_prop) leop.survey.cov$MoistDec.std<-(leop.survey.cov$Moist_dec_prop-mean(leop.survey.cov$Moist_dec_prop))/ sd(leop.survey.cov$Moist_dec_prop) leop.survey.cov$Orchard.std<-(leop.survey.cov$Orchard_prop-mean(leop.survey.cov$Orchard_prop))/ sd(leop.survey.cov$Orchard_prop) leop.survey.cov$Scrub.std<-(leop.survey.cov$Scrub_prop-mean(leop.survey.cov$Scrub_prop))/ sd(leop.survey.cov$Scrub_prop) leop.survey.cov$Rocky.std<-(leop.survey.cov$Rocky_prop-mean(leop.survey.cov$Rocky_prop))/ sd(leop.survey.cov$Rocky_prop) leop.survey.cov$Irrigated.std<-(leop.survey.cov$Irrigated_prop-mean(leop.survey.cov$Irrigated_prop))/ sd(leop.survey.cov$Irrigated_prop) leop.survey.cov$Natural.std<-(leop.survey.cov$Natural_prop-mean(leop.survey.cov$Natural_prop))/ sd(leop.survey.cov$Natural_prop) leop.survey.cov$Raifed.std<-(leop.survey.cov$Rainfed_prop-mean(leop.survey.cov$Rainfed_prop))/ sd(leop.survey.cov$Rainfed_prop) leop.survey.cov$Livestock.std<-(leop.survey.cov$Livestock_Prop-mean(leop.survey.cov$Livestock_Prop))/ sd(leop.survey.cov$Livestock_Prop) # Check headers and column names of the standardized survey covariates names(leop.site.cov) head(leop.site.cov) names(leop.survey.cov) head(leop.survey.cov) # create leop.data pao data object using leopard detection history, and site & survey covariates leop.data<-createPao(data= leop.hist10, unitcov=leop.site.cov, survcov = leop.survey.cov) save(leop.data,file="LeopData.RData") #Save data for future use without having to repeat the above steps # fit psi(.)p(.) or null model leop.null<-occMod(model = list(psi~1, p~1), data = leop.data, type = "so") summary(leop.null) #shows AIC value and number of parameters of model #fit null model with spatial correlation in detection leop.null.theta<-occMod(model = list(psi~1, p~1, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") summary(leop.null.theta) #shows AIC value and number of parameters of model # review estimated probabilities psi and p unique(fitted(leop.null,"psi")) # prints just the unique values unique(fitted(leop.null,"p")) # prints just the unique values unique(fitted(leop.null.theta,"psi")) # prints just the unique values unique(fitted(leop.null.theta,"p")) # prints just the unique values # review coefficients for psi and p coef(leop.null,"psi") #coefficients for occupancy coef(leop.null,"p") #coefficients for detection coef(leop.null.theta,"psi") #regression coefficients for occupancy coef(leop.null.theta,"p") #regression coefficients for detection # fit model with constant psi and p as a function of livestock leop2 <- occMod(model = list(psi~1,p~Livestock.std), data = leop.data, type = "so") #fit model with spatial correlation in detection and livestock as survey covariate leop2.theta<-occMod(model = list(psi~1, p~Livestock.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") # review coefficients for psi and p coef(leop2,"psi") #coefficients for occupancy coef(leop2,"p") #coefficients for detection # review estimated probabilities unique(fitted(leop2,"psi")) # prints just the unique values unique(fitted(leop2.theta,"psi")) # prints just the unique values unique(fitted(leop2,"p")) # prints just the unique estimated values unique(fitted(leop2.theta,"p")) # prints just the unique estimated values # Compile and compare all models run so far models1<-list(leop.null, leop2, leop2.theta, leop.null.theta) results<-createAicTable(models1) # Put all models and their AIC values together summary(results) # Create comparable AIC table from the above list # Create models using various combinations of survey covariates and testing them against models # with spatial correlation in detections leop3 <- occMod(model = list(psi~1,p~Dogs.std), data = leop.data, type = "so") leop3.theta<-occMod(model=list(psi~1, p~Dogs.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop4 <- occMod(model = list(psi~1,p~Mining.std), data = leop.data, type = "so") leop4.theta<-occMod(model=list(psi~1, p~Mining.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop5<- occMod(model = list(psi~1,p~DryDec.std), data = leop.data, type = "so") leop5.theta<-occMod(model=list(psi~1, p~DryDec.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop6 <- occMod(model = list(psi~1,p~Eucaly.std), data = leop.data, type = "so") leop6.theta<-occMod(model=list(psi~1, p~Eucaly.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop7 <- occMod(model = list(psi~1,p~Grassland.std), data = leop.data, type = "so") leop7.theta<-occMod(model=list(psi~1, p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop8 <- occMod(model = list(psi~1,p~MoistDec.std), data = leop.data, type = "so") leop8.theta<-occMod(model=list(psi~1, p~MoistDec.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop9 <- occMod(model = list(psi~1,p~Orchard.std), data = leop.data, type = "so") leop9.theta<-occMod(model=list(psi~1, p~Orchard.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop10 <- occMod(model = list(psi~1,p~Scrub.std), data = leop.data, type = "so") leop10.theta <- occMod(model = list(psi~1,p~Scrub.std, theta~PRIME,th0pi~1), data = leop.data, type = "so.cd") leop11 <- occMod(model = list(psi~1,p~Rocky.std), data = leop.data, type = "so") leop11.theta <- occMod(model = list(psi~1,p~Rocky.std, theta~PRIME,th0pi~1), data = leop.data, type = "so.cd") leop12 <- occMod(model = list(psi~1,p~Irrigated.std), data = leop.data, type = "so") leop12.theta<-occMod(model=list(psi~1, p~Irrigated.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop13 <- occMod(model = list(psi~1,p~Natural.std), data = leop.data, type = "so") leop13.theta <- occMod(model = list(psi~1,p~Natural.std, theta~PRIME,th0pi~1), data = leop.data, type = "so.cd") leop14 <- occMod(model = list(psi~1,p~Raifed.std), data = leop.data, type = "so") leop14.theta <- occMod(model = list(psi~1,p~Raifed.std, theta~PRIME,th0pi~1), data = leop.data, type = "so.cd") leop15 <- occMod(model = list(psi~1,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop15.theta <- occMod(model = list(psi~1,p~MoistDec.std+Grassland.std+Irrigated.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop16 <- occMod(model = list(psi~1,p~MoistDec.std+Grassland.std+Irrigated.std+Raifed.std+Natural.std), data = leop.data, type = "so") leop16.theta <- occMod(model = list(psi~1,p~MoistDec.std+Grassland.std+Irrigated.std+Raifed.std+Natural.std, theta~PRIME,th0pi~1),data = leop.data, type = "so.cd") leop28.theta <- occMod(model = list(psi~1,p~Grassland.std+Rocky.std+MoistDec.std,theta~PRIME,th0pi~1), data = leop.data, type = "so.cd") leop29.theta <- occMod(model = list(psi~1,p~Grassland.std+Rocky.std,theta~PRIME,th0pi~1), data = leop.data, type = "so.cd") #Run Goodness of fit with most parameterized model by bootstrapping 1000 times Mod.gof1<-occMod(model=list(psi~std.habitat+std.irrigated+std.rainfed+Large.Prey+Small.Prey+MinesQuaries, p~Grassland.std, theta ~ PRIME, th0pi~1), data=leop.data,type="so.cd", modfitboot = 1000) Mod.gof1$gof #Check the value of c-hat # Run additional models with site covariates leop17 <- occMod(model = list(psi~std.habitat,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop17.theta <- occMod(model = list(psi~std.habitat, p~1, theta ~ PRIME, th0pi~1), data = leop.data, type = "so.cd") leop18 <- occMod(model = list(psi~std.irrigated,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop18.theta <- occMod(model = list(psi~std.irrigated,p~MoistDec.std+Grassland.std+Irrigated.std, theta~PRIME,th0pi~1) , data = leop.data, type = "so.cd") leop18g.theta <- occMod(model = list(psi~std.irrigated,p~Grassland.std, theta~PRIME,th0pi~1) , data = leop.data, type = "so.cd") leop19 <- occMod(model = list(psi~std.rainfed,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop19g.theta <- occMod(model = list(psi~std.rainfed,p~Grassland.std,theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop20 <- occMod(model = list(psi~Category1,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop20g.theta <- occMod(model = list(psi~Category1,p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop21 <- occMod(model = list(psi~Large.Prey,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop21g.theta <- occMod(model = list(psi~Large.Prey,p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop22 <- occMod(model = list(psi~Small.Prey,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop22g.theta <- occMod(model = list(psi~Small.Prey,p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop23 <- occMod(model = list(psi~MinesQuaries,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop23g.theta <- occMod(model = list(psi~MinesQuaries,p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd") leop24 <- occMod(model = list(psi~std.habitat+std.rainfed+Large.Prey, p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop24g.theta <- occMod(model = list(psi~std.habitat+std.rainfed+Large.Prey, p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type = "so.cd", initvals = 2) leop25 <- occMod(model = list(psi~std.rainfed+Large.Prey,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop26 <- occMod(model = list(psi~std.habitat+Large.Prey,p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop30 <- occMod(model = list(psi~std.habitat+std.rainfed+Large.Prey+std.irrigated, p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop31 <- occMod(model = list(psi~Category1+std.rainfed+Large.Prey+std.irrigated, p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop32 <- occMod(model = list(psi~std.habitat+Large.Prey+std.irrigated, p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop33 <- occMod(model = list(psi~std.habitat+std.rainfed+std.irrigated, p~MoistDec.std+Grassland.std+Irrigated.std), data = leop.data, type = "so") leop34.theta<-occMod(model=list(psi~std.habitat+std.rainfed, p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type="so.cd") leop34.theta.1<-occMod(model=list(psi~std.habitat+std.rainfed, p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type="so.cd", initvals=2) leop35.theta<-occMod(model=list(psi~std.habitat, p~Grassland.std, theta~PRIME, th0pi~1), data = leop.data, type="so.cd", randinit =3, initvals=2) leop36.theta<-occMod(model=list(psi~std.habitat+std.rainfed, p~1, theta~PRIME, th0pi~1), data = leop.data, type="so.cd") # Put together all models run so far and compare using AIC table models3<-list(leop.null, leop.null.theta, leop2, leop2.theta, leop3, leop3.theta, leop4, leop4.theta, leop5, leop5.theta, leop6, leop6.theta, leop7, leop7.theta, leop8, leop8.theta, leop9, leop9.theta, leop10, leop10.theta, leop11, leop11.theta, leop12, leop12.theta, leop13, leop13.theta, leop14, leop14.theta, leop15, leop15.theta, leop16, leop16.theta, leop28.theta, leop29.theta,leop17, leop17.theta, leop18, leop18.theta, leop19, leop19g.theta, leop20, leop20g.theta, leop21, leop21g.theta, leop22, leop22g.theta, leop23, leop23g.theta, leop24, leop24g.theta, leop25, leop26, leop30, leop31, leop32, leop33, leop34.theta.1, leop35.theta) results3<-createAicTable(models3) summary(results3) # Save the AIC table as a csv file write.csv(summary(results3),file = "AICTable.csv") # List the top 3 models top.model3<-results3$models[[1]] summary(top.model3) top.model3.2<-results3$models[[2]] summary(top.model3.2) top.model3.3<-results3$models[[3]] summary(top.model3.3) # Extract coefficients of top 3 models coef(top.model3,"psi") coef(top.model3,"p") coef(top.model3,"theta") coef(top.model3.2,"psi") coef(top.model3.2,"p") coef(top.model3.2,"theta") coef(top.model3.3,"psi") coef(top.model3.3,"p") coef(top.model3,"theta") unique(fitted(top.model3,"psi")) # prints just the unique values of psi for the top model # Save the modeled ouptputs to be used later without having to go through the above steps every time save(leop.null, leop.null.theta, leop2, leop2.theta, leop3, leop3.theta, leop4, leop4.theta, leop5, leop5.theta, leop6, leop6.theta, leop7, leop7.theta, leop8, leop8.theta, leop9, leop9.theta, leop10, leop10.theta, leop11, leop11.theta, leop12, leop12.theta, leop13, leop13.theta, leop14, leop14.theta, leop15, leop15.theta, leop16, leop16.theta, leop28.theta, leop29.theta,leop17, leop17.theta, leop18, leop18.theta, leop19, leop19g.theta, leop20, leop20g.theta, leop21, leop21g.theta, leop22, leop22g.theta, leop23, leop23g.theta, leop24, leop24g.theta, leop25, leop26, leop30, leop31, leop32, leop33, leop34.theta, leop35.theta, leop34.theta.1, file="Sanjaymodels2.RData") #Saves model outputs # Estimate bias using naive estimate and model with constant psi NaiveBias<-(unique(fitted(leop7.theta, "psi")$est)-0.322097)/unique(fitted(leop7.theta, "psi")$est) NaiveBias ################################################################################################# ## Create species distribution map from top models ################################################################################################# #Read shapefile with all grid centers (point feature file) AllGrids <- readOGR("C:/Users/koust/Dropbox (Snow Leopard Trust)/Papers in Prep/Sanjay Gubbi/Shapes/All_Grids_pts.shp", stringsAsFactors = F) head(AllGrids) plot(AllGrids) # Convert the shapefile to a dataframe for use in ggplot2 AllGrids_df <- as(AllGrids, "data.frame") head(AllGrids_df) names(AllGrids) names(AllGrids_df) names(AllGrids_df)[27]="x" #Rename coords.x1 to Longitude names(AllGrids_df)[28]="y" #Rename coords.x2 to Latitude #Create variables for easy review of categorical analysis AllGrids_df$CatVal<-ifelse(AllGrids_df$Had_cat=="PREDOMINANTLY ANTHROPOGENIC",1,ifelse(AllGrids_df$Had_cat=="INTERMEDIATE", 2,3)) AllGrids_df$Had_cat<-ifelse(AllGrids_df$Had_cat=="PREDOMINANTLY ANTHROPOGENIC","Anthropogenic", ifelse(AllGrids_df$Had_cat=="INTERMEDIATE","Intermediate","Natural")) # Plot shapefile as geom_point map1 <- ggplot() + geom_point(data = AllGrids_df, aes(x = x, y = y, color=Nat_hab)) + geom_point(data=AllGrids_df[which(AllGrids_df$Sampled==1),], aes(x=x, y = y), shape = 0, color = "red", size = 4) map1 # read landscape data for predicting distibution using occupancy models AllGrids_df$std.habitat<-(AllGrids_df$Nat_hab-mean(AllGrids_df$Nat_hab))/sd(AllGrids_df$Nat_hab) #Standardize natural habitat AllGrids_df$std.rainfed<-(AllGrids_df$Rainfed-mean(AllGrids_df$Rainfed))/sd(AllGrids_df$Rainfed) #Standardize rainfed AllGrids_df$std.irrigated<-(AllGrids_df$Irrigated-mean(AllGrids_df$Irrigated))/sd(AllGrids_df$Irrigated) #Standardize irrigated # model predicted distributions using second best model as per AIC table SpecDist<-predict(top.model3.2, newdata=AllGrids_df, param="psi") # merge predicted distribution with shapefile table from all grids SpecDist<-data.frame(SpecDist,AllGrids_df) # Save the Species Distribution model output as a csv file write.csv(SpecDist,file = "SpecDist.csv") # Creae a map frame using coordinates of the data frame SanjayBox <- make_bbox(lon = SpecDist$x, lat = SpecDist$y, f = .1) # Download terrain data from Google for base map SanjayMap <- get_map(location = SanjayBox, maptype = "terrain", source = "google") # Create variations of the map for display map2 <- ggmap(SanjayMap) + scale_fill_gradient(low="white", high= "red")+ geom_point(data=SpecDist,aes(x=x,y=y,alpha=est), shape = 15,size=3.5) map2 map3 <- ggmap(SanjayMap) + geom_point(data=SpecDist,aes(x=x,y=y,colour=est), alpha = 0.5, shape = 15,size=3.5)+ scale_colour_gradientn(colours = terrain.colors(10)) map3 map4 <- ggmap(SanjayMap) + geom_point(data=SpecDist,aes(x=x,y=y,size=est/10), alpha = 0.3) map4 map5 <- ggmap(SanjayMap) + geom_point(data=SpecDist,aes(x=x,y=y,colour=est), alpha = 0.5, shape = 15,size=3.5)+ scale_colour_gradient(low = "white", high = "maroon") map5 map6 <- ggmap(SanjayMap) + geom_point(data=SpecDist,aes(x=x,y=y,colour=est), alpha = 0.5, shape = 15,size=3.5)+ scale_colour_gradient("Psi", low = "white", high = "blue4") + geom_point(data=AllGrids_df[which(AllGrids_df$Sampled==1),], aes(x=x, y = y), shape = 0, color = "red", size = 3.8, alpha=0.3) + xlab("Longitude")+ ylab("Latitude") map6 map6a <- ggmap(SanjayMap) + geom_point(data=SpecDist,aes(x=x,y=y,colour=est), alpha = 0.5, shape = 15,size=3.5)+ scale_colour_gradient("Psi", low = "white", high = "blue4") + geom_point(data=AllGrids_df[which(AllGrids_df$Sampled==1),], aes(x=x, y = y), shape = 4, color = "red", size = 3.5, alpha=0.3) + xlab("Longitude")+ ylab("Latitude") map6a ggsave("PsiMap.png",plot = map6a, scale = 1, width = 10, height = 7, units="cm", dpi = 300) #Create plot of occupancy as a function of habitat plot1<-ggplot()+geom_point(data=SpecDist, aes(x=std.habitat, y=est, colour= Had_cat), shape = 1)+ xlab("Habitat")+ylab("Psi")+ scale_color_manual("Category", values = c("dark blue", "dark green", "maroon"))+ geom_segment(aes(x = AllGrids_df[which(AllGrids_df$Sampled==1),"std.habitat"], y = 0, xend = AllGrids_df[which(AllGrids_df$Sampled==1),"std.habitat"], yend = -.05))+ geom_vline(xintercept=min(SpecDist$std.habitat))+geom_hline(yintercept=0) plot1 #Add 95% CI lines plot2<-plot1+geom_segment(data=SpecDist, aes(x=std.habitat, y = upper_0.95, xend = std.habitat, yend=lower_0.95, colour=Had_cat), linetype="solid", alpha=0.3)+theme_bw() plot2 ggsave("Psi_hab_category.png",plot = plot2, scale = 1, width = 15, height = 10, units="cm", dpi = 300) SpecDist[which(SpecDist$est==min(SpecDist$est)),c("est","lower_0.95", "upper_0.95")] SpecDist[which(SpecDist$est==max(SpecDist$est)),c("est","lower_0.95", "upper_0.95")] # Review coefficients of the top 3 models coef(top.model3,"psi") coef(top.model3.2,"psi") coef(top.model3.3,"psi") coef(top.model3,"p") coef(top.model3.2,"p") coef(top.model3.3,"p") #Create plot of occupancy as a function of habitat and large prey for sampled units only plot3<-ggplot()+geom_point(aes(x=leop.data$unitcov$std.habitat, y=fitted(top.model3,"psi")$est, colour= factor(leop.data$unitcov$Large.Prey)),shape = 1)+ xlab("Habitat")+ylab("Psi")+ scale_color_manual("Large Prey", values = c("black", "red"))+ geom_vline(xintercept=min(SpecDist$std.habitat))+geom_hline(yintercept=0) plot3 #Add 95% CI lines to occupancy output (not predicted) plot4<-plot3+geom_segment(aes(x=leop.data$unitcov$std.habitat, y = fitted(top.model3,"psi")$lower_0.95, xend = leop.data$unitcov$std.habitat, yend=fitted(top.model3,"psi")$upper_0.95, colour=factor(leop.data$unitcov$Large.Prey)), linetype="solid", alpha=0.3)+theme_bw() plot4 ggsave("Psi_hab_LPrey.png",plot = plot4, scale = 1, width = 15, height = 10, units="cm", dpi = 300) #Create plot of occupancy as a function of habitat and large prey for sampled units only plot5<-ggplot()+geom_point(aes(x=leop.data$unitcov$std.rainfed, y=fitted(top.model3,"psi")$est, colour= factor(leop.data$unitcov$Large.Prey)),shape = 1)+ xlab("Rainfed")+ylab("Psi")+ scale_color_manual("Large Prey", values = c("black", "red"))+ geom_vline(xintercept=min(SpecDist$std.rainfed))+geom_hline(yintercept=0) plot5 #Add 95% CI lines to occupancy output (not predicted) plot6<-plot5+geom_segment(aes(x=leop.data$unitcov$std.rainfed, y = fitted(top.model3,"psi")$lower_0.95, xend = leop.data$unitcov$std.rainfed, yend=fitted(top.model3,"psi")$upper_0.95, colour=factor(leop.data$unitcov$Large.Prey)), linetype="solid", alpha=0.3)+theme_bw() plot6 ggsave("Psi_Rainfed_LPrey.png",plot = plot6, scale = 1, width = 15, height = 10, units="cm", dpi = 300) ############################################################################### ########################## End of Code ######################################## ###############################################################################