# Multi state single season
# Nichols et al (2007) Ecology 88:1395-1400 looked at breeding and
# non-breeding California Spotted Owls in s = 54 sites over k = 5
# surveys.
# Hand-fed mice con rmed occupancy; breeding status only
# confrmed when owl took mouse to nest.
# . = site not visited.
# 0 = owl not detected.
# 1 = owl detected, but no detection of young.
# 2 = owl detected, along with evidence of young.
# Territories that were expected to be occupied were selectively monitored
# and so psi is expected to be close to 1 and not really of interest.
# See MacKenzie et al, Section 5.7.1 for more details
# 2018-10-01 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
#
library(ggplot2)
library(plyr)
library(readxl)
library(reshape2)
library(RPresence)
# get the data
occ.data <- readxl::read_excel(file.path("..","CalSpottedOwl.xls"),
sheet="RawData",
skip=10, na='.')
head(occ.data)
## # A tibble: 6 x 6
## Site `1` `2` `3` `4` `5`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 NA 2 2 1
## 2 2 0 1 1 NA NA
## 3 3 1 NA 1 NA NA
## 4 4 1 NA 1 NA NA
## 5 5 NA 1 2 1 2
## 6 6 1 1 1 1 NA
input.history <- occ.data[,-1]
input.history[1:5,]
## # A tibble: 5 x 5
## `1` `2` `3` `4` `5`
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 NA 2 2 1
## 2 0 1 1 NA NA
## 3 1 NA 1 NA NA
## 4 1 NA 1 NA NA
## 5 NA 1 2 1 2
Nsites <- nrow(input.history)
Nvisits<- ncol(input.history)
# For biological reasons, we believe that errors in identifying states
# varies between the first 2 visits and the last 3 visit (early bs late
# breeding season). We create a site-time covariates.
# The format is the "long" format ordered by visit major order,i.e.
# stack the columns of the matrix of period covariates
#
# The covariate matrix would be arranged originally as
# 1 2 3 4 5
# site1 [E E L L L]
# site2 [E E L L L]
# : [: : : : :]
# siteN [E E L L L]
#
# and this is stacked by columns
period <- data.frame(Site=rep(1:Nsites, Nvisits),
Visit=rep(1:Nvisits, each=Nsites),
Per =rep(c("E","E","L","L","L"), each=Nsites))
period[1:10,]
## Site Visit Per
## 1 1 1 E
## 2 2 1 E
## 3 3 1 E
## 4 4 1 E
## 5 5 1 E
## 6 6 1 E
## 7 7 1 E
## 8 8 1 E
## 9 9 1 E
## 10 10 1 E
# create the pao file
owl.pao <- createPao(input.history,
survcov=period,
title="owl multi state single season")
summary(owl.pao)
## paoname=pres.pao
## title=owl multi state single season
## Naive occ=0.8703704
## naiveR =0.3518519
## nunits nsurveys nseasons nsurveyseason nmethods
## "54" "5" "1" "5" "1"
## nunitcov nsurvcov
## "1" "4"
## unit covariates : TEMP
## survey covariates: Site Visit Per SURVEY
# Use the multi-state model, but only for a single season
#
# psi occupancy probability,
# r probability of being in the second state, conditional on the unit being occupied.
# p detection probability. .
# delta probability of detecting the second state in a survey, conditional on the species being detected in the survey.
mod1 <- occMod(model=list(
psi~1, # occupancy regardless of state
r~1, # occupancy in state 2 | occupied
p~1, # detection in state i
delta~1),# identified as state 2 if detected (and in state 2)
data=owl.pao,
type="do.ms.2")
## PRESENCE Version 2.12.21.
summary(mod1)
## Model name=psi()r()p()delta()
## AIC=335.3784
## -2*log-likelihood=327.3784
## num. par=4
names(mod1$real)
## [1] "psi" "Cpsi0" "Cpsi1" "Cpsi2" "r" "CR0" "CR1" "CR2"
## [9] "delta" "p1" "p2"
mod1$real$psi[1,] # prob site is occupied
## est se lower upper
## unit1 0.9752 0.0378 0.6475 0.9988
mod1$real$r[1,] # prob in state 2 | occupied
## est se lower upper
## unit1 0.5543 0.1124 0.3376 0.7521
mod1$real$p1[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 1
## est se lower upper
## p1(1-1)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-2)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-3)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-4)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-5)1unit1 0.7448 0.0363 0.6674 0.8093
mod1$real$p2[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 2
## est se lower upper
## p2(1-1)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-2)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-3)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-4)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-5)1unit1 0.7448 0.0363 0.6674 0.8093
mod1$real$delta[seq(1,by=Nsites, length.out=Nvisits),] # prob identify in in state 2 if detect in state 2
## est se lower upper
## delta(1-1)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-2)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-3)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-4)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-5)1unit1 0.423 0.0843 0.2713 0.5907
# Need a sutiable model for the detection probabilities
mod2 <- occMod(model=list(
psi~1, # occupancy regardless of state
r~1, # occupancy in state 2 | occupied
p~SURVEY, # detection in each state vary by survey
delta~1),# identified as state 2 if detected (and in state 2)
data=owl.pao,
type="do.ms.2")
## PRESENCE Version 2.12.21.
summary(mod2)
## Model name=psi()r()p(SURVEY)delta()
## AIC=326.3843
## -2*log-likelihood=310.3843
## num. par=8
names(mod2$real)
## [1] "psi" "Cpsi0" "Cpsi1" "Cpsi2" "r" "CR0" "CR1" "CR2"
## [9] "delta" "p1" "p2"
mod2$real$psi[1,] # prob site is occupied
## est se lower upper
## unit1 0.9789 0.0324 0.6817 0.999
mod2$real$r[1,] # prob in state 2 | occupied
## est se lower upper
## unit1 0.5543 0.1124 0.3376 0.7521
mod2$real$p1[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 1
## est se lower upper
## p1(1-1)1unit1 0.5790 0.0821 0.4154 0.7269
## p1(1-2)1unit1 0.8238 0.0678 0.6518 0.9211
## p1(1-3)1unit1 0.9224 0.0429 0.7859 0.9747
## p1(1-4)1unit1 0.6970 0.0771 0.5293 0.8248
## p1(1-5)1unit1 0.6196 0.1123 0.3904 0.8056
mod2$real$p2[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 2
## est se lower upper
## p2(1-1)1unit1 0.5790 0.0821 0.4154 0.7269
## p2(1-2)1unit1 0.8238 0.0678 0.6518 0.9211
## p2(1-3)1unit1 0.9224 0.0429 0.7859 0.9747
## p2(1-4)1unit1 0.6970 0.0771 0.5293 0.8248
## p2(1-5)1unit1 0.6196 0.1123 0.3904 0.8056
mod2$real$delta[seq(1,by=Nsites, length.out=Nvisits),] # prob identify in in state 2 if detect in state 2
## est se lower upper
## delta(1-1)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-2)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-3)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-4)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-5)1unit1 0.423 0.0843 0.2713 0.5907
# Need a sutiable model for the detection probabilities
mod3 <- occMod(model=list(
psi~1, # occupancy regardless of state
r~1, # occupancy in state 2 | occupied
p~STATE, # detection in each state vary by state
delta~1),# identified as state 2 if detected (and in state 2)
data=owl.pao,
type="do.ms.2")
## PRESENCE Version 2.12.21.
summary(mod3)
## Model name=psi()r()p(STATE)delta()
## AIC=337.3627
## -2*log-likelihood=327.3627
## num. par=5
names(mod3$real)
## [1] "psi" "Cpsi0" "Cpsi1" "Cpsi2" "r" "CR0" "CR1" "CR2"
## [9] "delta" "p1" "p2"
mod3$real$psi[1,] # prob site is occupied
## est se lower upper
## unit1 0.9753 0.0378 0.6451 0.9988
mod3$real$r[1,] # prob in state 2 | occupied
## est se lower upper
## unit1 0.5623 0.1321 0.3097 0.7862
mod3$real$p1[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 1
## est se lower upper
## p1(1-1)1unit1 0.7531 0.0767 0.5761 0.8726
## p1(1-2)1unit1 0.7531 0.0767 0.5761 0.8726
## p1(1-3)1unit1 0.7531 0.0767 0.5761 0.8726
## p1(1-4)1unit1 0.7531 0.0767 0.5761 0.8726
## p1(1-5)1unit1 0.7531 0.0767 0.5761 0.8726
mod3$real$p2[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 2
## est se lower upper
## p2(1-1)1unit1 0.738 0.0659 0.5908 0.846
## p2(1-2)1unit1 0.738 0.0659 0.5908 0.846
## p2(1-3)1unit1 0.738 0.0659 0.5908 0.846
## p2(1-4)1unit1 0.738 0.0659 0.5908 0.846
## p2(1-5)1unit1 0.738 0.0659 0.5908 0.846
mod3$real$delta[seq(1,by=Nsites, length.out=Nvisits),] # prob identify in in state 2 if detect in state 2
## est se lower upper
## delta(1-1)1unit1 0.4203 0.0874 0.2641 0.5943
## delta(1-2)1unit1 0.4203 0.0874 0.2641 0.5943
## delta(1-3)1unit1 0.4203 0.0874 0.2641 0.5943
## delta(1-4)1unit1 0.4203 0.0874 0.2641 0.5943
## delta(1-5)1unit1 0.4203 0.0874 0.2641 0.5943
# Need a sutiable model for the detection probabilities
# On the surface, this model does not appear to converge, but the problem is
# tht delta is estimated to be 0 in the first period which causes numerical
# problems (logit(0) is -infinity.)
mod4 <- occMod(model=list(
psi~1, # occupancy regardless of state
r~1, # occupancy in state 2 | occupied
p~1, # detection in each state vary by state
delta~Per),# identified as state 2 if detected (and in state 2)
data=owl.pao,
type="do.ms.2")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 2.62 signifcant digits.
summary(mod4)
## Model name=psi()r()p()delta(Per)
## AIC=288.3492
## -2*log-likelihood=278.3492
## num. par=5
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 2.62 signifcant digits.
## Warning: Problem with estimation of VC covariance. Ignore all SE's and values in VC matrix.
names(mod4$real)
## [1] "psi" "Cpsi0" "Cpsi1" "Cpsi2" "r" "CR0" "CR1" "CR2"
## [9] "delta" "p1" "p2"
mod4$beta$delta
## est se
## C1_delta(1-1).delta.int -30.25791 NaN
## C2_delta(1-1).delta.PerL 32.14419 NaN
mod4$real$psi[1,] # prob site is occupied
## est se lower upper
## unit1 0.9752 0.0378 0.6475 0.9988
mod4$real$r[1,] # prob in state 2 | occupied
## est se lower upper
## unit1 0.4495 0.0802 0.302 0.6065
mod4$real$p1[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 1
## est se lower upper
## p1(1-1)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-2)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-3)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-4)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-5)1unit1 0.7448 0.0363 0.6674 0.8093
mod4$real$p2[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 2
## est se lower upper
## p2(1-1)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-2)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-3)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-4)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-5)1unit1 0.7448 0.0363 0.6674 0.8093
mod4$real$delta[seq(1,by=Nsites, length.out=Nvisits),] # prob identify in in state 2 if detect in state 2
## est se lower upper
## delta(1-1)1unit1 0.0000 NaN NaN NaN
## delta(1-2)1unit1 0.0000 NaN NaN NaN
## delta(1-3)1unit1 0.8683 0.0727 0.6547 0.9582
## delta(1-4)1unit1 0.8683 0.0727 0.6547 0.9582
## delta(1-5)1unit1 0.8683 0.0727 0.6547 0.9582
# Need a sutiable model for the detection probabilities
mod5 <- occMod(model=list(
psi~1, # occupancy regardless of state
r~1, # occupancy in state 2 | occupied
p~SURVEY, # detection in each state vary by state
delta~Per),# identified as state 2 if detected (and in state 2)
data=owl.pao,
type="do.ms.2")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 0.39 signifcant digits.
summary(mod5)
## Model name=psi()r()p(SURVEY)delta(Per)
## AIC=279.355
## -2*log-likelihood=261.355
## num. par=9
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 0.39 signifcant digits.
## Warning: Problem with estimation of VC covariance. Ignore all SE's and values in VC matrix.
names(mod5$real)
## [1] "psi" "Cpsi0" "Cpsi1" "Cpsi2" "r" "CR0" "CR1" "CR2"
## [9] "delta" "p1" "p2"
mod5$real$psi[1,] # prob site is occupied
## est se lower upper
## unit1 0.9789 0.0324 0.6817 0.999
mod5$real$r[1,] # prob in state 2 | occupied
## est se lower upper
## unit1 0.4495 0.0802 0.302 0.6065
mod5$real$p1[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 1
## est se lower upper
## p1(1-1)1unit1 0.5790 0.0821 0.4154 0.7269
## p1(1-2)1unit1 0.8238 0.0678 0.6518 0.9211
## p1(1-3)1unit1 0.9224 0.0429 0.7859 0.9747
## p1(1-4)1unit1 0.6970 0.0771 0.5293 0.8248
## p1(1-5)1unit1 0.6196 0.1123 0.3904 0.8056
mod5$real$p2[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 2
## est se lower upper
## p2(1-1)1unit1 0.5790 0.0821 0.4154 0.7269
## p2(1-2)1unit1 0.8238 0.0678 0.6518 0.9211
## p2(1-3)1unit1 0.9224 0.0429 0.7859 0.9747
## p2(1-4)1unit1 0.6970 0.0771 0.5293 0.8248
## p2(1-5)1unit1 0.6196 0.1123 0.3904 0.8056
mod5$real$delta[seq(1,by=Nsites, length.out=Nvisits),] # prob identify in in state 2 if detect in state 2
## est se lower upper
## delta(1-1)1unit1 0.0000 NaN NaN NaN
## delta(1-2)1unit1 0.0000 NaN NaN NaN
## delta(1-3)1unit1 0.8683 0.0727 0.6547 0.9582
## delta(1-4)1unit1 0.8683 0.0727 0.6547 0.9582
## delta(1-5)1unit1 0.8683 0.0727 0.6547 0.9582
# Need a sutiable model for the detection probabilities
mod6 <- occMod(model=list(
psi~1, # occupancy regardless of state
r~1, # occupancy in state 2 | occupied
p~STATE, # detection in each state vary by state
delta~Per),# identified as state 2 if detected (and in state 2)
data=owl.pao,
type="do.ms.2")
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 3.25 signifcant digits.
summary(mod6)
## Model name=psi()r()p(STATE)delta(Per)
## AIC=290.24
## -2*log-likelihood=278.24
## num. par=6
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 3.25 signifcant digits.
## Warning: Problem with estimation of VC covariance. Ignore all SE's and values in VC matrix.
names(mod6$real)
## [1] "psi" "Cpsi0" "Cpsi1" "Cpsi2" "r" "CR0" "CR1" "CR2"
## [9] "delta" "p1" "p2"
mod6$real$psi[1,] # prob site is occupied
## est se lower upper
## unit1 0.9754 0.0378 0.6428 0.9989
mod6$real$r[1,] # prob in state 2 | occupied
## est se lower upper
## unit1 0.4421 0.0823 0.2918 0.6037
mod6$real$p1[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 1
## est se lower upper
## p1(1-1)1unit1 0.7322 0.053 0.6168 0.8229
## p1(1-2)1unit1 0.7322 0.053 0.6168 0.8229
## p1(1-3)1unit1 0.7322 0.053 0.6168 0.8229
## p1(1-4)1unit1 0.7322 0.053 0.6168 0.8229
## p1(1-5)1unit1 0.7322 0.053 0.6168 0.8229
mod6$real$p2[seq(1,by=Nsites, length.out=Nvisits),] # prob detect if in state 2
## est se lower upper
## p2(1-1)1unit1 0.7601 0.0578 0.6299 0.8551
## p2(1-2)1unit1 0.7601 0.0578 0.6299 0.8551
## p2(1-3)1unit1 0.7601 0.0578 0.6299 0.8551
## p2(1-4)1unit1 0.7601 0.0578 0.6299 0.8551
## p2(1-5)1unit1 0.7601 0.0578 0.6299 0.8551
mod6$real$delta[seq(1,by=Nsites, length.out=Nvisits),] # prob identify in in state 2 if detect in state 2
## est se lower upper
## delta(1-1)1unit1 0.0000 NaN NaN NaN
## delta(1-2)1unit1 0.0000 NaN NaN NaN
## delta(1-3)1unit1 0.8708 0.0718 0.6587 0.9593
## delta(1-4)1unit1 0.8708 0.0718 0.6587 0.9593
## delta(1-5)1unit1 0.8708 0.0718 0.6587 0.9593
# Model averaging
models<-list(mod1, mod2, mod3,
mod4, mod5, mod6
)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll warn.conv warn.VC
## 1 psi()r()p(SURVEY)delta(Per) 0.00 0.9848 9 261.36 0.39 1
## 2 psi()r()p()delta(Per) 8.99 0.0110 5 278.35 2.62 1
## 3 psi()r()p(STATE)delta(Per) 10.88 0.0043 6 278.24 3.25 1
## 4 psi()r()p(SURVEY)delta() 47.03 0.0000 8 310.38 0.00 0
## 5 psi()r()p()delta() 56.02 0.0000 4 327.38 0.00 0
## 6 psi()r()p(STATE)delta() 58.01 0.0000 5 327.36 0.00 0
RPresence::modAvg(results, param="psi")[1,]
## est se lower_0.95 upper_0.95
## unit1 0.9788445 0.03249207 0.681173 0.999003
RPresence::modAvg(results, param="r") [1,]
## est se lower_0.95 upper_0.95
## unit1 0.4494685 0.08021052 0.3019253 0.6064727
RPresence::modAvg(results, param="p1")[seq(1,by=Nsites, length.out=Nvisits),]
## est se lower_0.95 upper_0.95
## p1(1-1)1unit1 0.5814721 0.08402199 0.4138997 0.7321398
## p1(1-2)1unit1 0.8225428 0.06823416 0.6496395 0.9205532
## p1(1-3)1unit1 0.9196407 0.04828519 0.7607676 0.9762945
## p1(1-4)1unit1 0.6976745 0.07687768 0.5304329 0.8250024
## p1(1-5)1unit1 0.6214536 0.11255304 0.3912629 0.8074386
RPresence::modAvg(results, param="p2")[seq(1,by=Nsites, length.out=Nvisits),]
## est se lower_0.95 upper_0.95
## p2(1-1)1unit1 0.5815911 0.08426813 0.4135243 0.7326347
## p2(1-2)1unit1 0.8226617 0.06811742 0.6500722 0.9205333
## p2(1-3)1unit1 0.9197596 0.04787949 0.7627009 0.9761219
## p2(1-4)1unit1 0.6977934 0.07696727 0.5303263 0.8252270
## p2(1-5)1unit1 0.6215725 0.11269472 0.3910747 0.8077185
RPresence::modAvg(results, param="delta")[seq(1,by=Nsites, length.out=Nvisits),]
## est se lower_0.95 upper_0.95
## delta(1-1)1unit1 2.593953e-11 NaN NaN NaN
## delta(1-2)1unit1 2.593953e-11 NaN NaN NaN
## delta(1-3)1unit1 8.683107e-01 0.07269637 0.6547659 0.9581997
## delta(1-4)1unit1 8.683107e-01 0.07269637 0.6547659 0.9581997
## delta(1-5)1unit1 8.683107e-01 0.07269637 0.6547659 0.9581997