# load required R libraries
library(RPresence) # not available on CRAN. https://www.mbr-pwrc.usgs.gov/software/presence.shtml
## Warning: package 'RPresence' was built under R version 3.5.0
library(ggplot2)
library(AICcmodavg)
# load in data file
# open this file in excel to examine it's structure if you're not familiar with R
pronghorn<-read.csv(file.path("..","pronghorn.csv"), header=TRUE, as.is=TRUE, strip.white=TRUE)
# slope covariate in degrees. Convert to proportion slope
pronghorn$slope<-pronghorn$slope/90
# distance to water in m. Convert to km
pronghorn$DW<-pronghorn$DW/1000
## create pao data object. Survey.1 and Survey.2 is the detection/nondetection data
prong<-createPao(data=pronghorn[,c("Survey.1","Survey.2")],
unitcov = pronghorn[,4:7])
###############################
# prepare some data frames that will be used later for creating graphs
sagebrush<-expand.grid(sagebrush=seq(0,20,length.out=100),slope=median(prong$unitcov$slope),DW=median(prong$unitcov$DW),aspect=c("N","E","S","W"))
slope<-expand.grid(sagebrush=median(prong$unitcov$sagebrush),slope=seq(0,1,length.out=100),DW=median(prong$unitcov$DW),aspect=c("N","E","S","W"))
DW<-expand.grid(sagebrush=median(prong$unitcov$sagebrush),slope=median(prong$unitcov$slope),DW=seq(0,3,length.out=100),aspect=c("N","E","S","W"))
aspect<-expand.grid(sagebrush=median(prong$unitcov$sagebrush),slope=median(prong$unitcov$slope),DW=median(prong$unitcov$DW),aspect=c("N","E","S","W"))
#################################################################################################
### logistic regression using GLM
### 'naive' analysis that doesn't account for detection probability
#################################################################################################
glm.data<-prong$unitcov
# y =1 if ever detected at plot, =0 if never detected
glm.data$y<-as.numeric(rowSums(prong$det)>0)
glm.mod1<-glm(y~1,data=glm.data,family=binomial)
glm.mod2<-glm(y~sagebrush,data=glm.data,family=binomial)
glm.mod3<-glm(y~slope,data=glm.data,family=binomial)
glm.mod4<-glm(y~sagebrush+slope,data=glm.data,family=binomial)
glm.mod5<-glm(y~DW,data=glm.data,family=binomial)
glm.mod6<-glm(y~sagebrush+DW,data=glm.data,family=binomial)
glm.mod7<-glm(y~slope+DW,data=glm.data,family=binomial)
glm.mod8<-glm(y~sagebrush+slope+DW,data=glm.data,family=binomial)
glm.mod9<-glm(y~aspect,data=glm.data,family=binomial)
glm.mod10<-glm(y~sagebrush+aspect,data=glm.data,family=binomial)
glm.mod11<-glm(y~slope+aspect,data=glm.data,family=binomial)
glm.mod12<-glm(y~sagebrush+slope+aspect,data=glm.data,family=binomial)
glm.mod13<-glm(y~DW+aspect,data=glm.data,family=binomial)
glm.mod14<-glm(y~sagebrush+DW+aspect,data=glm.data,family=binomial)
glm.mod15<-glm(y~slope+DW+aspect,data=glm.data,family=binomial)
glm.mod16<-glm(y~sagebrush+slope+DW+aspect,data=glm.data,family=binomial)
mods<-list(glm.mod1,glm.mod2,glm.mod3,glm.mod4,
glm.mod5,glm.mod6,glm.mod7,glm.mod8,
glm.mod9,glm.mod10,glm.mod11,glm.mod12,
glm.mod13,glm.mod14,glm.mod15,glm.mod16)
## AIC comparison of GLM models
aictab(mods,second.ord=FALSE)
## Warning in aictab.AICglm.lm(mods, second.ord = FALSE):
## Model names have been supplied automatically in the table
##
## Model selection based on AIC:
##
## K AIC Delta_AIC AICWt Cum.Wt LL
## Mod7 3 351.26 0.00 0.23 0.23 -172.63
## Mod5 2 351.48 0.22 0.21 0.44 -173.74
## Mod8 4 352.08 0.82 0.16 0.60 -172.04
## Mod6 3 352.44 1.18 0.13 0.73 -173.22
## Mod15 6 354.05 2.79 0.06 0.79 -171.03
## Mod13 5 354.34 3.08 0.05 0.84 -172.17
## Mod3 2 355.07 3.81 0.03 0.87 -175.54
## Mod16 7 355.31 4.05 0.03 0.91 -170.65
## Mod14 6 355.71 4.45 0.03 0.93 -171.85
## Mod4 3 355.93 4.67 0.02 0.95 -174.97
## Mod1 1 356.89 5.63 0.01 0.97 -177.45
## Mod11 5 357.37 6.11 0.01 0.98 -173.69
## Mod2 2 357.91 6.65 0.01 0.99 -176.95
## Mod12 6 358.71 7.45 0.01 0.99 -173.36
## Mod9 4 358.93 7.67 0.01 1.00 -175.47
## Mod10 5 360.39 9.13 0.00 1.00 -175.19
summary(mods[[7]]) # top aic-ranked model
##
## Call:
## glm(formula = y ~ slope + DW, family = binomial, data = glm.data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.44059 -1.12575 0.08584 1.11069 1.70310
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6088 0.2369 2.570 0.0102 *
## slope -0.8988 0.6097 -1.474 0.1405
## DW -0.3429 0.1437 -2.387 0.0170 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 354.89 on 255 degrees of freedom
## Residual deviance: 345.26 on 253 degrees of freedom
## AIC: 351.26
##
## Number of Fisher Scoring iterations: 4
#####################################################
# produce some plots of model-averaged estimates to
# illustrate estimated effects
#####################################################
# predict presence probability for different sagebrush values, while fixing other values
glm.ma.sagebrush <- modavgPred(mods,type="response",newdata=sagebrush,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = sagebrush, :
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.sagebrush$mod.avg.pred)
l.se<-glm.ma.sagebrush$uncond.se/(glm.ma.sagebrush$mod.avg.pred*(1-glm.ma.sagebrush$mod.avg.pred))
glm.ma.sagebrush$lower<-plogis(logit-1.96*l.se)
glm.ma.sagebrush$upper<-plogis(logit+1.96*l.se)
#jpeg("PronghornNaiveSagebrush.jpg",width=800,height=800,res=144)
#Convert model averaged results to a data frame
glm.ma.sagebrush = as.data.frame(glm.ma.sagebrush)
#Merge model averaged results with sagebrush data set so that plot can be generated
glm.ma.sagebrush = cbind(glm.ma.sagebrush, sagebrush)
# plot only the esimates for northern aspect
glm.ma.sagebrush = glm.ma.sagebrush[glm.ma.sagebrush$aspect == "N",]
ggplot(data = glm.ma.sagebrush, aes(x = sagebrush, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Sagebrush Density")+
geom_line()+
geom_ribbon(aes(ymin = lower.CL, ymax=upper.CL), alpha = .2)+
xlab("Sagebrush density")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict presence probability for different slope values, while fixing other values
glm.ma.slope <- modavgPred(mods,type="response",newdata=slope,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = slope, :
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.slope$mod.avg.pred)
l.se<-glm.ma.slope$uncond.se/(glm.ma.slope$mod.avg.pred*(1-glm.ma.slope$mod.avg.pred))
glm.ma.slope$lower<-plogis(logit-1.96*l.se)
glm.ma.slope$upper<-plogis(logit+1.96*l.se)
#jpeg("PronghornNaiveSlope.jpg",width=800,height=800,res=144)
#Convert model averaged results to a data frame
glm.ma.slope = as.data.frame(glm.ma.slope)
#Merge model averaged results with slope data set so that plot can be generated
glm.ma.slope = cbind(glm.ma.slope, slope)
# plot only the esimates for northern aspect
glm.ma.slope = glm.ma.slope[glm.ma.slope$aspect == "N",]
ggplot(data = glm.ma.slope, aes(x = slope, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Percent Slope")+
geom_line()+
geom_ribbon(aes(ymin = lower.CL, ymax=upper.CL), alpha = .2)+
xlab("Percent Slope")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict presence probability for different sagebrush values, while fixing other values
glm.ma.DW <- modavgPred(mods,type="response",newdata=DW,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = DW, uncond.se = "revised"):
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.DW$mod.avg.pred)
l.se<-glm.ma.DW$uncond.se/(glm.ma.DW$mod.avg.pred*(1-glm.ma.DW$mod.avg.pred))
glm.ma.DW$lower<-plogis(logit-1.96*l.se)
glm.ma.DW$upper<-plogis(logit+1.96*l.se)
#Convert model averaged results to a data frame
glm.ma.DW = as.data.frame(glm.ma.DW)
#Merge model averaged results with distance from water data set so that plot can be generated
glm.ma.DW = cbind(glm.ma.DW, DW)
# plot only the esimates for northern aspect
glm.ma.DW = glm.ma.DW[glm.ma.DW$aspect == "N",]
#jpeg("PronghornNaiveDW.jpg",width=800,height=800,res=144)
ggplot(data = glm.ma.DW, aes(x = DW, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Distance from Water")+
geom_line()+
geom_ribbon(aes(ymin = lower.CL, ymax=upper.CL), alpha = .2)+
xlab("Distance from Water")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict presence probability for different sagebrush values, while fixing other values
glm.ma.aspect <- modavgPred(mods,type="response",newdata=aspect,uncond.se="revised")
## Warning in modavgPred.AICglm.lm(mods, type = "response", newdata = aspect, :
## Model names have been supplied automatically in the table
## manually calculate confidence intervals
logit<-qlogis(glm.ma.aspect$mod.avg.pred)
l.se<-glm.ma.aspect$uncond.se/(glm.ma.aspect$mod.avg.pred*(1-glm.ma.aspect$mod.avg.pred))
glm.ma.aspect$lower<-plogis(logit-1.96*l.se)
glm.ma.aspect$upper<-plogis(logit+1.96*l.se)
#Convert model averaged results to a data frame
glm.ma.aspect = as.data.frame(glm.ma.aspect)
#Merge model averaged results with distance from water data set so that plot can be generated
glm.ma.aspect = cbind(glm.ma.aspect, aspect)
#jpeg("PronghornNaiveAspect.jpg",width=800,height=800,res=144)
ggplot(data = glm.ma.aspect, aes(x = aspect, y = mod.avg.pred)) +
ggtitle("Naive Occupancy Probability as a function of Aspect")+
geom_point()+
geom_errorbar(aes(ymin = lower.CL, ymax=upper.CL), width=.1)+
xlab("Aspect")+
ylab("Model averaged naive occupancy probability")+
ylim(c(0,1))+
scale_x_discrete(labels=c("North","East","South","West"))
#dev.off()
#######################################################
### use occupancy models to fit 16 models for occupancy component
### with general model for detection component
#######################################################
psi.mod1<-occMod(model=list(psi~1,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod2<-occMod(model=list(psi~sagebrush,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
## Warning in file.remove(outfile): cannot remove file 'del_20016_50989.out',
## reason 'Permission denied'
psi.mod3<-occMod(model=list(psi~slope,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.82 signifcant digits.
psi.mod4<-occMod(model=list(psi~sagebrush+slope,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod5<-occMod(model=list(psi~DW,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod6<-occMod(model=list(psi~sagebrush+DW,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod7<-occMod(model=list(psi~slope+DW,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod8<-occMod(model=list(psi~sagebrush+slope+DW,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.74 signifcant digits.
psi.mod9<-occMod(model=list(psi~aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod10<-occMod(model=list(psi~sagebrush+aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod11<-occMod(model=list(psi~slope+aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.19 signifcant digits.
psi.mod12<-occMod(model=list(psi~sagebrush+slope+aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
## Warning in file.remove(outfile): cannot remove file 'del_20016_63188.out',
## reason 'Permission denied'
psi.mod13<-occMod(model=list(psi~DW+aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod14<-occMod(model=list(psi~sagebrush+DW+aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
psi.mod15<-occMod(model=list(psi~slope+DW+aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.7 signifcant digits.
psi.mod16<-occMod(model=list(psi~sagebrush+slope+DW+aspect,p~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.49 signifcant digits.
# compare models by AIC
psi.results<-createAicTable(list(psi.mod1,psi.mod2,psi.mod3,psi.mod4,
psi.mod5,psi.mod6,psi.mod7,psi.mod8,
psi.mod9,psi.mod10,psi.mod11,psi.mod12,
psi.mod13,psi.mod14,psi.mod15,psi.mod16))
# look at AIC table
summary(psi.results)
## Model
## 1 psi(slope)p(sagebrush P slope P DW P aspect)
## 2 psi()p(sagebrush P slope P DW P aspect)
## 3 psi(sagebrush P slope)p(sagebrush P slope P DW P aspect)
## 4 psi(sagebrush)p(sagebrush P slope P DW P aspect)
## 5 psi(slope P DW)p(sagebrush P slope P DW P aspect)
## 6 psi(DW)p(sagebrush P slope P DW P aspect)
## 7 psi(sagebrush P slope P DW)p(sagebrush P slope P DW P aspect)
## 8 psi(sagebrush P DW)p(sagebrush P slope P DW P aspect)
## 9 psi(slope P aspect)p(sagebrush P slope P DW P aspect)
## 10 psi(aspect)p(sagebrush P slope P DW P aspect)
## 11 psi(sagebrush P slope P aspect)p(sagebrush P slope P DW P aspect)
## 12 psi(DW P aspect)p(sagebrush P slope P DW P aspect)
## 13 psi(slope P DW P aspect)p(sagebrush P slope P DW P aspect)
## 14 psi(sagebrush P aspect)p(sagebrush P slope P DW P aspect)
## 15 psi(sagebrush P slope P DW P aspect)p(sagebrush P slope P DW P aspect)
## 16 psi(sagebrush P DW P aspect)p(sagebrush P slope P DW P aspect)
## DAIC wgt npar neg2ll warn.conv warn.VC
## 1 0.00 0.2287 9 615.48 5.82 0
## 2 0.72 0.1593 8 618.20 0.00 0
## 3 0.85 0.1497 10 614.33 0.00 0
## 4 1.12 0.1307 9 616.60 0.00 0
## 5 1.95 0.0861 10 615.43 0.00 0
## 6 2.25 0.0743 9 617.73 0.00 0
## 7 2.85 0.0551 11 614.33 5.74 0
## 8 2.99 0.0513 10 616.47 0.00 0
## 9 5.53 0.0144 12 615.01 5.19 0
## 10 6.05 0.0111 11 617.53 0.00 0
## 11 6.41 0.0093 13 613.89 0.00 0
## 12 6.60 0.0084 12 616.08 0.00 0
## 13 6.87 0.0074 13 614.35 5.70 0
## 14 7.03 0.0068 12 616.51 0.00 0
## 15 8.21 0.0038 14 613.69 5.49 0
## 16 8.34 0.0035 13 615.83 0.00 0
# regression coefficients for occupancy (psi) and detection (p) from top-ranked model
coef(psi.results$models[[1]], "psi")
## est se
## A1_psi 1.661703 0.561173
## A2_psi.psi.slope -2.442779 1.230053
coef(psi.results$models[[1]], "p")
## est se
## B1_p1 -0.010720 0.462496
## B2_p1.p.sagebrush -0.008967 0.030615
## B3_p1.p.slope 0.759702 0.915916
## B4_p1.p.DW -0.418313 0.138760
## B5_p1.p.aspectN 0.954721 0.542644
## B6_p1.p.aspectS 0.271184 0.476801
## B7_p1.p.aspectW -0.007080 0.400031
#####################################################
# produce some plots of model-averaged estimates to
# illustrate estimated effects
#####################################################
# predict occupancy for different sagebrush values, while fixing other values
# plot only the esimates for northern aspect
psi.ma.sagebrush <- modAvg(psi.results,param="psi",predict=TRUE,newdata=sagebrush[sagebrush$aspect == "N",])
#Merge data sets together for plotting
psi.ma.sagebrush = cbind(psi.ma.sagebrush, sagebrush[sagebrush$aspect=="N",])
#jpeg("PronghornPsiSagebrush.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.sagebrush, aes(x = sagebrush, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Sagebrush Density")+
geom_line()+
geom_ribbon(aes(ymin = lower_0.95, ymax=upper_0.95), alpha = .2)+
xlab("Sagebrush density")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict occupancy for different percent slope values, while fixing other values
# plot only the esimates for northern aspect
psi.ma.slope <- modAvg(psi.results,param="psi",predict=TRUE,newdata=slope[slope$aspect == "N",])
#Merge data sets together for plotting
psi.ma.slope = cbind(psi.ma.slope, slope[slope$aspect=="N",])
#jpeg("PronghornPsislope.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.slope, aes(x = slope, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Percent Slope")+
geom_line()+
geom_ribbon(aes(ymin = lower_0.95, ymax=upper_0.95), alpha = .2)+
xlab("Percent Slope")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict occupancy for different distance to water values, while fixing other values
# plot only the esimates for northern aspect
psi.ma.DW <- modAvg(psi.results,param="psi",predict=TRUE,newdata=DW[DW$aspect == "N",])
#Merge data sets together for plotting
psi.ma.DW = cbind(psi.ma.DW, DW[DW$aspect=="N",])
#jpeg("PronghornPsiDW.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.DW, aes(x = DW, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Distance to Water (km)")+
geom_line()+
geom_ribbon(aes(ymin = lower_0.95, ymax=upper_0.95), alpha = .2)+
xlab("Distance to Water (km)")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict occupancy for different aspects, while fixing other values
psi.ma.aspect <- modAvg(psi.results,param="psi",predict=TRUE,newdata=aspect)
#Create data set for plotting
asp = c("North","East","South","West")
psi.ma.aspect = cbind(psi.ma.aspect, asp)
#jpeg("PronghornPsiaspect.jpg",width=800,height=800,res=144)
ggplot(data = psi.ma.aspect, aes(x = asp, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Aspect")+
geom_point()+
geom_errorbar(aes(ymin = lower_0.95, ymax=upper_0.95), width = .1)+
xlab("Aspect")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################################
### use occupancy models to fit 16 different models for detection component
### with general model for occupancy component
##############################################################
p.mod1<-occMod(model=list(p~1,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod2<-occMod(model=list(p~sagebrush,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod3<-occMod(model=list(p~slope,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod4<-occMod(model=list(p~sagebrush+slope,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
## Warning in file.remove(outfile): cannot remove file 'del_20016_57471.out',
## reason 'Permission denied'
p.mod5<-occMod(model=list(p~DW,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.39 signifcant digits.
p.mod6<-occMod(model=list(p~sagebrush+DW,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod7<-occMod(model=list(p~slope+DW,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod8<-occMod(model=list(p~sagebrush+slope+DW,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod9<-occMod(model=list(p~aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod10<-occMod(model=list(p~sagebrush+aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.79 signifcant digits.
p.mod11<-occMod(model=list(p~slope+aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
## Warning in file.remove(outfile): cannot remove file 'del_20016_56263.out',
## reason 'Permission denied'
p.mod12<-occMod(model=list(p~sagebrush+slope+aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.98 signifcant digits.
p.mod13<-occMod(model=list(p~DW+aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod14<-occMod(model=list(p~sagebrush+DW+aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
## Warning in file.remove(outfile): cannot remove file 'del_20016_84746.out',
## reason 'Permission denied'
p.mod15<-occMod(model=list(p~slope+DW+aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
p.mod16<-occMod(model=list(p~sagebrush+slope+DW+aspect,psi~sagebrush+slope+DW+aspect),data=prong,type="so")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.49 signifcant digits.
# compare models by AIC
p.results<-createAicTable(list(p.mod1,p.mod2,p.mod3,p.mod4,
p.mod5,p.mod6,p.mod7,p.mod8,
p.mod9,p.mod10,p.mod11,p.mod12,
p.mod13,p.mod14,p.mod15,p.mod16))
# look at AIC table
summary(p.results)
## Model
## 1 p(DW)psi(sagebrush P slope P DW P aspect)
## 2 p(slope P DW)psi(sagebrush P slope P DW P aspect)
## 3 p(aspect)psi(sagebrush P slope P DW P aspect)
## 4 p(sagebrush P DW)psi(sagebrush P slope P DW P aspect)
## 5 p(DW P aspect)psi(sagebrush P slope P DW P aspect)
## 6 p()psi(sagebrush P slope P DW P aspect)
## 7 p(sagebrush P slope P DW)psi(sagebrush P slope P DW P aspect)
## 8 p(slope P DW P aspect)psi(sagebrush P slope P DW P aspect)
## 9 p(sagebrush P DW P aspect)psi(sagebrush P slope P DW P aspect)
## 10 p(sagebrush P aspect)psi(sagebrush P slope P DW P aspect)
## 11 p(slope P aspect)psi(sagebrush P slope P DW P aspect)
## 12 p(slope)psi(sagebrush P slope P DW P aspect)
## 13 p(sagebrush)psi(sagebrush P slope P DW P aspect)
## 14 p(sagebrush P slope P DW P aspect)psi(sagebrush P slope P DW P aspect)
## 15 p(sagebrush P slope P aspect)psi(sagebrush P slope P DW P aspect)
## 16 p(sagebrush P slope)psi(sagebrush P slope P DW P aspect)
## DAIC wgt npar neg2ll warn.conv warn.VC
## 1 0.00 0.2385 9 618.54 5.39 0
## 2 1.37 0.1200 10 617.91 0.00 0
## 3 1.75 0.0996 11 616.29 0.00 0
## 4 1.79 0.0977 10 618.32 0.00 0
## 5 1.82 0.0960 12 614.36 0.00 0
## 6 2.81 0.0587 8 623.34 0.00 0
## 7 3.17 0.0490 11 617.71 0.00 0
## 8 3.56 0.0402 13 614.10 0.00 0
## 9 3.58 0.0399 13 614.12 0.00 0
## 10 3.62 0.0391 12 616.16 5.79 0
## 11 3.73 0.0369 12 616.27 0.00 0
## 12 4.76 0.0220 9 623.30 0.00 0
## 13 4.79 0.0217 9 623.33 0.00 0
## 14 5.15 0.0182 14 613.69 5.49 0
## 15 5.61 0.0144 13 616.15 5.98 0
## 16 6.75 0.0081 10 623.29 0.00 0
# regression coefficients for occupancy (psi) and detection (p) from top-ranked model
coef(p.results$models[[1]], "psi")
## est se
## A1_psi 0.834604 0.848578
## A2_psi.psi.sagebrush 0.046531 0.065731
## A3_psi.psi.slope -1.754836 1.189384
## A4_psi.psi.DW 0.216324 0.687031
## A5_psi.psi.aspectN 2.686885 5.364325
## A6_psi.psi.aspectS 0.349786 0.960202
## A7_psi.psi.aspectW 0.018775 0.851608
coef(p.results$models[[1]], "p")
## est se
## B1_p1 0.311566 0.264343
## B2_p1.p.DW -0.478974 0.220606
#####################################################
# produce some plots of model-averaged estimates to
# illustrate estimated effects
#####################################################
# predict detection for different sagebrush values, while fixing other values
# plot only the esimates for northern aspect
p.ma.sagebrush <- modAvg(p.results,param="p",predict=TRUE,newdata=sagebrush[sagebrush$aspect == "N",])
#Merge data sets together for plotting
p.ma.sagebrush = cbind(p.ma.sagebrush, sagebrush[sagebrush$aspect=="N",])
#jpeg("PronghornPSagebrush.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.sagebrush, aes(x = sagebrush, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Sagebrush Density")+
geom_line()+
geom_ribbon(aes(ymin = lower_0.95, ymax=upper_0.95), alpha = .2)+
xlab("Sagebrush density")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict detection for different slope values, while fixing other values
# plot only the esimates for northern aspect
p.ma.slope <- modAvg(p.results,param="p",predict=TRUE,newdata=slope[slope$aspect == "N",])
#Merge data sets together for plotting
p.ma.slope = cbind(p.ma.slope, slope[slope$aspect=="N",])
#jpeg("PronghornPslope.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.slope, aes(x = slope, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Percent Slope")+
geom_line()+
geom_ribbon(aes(ymin = lower_0.95, ymax=upper_0.95), alpha = .2)+
xlab("Percent Slope")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict detection for different distance to water values, while fixing other values
# plot only the esimates for northern aspect
p.ma.DW <- modAvg(p.results,param="p",predict=TRUE,newdata=DW[DW$aspect == "N",])
#Merge data sets together for plotting
p.ma.DW = cbind(p.ma.DW, DW[DW$aspect=="N",])
#jpeg("PronghornPDW.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.DW, aes(x = DW, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Distance to Water (km)")+
geom_line()+
geom_ribbon(aes(ymin = lower_0.95, ymax=upper_0.95), alpha = .2)+
xlab("Distance to Water (km)")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()
##############################################
# predict detection for different aspects, while fixing other values
p.ma.aspect <- modAvg(p.results,param="p",predict=TRUE,newdata=aspect)
#Create data set for plotting
asp = c("North","East","South","West")
p.ma.aspect = cbind(p.ma.aspect, asp)
#jpeg("PronghornPaspect.jpg",width=800,height=800,res=144)
ggplot(data = p.ma.aspect, aes(x = asp, y = est)) +
ggtitle("Model Averaged Occupancy Probability as a function of Aspect")+
geom_point()+
geom_errorbar(aes(ymin = lower_0.95, ymax=upper_0.95), width = .1)+
xlab("Aspect")+
ylab("Model averaged occupancy probability")+
ylim(c(0,1))
#dev.off()