# 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
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="MissingData",
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 NA 0 0
## 3 0 NA NA 0 0
## 4 1 NA NA 1 0
## 5 0 NA NA 0 0
## 6 0 0 NA 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 01NA00
## 3 1 0NANA00
## 4 1 1NANA10
## 5 1 0NANA00
## 6 1 00NA00
# 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 01.00
## 3 1 0..00
## 4 1 1..10
## 5 1 0..00
## 6 1 00.00
# 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
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")
#browser()
#fit <- myoccMod(model=list(as.formula(paste("psi",x$psi)),
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: 137.613
## AICc : 141.9463
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.9444716 0.3316785 -1.5945615 -0.2943816
## Psi:(Intercept) -0.0226465 0.4636926 -0.9314841 0.8861911
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.279998 0.279998 0.279998 0.279998 0.279998
##
##
## Real Parameter Psi
## 1
## 0.4943386
##
##
## ***** Starting ~time ~1 2
##
## Output summary for Occupancy model
## Name : p(~time)Psi(~1)
##
## Npar : 6
## -2lnL: 128.2942
## AICc : 142.9192
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.2900290 0.5981794 -2.4624605 -0.1175974
## p:time2 -1.4899440 1.1761765 -3.7952500 0.8153621
## p:time3 0.9155972 0.7818234 -0.6167766 2.4479711
## p:time4 1.0347380 0.7392732 -0.4142374 2.4837135
## p:time5 0.6688581 0.7605476 -0.8218153 2.1595314
## Psi:(Intercept) -0.0994084 0.4357577 -0.9534935 0.7546767
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.2158479 0.058416 0.4074706 0.4365216 0.3495152
##
##
## Real Parameter Psi
## 1
## 0.4751684
##
##
## ***** Starting ~Effort ~1 3
##
## Output summary for Occupancy model
## Name : p(~Effort)Psi(~1)
##
## Npar : 3
## -2lnL: 130.5309
## AICc : 137.2166
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.820028 0.5149302 -2.8292912 -0.8107647
## p:Effortd2 1.397350 0.5669018 0.2862225 2.5084776
## Psi:(Intercept) -0.086860 0.4397341 -0.9487388 0.7750187
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.1394305 0.1394305 0.3958761 0.3958761 0.3958761
##
##
## Real Parameter Psi
## 1
## 0.4782986
# 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: 128.2942
## AICc : 142.9192
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.2900290 0.5981794 -2.4624605 -0.1175974
## p:time2 -1.4899440 1.1761765 -3.7952500 0.8153621
## p:time3 0.9155972 0.7818234 -0.6167766 2.4479711
## p:time4 1.0347380 0.7392732 -0.4142374 2.4837135
## p:time5 0.6688581 0.7605476 -0.8218153 2.1595314
## Psi:(Intercept) -0.0994084 0.4357577 -0.9534935 0.7546767
##
##
## Real Parameter p
## 1 2 3 4 5
## 0.2158479 0.058416 0.4074706 0.4365216 0.3495152
##
##
## Real Parameter Psi
## 1
## 0.4751684
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## p g1 a0 t1 0.2158479 0.1012464 0.0785321 0.4706345
## p g1 a1 t2 0.0584160 0.0573941 0.0079614 0.3241443
## p g1 a2 t3 0.4074706 0.1443567 0.1756182 0.6894305
## p g1 a3 t4 0.4365216 0.1331347 0.2114569 0.6911677
## p g1 a4 t5 0.3495152 0.1270821 0.1522924 0.6164226
## Psi g1 a0 t1 0.4751684 0.1086707 0.2781828 0.6801969
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## p:(Intercept) -1.2900290 0.5981794 -2.4624605 -0.1175974
## p:time2 -1.4899440 1.1761765 -3.7952500 0.8153621
## p:time3 0.9155972 0.7818234 -0.6167766 2.4479711
## p:time4 1.0347380 0.7392732 -0.4142374 2.4837135
## p:time5 0.6688581 0.7605476 -0.8218153 2.1595314
## Psi:(Intercept) -0.0994084 0.4357577 -0.9534935 0.7546767
model.fits[[model.number]]$results$derived
## $Occupancy
## estimate se lcl ucl
## 1 0.4751684 0.1086707 0.2781828 0.6801969
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 6 6 0.4751684 0.1086707 0.2781828
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.6801969 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.2158479 0.1012464 0.0785321
## p g1 a1 t2 2 2 0.0584160 0.0573941 0.0079614
## p g1 a2 t3 3 3 0.4074706 0.1443567 0.1756182
## p g1 a3 t4 4 4 0.4365216 0.1331347 0.2114569
## p g1 a4 t5 5 5 0.3495152 0.1270821 0.1522924
## ucl fixed note group age time Age Time
## p g1 a0 t1 0.4706345 1 0 1 0 0
## p g1 a1 t2 0.3241443 1 1 2 1 1
## p g1 a2 t3 0.6894305 1 2 3 2 2
## p g1 a3 t4 0.6911677 1 3 4 3 3
## p g1 a4 t5 0.6164226 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 137.2166 0.000000 0.86825465 -39.48951
## 1 p(~1)Psi(~1) 2 141.9463 4.729639 0.08158664 -32.40749
## 2 p(~time)Psi(~1) 6 142.9192 5.702586 0.05015871 -41.72621
names(model.set)
## [1] "m...001" "m...002" "m...003" "model.table"
model.set$model.table
## p Psi model npar AICc DeltaAICc weight
## 3 ~Effort ~1 p(~Effort)Psi(~1) 3 137.2166 0.000000 0.86825465
## 1 ~1 ~1 p(~1)Psi(~1) 2 141.9463 4.729639 0.08158664
## 2 ~time ~1 p(~time)Psi(~1) 6 142.9192 5.702586 0.05015871
## Deviance
## 3 -39.48951
## 1 -32.40749
## 2 -41.72621
# 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.4943386 0.1159083 0.2826237
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.7081035 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.4751684 0.1086707 0.2781828
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.6801969 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.4782986 0.1097264 0.2791385
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.6846055 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.4794503 0.1102827 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 ucl
## p g1 a0 t1 1 1 0.279998 0.0668661 0.1687431 0.4269315
## p g1 a1 t2 2 1 0.279998 0.0668661 0.1687431 0.4269315
## p g1 a2 t3 3 1 0.279998 0.0668661 0.1687431 0.4269315
## p g1 a3 t4 4 1 0.279998 0.0668661 0.1687431 0.4269315
## p g1 a4 t5 5 1 0.279998 0.0668661 0.1687431 0.4269315
## fixed note group age time Age Time
## p g1 a0 t1 1 0 1 0 0
## p g1 a1 t2 1 1 2 1 1
## p g1 a2 t3 1 2 3 2 2
## p g1 a3 t4 1 3 4 3 3
## p g1 a4 t5 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.2158479 0.1012464 0.0785321
## p g1 a1 t2 2 2 0.0584160 0.0573941 0.0079614
## p g1 a2 t3 3 3 0.4074706 0.1443567 0.1756182
## p g1 a3 t4 4 4 0.4365216 0.1331347 0.2114569
## p g1 a4 t5 5 5 0.3495152 0.1270821 0.1522924
## ucl fixed note group age time Age Time
## p g1 a0 t1 0.4706345 1 0 1 0 0
## p g1 a1 t2 0.3241443 1 1 2 1 1
## p g1 a2 t3 0.6894305 1 2 3 2 2
## p g1 a3 t4 0.6911677 1 3 4 3 3
## p g1 a4 t5 0.6164226 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.1394305 0.0617863 0.0557617
## p g1 a1 t2 2 1 0.1394305 0.0617863 0.0557617
## p g1 a2 t3 3 2 0.3958761 0.0921396 0.2354481
## p g1 a3 t4 4 2 0.3958761 0.0921396 0.2354481
## p g1 a4 t5 5 2 0.3958761 0.0921396 0.2354481
## ucl fixed note group age time Age Time
## p g1 a0 t1 0.3077276 1 0 1 0 0
## p g1 a1 t2 0.3077276 1 1 2 1 1
## p g1 a2 t3 0.5823539 1 2 3 2 2
## p g1 a3 t4 0.5823539 1 3 4 3 3
## p g1 a4 t5 0.5823539 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.1547319 0.07657945 1 0 1 0
## p g1 a1 t2 2 0.1468354 0.07570386 1 1 2 1
## p g1 a2 t3 3 0.3870036 0.09901919 1 2 3 2
## p g1 a3 t4 4 0.3884608 0.09873430 1 3 4 3
## p g1 a4 t5 5 0.3840966 0.09803385 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)