-
Notifications
You must be signed in to change notification settings - Fork 1
/
woodp_Occ.R
82 lines (66 loc) · 2.61 KB
/
woodp_Occ.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
library(unmarked)
library(AICcmodavg)
library(RCurl)
#Choose a directory where the data are downloaded to and load the data
x <- getURL("https://raw.githubusercontent.com/chrissuthy/Occupancy_Chapter/master/woodp_Occ.csv")
woodpecker <- read.csv(text = x)
head(woodpecker)
n.sites <- nrow(woodpecker)
woodp.y <- woodpecker[,c("y.1","y.2","y.3")]
woodp.sitecovs <- woodpecker[,c("snags"),drop=F]
woodp.obsCovs <- woodpecker[,c("date.1","date.2","date.3")]
# create unmarked frame with data
umf <- unmarkedFrameOccu(y = woodp.y,
siteCovs = woodp.sitecovs,
obsCovs = list(date=woodp.obsCovs))
umf
summary(umf)
# fit a collection of models
modList <- fitList(
#~p ~psi
"p(.)psi(.)" = occu(~1 ~1, umf),
"p(date)psi(.)" = occu(~scale(date) ~1, umf),
"p(.)psi(snags)" = occu(~1 ~scale(log(snags)), umf),
"p(date)psi(snags)" = occu(~scale(date) ~scale(log(snags)), umf)
)
# AIC model selection
modSel(modList)
#mean detection (95% CI)
plogis(coef(modList@fits[[4]])[3])
plogis(confint(modList@fits[[4]],type="det")[1,])
#mean occupancy (95% CI)
plogis(coef(modList@fits[[4]])[1])
plogis(confint(modList@fits[[4]],type="state")[1,])
#finite sample occupancy
sum(bup(ranef(modList@fits[[4]]),stat="mode"))/n.sites
apply(confint(ranef(modList@fits[[4]]),level=.95),2,sum)/n.sites
#------------------
# model predictions
#------------------
#png("woodp_Occ.png",width = 6,height=3.5,units="in",res=600)
#postscript("Figure7.2.eps",width = 6,height=3.5)
#par(mfrow=c(1,2),cex=.8)
postscript("Figure 7.2a.eps",width = 3,height=3.5)
par(cex=0.8)
# predict/plot date effect on detection
date.dat <- data.frame(date=150:190)
date.fit <- predict(modList@fits[[4]],type="det",newdata=date.dat,appendData=T)
with(date.fit, {
plot(date, Predicted, type="l",xlim=c(150,190),ylim=c(0,1),las=1,bty="n",lwd=2,
xlab="Ordinal date", ylab="Detection probability")
lines(date, Predicted+1.96*SE, lty=2, col="darkgray")
lines(date, Predicted-1.96*SE, lty=2, col="darkgray")
})
dev.off()
postscript("Figure 7.2b.eps",width = 3,height=3.5)
par(cex=0.8)
# predict/plot snag density effect on occupancy
snags.dat <- data.frame(snags=seq(0.25,6,by=0.25))
snags.fit <- predict(modList@fits[[4]],type="state",newdata=snags.dat,appendData=T)
with(snags.fit, {
plot(snags, Predicted, type="l",xlim=c(0,6),ylim=c(0,1),las=1,bty="n",lwd=2,
xlab="Snag density (#/ha)", ylab="Nest occupancy probability")
lines(snags, Predicted+1.96*SE, lty=2, col="darkgray")
lines(snags, Predicted-1.96*SE, lty=2, col="darkgray")
})
dev.off()