# Blue Ridge Salamander
# Single Species Single Season Occupancy
# Fitting a single model
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
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 <- readxl::read_excel("../salamander.xls",
sheet="CompleteData",
na="-",
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 and create a capture history
input.history <- data.frame(freq=1,
ch=apply(input.data[,1:5],1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 00011
## 2 1 01000
## 3 1 01000
## 4 1 11110
## 5 1 00100
## 6 1 00100
# 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 00011
## 2 1 01000
## 3 1 01000
## 4 1 11110
## 5 1 00100
## 6 1 00100
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 39
sal.data <- process.data(data=input.history,
model="Occupancy")
summary(sal.data)
## Length Class Mode
## data 2 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 39 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 5 -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("Occupancy", check=TRUE)
## [1] "p" "Psi"
# Fit a model
# Note that formula HAVE AN = SIGN
mod.fit <- RMark::mark(sal.data,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~1)
)
)
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 161.7586
## AICc : 166.0919
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.0525847 0.3008626 -1.642275 -0.462894
## Psi:(Intercept) 0.3831084 0.5086094 -0.613766 1.379983
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.2587291 0.2587291 0.2587291 0.2587291 0.2587291
##
##
## Real Parameter Psi
## 1
## 0.5946226
summary(mod.fit)
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 161.7586
## AICc : 166.0919
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.0525847 0.3008626 -1.642275 -0.462894
## Psi:(Intercept) 0.3831084 0.5086094 -0.613766 1.379983
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.2587291 0.2587291 0.2587291 0.2587291 0.2587291
##
##
## Real Parameter Psi
## 1
## 0.5946226
# Look the objects returned in more details
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" "profile.int"
## [22] "simplify" "model.parameters" "results"
## [25] "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"
# look at estimates on beta and original scale
mod.fit$results$beta # on the logit scale
## estimate se lcl ucl
## p:(Intercept) -1.0525847 0.3008626 -1.642275 -0.462894
## Psi:(Intercept) 0.3831084 0.5086094 -0.613766 1.379983
mod.fit$results$real# on the regular 0-1 scale for each site
## estimate se lcl ucl fixed note
## p g1 a0 t1 0.2587291 0.0577019 0.1621557 0.3862995
## Psi g1 a0 t1 0.5946226 0.1225985 0.3512006 0.7989882
# derived variabldes is the occupancy probability
names(mod.fit$results$derived)
## [1] "Occupancy"
mod.fit$results$derived$Occupancy
## estimate se lcl ucl
## 1 0.5946226 0.1225985 0.3512006 0.7989882
# alternatively
get.real(mod.fit, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 6 2 0.5946226 0.1225985 0.3512006
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.7989882 1 0 1 0 0
get.real(mod.fit, "Psi", pim=TRUE)
## 1
## 0.5946226
get.real(mod.fit, "p", se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 a0 t1 1 1 0.2587291 0.0577019 0.1621557
## p g1 a1 t2 2 1 0.2587291 0.0577019 0.1621557
## p g1 a2 t3 3 1 0.2587291 0.0577019 0.1621557
## p g1 a3 t4 4 1 0.2587291 0.0577019 0.1621557
## p g1 a4 t5 5 1 0.2587291 0.0577019 0.1621557
## ucl fixed note group age time Age Time
## p g1 a0 t1 0.3862995 1 0 1 0 0
## p g1 a1 t2 0.3862995 1 1 2 1 1
## p g1 a2 t3 0.3862995 1 2 3 2 2
## p g1 a3 t4 0.3862995 1 3 4 3 3
## p g1 a4 t5 0.3862995 1 4 5 4 4
#######################################
# Fit a model with p constant over time
mod.fit2 <- RMark::mark(sal.data,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~time)
)
)
##
## Output summary for Occupancy model
## Name : p(~time)Psi(~1)
##
## Npar : 6
## -2lnL: 155.7144
## AICc : 170.3394
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.5376880 0.5804914 -2.6754511 -0.3999249
## p:time2 -0.3400080 0.8295311 -1.9658889 1.2858730
## p:time3 1.1237209 0.7019705 -0.2521412 2.4995830
## p:time4 0.9350629 0.7068468 -0.4503569 2.3204827
## p:time5 0.5191254 0.7287348 -0.9091947 1.9474456
## Psi:(Intercept) 0.3222754 0.4825851 -0.6235913 1.2681422
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.1768716 0.1326537 0.3979613 0.3537433 0.2653075
##
##
## Real Parameter Psi
## 1
## 0.5798787
summary(mod.fit2)
## Output summary for Occupancy model
## Name : p(~time)Psi(~1)
##
## Npar : 6
## -2lnL: 155.7144
## AICc : 170.3394
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.5376880 0.5804914 -2.6754511 -0.3999249
## p:time2 -0.3400080 0.8295311 -1.9658889 1.2858730
## p:time3 1.1237209 0.7019705 -0.2521412 2.4995830
## p:time4 0.9350629 0.7068468 -0.4503569 2.3204827
## p:time5 0.5191254 0.7287348 -0.9091947 1.9474456
## Psi:(Intercept) 0.3222754 0.4825851 -0.6235913 1.2681422
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.1768716 0.1326537 0.3979613 0.3537433 0.2653075
##
##
## Real Parameter Psi
## 1
## 0.5798787
get.real(mod.fit2, "p", se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 a0 t1 1 1 0.1768716 0.0845126 0.0644376
## p g1 a1 t2 2 2 0.1326537 0.0740541 0.0415184
## p g1 a2 t3 3 3 0.3979613 0.1190041 0.1998064
## p g1 a3 t4 4 4 0.3537433 0.1137000 0.1711580
## p g1 a4 t5 5 5 0.2653075 0.1010181 0.1156439
## ucl fixed note group age time Age Time
## p g1 a0 t1 0.4013304 1 0 1 0 0
## p g1 a1 t2 0.3506509 1 1 2 1 1
## p g1 a2 t3 0.6363532 1 2 3 2 2
## p g1 a3 t4 0.5919885 1 3 4 3 3
## p g1 a4 t5 0.4993047 1 4 5 4 4
####################################################
# fit a model with p equal in first two visits and last 3 visits
# This is a survey-level covariate and so you need to modify the ddl
#
# add a survey level covariate. by adding a column to the design matrix
sal.ddl <- make.design.data(sal.data)
sal.ddl
## $p
## par.index model.index group age time Age Time
## 1 1 1 1 0 1 0 0
## 2 2 2 1 1 2 1 1
## 3 3 3 1 2 3 2 2
## 4 4 4 1 3 4 3 3
## 5 5 5 1 4 5 4 4
##
## $Psi
## par.index model.index group age time Age Time
## 1 1 6 1 0 1 0 0
##
## $pimtypes
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
##
##
## $pimtypes$Psi
## $pimtypes$Psi$pim.type
## [1] "all"
sal.ddl$p$Effort <- factor(c("d1","d1","d2","d2","d2"))
sal.ddl
## $p
## par.index model.index group age time Age Time Effort
## 1 1 1 1 0 1 0 0 d1
## 2 2 2 1 1 2 1 1 d1
## 3 3 3 1 2 3 2 2 d2
## 4 4 4 1 3 4 3 3 d2
## 5 5 5 1 4 5 4 4 d2
##
## $Psi
## par.index model.index group age time Age Time
## 1 1 6 1 0 1 0 0
##
## $pimtypes
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
##
##
## $pimtypes$Psi
## $pimtypes$Psi$pim.type
## [1] "all"
mod.fit3 <- RMark::mark(sal.data, ddl=sal.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~Effort)
)
)
##
## Output summary for Occupancy model
## Name : p(~Effort)Psi(~1)
##
## Npar : 3
## -2lnL: 156.8176
## AICc : 163.5033
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.7043681 0.4484046 -2.5832411 -0.8254952
## p:Effortd2 1.0281446 0.4867145 0.0741842 1.9821050
## Psi:(Intercept) 0.3357008 0.4881036 -0.6209822 1.2923839
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.1538956 0.1538956 0.3371047 0.3371047 0.3371047
##
##
## Real Parameter Psi
## 1
## 0.5831458
summary(mod.fit3)
## Output summary for Occupancy model
## Name : p(~Effort)Psi(~1)
##
## Npar : 3
## -2lnL: 156.8176
## AICc : 163.5033
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.7043681 0.4484046 -2.5832411 -0.8254952
## p:Effortd2 1.0281446 0.4867145 0.0741842 1.9821050
## Psi:(Intercept) 0.3357008 0.4881036 -0.6209822 1.2923839
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.1538956 0.1538956 0.3371047 0.3371047 0.3371047
##
##
## Real Parameter Psi
## 1
## 0.5831458
get.real(mod.fit3, "p", se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 a0 t1 1 1 0.1538956 0.0583875 0.0702248
## p g1 a1 t2 2 1 0.1538956 0.0583875 0.0702248
## p g1 a2 t3 3 2 0.3371047 0.0767915 0.2059100
## p g1 a3 t4 4 2 0.3371047 0.0767915 0.2059100
## p g1 a4 t5 5 2 0.3371047 0.0767915 0.2059100
## ucl fixed note group age time Age Time
## p g1 a0 t1 0.3045984 1 0 1 0 0
## p g1 a1 t2 0.3045984 1 1 2 1 1
## p g1 a2 t3 0.4993277 1 2 3 2 2
## p g1 a3 t4 0.4993277 1 3 4 3 3
## p g1 a4 t5 0.4993277 1 4 5 4 4
##############################################
# Collect models
model.set <- RMark::collect.models( type="Occupancy")
model.set
## model npar AICc DeltaAICc weight Deviance
## 3 p(~Effort)Psi(~1) 3 163.5033 0.000000 0.7651935 26.11443
## 1 p(~1)Psi(~1) 2 166.0919 2.588649 0.2097265 31.05546
## 2 p(~time)Psi(~1) 6 170.3394 6.836116 0.0250800 25.01126
model.set <- RMark::collect.models( type="Occupancy")
model.set
## model npar AICc DeltaAICc weight Deviance
## 3 p(~Effort)Psi(~1) 3 163.5033 0.000000 0.7651935 26.11443
## 1 p(~1)Psi(~1) 2 166.0919 2.588649 0.2097265 31.05546
## 2 p(~time)Psi(~1) 6 170.3394 6.836116 0.0250800 25.01126
names(model.set)
## [1] "mod.fit" "mod.fit2" "mod.fit3" "model.table"
model.set$model.table
## p Psi model npar AICc DeltaAICc weight Deviance
## 3 ~Effort ~1 p(~Effort)Psi(~1) 3 163.5033 0.000000 0.7651935 26.11443
## 1 ~1 ~1 p(~1)Psi(~1) 2 166.0919 2.588649 0.2097265 31.05546
## 2 ~time ~1 p(~time)Psi(~1) 6 170.3394 6.836116 0.0250800 25.01126
# model averaged values
get.real(mod.fit , "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 6 2 0.5946226 0.1225985 0.3512006
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.7989882 1 0 1 0 0
get.real(mod.fit2, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 6 6 0.5798787 0.1175671 0.3489651
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.7804245 1 0 1 0 0
get.real(mod.fit3, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 6 3 0.5831458 0.1186515 0.3495581
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.7845504 1 0 1 0 0
Psi.ma <- RMark::model.average(model.set, param="Psi")
Psi.ma
## par.index estimate se fixed note group age time
## Psi g1 a0 t1 6 0.5854709 0.1195573 1 0 1
## Age Time
## Psi g1 a0 t1 0 0
# cleanup
cleanup(ask=FALSE)