# Blue Ridge Salamander with some added missing values
# Notice how we change the NA to . in the capture history
# Single Species Single Season Occupancy
# Fitting multiple models and model averaging with a general structure
# 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",
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
# 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"
# What are the parameter names for Single Season Single Species models
setup.parameters("Occupancy", check=TRUE)
## [1] "p" "Psi"
# Get the list of models. NOtice NO equal signs here
model.list.csv <- textConnection("
p, Psi
~1, ~1
~time, ~1
~Effort, ~1")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list$model.number <- 1:nrow(model.list)
model.list
## p Psi model.number
## 1 ~1 ~1 1
## 2 ~time ~1 2
## 3 ~Effort ~1 3
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing
rm(list=myobs)
model.fits <- plyr::dlply(model.list, "model.number", function(x,input.data, input.ddl){
cat("\n\n***** Starting ", unlist(x), "\n")
fit <- RMark::mark(input.data, ddl=input.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=as.formula(eval(x$Psi))),
p =list(formula=as.formula(eval(x$p)))
)
#,brief=TRUE,output=FALSE, delete=TRUE
#,invisible=TRUE,output=TRUE # set for debugging
)
mnumber <- paste("m...",formatC(x$model.number, width = 3, format = "d", flag = "0"),sep="")
assign( mnumber, fit, envir=.GlobalEnv)
#browser()
fit
},input.data=sal.data, input.ddl=sal.ddl)
##
##
## ***** Starting ~1 ~1 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
##
##
## ***** Starting ~time ~1 2
##
## 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
##
##
## ***** Starting ~Effort ~1 3
##
## 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
# examine individula model results
model.number <-2
summary(model.fits[[model.number]])
## 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
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## p g1 a0 t1 0.1768716 0.0845126 0.0644376 0.4013304
## p g1 a1 t2 0.1326537 0.0740541 0.0415184 0.3506509
## p g1 a2 t3 0.3979613 0.1190041 0.1998064 0.6363532
## p g1 a3 t4 0.3537433 0.1137000 0.1711580 0.5919885
## p g1 a4 t5 0.2653075 0.1010181 0.1156439 0.4993047
## Psi g1 a0 t1 0.5798787 0.1175671 0.3489651 0.7804245
model.fits[[model.number]]$results$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
model.fits[[model.number]]$results$derived
## $Occupancy
## estimate se lcl ucl
## 1 0.5798787 0.1175671 0.3489651 0.7804245
get.real(model.fits[[model.number]], "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(model.fits[[model.number]], "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
# collect models and make AICc table
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] "m...001" "m...002" "m...003" "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(model.set[[1]], "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(model.set[[2]], "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(model.set[[3]], "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
get.real(model.set[[1]], "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
get.real(model.set[[2]], "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
get.real(model.set[[3]], "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
p.ma <- RMark::model.average(model.set, param="p")
p.ma
## par.index estimate se fixed note group age time Age
## p g1 a0 t1 1 0.1764582 0.07276756 1 0 1 0
## p g1 a1 t2 2 0.1753492 0.07280396 1 1 2 1
## p g1 a2 t3 3 0.3221935 0.08200370 1 2 3 2
## p g1 a3 t4 4 0.3210846 0.08107655 1 3 4 3
## p g1 a4 t5 5 0.3188666 0.08094115 1 4 5 4
## Time
## p g1 a0 t1 0
## p g1 a1 t2 1
## p g1 a2 t3 2
## p g1 a3 t4 3
## p g1 a4 t5 4
# cleanup
cleanup(ask=FALSE)