# Single Species, Multi Season Occupancy analyais
# Northern Spotted Owl (Strix occidentalis caurina) in California.
# s=55 sites visited up to K=5 times per season between 1997 and 2001 (Y=5).
# Detection probabilities relatively constant within years, but likely different among years.
# Using RMark
# RPresence package
library(readxl)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
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 <- read.csv(file.path("..","NSO.csv"), header=FALSE, skip=2, na.strings="-")
input.data$V1 <- NULL # drop the site number
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.data)
## [1] 55
ncol(input.data)
## [1] 40
range(input.data, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.data))
## [1] 752
# Extract the history records and create a capture history
input.history <- data.frame(freq=1,
ch=apply(input.data,1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 0111NANANANA010NANANANANA11NANANANANANA011NANANANANA0101NANANANA
## 2 1 00NANANANANANA000NANANANANA000NANANANANA000110NANA0011NANANANA
## 3 1 1000111NA0100NANANANA0001NANANANA000000NANA00000NANANA
## 4 1 111NANANANANA10NANANANANANA011NANANANANA00000100000011NANA
## 5 1 000000NANA0000000NA000000NANA000000NANA0111NANANANA
## 6 1 1NANANANANANANA000111NANA0111NANANANA01111NANANA011111NANA
# Change any NA to . in the chapter history
input.history$ch <- gsub("NA",".", input.history$ch, fixed=TRUE)
head(input.history)
## freq ch
## 1 1 0111....010.....11......011.....0101....
## 2 1 00......000.....000.....000110..0011....
## 3 1 1000111.0100....0001....000000..00000...
## 4 1 111.....10......011.....00000100000011..
## 5 1 000000..0000000.000000..000000..0111....
## 6 1 1.......000111..0111....01111...011111..
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 55
# Create the data structure
max.visit.per.year <- 8
n.season <- 5
nso.data <- process.data(data=input.history, model="RDOccupEG",
time.intervals=c( rep( c(rep(0,max.visit.per.year-1),1),n.season-1),
rep(0,max.visit.per.year-1)))
summary(nso.data)
## Length Class Mode
## data 2 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 55 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 5 -none- numeric
## time.intervals 39 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 1 -none- numeric
## group.covariates 0 -none- NULL
## nstrata 1 -none- numeric
## strata.labels 0 -none- NULL
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
# What are the parameter names for Single Season Single Species models
setup.parameters("RDOccupEG", check=TRUE)
## [1] "Psi" "Epsilon" "Gamma" "p"
#there are other parameterizations available
setup.parameters("RDOccupPE", check=TRUE) # psi, epsilon, p
## [1] "Psi" "Epsilon" "p"
setup.parameters("RDOccupPG", check=TRUE) # psi, gamma, p
## [1] "Psi" "Gamma" "p"
# Fit a models.
# We start with the first formulation in terms of psi, gamma, epsilon, p (do.1_)
# Note that formula DO NOT HAVE AN = SIGN
mod.fit <- RMark::mark(nso.data,
model="RDOccupEG",
model.parameters=list(
Psi =list(formula=~1), # initial occupancy
p =list(formula=~1),
Epsilon=list(formula=~1),
Gamma =list(formula=~1)))
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~1)
##
## Npar : 4
## -2lnL: 1355.315
## AICc : 1363.464
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5372049 0.2889444 -0.0291261 1.1035359
## Epsilon:(Intercept) -1.7288379 0.2595729 -2.2376007 -1.2200751
## Gamma:(Intercept) -1.4883083 0.2843058 -2.0455478 -0.9310689
## p:(Intercept) -0.0210406 0.0743374 -0.1667418 0.1246607
##
##
## Real Parameter Psi
##
## 1
## 0.631162
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.1507363 0.1507363 0.1507363 0.1507363
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1841758 0.1841758 0.1841758 0.1841758
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:2
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:4
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:5
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
names(mod.fit)
## [1] "data" "model" "title"
## [4] "model.name" "links" "mixtures"
## [7] "call" "parameters" "time.intervals"
## [10] "number.of.groups" "group.labels" "nocc"
## [13] "begin.time" "covariates" "fixed"
## [16] "design.matrix" "pims" "design.data"
## [19] "strata.labels" "mlogit.list" "nocc.secondary"
## [22] "profile.int" "simplify" "model.parameters"
## [25] "results" "output"
names(mod.fit$results)
## [1] "lnl" "deviance" "deviance.df"
## [4] "npar" "n" "AICc"
## [7] "beta" "real" "beta.vcv"
## [10] "derived" "derived.vcv" "covariate.values"
## [13] "singular" "real.vcv"
# Estimate of initial occupancy
names(mod.fit$results$real)
## [1] "estimate" "se" "lcl" "ucl" "fixed" "note"
mod.fit$results$real
## estimate se lcl ucl fixed note
## Psi g1 a0 t1 0.6311620 0.0672653 0.4927190 0.7509220
## Epsilon g1 a0 t1 0.1507363 0.0332292 0.0964244 0.2279232
## Gamma g1 a0 t1 0.1841758 0.0427184 0.1145030 0.2827079
## p g1 s1 t1 0.4947401 0.0185823 0.4584109 0.5311249
get.real(mod.fit, parameter="Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 1 1 0.631162 0.0672653 0.492719 0.750922
## fixed note group age time Age Time
## Psi g1 a0 t1 1 0 1 0 0
get.real(mod.fit, parameter="p" , se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 s1 t1 10 4 0.4947401 0.0185823 0.4584109
## p g1 s1 t2 11 4 0.4947401 0.0185823 0.4584109
## p g1 s1 t3 12 4 0.4947401 0.0185823 0.4584109
## p g1 s1 t4 13 4 0.4947401 0.0185823 0.4584109
## p g1 s1 t5 14 4 0.4947401 0.0185823 0.4584109
## p g1 s1 t6 15 4 0.4947401 0.0185823 0.4584109
## p g1 s1 t7 16 4 0.4947401 0.0185823 0.4584109
## p g1 s1 t8 17 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t1 18 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t2 19 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t3 20 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t4 21 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t5 22 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t6 23 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t7 24 4 0.4947401 0.0185823 0.4584109
## p g1 s2 t8 25 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t1 26 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t2 27 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t3 28 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t4 29 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t5 30 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t6 31 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t7 32 4 0.4947401 0.0185823 0.4584109
## p g1 s3 t8 33 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t1 34 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t2 35 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t3 36 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t4 37 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t5 38 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t6 39 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t7 40 4 0.4947401 0.0185823 0.4584109
## p g1 s4 t8 41 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t1 42 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t2 43 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t3 44 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t4 45 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t5 46 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t6 47 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t7 48 4 0.4947401 0.0185823 0.4584109
## p g1 s5 t8 49 4 0.4947401 0.0185823 0.4584109
## ucl fixed note group time session Time
## p g1 s1 t1 0.5311249 1 1 1 0
## p g1 s1 t2 0.5311249 1 2 1 1
## p g1 s1 t3 0.5311249 1 3 1 2
## p g1 s1 t4 0.5311249 1 4 1 3
## p g1 s1 t5 0.5311249 1 5 1 4
## p g1 s1 t6 0.5311249 1 6 1 5
## p g1 s1 t7 0.5311249 1 7 1 6
## p g1 s1 t8 0.5311249 1 8 1 7
## p g1 s2 t1 0.5311249 1 1 2 0
## p g1 s2 t2 0.5311249 1 2 2 1
## p g1 s2 t3 0.5311249 1 3 2 2
## p g1 s2 t4 0.5311249 1 4 2 3
## p g1 s2 t5 0.5311249 1 5 2 4
## p g1 s2 t6 0.5311249 1 6 2 5
## p g1 s2 t7 0.5311249 1 7 2 6
## p g1 s2 t8 0.5311249 1 8 2 7
## p g1 s3 t1 0.5311249 1 1 3 0
## p g1 s3 t2 0.5311249 1 2 3 1
## p g1 s3 t3 0.5311249 1 3 3 2
## p g1 s3 t4 0.5311249 1 4 3 3
## p g1 s3 t5 0.5311249 1 5 3 4
## p g1 s3 t6 0.5311249 1 6 3 5
## p g1 s3 t7 0.5311249 1 7 3 6
## p g1 s3 t8 0.5311249 1 8 3 7
## p g1 s4 t1 0.5311249 1 1 4 0
## p g1 s4 t2 0.5311249 1 2 4 1
## p g1 s4 t3 0.5311249 1 3 4 2
## p g1 s4 t4 0.5311249 1 4 4 3
## p g1 s4 t5 0.5311249 1 5 4 4
## p g1 s4 t6 0.5311249 1 6 4 5
## p g1 s4 t7 0.5311249 1 7 4 6
## p g1 s4 t8 0.5311249 1 8 4 7
## p g1 s5 t1 0.5311249 1 1 5 0
## p g1 s5 t2 0.5311249 1 2 5 1
## p g1 s5 t3 0.5311249 1 3 5 2
## p g1 s5 t4 0.5311249 1 4 5 3
## p g1 s5 t5 0.5311249 1 5 5 4
## p g1 s5 t6 0.5311249 1 6 5 5
## p g1 s5 t7 0.5311249 1 7 5 6
## p g1 s5 t8 0.5311249 1 8 5 7
get.real(mod.fit, parameter="Epsilon", se=TRUE)
## all.diff.index par.index estimate se lcl
## Epsilon g1 a0 t1 2 2 0.1507363 0.0332292 0.0964244
## Epsilon g1 a1 t2 3 2 0.1507363 0.0332292 0.0964244
## Epsilon g1 a2 t3 4 2 0.1507363 0.0332292 0.0964244
## Epsilon g1 a3 t4 5 2 0.1507363 0.0332292 0.0964244
## ucl fixed note group age time Age Time
## Epsilon g1 a0 t1 0.2279232 1 0 1 0 0
## Epsilon g1 a1 t2 0.2279232 1 1 2 1 1
## Epsilon g1 a2 t3 0.2279232 1 2 3 2 2
## Epsilon g1 a3 t4 0.2279232 1 3 4 3 3
get.real(mod.fit, parameter="Gamma", se=TRUE)
## all.diff.index par.index estimate se lcl
## Gamma g1 a0 t1 6 3 0.1841758 0.0427184 0.114503
## Gamma g1 a1 t2 7 3 0.1841758 0.0427184 0.114503
## Gamma g1 a2 t3 8 3 0.1841758 0.0427184 0.114503
## Gamma g1 a3 t4 9 3 0.1841758 0.0427184 0.114503
## ucl fixed note group age time Age Time
## Gamma g1 a0 t1 0.2827079 1 0 1 0 0
## Gamma g1 a1 t2 0.2827079 1 1 2 1 1
## Gamma g1 a2 t3 0.2827079 1 2 3 2 2
## Gamma g1 a3 t4 0.2827079 1 3 4 3 3
# Derived parameters - estimated occupancy for each unit in years 2....
names(mod.fit$results$derived)
## [1] "psi Probability Occupied" "lambda Rate of Change"
## [3] "log odds lambda"
mod.fit$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.6311620 0.06726525 0.4993221 0.7630018
## 2 0.6039540 0.05095844 0.5040754 0.7038325
## 3 0.5858583 0.05112760 0.4856482 0.6860684
## 4 0.5738230 0.05661538 0.4628569 0.6847892
## 5 0.5658186 0.06213794 0.4440282 0.6876089
mod.fit$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 0.9568922 0.05137574 0.8561958 1.057589
## 2 0.9700379 0.03736524 0.8968021 1.043274
## 3 0.9794571 0.02655675 0.9274059 1.031508
## 4 0.9860506 0.01858500 0.9496240 1.022477
mod.fit$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 0.8911547 0.13254333 0.6313697 1.150940
## 2 0.9276527 0.08800695 0.7551590 1.100146
## 3 0.9517972 0.05914705 0.8358690 1.067725
## 4 0.9678720 0.03998058 0.8895100 1.046234