# Single Species, MultiSeason Occupancy analyais
# Grand Skinks
# Data has been collected on 352 tors over a 5 year (the ÂseasonsÂ) period,^
# although not all tors (rock piles) were surveyed each year, with up to 3 surveys^
# of each tor per year.
# The 15 columns are in 5 blocks of 3.
# There is also a site-specific covariate Pasture indicating whether the surrounding^
# matrix is either predominately the modified habitat (farm pasture, Pasture =1) or
# "native" grassland (tussock, Pasture = 0).
# 2018-11-09 Code submitted by Neil Faught
# Using the RMark package
library(car)
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 RMark additional functions
source(file.path("..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
# 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(file.path("..","GrandSkinks.xls"),
sheet="DetectionHistory",
skip=3, na="-",
col_names=FALSE)
head(input.data)
## # A tibble: 6 x 18
## X__1 X__2 X__3 X__4 X__5 X__6 X__7 X__8 X__9 X__10 X__11 X__12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00 0 0 0 0 0 0 0 0 NA 0 0
## 2 2.00 0 0 0 0 0 NA 0 0 0 0 0
## 3 3.00 0 0 0 0 0 NA 0 0 0 0 0
## 4 4.00 0 NA NA 0 NA NA 0 NA NA 0 0
## 5 5.00 0 NA NA NA NA NA NA NA NA NA NA
## 6 6.00 0 NA NA NA NA NA NA NA NA NA NA
## # ... with 6 more variables: X__13 <dbl>, X__14 <dbl>, X__15 <dbl>, X__16
## # <dbl>, X__17 <lgl>, X__18 <dbl>
input.history <- input.data[, 2:16]
head(input.history)
## # A tibble: 6 x 15
## X__2 X__3 X__4 X__5 X__6 X__7 X__8 X__9 X__10 X__11 X__12 X__13
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 0 0 0 0 NA 0 0 0
## 2 0 0 0 0 0 NA 0 0 0 0 0 0
## 3 0 0 0 0 0 NA 0 0 0 0 0 0
## 4 0 NA NA 0 NA NA 0 NA NA 0 0 NA
## 5 0 NA NA NA NA NA NA NA NA NA NA NA
## 6 0 NA NA NA NA NA NA NA NA NA NA NA
## # ... with 3 more variables: X__14 <dbl>, X__15 <dbl>, X__16 <dbl>
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 352
ncol(input.history)
## [1] 15
range(input.history, na.rm=TRUE)
## [1] 0 1
#Looks like there are some missing values
sum(is.na(input.history))
## [1] 2969
site.covar <- data.frame(Site=1:nrow(input.data),
Pasture=car::recode(input.data$X__18,
"1='Modified'; 0='Native';"))
head(site.covar)
## Site Pasture
## 1 1 Modified
## 2 2 Modified
## 3 3 Modified
## 4 4 Modified
## 5 5 Modified
## 6 6 Modified
#Format the capture history to be used by RMark
input.history <- data.frame(freq=1,
ch=apply(input.history,1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 00000000NA000000
## 2 1 00000NA000000000
## 3 1 00000NA000000000
## 4 1 0NANA0NANA0NANA00NA0NANA
## 5 1 0NANANANANANANANANANANA0NANA
## 6 1 0NANANANANANANANANANANA0NANA
# 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 00000000.000000
## 2 1 00000.000000000
## 3 1 00000.000000000
## 4 1 0..0..0..00.0..
## 5 1 0...........0..
## 6 1 0...........0..
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## freq ch Site Pasture
## 1 1 00000000.000000 1 Modified
## 2 1 00000.000000000 2 Modified
## 3 1 00000.000000000 3 Modified
## 4 1 0..0..0..00.0.. 4 Modified
## 5 1 0...........0.. 5 Modified
## 6 1 0...........0.. 6 Modified
#Create the RMark data structure
#Five Seasons, with 3 visits per season
max.visit.per.year <- 3
n.season <- 5
skink.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(skink.data)
## Length Class Mode
## data 4 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 352 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 5 -none- numeric
## time.intervals 14 -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
# add visit level covariates
# In this case none
skink.ddl <- make.design.data(skink.data)
skink.ddl
## $Psi
## par.index model.index group age time Age Time
## 1 1 1 1 0 1 0 0
##
## $Epsilon
## par.index model.index group age time Age Time
## 1 1 2 1 0 1 0 0
## 2 2 3 1 1 2 1 1
## 3 3 4 1 2 3 2 2
## 4 4 5 1 3 4 3 3
##
## $Gamma
## par.index model.index group age time Age Time
## 1 1 6 1 0 1 0 0
## 2 2 7 1 1 2 1 1
## 3 3 8 1 2 3 2 2
## 4 4 9 1 3 4 3 3
##
## $p
## par.index model.index group time session Time
## 1 1 10 1 1 1 0
## 2 2 11 1 2 1 1
## 3 3 12 1 3 1 2
## 4 4 13 1 1 2 0
## 5 5 14 1 2 2 1
## 6 6 15 1 3 2 2
## 7 7 16 1 1 3 0
## 8 8 17 1 2 3 1
## 9 9 18 1 3 3 2
## 10 10 19 1 1 4 0
## 11 11 20 1 2 4 1
## 12 12 21 1 3 4 2
## 13 13 22 1 1 5 0
## 14 14 23 1 2 5 1
## 15 15 24 1 3 5 2
##
## $pimtypes
## $pimtypes$Psi
## $pimtypes$Psi$pim.type
## [1] "all"
##
##
## $pimtypes$Epsilon
## $pimtypes$Epsilon$pim.type
## [1] "all"
##
##
## $pimtypes$Gamma
## $pimtypes$Gamma$pim.type
## [1] "all"
##
##
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
# Fit some models.
# We start with the first formulation in terms of psi, gamma, epsilon, p (RDOccupEG)
# Note that formula DO NOT HAVE AN = SIGN
mod.fit <- RMark::mark(skink.data, ddl = skink.ddl,
model="RDOccupEG",
model.parameters=list(
Psi =list(formula=~1),
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: 1767.015
## AICc : 1775.043
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.4284161 0.1318158 -0.6867751 -0.170057
## Epsilon:(Intercept) -2.1970569 0.2028714 -2.5946848 -1.799429
## Gamma:(Intercept) -2.6132991 0.1864902 -2.9788199 -2.247778
## p:(Intercept) 0.7961691 0.1077378 0.5850030 1.007335
##
##
## Real Parameter Psi
##
## 1
## 0.3945046
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.1000151 0.1000151 0.1000151 0.1000151
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.0682874 0.0682874 0.0682874 0.0682874
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:2
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:3
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:4
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:5
## 1 2 3
## 0.6891544 0.6891544 0.6891544
summary(mod.fit)
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~1)
##
## Npar : 4
## -2lnL: 1767.015
## AICc : 1775.043
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.4284161 0.1318158 -0.6867751 -0.170057
## Epsilon:(Intercept) -2.1970569 0.2028714 -2.5946848 -1.799429
## Gamma:(Intercept) -2.6132991 0.1864902 -2.9788199 -2.247778
## p:(Intercept) 0.7961691 0.1077378 0.5850030 1.007335
##
##
## Real Parameter Psi
##
## 1
## 0.3945046
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.1000151 0.1000151 0.1000151 0.1000151
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.0682874 0.0682874 0.0682874 0.0682874
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:2
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:3
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:4
## 1 2 3
## 0.6891544 0.6891544 0.6891544
##
## Session:5
## 1 2 3
## 0.6891544 0.6891544 0.6891544
# Estimates of real initial occupancy, detection, local exctinction,
# and local colonization probabilities
mod.fit$results$real
## estimate se lcl ucl fixed note
## Psi g1 a0 t1 0.3945046 0.0314869 0.3347509 0.4575879
## Epsilon g1 a0 t1 0.1000151 0.0182609 0.0694813 0.1419206
## Gamma g1 a0 t1 0.0682874 0.0118653 0.0483919 0.0955413
## p g1 s1 t1 0.6891544 0.0230797 0.6422178 0.7324983
# Estimates of regression coefficients
mod.fit$results$beta
## estimate se lcl ucl
## Psi:(Intercept) -0.4284161 0.1318158 -0.6867751 -0.170057
## Epsilon:(Intercept) -2.1970569 0.2028714 -2.5946848 -1.799429
## Gamma:(Intercept) -2.6132991 0.1864902 -2.9788199 -2.247778
## p:(Intercept) 0.7961691 0.1077378 0.5850030 1.007335
# Estimates of various derived parameters for estimated occpancy in each year,
# estimated ratio of in occupancy from year to year, and estimated log
# of the ratio of odds in occupancy from year to year
mod.fit$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.3945046 0.03148694 0.3327902 0.4562190
## 2 0.3963959 0.02657622 0.3443065 0.4484853
## 3 0.3979689 0.02558104 0.3478301 0.4481077
## 4 0.3992771 0.02719514 0.3459747 0.4525796
## 5 0.4003652 0.02998981 0.3415852 0.4591452
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.004794 0.02813991 0.9496399 1.059948
## 2 1.003968 0.02319011 0.9585156 1.049421
## 3 1.003287 0.01914292 0.9657672 1.040807
## 4 1.002725 0.01582392 0.9717102 1.033740
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 1.007942 0.04651943 0.9167644 1.099121
## 2 1.006591 0.03853687 0.9310591 1.082124
## 3 1.005472 0.03195039 0.9428495 1.068095
## 4 1.004545 0.02650735 0.9525902 1.056499
mod.fit$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.3945046 0.03148694 0.3327902 0.4562190
## 2 0.3963959 0.02657622 0.3443065 0.4484853
## 3 0.3979689 0.02558104 0.3478301 0.4481077
## 4 0.3992771 0.02719514 0.3459747 0.4525796
## 5 0.4003652 0.02998981 0.3415852 0.4591452
mod.fit$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 1.004794 0.02813991 0.9496399 1.059948
## 2 1.003968 0.02319011 0.9585156 1.049421
## 3 1.003287 0.01914292 0.9657672 1.040807
## 4 1.002725 0.01582392 0.9717102 1.033740
mod.fit$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 1.007942 0.04651943 0.9167644 1.099121
## 2 1.006591 0.03853687 0.9310591 1.082124
## 3 1.005472 0.03195039 0.9428495 1.068095
## 4 1.004545 0.02650735 0.9525902 1.056499
# Plot derived occupancy probabilities
plotdata = mod.fit$results$derived$"psi Probability Occupied"
plotdata$Season = c(1:5)
head(plotdata)
## estimate se lcl ucl Season
## 1 0.3945046 0.03148694 0.3327902 0.4562190 1
## 2 0.3963959 0.02657622 0.3443065 0.4484853 2
## 3 0.3979689 0.02558104 0.3478301 0.4481077 3
## 4 0.3992771 0.02719514 0.3459747 0.4525796 4
## 5 0.4003652 0.02998981 0.3415852 0.4591452 5
ggplot(data=plotdata, aes(x=Season,y=estimate))+
ggtitle("Estimated occupancy over time")+
geom_point(position=position_dodge(w=0.2))+
geom_line(position=position_dodge(w=0.2))+
ylim(0,1)+
geom_errorbar(aes(ymin=lcl, ymax=ucl), width=.1,position=position_dodge(w=0.2))+
scale_x_continuous(breaks=1:10)
cleanup(ask=FALSE)