# Single Species, Multi Season Occupancy analyais
# Salamanders over four seasons with covariates on elevation and stream type.
# We will use only the last 2 years (unsure why)
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
# 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 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("..","SalamanderMultiSeason.xls"),
sheet="Detections",
skip=2, na="-",
col_names=FALSE)
head(input.data)
## # A tibble: 6 x 22
## 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 1.00 0 1.00 0 1.00 1.00 0 1.00 1.00 1.00 0
## 2 2.00 0 0 0 0 0 0 0 0 0 0 0
## 3 3.00 0 1.00 0 1.00 1.00 1.00 0 0 0 1.00 1.00
## 4 4.00 1.00 1.00 1.00 1.00 0 0 0 1.00 1.00 1.00 1.00
## 5 5.00 1.00 1.00 1.00 1.00 0 1.00 0 1.00 0 0 0
## 6 6.00 0 0 0 0 0 0 0 0 1.00 0 0
## # ... with 10 more variables: X__13 <dbl>, X__14 <dbl>, X__15 <dbl>, X__16
## # <dbl>, X__17 <dbl>, X__18 <dbl>, X__19 <dbl>, X__20 <lgl>, X__21
## # <dbl>, X__22 <dbl>
input.history <- input.data[, 10:19]
head(input.history)
## # A tibble: 6 x 10
## X__10 X__11 X__12 X__13 X__14 X__15 X__16 X__17 X__18 X__19
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00 1.00 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 1.00 1.00 0 0 0 0 1.00 0 0
## 4 1.00 1.00 1.00 1.00 0 1.00 1.00 1.00 1.00 0
## 5 0 0 0 0 1.00 0 0 0 0 0
## 6 1.00 0 0 0 0 0 0 1.00 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] 10
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 0
site.covar <- data.frame(Site=1:nrow(input.data),
Elevation=unlist(input.data[,21]),
Prox=car::recode(unlist(input.data[,22]),
"1='Yes'; 0='No';"))
row.names(site.covar) <- NULL
head(site.covar)
## Site Elevation Prox
## 1 1 841.248 Yes
## 2 2 853.440 No
## 3 3 780.288 Yes
## 4 4 707.136 Yes
## 5 5 670.560 No
## 6 6 670.560 Yes
#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 1100000000
## 2 1 0000000000
## 3 1 0110000100
## 4 1 1111011110
## 5 1 0000100000
## 6 1 1000000100
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## freq ch Site Elevation Prox
## 1 1 1100000000 1 841.248 Yes
## 2 1 0000000000 2 853.440 No
## 3 1 0110000100 3 780.288 Yes
## 4 1 1111011110 4 707.136 Yes
## 5 1 0000100000 5 670.560 No
## 6 1 1000000100 6 670.560 Yes
#Create the RMark data structure
#Two Seasons, with 5 visits per season
max.visit.per.year <- 5
n.season <- 2
salamander.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(salamander.data)
## Length Class Mode
## data 5 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 39 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 2 -none- numeric
## time.intervals 9 -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
# Fit a models.
# We start with the first formulation in terms of psi, gamma, epsilon, p
mod.fit <- RMark::mark(salamander.data,
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: 316.773
## AICc : 325.321
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.1254139 0.4123723 -0.6828358 0.9336637
## Epsilon:(Intercept) -1.8185768 1.2131230 -4.1962979 0.5591442
## Gamma:(Intercept) -1.2370967 0.8684573 -2.9392731 0.4650796
## p:(Intercept) -0.8921721 0.2030460 -1.2901422 -0.4942019
##
##
## Real Parameter Psi
##
## 1
## 0.5313125
##
##
## Real Parameter Epsilon
##
## 1
## 0.1396047
##
##
## Real Parameter Gamma
##
## 1
## 0.2249417
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
##
## Session:2
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
summary(mod.fit)
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~1)
##
## Npar : 4
## -2lnL: 316.773
## AICc : 325.321
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.1254139 0.4123723 -0.6828358 0.9336637
## Epsilon:(Intercept) -1.8185768 1.2131230 -4.1962979 0.5591442
## Gamma:(Intercept) -1.2370967 0.8684573 -2.9392731 0.4650796
## p:(Intercept) -0.8921721 0.2030460 -1.2901422 -0.4942019
##
##
## Real Parameter Psi
##
## 1
## 0.5313125
##
##
## Real Parameter Epsilon
##
## 1
## 0.1396047
##
##
## Real Parameter Gamma
##
## 1
## 0.2249417
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
##
## Session:2
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
# 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.5313125 0.1026888 0.3356287 0.7178180
## Epsilon g1 a0 t1 0.1396047 0.1457146 0.0148280 0.6362545
## Gamma g1 a0 t1 0.2249417 0.1514094 0.0502460 0.6142185
## p g1 s1 t1 0.2906618 0.0418635 0.2158287 0.3789042
# Estimates of regression coefficients
mod.fit$results$beta
## estimate se lcl ucl
## Psi:(Intercept) 0.1254139 0.4123723 -0.6828358 0.9336637
## Epsilon:(Intercept) -1.8185768 1.2131230 -4.1962979 0.5591442
## Gamma:(Intercept) -1.2370967 0.8684573 -2.9392731 0.4650796
## p:(Intercept) -0.8921721 0.2030460 -1.2901422 -0.4942019
# 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
# of the ratio of odds of occupancy from year to year
mod.fit$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.5313125 0.1026888 0.3300425 0.7325824
## 2 0.5625661 0.1038643 0.3589920 0.7661402
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.058824 0.2007517 0.6653502 1.452297
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 1.134474 0.4746753 0.2041103 2.064838
mod.fit$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.5313125 0.1026888 0.3300425 0.7325824
## 2 0.5625661 0.1038643 0.3589920 0.7661402
mod.fit$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 1.058824 0.2007517 0.6653502 1.452297
mod.fit$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 1.134474 0.4746753 0.2041103 2.064838
cleanup(ask = FALSE)