# Single Species Single Season Occupancy models using unmarked packages
# unmarked package
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.data <- readxl::read_excel("../salamander.xls",
sheet="CompleteData",
col_names=FALSE) # notice no column names in row 1 of data file.
head(input.data)
## # A tibble: 6 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 0 0 0
## 3 0 1 0 0 0
## 4 1 1 1 1 0
## 5 0 0 1 0 0
## 6 0 0 1 0 0
# Extract the history records
input.history <- input.data[, 1:5] # the history extracted
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 1 1
## 2 0 1 0 0 0
## 3 0 1 0 0 0
## 4 1 1 1 1 0
## 5 0 0 1 0 0
## 6 0 0 1 0 0
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 39
ncol(input.history)
## [1] 5
range(input.history, na.rm=TRUE) # check that all values are either 0 or 1
## [1] 0 1
sum(is.na(input.history)) # are there any missing values?
## [1] 0
# Create the data structore for the occupancy model
# 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.
# 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
obs.covar <- 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)),
D =rep(c("D1","D1","D2","D2","D2"), nrow(input.history)))
obs.covar[1:10,]
## Site Visit Time D
## 1 1 1 T1 D1
## 2 1 2 T2 D1
## 3 1 3 T3 D2
## 4 1 4 T4 D2
## 5 1 5 T5 D2
## 6 2 1 T1 D1
## 7 2 2 T2 D1
## 8 2 3 T3 D2
## 9 2 4 T4 D2
## 10 2 5 T5 D2
salamander.UMF <- unmarked::unmarkedFrameOccu(
y = input.history,
obsCovs=obs.covar)
summary(salamander.UMF)
## unmarkedFrame Object
##
## 39 sites
## Maximum number of observations per site: 5
## Mean number of observations per site: 5
## Sites with at least one detection: 18
##
## Tabulation of y observations:
## 0 1
## 165 30
##
## Observation-level covariates:
## Site Visit Time D
## Min. : 1 Min. :1 T1:39 D1: 78
## 1st Qu.:10 1st Qu.:2 T2:39 D2:117
## Median :20 Median :3 T3:39
## Mean :20 Mean :3 T4:39
## 3rd Qu.:30 3rd Qu.:4 T5:39
## Max. :39 Max. :5
plot(salamander.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=salamander.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"
# These functions return estimates on the LOGIT scale - not very useful for most people
coef(mod.pdot.u, type="state")
## psi(Int)
## 0.383108
SE(mod.pdot.u, type="state")
## psi(Int)
## 0.5086094
confint(mod.pdot.u, type="state")
## 0.025 0.975
## psi(Int) -0.6137481 1.379964
# 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.595 0.123 0.383 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.5946225 0.1225986 0.3512047 0.7989852
# Ditto for the detection probability
coef(mod.pdot.u, type="det")
## p(Int)
## -1.052585
SE(mod.pdot.u, type="det")
## p(Int)
## 0.3008626
confint(mod.pdot.u, type="det")
## 0.025 0.975
## p(Int) -1.642265 -0.4629054
backTransform(mod.pdot.u, type="det")
## Backtransformed linear combination(s) of Detection estimate(s)
##
## Estimate SE LinComb (Intercept)
## 0.259 0.0577 -1.05 1
##
## Transformation: logistic
predict(mod.pdot.u, type='det', newdata=newdata)
## Predicted SE lower upper
## 1 0.258729 0.05770192 0.1621571 0.3862968
# Look at the posterior probability of detection
# The best unbiased predictor of the posterior probability of detection
bup(ranef(mod.pdot.u))
## [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [8] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [15] 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563
## [22] 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563
## [29] 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563 0.2471563
## [36] 1.0000000 1.0000000 1.0000000 1.0000000
# Model where p(t) varies across survey occasions
#
# Note that formula DO NOT HAVE AN = SIGN
# The two formula are for detecton and occupancy in that order
mod.pt.u <- unmarked::occu(~Time ~1 , data=salamander.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
# These functions return estimates on the LOGIT scale - not very useful for most people
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(Time=c("T1","T2","T3","T4","T5"))
predict(mod.pt.u, type='state', newdata=newdata)[1,]
## Predicted SE lower upper
## 1 0.5798805 0.1175673 0.34897 0.7804233
# Ditto for the detection probability
cbind(newdata, predict(mod.pt.u, type='det', newdata=newdata))
## Time Predicted SE lower upper
## 1 T1 0.1768709 0.08451237 0.06443848 0.4013243
## 2 T2 0.1326551 0.07405439 0.04152004 0.3506471
## 3 T3 0.3979631 0.11900439 0.19981035 0.6363510
## 4 T4 0.3537424 0.11369993 0.17115987 0.5919834
## 5 T5 0.2653093 0.10101845 0.11564692 0.4993021
# Look at the posterior probability of detection
# The best unbiased predictor of the posterior probability of detection
bup(ranef(mod.pt.u))
## [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [8] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [15] 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751
## [22] 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751
## [29] 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751 0.2197751
## [36] 1.0000000 1.0000000 1.0000000 1.0000000
sum(bup(ranef(mod.pt.u)))
## [1] 22.61528
# Note that formula DO NOT HAVE AN = SIGN
# The two formula are for detecton and occupancy in that order
mod.pcustom.u <- unmarked::occu(~D ~1 , data=salamander.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
# These functions return estimates on the LOGIT scale - not very useful for most people
# better to use predict() to do the estimation.
# Because there are no covariates here, the data.frame is empty
newdata <- data.frame(D=c("D1","D1","D2","D2","D2"))
predict(mod.pcustom.u, type='state', newdata=newdata)[1,]
## Predicted SE lower upper
## 1 0.583171 0.1186577 0.3495705 0.7845761
# Ditto for the detection probability
cbind(newdata, predict(mod.pcustom.u, type='det', newdata=newdata))
## D Predicted SE lower upper
## 1 D1 0.1539190 0.05839225 0.07024011 0.3046249
## 2 D1 0.1539190 0.05839225 0.07024011 0.3046249
## 3 D2 0.3370867 0.07679144 0.20589609 0.4993087
## 4 D2 0.3370867 0.07679144 0.20589609 0.4993087
## 5 D2 0.3370867 0.07679144 0.20589609 0.4993087
# Look at the posterior probability of detection
# The best unbiased predictor of the posterior probability of detection
bup(ranef(mod.pcustom.u))
## [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [8] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [15] 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649
## [22] 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649
## [29] 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649 0.2258649
## [36] 1.0000000 1.0000000 1.0000000 1.0000000
sum(bup(ranef(mod.pcustom.u)))
## [1] 22.74316
# Model averaging
models.u <-unmarked::fitList(mod.pdot.u, mod.pt.u, mod.pcustom.u)
## Warning in unmarked::fitList(mod.pdot.u, mod.pt.u, mod.pcustom.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.pcustom.u 3 162.82 0.00 0.760 0.76
## mod.pdot.u 2 165.76 2.94 0.175 0.93
## mod.pt.u 6 167.71 4.90 0.066 1.00
# Get model averaged estimates of occupancy
predict(models.u, type="state")[1,]
## Predicted SE lower upper
## 1 0.5849543 0.119358 0.3498164 0.7868192
# 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(obs.covar, predict(models.u, type="det"))[1:10,]
## Site Visit Time D Predicted SE lower upper
## 1 1 1 T1 D1 0.1737253 0.07033638 0.08590731 0.3252341
## 2 1 2 T2 D1 0.1708219 0.07001399 0.08440237 0.3219064
## 3 1 3 T3 D2 0.3274033 0.08354904 0.19785991 0.4885763
## 4 1 4 T4 D2 0.3244996 0.08210524 0.19597858 0.4856629
## 5 1 5 T5 D2 0.3186927 0.08202387 0.19233333 0.4795770
## 6 2 1 T1 D1 0.1737253 0.07033638 0.08590731 0.3252341
## 7 2 2 T2 D1 0.1708219 0.07001399 0.08440237 0.3219064
## 8 2 3 T3 D2 0.3274033 0.08354904 0.19785991 0.4885763
## 9 2 4 T4 D2 0.3244996 0.08210524 0.19597858 0.4856629
## 10 2 5 T5 D2 0.3186927 0.08202387 0.19233333 0.4795770