# Bullfrogs
# Single Species Single Season Occupancy - Mixture model
# Fitting a single model
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 <- read.table(file.path("..","MARK","bullfrog.inp"), as.is=TRUE, colClasses="character")
head(input.data)
## V1 V2
## 1 0001 1;
## 2 0011 1;
## 3 0000 1;
## 4 0000 1;
## 5 0000 1;
## 6 0000 1;
# Extract the history records and create a capture history
input.history <- data.frame(freq=1,
ch=input.data$V1, stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 0001
## 2 1 0011
## 3 1 0000
## 4 1 0000
## 5 1 0000
## 6 1 0000
# 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 0001
## 2 1 0011
## 3 1 0000
## 4 1 0000
## 5 1 0000
## 6 1 0000
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 500
bullfrog.data <- process.data(data=input.history,
model="Occupancy")
summary(bullfrog.data)
## Length Class Mode
## data 2 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 500 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 4 -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("Occupamcut", check=TRUE)
## character(0)
# Fit a model without heterogeneity
mod.fit <- RMark::mark(bullfrog.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: 1887.268
## AICc : 1891.292
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2084764 0.0740655 0.0633079 0.3536449
## Psi:(Intercept) -0.1413071 0.0939972 -0.3255417 0.0429275
##
##
## Real Parameter p
## 1 2 3 4
## 0.5519311 0.5519311 0.5519311 0.5519311
##
##
## Real Parameter Psi
## 1
## 0.4647319
summary(mod.fit)
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 1887.268
## AICc : 1891.292
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2084764 0.0740655 0.0633079 0.3536449
## Psi:(Intercept) -0.1413071 0.0939972 -0.3255417 0.0429275
##
##
## Real Parameter p
## 1 2 3 4
## 0.5519311 0.5519311 0.5519311 0.5519311
##
##
## Real Parameter Psi
## 1
## 0.4647319
# 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) 0.2084764 0.0740655 0.0633079 0.3536449
## Psi:(Intercept) -0.1413071 0.0939972 -0.3255417 0.0429275
mod.fit$results$real# on the regular 0-1 scale for each site
## estimate se lcl ucl fixed note
## p g1 a0 t1 0.5519311 0.0183166 0.5158217 0.5875012
## Psi g1 a0 t1 0.4647319 0.0233824 0.4193258 0.5107302
# derived variabldes is the occupancy probability
names(mod.fit$results$derived)
## [1] "Occupancy"
mod.fit$results$derived$Occupancy
## estimate se lcl ucl
## 1 0.4647319 0.02338239 0.4193258 0.5107302
# alternatively
get.real(mod.fit, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 5 2 0.4647319 0.0233824 0.4193258
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.5107302 1 0 1 0 0
get.real(mod.fit, "Psi", pim=TRUE)
## 1
## 0.4647319
get.real(mod.fit, "p", se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 a0 t1 1 1 0.5519311 0.0183166 0.5158217
## p g1 a1 t2 2 1 0.5519311 0.0183166 0.5158217
## p g1 a2 t3 3 1 0.5519311 0.0183166 0.5158217
## p g1 a3 t4 4 1 0.5519311 0.0183166 0.5158217
## ucl fixed note group age time Age Time
## p g1 a0 t1 0.5875012 1 0 1 0 0
## p g1 a1 t2 0.5875012 1 1 2 1 1
## p g1 a2 t3 0.5875012 1 2 3 2 2
## p g1 a3 t4 0.5875012 1 3 4 3 3
# Fit a model allowing for heterogeneity
# Note that formula HAVE AN = SIGN
bullfrog.data <- process.data(data=input.history,
model="OccupHet")
summary(bullfrog.data)
## Length Class Mode
## data 2 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 500 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 4 -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("OccupHet", check=TRUE)
## [1] "pi" "p" "Psi"
mod.fit.het <- RMark::mark(bullfrog.data,
model="OccupHety", # notice change of modelnames
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~mixture),
pi =list(formula=~1)
)
)
##
## Output summary for OccupHet model
## Name : pi(~1)p(~mixture)Psi(~1)
##
## Npar : 4
## -2lnL: 1862.765
## AICc : 1870.846
##
## Beta
## estimate se lcl ucl
## pi:(Intercept) 1.2149812 0.6648739 -0.0881716 2.5181340
## p:(Intercept) -0.3203396 0.2522127 -0.8146764 0.1739973
## p:mixture2 2.2874400 0.6917540 0.9316020 3.6432779
## Psi:(Intercept) -0.0460831 0.1124190 -0.2664243 0.1742581
##
##
## Real Parameter pi
##
## 1
## mixture:1 0.7711791
##
##
## Real Parameter p
##
## 1 2 3 4
## mixture:1 0.4205930 0.4205930 0.4205930 0.4205930
## mixture:2 0.8772993 0.8772993 0.8772993 0.8772993
##
##
## Real Parameter Psi
##
## 1
## 0.4884813
summary(mod.fit.het)
## Output summary for OccupHet model
## Name : pi(~1)p(~mixture)Psi(~1)
##
## Npar : 4
## -2lnL: 1862.765
## AICc : 1870.846
##
## Beta
## estimate se lcl ucl
## pi:(Intercept) 1.2149812 0.6648739 -0.0881716 2.5181340
## p:(Intercept) -0.3203396 0.2522127 -0.8146764 0.1739973
## p:mixture2 2.2874400 0.6917540 0.9316020 3.6432779
## Psi:(Intercept) -0.0460831 0.1124190 -0.2664243 0.1742581
##
##
## Real Parameter pi
##
## 1
## mixture:1 0.7711791
##
##
## Real Parameter p
##
## 1 2 3 4
## mixture:1 0.4205930 0.4205930 0.4205930 0.4205930
## mixture:2 0.8772993 0.8772993 0.8772993 0.8772993
##
##
## Real Parameter Psi
##
## 1
## 0.4884813
# Look the objects returned in more details
names(mod.fit.het)
## [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.het$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.het$results$beta # on the logit scale
## estimate se lcl ucl
## pi:(Intercept) 1.2149812 0.6648739 -0.0881716 2.5181340
## p:(Intercept) -0.3203396 0.2522127 -0.8146764 0.1739973
## p:mixture2 2.2874400 0.6917540 0.9316020 3.6432779
## Psi:(Intercept) -0.0460831 0.1124190 -0.2664243 0.1742581
mod.fit.het$results$real# on the regular 0-1 scale for each site
## estimate se lcl ucl fixed note
## pi g1 a0 t1 m1 0.7711791 0.1173249 0.4779714 0.9254033
## p g1 a0 t1 m1 0.4205930 0.0614628 0.3068949 0.5433899
## p g1 a0 t1 m2 0.8772993 0.0956976 0.5559254 0.9760970
## Psi g1 a0 t1 0.4884813 0.0280898 0.4337851 0.5434546
# derived variabldes is the occupancy probability
names(mod.fit.het$results$derived)
## [1] "Occupancy"
mod.fit.het$results$derived$Occupancy
## estimate se lcl ucl
## 1 0.4884813 0.02808983 0.4337851 0.5434546
# alternatively
get.real(mod.fit.het, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 10 4 0.4884813 0.0280898 0.4337851
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.5434546 1 0 1 0 0
get.real(mod.fit.het, "Psi", pim=TRUE)
## [[1]]
## 1
## 0.4884813
get.real(mod.fit.het, "p", se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 a0 t1 m1 2 2 0.4205930 0.0614628 0.3068949
## p g1 a1 t2 m1 3 2 0.4205930 0.0614628 0.3068949
## p g1 a2 t3 m1 4 2 0.4205930 0.0614628 0.3068949
## p g1 a3 t4 m1 5 2 0.4205930 0.0614628 0.3068949
## p g1 a0 t1 m2 6 3 0.8772993 0.0956976 0.5559254
## p g1 a1 t2 m2 7 3 0.8772993 0.0956976 0.5559254
## p g1 a2 t3 m2 8 3 0.8772993 0.0956976 0.5559254
## p g1 a3 t4 m2 9 3 0.8772993 0.0956976 0.5559254
## ucl fixed note group age time mixture Age Time
## p g1 a0 t1 m1 0.5433899 1 0 1 1 0 0
## p g1 a1 t2 m1 0.5433899 1 1 2 1 1 1
## p g1 a2 t3 m1 0.5433899 1 2 3 1 2 2
## p g1 a3 t4 m1 0.5433899 1 3 4 1 3 3
## p g1 a0 t1 m2 0.9760970 1 0 1 2 0 0
## p g1 a1 t2 m2 0.9760970 1 1 2 2 1 1
## p g1 a2 t3 m2 0.9760970 1 2 3 2 2 2
## p g1 a3 t4 m2 0.9760970 1 3 4 2 3 3
# Fit a Royle Nichols Poisson model
bullfrog.data <- process.data(data=input.history,
model="OccupRNPoisson")
summary(bullfrog.data)
## Length Class Mode
## data 2 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 500 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 4 -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("OccupRNPoisson", check=TRUE)
## [1] "r" "Lambda"
mod.fit.rp <- RMark::mark(bullfrog.data,
model="OccupRNPoisson", # notice change of modelnames
model.parameters=list(
r =list(formula=~1),
Lambda =list(formula=~1))
)
##
## Output summary for OccupRNPoisson model
## Name : r(~1)Lambda(~1)
##
## Npar : 2
## -2lnL: 1870.52
## AICc : 1874.544
##
## Beta
## estimate se lcl ucl
## r:(Intercept) -0.1759362 0.0929685 -0.3581544 0.0062819
## Lambda:(Intercept) -0.4168280 0.0732902 -0.5604768 -0.2731791
##
##
## Real Parameter r
## 1
## 0.456129
##
##
## Real Parameter Lambda
## 1
## 0.6591343
summary(mod.fit.rp)
## Output summary for OccupRNPoisson model
## Name : r(~1)Lambda(~1)
##
## Npar : 2
## -2lnL: 1870.52
## AICc : 1874.544
##
## Beta
## estimate se lcl ucl
## r:(Intercept) -0.1759362 0.0929685 -0.3581544 0.0062819
## Lambda:(Intercept) -0.4168280 0.0732902 -0.5604768 -0.2731791
##
##
## Real Parameter r
## 1
## 0.456129
##
##
## Real Parameter Lambda
## 1
## 0.6591343
# Look the objects returned in more details
names(mod.fit.rp)
## [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.rp$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.rp$results$beta # on the logit scale
## estimate se lcl ucl
## r:(Intercept) -0.1759362 0.0929685 -0.3581544 0.0062819
## Lambda:(Intercept) -0.4168280 0.0732902 -0.5604768 -0.2731791
mod.fit.rp$results$real# on the regular 0-1 scale for each site
## estimate se lcl ucl fixed note
## r g1 a0 t1 0.4561290 0.0230632 0.4114064 0.5015705
## Lambda g1 a0 t1 0.6591343 0.0483081 0.5710466 0.7608101
# derived variabldes is the occupancy probability
names(mod.fit.rp$results$derived)
## [1] "Occupancy" "Expected p"
mod.fit.rp$results$derived$Occupancy
## estimate se lcl ucl
## 1 0.482701 0.02498972 0.4340422 0.5316902
collect.models()
## Warning in model.table(x, type, pf = 2, adjust = adjust): Model list contains models of differing types
## model npar AICc DeltaAICc weight Deviance
## 2 pi(~1)p(~mixture)Psi(~1) 4 1870.846 0.000000 8.639553e-01 8.516866
## 3 r(~1)Lambda(~1) 2 1874.544 3.697537 1.360133e-01 16.270996
## 1 p(~1)Psi(~1) 2 1891.292 20.445737 3.138744e-05 33.019267
# cleanup
cleanup(ask=FALSE)