# Single Species Single Season Occupancy models
# Weta Example
# Using unmarked
library(readxl)
library(unmarked)
## Loading required package: reshape
## Loading required package: lattice
## Loading required package: parallel
## Loading required package: Rcpp
library(ggplot2)
# get the data read in
# Data for detections should be a data frame with rows corresponding to sites
# and columns to visits.
# The usual 1=detected; 0=not detected; NA=not visited is used.
input.history <- readxl::read_excel(file.path("..","weta.xls"),
sheet="detection_histories",
na="-",
col_names=FALSE) # notice no column names in row 1 of data file.
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 72
ncol(input.history)
## [1] 5
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 98
head(input.history)
## # A tibble: 6 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 NA
## 2 0 0 0 0 NA
## 3 0 0 0 1 NA
## 4 0 0 0 0 NA
## 5 0 0 0 0 NA
## 6 0 0 0 0 NA
# Get the site level covariates
site_covar <- readxl::read_excel(file.path("..","weta.xls"),
sheet="site_covar",
na="-",
col_names=TRUE) # notice col_names in row 1 of table.
# Create an alternate site level covariate that is a categorical variable rather
# than indicator variables
site_covar$BrowCat <- paste(c("","B")[1+unlist(site_covar[,1])], c("","N")[1+unlist(site_covar[,2])], sep="")
xtabs(~BrowCat, data=site_covar,exclude=NULL, na.action=na.pass)
## BrowCat
## B N
## 35 37
colSums(site_covar[,1:2])
## Browsed Unbrowsed
## 35 37
head(site_covar)
## # A tibble: 6 x 3
## Browsed Unbrowsed BrowCat
## <dbl> <dbl> <chr>
## 1 1 0 B
## 2 1 0 B
## 3 1 0 B
## 4 0 1 N
## 5 1 0 B
## 6 0 1 N
# Visit Covariates
# unmarked does not have any reserved keywords like RPresence so
# we need to construct a covariate data frame with the covariates
# again stacked. The ordering is different in unmarked.
# Here you want for site 1, the visit specific covariates, then
# for site 2, etc.
# This needs to be added to the unmarked dataframe.
# Get the individual covariates.
obs1 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs1",
na="-",
col_names=FALSE)
obs2 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs2",
na="-",
col_names=FALSE)
obs3 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs3",
na="-",
col_names=FALSE)
obs <- obs1*1 + obs2*2 + obs3*3
head(obs)
## X__1 X__2 X__3 X__4 X__5
## 1 1 3 2 3 NA
## 2 1 3 2 3 NA
## 3 1 3 2 3 NA
## 4 1 3 2 3 NA
## 5 1 3 2 3 NA
## 6 1 3 2 3 NA
# We also want to set up covariate values where we
# set the detection probability equal in first 2 occasions and last 3 occasions.
# We need to define a survey covariate that has two levels
survey.covar.u <- data.frame( Site=rep(1:nrow(input.history), each=ncol(input.history)),
Visit=rep(1:ncol(input.history), nrow(input.history)),
Time =rep(c("T1","T2","T3","T4","T5"), nrow(input.history)),
obs1 =as.vector(t(as.matrix(obs1))),
obs2 =as.vector(t(as.matrix(obs2))),
obs3 =as.vector(t(as.matrix(obs3))),
obs =as.character(as.vector(t(as.matrix(obs)))), stringsAsFactors=FALSE)
survey.covar.u[1:10,]
## Site Visit Time obs1 obs2 obs3 obs
## 1 1 1 T1 1 0 0 1
## 2 1 2 T2 0 0 1 3
## 3 1 3 T3 0 1 0 2
## 4 1 4 T4 0 0 1 3
## 5 1 5 T5 NA NA NA <NA>
## 6 2 1 T1 1 0 0 1
## 7 2 2 T2 0 0 1 3
## 8 2 3 T3 0 1 0 2
## 9 2 4 T4 0 0 1 3
## 10 2 5 T5 NA NA NA <NA>
# check that missing values in history and observer covariates align
select <- is.na(as.vector(t(as.matrix(input.history))))
survey.covar.u[select,]
## Site Visit Time obs1 obs2 obs3 obs
## 5 1 5 T5 NA NA NA <NA>
## 10 2 5 T5 NA NA NA <NA>
## 15 3 5 T5 NA NA NA <NA>
## 20 4 5 T5 NA NA NA <NA>
## 25 5 5 T5 NA NA NA <NA>
## 30 6 5 T5 NA NA NA <NA>
## 35 7 5 T5 NA NA NA <NA>
## 40 8 5 T5 NA NA NA <NA>
## 45 9 5 T5 NA NA NA <NA>
## 50 10 5 T5 NA NA NA <NA>
## 53 11 3 T3 NA NA NA <NA>
## 58 12 3 T3 NA NA NA <NA>
## 63 13 3 T3 NA NA NA <NA>
## 68 14 3 T3 NA NA NA <NA>
## 73 15 3 T3 NA NA NA <NA>
## 76 16 1 T1 NA NA NA <NA>
## 79 16 4 T4 NA NA NA <NA>
## 81 17 1 T1 NA NA NA <NA>
## 84 17 4 T4 NA NA NA <NA>
## 86 18 1 T1 NA NA NA <NA>
## 89 18 4 T4 NA NA NA <NA>
## 91 19 1 T1 NA NA NA <NA>
## 94 19 4 T4 NA NA NA <NA>
## 99 20 4 T4 NA NA NA <NA>
## 102 21 2 T2 NA NA NA <NA>
## 104 21 4 T4 NA NA NA <NA>
## 107 22 2 T2 NA NA NA <NA>
## 109 22 4 T4 NA NA NA <NA>
## 112 23 2 T2 NA NA NA <NA>
## 114 23 4 T4 NA NA NA <NA>
## 117 24 2 T2 NA NA NA <NA>
## 119 24 4 T4 NA NA NA <NA>
## 122 25 2 T2 NA NA NA <NA>
## 124 25 4 T4 NA NA NA <NA>
## 129 26 4 T4 NA NA NA <NA>
## 130 26 5 T5 NA NA NA <NA>
## 134 27 4 T4 NA NA NA <NA>
## 135 27 5 T5 NA NA NA <NA>
## 139 28 4 T4 NA NA NA <NA>
## 140 28 5 T5 NA NA NA <NA>
## 144 29 4 T4 NA NA NA <NA>
## 145 29 5 T5 NA NA NA <NA>
## 150 30 5 T5 NA NA NA <NA>
## 155 31 5 T5 NA NA NA <NA>
## 160 32 5 T5 NA NA NA <NA>
## 165 33 5 T5 NA NA NA <NA>
## 186 38 1 T1 NA NA NA <NA>
## 188 38 3 T3 NA NA NA <NA>
## 191 39 1 T1 NA NA NA <NA>
## 193 39 3 T3 NA NA NA <NA>
## 196 40 1 T1 NA NA NA <NA>
## 201 41 1 T1 NA NA NA <NA>
## 206 42 1 T1 NA NA NA <NA>
## 211 43 1 T1 NA NA NA <NA>
## 212 43 2 T2 NA NA NA <NA>
## 216 44 1 T1 NA NA NA <NA>
## 217 44 2 T2 NA NA NA <NA>
## 221 45 1 T1 NA NA NA <NA>
## 222 45 2 T2 NA NA NA <NA>
## 226 46 1 T1 NA NA NA <NA>
## 227 46 2 T2 NA NA NA <NA>
## 231 47 1 T1 NA NA NA <NA>
## 232 47 2 T2 NA NA NA <NA>
## 240 48 5 T5 NA NA NA <NA>
## 245 49 5 T5 NA NA NA <NA>
## 250 50 5 T5 NA NA NA <NA>
## 255 51 5 T5 NA NA NA <NA>
## 260 52 5 T5 NA NA NA <NA>
## 265 53 5 T5 NA NA NA <NA>
## 270 54 5 T5 NA NA NA <NA>
## 275 55 5 T5 NA NA NA <NA>
## 280 56 5 T5 NA NA NA <NA>
## 285 57 5 T5 NA NA NA <NA>
## 288 58 3 T3 NA NA NA <NA>
## 293 59 3 T3 NA NA NA <NA>
## 298 60 3 T3 NA NA NA <NA>
## 303 61 3 T3 NA NA NA <NA>
## 308 62 3 T3 NA NA NA <NA>
## 311 63 1 T1 NA NA NA <NA>
## 314 63 4 T4 NA NA NA <NA>
## 316 64 1 T1 NA NA NA <NA>
## 319 64 4 T4 NA NA NA <NA>
## 321 65 1 T1 NA NA NA <NA>
## 324 65 4 T4 NA NA NA <NA>
## 326 66 1 T1 NA NA NA <NA>
## 329 66 4 T4 NA NA NA <NA>
## 331 67 1 T1 NA NA NA <NA>
## 334 67 4 T4 NA NA NA <NA>
## 337 68 2 T2 NA NA NA <NA>
## 339 68 4 T4 NA NA NA <NA>
## 342 69 2 T2 NA NA NA <NA>
## 344 69 4 T4 NA NA NA <NA>
## 347 70 2 T2 NA NA NA <NA>
## 349 70 4 T4 NA NA NA <NA>
## 352 71 2 T2 NA NA NA <NA>
## 354 71 4 T4 NA NA NA <NA>
## 357 72 2 T2 NA NA NA <NA>
## 359 72 4 T4 NA NA NA <NA>
sum(is.na(survey.covar.u[!select,]))
## [1] 0
weta.UMF <- unmarked::unmarkedFrameOccu(
y = input.history,
siteCovs=site_covar,
obsCovs=survey.covar.u
)
summary(weta.UMF)
## unmarkedFrame Object
##
## 72 sites
## Maximum number of observations per site: 5
## Mean number of observations per site: 3.64
## Sites with at least one detection: 35
##
## Tabulation of y observations:
## 0 1 <NA>
## 206 56 98
##
## Site-level covariates:
## Browsed Unbrowsed BrowCat
## Min. :0.0000 Min. :0.0000 Length:72
## 1st Qu.:0.0000 1st Qu.:0.0000 Class :character
## Median :0.0000 Median :1.0000 Mode :character
## Mean :0.4861 Mean :0.5139
## 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000
##
## Observation-level covariates:
## Site Visit Time obs1
## Min. : 1.00 Min. :1 Length:360 Min. :0.0000
## 1st Qu.:18.75 1st Qu.:2 Class :character 1st Qu.:0.0000
## Median :36.50 Median :3 Mode :character Median :0.0000
## Mean :36.50 Mean :3 Mean :0.3473
## 3rd Qu.:54.25 3rd Qu.:4 3rd Qu.:1.0000
## Max. :72.00 Max. :5 Max. :1.0000
## NA's :98
## obs2 obs3 obs
## Min. :0.0000 Min. :0.0000 Length:360
## 1st Qu.:0.0000 1st Qu.:0.0000 Class :character
## Median :0.0000 Median :0.0000 Mode :character
## Mean :0.3092 Mean :0.3435
## 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000
## NA's :98 NA's :98
plot(weta.UMF)
# Fit some models.
# Note that formula DO NOT HAVE AN = SIGN
# The two formula are for detecton and occupancy in that order
mod.pdot.u <- unmarked::occu(~1 ~1 , data=weta.UMF, se=TRUE)
# This creates a complicated S4 object with SLOTS that can be accessed using various
# functions. This is more complcated than with RPresence
slotNames(mod.pdot.u)
## [1] "knownOcc" "fitType" "call"
## [4] "formula" "data" "sitesRemoved"
## [7] "estimates" "AIC" "opt"
## [10] "negLogLike" "nllFun" "bootstrapSamples"
## [13] "covMatBS"
mod.pdot.u@estimates
## Occupancy:
## Estimate SE z P(>|z|)
## 0.475 0.375 1.27 0.205
##
## Detection:
## Estimate SE z P(>|z|)
## -0.622 0.236 -2.63 0.00853
# These functions return estimates on the LOGIT scale - not very useful for most people
coef(mod.pdot.u, type="state")
## psi(Int)
## 0.4752421
SE(mod.pdot.u, type="state")
## psi(Int)
## 0.3745705
confint(mod.pdot.u, type="state")
## 0.025 0.975
## psi(Int) -0.2589026 1.209387
# It is possible to use backTransform() but only if NO covariates
backTransform(mod.pdot.u, type='state')
## Backtransformed linear combination(s) of Occupancy estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.617 0.0885 0.475 1
##
## Transformation: logistic
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(factor=1:1)
predict(mod.pdot.u, type='state', newdata=newdata)
## Predicted SE lower upper
## 1 0.6166237 0.08854806 0.4356335 0.7701904
# Ditto for the detection probability
coef(mod.pdot.u, type="det")
## p(Int)
## -0.6217868
SE(mod.pdot.u, type="det")
## p(Int)
## 0.2363836
confint(mod.pdot.u, type="det")
## 0.025 0.975
## p(Int) -1.08509 -0.1584834
backTransform(mod.pdot.u, type="det")
## Backtransformed linear combination(s) of Detection estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.349 0.0537 -0.622 1
##
## Transformation: logistic
predict(mod.pdot.u, type='det', newdata=newdata)
## Predicted SE lower upper
## 1 0.3493752 0.05373287 0.252544 0.4604619
# Look at the posterior probability of detection
# The best unbiased predictor of the posterior probability of detection
bup(ranef(mod.pdot.u))
## [1] 0.2237324 0.2237324 1.0000000 0.2237324 0.2237324 0.2237324 0.2237324
## [8] 0.2237324 0.2237324 1.0000000 1.0000000 1.0000000 1.0000000 0.2237324
## [15] 1.0000000 1.0000000 0.3069911 1.0000000 0.3069911 0.2237324 0.3069911
## [22] 0.3069911 1.0000000 0.3069911 0.3069911 0.3069911 0.3069911 1.0000000
## [29] 0.3069911 0.2237324 0.2237324 0.2237324 0.2237324 1.0000000 0.1579091
## [36] 0.1579091 1.0000000 0.3069911 1.0000000 0.2237324 1.0000000 1.0000000
## [43] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.2237324 1.0000000
## [50] 1.0000000 1.0000000 1.0000000 1.0000000 0.2237324 0.2237324 0.2237324
## [57] 1.0000000 0.2237324 1.0000000 1.0000000 1.0000000 1.0000000 0.3069911
## [64] 0.3069911 1.0000000 0.3069911 1.0000000 1.0000000 1.0000000 0.3069911
## [71] 0.3069911 1.0000000
# p(.) psi(browse)
mod.pdot.psiB.1.u <- unmarked::occu(~1 ~-1+Browsed+Unbrowsed, data=weta.UMF, se=TRUE)
names(mod.pdot.psiB.1.u)
## [1] "state" "det"
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(Browsed=c(1,0), Unbrowsed=c(0,1))
cbind(newdata,predict(mod.pdot.psiB.1.u, type='state', newdata=newdata))
## Browsed Unbrowsed Predicted SE lower upper
## 1 1 0 0.7593382 0.1198171 0.4660501 0.9193922
## 2 0 1 0.4809820 0.1078825 0.2843231 0.6837155
coef(mod.pdot.psiB.1.u)
## psi(Browsed) psi(Unbrowsed) p(Int)
## 1.14905452 -0.07610881 -0.62220870
vcov(mod.pdot.psiB.1.u)
## psi(Browsed) psi(Unbrowsed) p(Int)
## psi(Browsed) 0.42988569 0.04877514 -0.07657277
## psi(Unbrowsed) 0.04877514 0.18675814 -0.03566979
## p(Int) -0.07657277 -0.03566979 0.05599849
mod.pdot.psiB.1.u.oddsratio.browse <- exp( sum(c(1,-1,0)*coef(mod.pdot.psiB.1.u)))
mod.pdot.psiB.1.u.oddsratio.browse
## [1] 3.404722
est <- linearComb(mod.pdot.psiB.1.u, c(1,-1), type="state")
exp(est@estimate)
## [1] 3.404722
# using the BrowCat variable directly as cell effect models
mod.pdot.psiB.2.u <- unmarked::occu(~1 ~BrowCat, data=weta.UMF, se=TRUE)
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(BrowCat=c("B","N"))
cbind(newdata,predict(mod.pdot.psiB.2.u, type='state', newdata=newdata))
## BrowCat Predicted SE lower upper
## 1 B 0.7593680 0.1198254 0.4660415 0.9194190
## 2 N 0.4809969 0.1078859 0.2843301 0.6837338
coef(mod.pdot.psiB.2.u)
## psi(Int) psi(BrowCatN) p(Int)
## 1.1492178 -1.2252670 -0.6222472
vcov(mod.pdot.psiB.2.u)
## psi(Int) psi(BrowCatN) p(Int)
## psi(Int) 0.43001803 -0.38122288 -0.07659795
## psi(BrowCatN) -0.38122288 0.51919696 0.04092222
## p(Int) -0.07659795 0.04092222 0.05600327
mod.pdot.psiB.2.u.oddsratio.browse <- exp( sum(c(0,-1,0)*coef(mod.pdot.psiB.2.u)))
mod.pdot.psiB.2.u.oddsratio.browse
## [1] 3.405075
est <- linearComb(mod.pdot.psiB.2.u, c(0,-1), type="state")
est
## Linear combination(s) of Occupancy estimate(s)
##
## Estimate SE (Intercept) BrowCatN
## 1.23 0.721 0 -1
confint(est)
## 0.025 0.975
## -0.1869914 2.637525
exp(est@estimate)
## [1] 3.405075
exp(confint(est))
## 0.025 0.975
## 0.8294509 13.97857
# Impact of browse on detectability as well?.
mod.pB.psiB.u <- unmarked::occu(~BrowCat ~BrowCat, data=weta.UMF, se=TRUE)
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(BrowCat=c("B","N"))
cbind(newdata,predict(mod.pB.psiB.u, type='state', newdata=newdata))
## BrowCat Predicted SE lower upper
## 1 B 0.7563981 0.1277022 0.4439049 0.9235363
## 2 N 0.4840246 0.1191192 0.2691862 0.7049346
newdata <- data.frame(BrowCat=c("B","N"))
cbind(newdata,predict(mod.pB.psiB.u, type='det', newdata=newdata))
## BrowCat Predicted SE lower upper
## 1 B 0.3519078 0.06891004 0.2309670 0.4953830
## 2 N 0.3451077 0.08602485 0.1999464 0.5263258
coef(mod.pB.psiB.u)
## psi(Int) psi(BrowCatN) p(Int) p(BrowCatN)
## 1.13303258 -1.19695607 -0.61066377 -0.02995047
mod.pB.psiB.u.oddsratio.browse <- exp( sum(c(0,-1,0,0)*coef(mod.pB.psiB.u)))
mod.pB.psiB.u.oddsratio.browse
## [1] 3.310026
est <- linearComb(mod.pB.psiB.u, c(0,-1), type="state")
exp(est@estimate)
## [1] 3.310026
# Model where p varies by time psi(browse)
#
mod.pt.psiB.u <- unmarked::occu(~Time ~BrowCat, data=weta.UMF, se=TRUE)
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(BrowCat=c("B","N"))
cbind(newdata,predict(mod.pt.psiB.u, type='state', newdata=newdata))
## BrowCat Predicted SE lower upper
## 1 B 0.7699156 0.1232913 0.4610099 0.9290342
## 2 N 0.4931714 0.1116875 0.2884153 0.7002442
newdata <- expand.grid(Time=c("T1","T2","T3","T4","T5"))
cbind(newdata,predict(mod.pt.psiB.u, type='det', newdata=newdata))
## Time Predicted SE lower upper
## 1 T1 0.3520566 0.09838300 0.18918319 0.5585560
## 2 T2 0.3175290 0.08921528 0.17192419 0.5104358
## 3 T3 0.1694829 0.06707747 0.07424059 0.3417987
## 4 T4 0.3115818 0.08815495 0.16822723 0.5031935
## 5 T5 0.5917250 0.10433923 0.38334101 0.7716398
coef(mod.pt.psiB.u)
## psi(Int) psi(BrowCatN) p(Int) p(TimeT2) p(TimeT3)
## 1.2078348 -1.2351507 -0.6100113 -0.1551397 -0.9792848
## p(TimeT4) p(TimeT5)
## -0.1827234 0.9811124
mod.pt.psiB.u.oddsratio.browse <- exp( sum(c(0,-1,0,0,0,0,0,0,0)*coef(mod.pt.psiB.u)))
## Warning in c(0, -1, 0, 0, 0, 0, 0, 0, 0) * coef(mod.pt.psiB.u): longer
## object length is not a multiple of shorter object length
mod.pt.psiB.u.oddsratio.browse
## [1] 3.438897
est <- linearComb(mod.pt.psiB.u, c(0,-1), type="state")
exp(est@estimate)
## [1] 3.438897
# Model where p(obs) psi(browse)
#
mod.pO.psiB.u <- unmarked::occu(~obs ~BrowCat, data=weta.UMF, se=TRUE)
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(BrowCat=c("B","N"))
cbind(newdata,predict(mod.pO.psiB.u, type='state', newdata=newdata))
## BrowCat Predicted SE lower upper
## 1 B 0.7535927 0.1164150 0.4723052 0.9126656
## 2 N 0.4846725 0.1084303 0.2865485 0.6877343
newdata <- expand.grid(obs=c("1",'2','3'))
cbind(newdata,predict(mod.pO.psiB.u, type='det', newdata=newdata))
## obs Predicted SE lower upper
## 1 1 0.2235216 0.06273691 0.1241470 0.3689348
## 2 2 0.3786085 0.07998866 0.2383337 0.5426262
## 3 3 0.4462474 0.08081043 0.2980095 0.6047043
coef(mod.pO.psiB.u)
## psi(Int) psi(BrowCatN) p(Int) p(obs2) p(obs3)
## 1.1178664 -1.1791957 -1.2452609 0.7498024 1.0294164
mod.pO.psiB.u.oddsratio.browse <- exp( sum(c(0,-1,0,0,0)*coef(mod.pO.psiB.u)))
mod.pO.psiB.u.oddsratio.browse
## [1] 3.251758
est <- linearComb(mod.pO.psiB.u, c(0,-1), type="state")
exp(est@estimate)
## [1] 3.251758
# Model where p varies by observer + visit but constant over time psi(browse)
#
mod.pOpV.psiB.u <- unmarked::occu(~obs+Time ~BrowCat, data=weta.UMF, se=TRUE)
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(BrowCat=c("B","N"))
cbind(newdata,predict(mod.pOpV.psiB.u, type='state', newdata=newdata))
## BrowCat Predicted SE lower upper
## 1 B 0.7673286 0.1206506 0.4672431 0.9253797
## 2 N 0.5059973 0.1148964 0.2938184 0.7160377
newdata <- expand.grid(obs=c("1",'2','3'), Time=c("T1","T2","T3","T4","T5"))
cbind(newdata,predict(mod.pOpV.psiB.u, type='det', newdata=newdata))
## obs Time Predicted SE lower upper
## 1 1 T1 0.21421538 0.08831907 0.08882487 0.4325798
## 2 2 T1 0.36049060 0.12278182 0.16560616 0.6155310
## 3 3 T1 0.44282792 0.11677497 0.23915456 0.6677296
## 4 1 T2 0.18929552 0.08154684 0.07613362 0.3981664
## 5 2 T2 0.32560546 0.10240877 0.16216606 0.5463518
## 6 3 T2 0.40502120 0.11960469 0.20466779 0.6429518
## 7 1 T3 0.09596579 0.04999011 0.03317466 0.2472154
## 8 2 T3 0.17999076 0.08137649 0.06932898 0.3927476
## 9 3 T3 0.23633705 0.09558452 0.09877780 0.4663370
## 10 1 T4 0.20256237 0.08233531 0.08553713 0.4082208
## 11 2 T4 0.34436763 0.10851937 0.16995622 0.5739918
## 12 3 T4 0.42547244 0.12515845 0.21351480 0.6688907
## 13 1 T5 0.43431646 0.12607706 0.21924963 0.6773285
## 14 2 T5 0.61353688 0.11669323 0.37697892 0.8064021
## 15 3 T5 0.69120237 0.11020106 0.44863353 0.8602883
coef(mod.pOpV.psiB.u)
## psi(Int) psi(BrowCatN) p(Int) p(obs2) p(obs3)
## 1.19328816 -1.16929772 -1.29970076 0.72646532 1.07000786
## p(TimeT2) p(TimeT3) p(TimeT4) p(TimeT5)
## -0.15489346 -0.94317467 -0.07065501 1.03543938
mod.pOpV.psiB.u.oddsratio.browse <- exp( sum(c(0,-1,0,0,0,0,0,0,0)*coef(mod.pOpV.psiB.u)))
mod.pOpV.psiB.u.oddsratio.browse
## [1] 3.219731
est <- linearComb(mod.pOpV.psiB.u, c(0,-1), type="state")
exp(est@estimate)
## [1] 3.219731
# Other models
mod.pOpVpB.psiB.u <- unmarked::occu(~obs+Time+BrowCat ~BrowCat, data=weta.UMF, se=TRUE)
mod.pOpVpB.psi.u <- unmarked::occu(~obs+Time+BrowCat ~1, data=weta.UMF, se=TRUE)
mod.pOpVpB.psi.u <- unmarked::occu(~obs+Time ~1, data=weta.UMF, se=TRUE)
# Model averaging
models.u <-unmarked::fitList(
mod.pdot.u,
mod.pdot.psiB.1.u,
mod.pB.psiB.u,
mod.pt.psiB.u,
mod.pO.psiB.u,
mod.pOpV.psiB.u,
mod.pOpVpB.psiB.u,
mod.pOpVpB.psi.u
)
## Warning in unmarked::fitList(mod.pdot.u, mod.pdot.psiB.1.u,
## mod.pB.psiB.u, : Your list was unnamed, so model names were added as object
## names
aic.table.u <- unmarked::modSel(models.u)
aic.table.u
## nPars AIC delta AICwt cumltvWt
## mod.pOpV.psiB.u 9 257.60 0.00 0.3901 0.39
## mod.pOpVpB.psi.u 8 258.55 0.95 0.2429 0.63
## mod.pt.psiB.u 7 259.44 1.84 0.1556 0.79
## mod.pOpVpB.psiB.u 10 259.60 2.00 0.1435 0.93
## mod.pO.psiB.u 5 262.04 4.44 0.0423 0.97
## mod.pdot.psiB.1.u 3 264.26 6.66 0.0139 0.99
## mod.pdot.u 2 265.79 8.19 0.0065 0.99
## mod.pB.psiB.u 4 266.26 8.66 0.0051 1.00
# Get model averaged estimates of occupancy
predict(models.u, type="state")[1,]
## Predicted SE lower upper
## 1 0.7354181 0.1283259 0.4575714 0.8945918
# Get model averaged estimates of detection.
# Notice these are one big list of nsites x nvisits long
# with the detection probabilities for each site
cbind(survey.covar.u, predict(models.u, type="det"))[1:10,]
## Site Visit Time obs1 obs2 obs3 obs Predicted SE lower
## 1 1 1 T1 1 0 0 1 0.2377461 0.10117589 0.10841329
## 2 1 2 T2 0 0 1 3 0.3921885 0.11811420 0.20417431
## 3 1 3 T3 0 1 0 2 0.1895520 0.08835743 0.08021826
## 4 1 4 T4 0 0 1 3 0.4075166 0.12539126 0.21032236
## 5 1 5 T5 NA NA NA <NA> NA NA NA
## 6 2 1 T1 1 0 0 1 0.2377461 0.10117589 0.10841329
## 7 2 2 T2 0 0 1 3 0.3921885 0.11811420 0.20417431
## 8 2 3 T3 0 1 0 2 0.1895520 0.08835743 0.08021826
## 9 2 4 T4 0 0 1 3 0.4075166 0.12539126 0.21032236
## 10 2 5 T5 NA NA NA <NA> NA NA NA
## upper
## 1 0.4500927
## 2 0.6177367
## 3 0.3941259
## 4 0.6375424
## 5 NA
## 6 0.4500927
## 7 0.6177367
## 8 0.3941259
## 9 0.6375424
## 10 NA