# Single Species Single Season Occupancy
# Weta Example
# RMark -
# categorical variable for site level covariate can be modelled directly as a covariate
# or as a grouping variable with the same results.
# showing how to use a covariate that varies by sites x visit (e.g. observers)
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library(readxl)
library(RMark)
## This is RMark 2.2.7
## 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(file.path("..","weta.xls"),
sheet="detection_histories",
na="-",
col_names=FALSE) # notice no column names in row 1 of data file.
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.data)
## [1] 72
ncol(input.data)
## [1] 5
range(input.data, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.data))
## [1] 98
head(input.data)
## # A tibble: 6 x 5
## ...1 ...2 ...3 ...4 ...5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 NA
## 2 0 0 0 0 NA
## 3 0 0 0 1 NA
## 4 0 0 0 0 NA
## 5 0 0 0 0 NA
## 6 0 0 0 0 NA
# 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 0000NA
## 2 1 0000NA
## 3 1 0001NA
## 4 1 0000NA
## 5 1 0000NA
## 6 1 0000NA
# 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 0000.
## 2 1 0000.
## 3 1 0001.
## 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] 72
# Get the site level covariates
site_covar <- readxl::read_excel(file.path("..","weta.xls"),
sheet="site_covar",
na="-",
col_names=TRUE) # notice col_names in row 1 of table.
head(site_covar)
## # A tibble: 6 x 3
## Browsed Unbrowsed BrowseCat
## <dbl> <dbl> <chr>
## 1 1 0 b
## 2 1 0 b
## 3 1 0 b
## 4 0 1 n
## 5 1 0 b
## 6 0 1 n
# Create a site level covariate that is a categorical variable rather
# than indicator variables. Be sure that it is declared as a FACTOR
input.history$BrowCat <- factor(site_covar$BrowseCat)
head(input.history)
## freq ch BrowCat
## 1 1 0000. b
## 2 1 0000. b
## 3 1 0001. b
## 4 1 0000. n
## 5 1 0000. b
## 6 1 0000. n
# Get the individual covariates.
obs1 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs1",
na="-",
col_names=FALSE)
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
obs2 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs2",
na="-",
col_names=FALSE)
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
obs3 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs3",
na="-",
col_names=FALSE)
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
Obs <-obs1*1 + obs2*2 + obs3*3
Obs<- as.data.frame(lapply(Obs, as.character), stringsAsFactors=FALSE)
Obs[ is.na(Obs)] <- "." # set all missing values to .
head(Obs)
## ...1 ...2 ...3 ...4 ...5
## 1 1 3 2 3 .
## 2 1 3 2 3 .
## 3 1 3 2 3 .
## 4 1 3 2 3 .
## 5 1 3 2 3 .
## 6 1 3 2 3 .
# make each column a factor variable. Be sure to use
Obs[] <- lapply(Obs, as.factor)
str(Obs)
## 'data.frame': 72 obs. of 5 variables:
## $ ...1: Factor w/ 4 levels ".","1","2","3": 2 2 2 2 2 2 2 2 2 2 ...
## $ ...2: Factor w/ 4 levels ".","1","2","3": 4 4 4 4 4 4 4 4 4 4 ...
## $ ...3: Factor w/ 4 levels ".","1","2","3": 3 3 3 3 3 3 3 3 3 3 ...
## $ ...4: Factor w/ 4 levels ".","1","2","3": 4 4 4 4 4 4 4 4 4 4 ...
## $ ...5: Factor w/ 4 levels ".","1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
colnames(Obs) = paste("Observer", 1:5, sep="") # assign the observer at each time point
head(Obs)
## Observer1 Observer2 Observer3 Observer4 Observer5
## 1 1 3 2 3 .
## 2 1 3 2 3 .
## 3 1 3 2 3 .
## 4 1 3 2 3 .
## 5 1 3 2 3 .
## 6 1 3 2 3 .
Obs.df <- make.time.factor(Obs, var.name="Observer", times=1:5, intercept=1)
head(Obs.df) # Column names are indicator variables for Obsever x at time y
## Observer21 Observer31 Observer22 Observer32 Observer23 Observer33 Observer24
## 1 0 0 0 1 1 0 0
## 2 0 0 0 1 1 0 0
## 3 0 0 0 1 1 0 0
## 4 0 0 0 1 1 0 0
## 5 0 0 0 1 1 0 0
## 6 0 0 0 1 1 0 0
## Observer34 Observer25 Observer35
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 1 0 0
## 6 1 0 0
input.history <- cbind(input.history, Obs.df)
head(input.history)
## freq ch BrowCat Observer21 Observer31 Observer22 Observer32 Observer23
## 1 1 0000. b 0 0 0 1 1
## 2 1 0000. b 0 0 0 1 1
## 3 1 0001. b 0 0 0 1 1
## 4 1 0000. n 0 0 0 1 1
## 5 1 0000. b 0 0 0 1 1
## 6 1 0000. n 0 0 0 1 1
## Observer33 Observer24 Observer34 Observer25 Observer35
## 1 0 0 1 0 0
## 2 0 0 1 0 0
## 3 0 0 1 0 0
## 4 0 0 1 0 0
## 5 0 0 1 0 0
## 6 0 0 1 0 0
# Illustration of using categorical covariate in the modelling process by creating groups
weta.data <- process.data(data=input.history, group="BrowCat",
model="Occupancy")
summary(weta.data)
## Length Class Mode
## data 14 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 2 data.frame list
## 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 2 -none- numeric
## group.covariates 1 data.frame list
## 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"
# modify the ddl if needed (e.g. for visit level covariates) - none in this case
weta.ddl <- make.design.data(weta.data)
weta.ddl
## $p
## par.index model.index group age time Age Time BrowCat
## 1 1 1 b 0 1 0 0 b
## 2 2 2 b 1 2 1 1 b
## 3 3 3 b 2 3 2 2 b
## 4 4 4 b 3 4 3 3 b
## 5 5 5 b 4 5 4 4 b
## 6 6 6 n 0 1 0 0 n
## 7 7 7 n 1 2 1 1 n
## 8 8 8 n 2 3 2 2 n
## 9 9 9 n 3 4 3 3 n
## 10 10 10 n 4 5 4 4 n
##
## $Psi
## par.index model.index group age time Age Time BrowCat
## 1 1 11 b 0 1 0 0 b
## 2 2 12 n 0 1 0 0 n
##
## $pimtypes
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
##
##
## $pimtypes$Psi
## $pimtypes$Psi$pim.type
## [1] "all"
# some simple models
mod.fit1 <- RMark::mark(weta.data, ddl=weta.ddl,
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: 261.7871
## AICc : 265.9611
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6218311 0.2363737 -1.085123 -0.1585386
## Psi:(Intercept) 0.4751239 0.3745433 -0.258981 1.2092288
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3493651 0.3493651 0.3493651 0.3493651 0.3493651
## Group:BrowCatn 0.3493651 0.3493651 0.3493651 0.3493651 0.3493651
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.6165958
## Group:BrowCatn 0.6165958
summary(mod.fit1)
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 261.7871
## AICc : 265.9611
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6218311 0.2363737 -1.085123 -0.1585386
## Psi:(Intercept) 0.4751239 0.3745433 -0.258981 1.2092288
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3493651 0.3493651 0.3493651 0.3493651 0.3493651
## Group:BrowCatn 0.3493651 0.3493651 0.3493651 0.3493651 0.3493651
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.6165958
## Group:BrowCatn 0.6165958
get.real(mod.fit1, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gb a0 t1 11 2 0.6165958 0.0885441 0.4356142 0.7701625
## Psi gn a0 t1 12 2 0.6165958 0.0885441 0.4356142 0.7701625
## fixed note group age time Age Time BrowCat
## Psi gb a0 t1 b 0 1 0 0 b
## Psi gn a0 t1 n 0 1 0 0 n
get.real(mod.fit1, "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p gb a0 t1 1 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gb a1 t2 2 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gb a2 t3 3 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gb a3 t4 4 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gb a4 t5 5 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gn a0 t1 6 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gn a1 t2 7 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gn a2 t3 8 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gn a3 t4 9 1 0.3493651 0.0537299 0.2525377 0.4604481
## p gn a4 t5 10 1 0.3493651 0.0537299 0.2525377 0.4604481
## fixed note group age time Age Time BrowCat
## p gb a0 t1 b 0 1 0 0 b
## p gb a1 t2 b 1 2 1 1 b
## p gb a2 t3 b 2 3 2 2 b
## p gb a3 t4 b 3 4 3 3 b
## p gb a4 t5 b 4 5 4 4 b
## p gn a0 t1 n 0 1 0 0 n
## p gn a1 t2 n 1 2 1 1 n
## p gn a2 t3 n 2 3 2 2 n
## p gn a3 t4 n 3 4 3 3 n
## p gn a4 t5 n 4 5 4 4 n
# some simple models
mod.fit2 <- RMark::mark(weta.data, ddl = weta.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~BrowCat),
p =list(formula=~1)
)
)
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~BrowCat)
##
## Npar : 3
## -2lnL: 258.2643
## AICc : 264.6172
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6222546 0.2366509 -1.0860905 -0.1584188
## Psi:(Intercept) 1.1492333 0.6557672 -0.1360704 2.4345370
## Psi:BrowCatn -1.2252777 0.7205601 -2.6375756 0.1870202
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3492688 0.3492688 0.3492688 0.3492688 0.3492688
## Group:BrowCatn 0.3492688 0.3492688 0.3492688 0.3492688 0.3492688
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7593708
## Group:BrowCatn 0.4809980
summary(mod.fit2)
## Output summary for Occupancy model
## Name : p(~1)Psi(~BrowCat)
##
## Npar : 3
## -2lnL: 258.2643
## AICc : 264.6172
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6222546 0.2366509 -1.0860905 -0.1584188
## Psi:(Intercept) 1.1492333 0.6557672 -0.1360704 2.4345370
## Psi:BrowCatn -1.2252777 0.7205601 -2.6375756 0.1870202
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3492688 0.3492688 0.3492688 0.3492688 0.3492688
## Group:BrowCatn 0.3492688 0.3492688 0.3492688 0.3492688 0.3492688
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7593708
## Group:BrowCatn 0.4809980
get.real(mod.fit2, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gb a0 t1 11 2 0.7593708 0.1198262 0.4660348 0.9194233
## Psi gn a0 t1 12 3 0.4809980 0.1078863 0.2843273 0.6837389
## fixed note group age time Age Time BrowCat
## Psi gb a0 t1 b 0 1 0 0 b
## Psi gn a0 t1 n 0 1 0 0 n
get.real(mod.fit2, "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p gb a0 t1 1 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gb a1 t2 2 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gb a2 t3 3 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gb a3 t4 4 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gb a4 t5 5 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gn a0 t1 6 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gn a1 t2 7 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gn a2 t3 8 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gn a3 t4 9 1 0.3492688 0.0537861 0.2523552 0.4604779
## p gn a4 t5 10 1 0.3492688 0.0537861 0.2523552 0.4604779
## fixed note group age time Age Time BrowCat
## p gb a0 t1 b 0 1 0 0 b
## p gb a1 t2 b 1 2 1 1 b
## p gb a2 t3 b 2 3 2 2 b
## p gb a3 t4 b 3 4 3 3 b
## p gb a4 t5 b 4 5 4 4 b
## p gn a0 t1 n 0 1 0 0 n
## p gn a1 t2 n 1 2 1 1 n
## p gn a2 t3 n 2 3 2 2 n
## p gn a3 t4 n 3 4 3 3 n
## p gn a4 t5 n 4 5 4 4 n
# some simple models
mod.fit3 <- RMark::mark(weta.data, ddl=weta.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~BrowCat),
p =list(formula=~BrowCat)
)
)
##
## Output summary for Occupancy model
## Name : p(~BrowCat)Psi(~BrowCat)
##
## Npar : 4
## -2lnL: 258.2606
## AICc : 266.8576
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6107331 0.3022129 -1.2030704 -0.0183958
## p:BrowCatn -0.0296222 0.4859588 -0.9821014 0.9228570
## Psi:(Intercept) 1.1336534 0.6935113 -0.2256287 2.4929356
## Psi:BrowCatn -1.1980362 0.8416144 -2.8476005 0.4515280
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3518920 0.3518920 0.3518920 0.3518920 0.3518920
## Group:BrowCatn 0.3451662 0.3451662 0.3451662 0.3451662 0.3451662
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7565125
## Group:BrowCatn 0.4839099
summary(mod.fit3)
## Output summary for Occupancy model
## Name : p(~BrowCat)Psi(~BrowCat)
##
## Npar : 4
## -2lnL: 258.2606
## AICc : 266.8576
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6107331 0.3022129 -1.2030704 -0.0183958
## p:BrowCatn -0.0296222 0.4859588 -0.9821014 0.9228570
## Psi:(Intercept) 1.1336534 0.6935113 -0.2256287 2.4929356
## Psi:BrowCatn -1.1980362 0.8416144 -2.8476005 0.4515280
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3518920 0.3518920 0.3518920 0.3518920 0.3518920
## Group:BrowCatn 0.3451662 0.3451662 0.3451662 0.3451662 0.3451662
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7565125
## Group:BrowCatn 0.4839099
get.real(mod.fit3, "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gb a0 t1 11 3 0.7565125 0.1277457 0.4438309 0.9236451
## Psi gn a0 t1 12 4 0.4839099 0.1190814 0.2691481 0.7047838
## fixed note group age time Age Time BrowCat
## Psi gb a0 t1 b 0 1 0 0 b
## Psi gn a0 t1 n 0 1 0 0 n
get.real(mod.fit3, "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p gb a0 t1 1 1 0.3518920 0.0689239 0.2309295 0.4954012
## p gb a1 t2 2 1 0.3518920 0.0689239 0.2309295 0.4954012
## p gb a2 t3 3 1 0.3518920 0.0689239 0.2309295 0.4954012
## p gb a3 t4 4 1 0.3518920 0.0689239 0.2309295 0.4954012
## p gb a4 t5 5 1 0.3518920 0.0689239 0.2309295 0.4954012
## p gn a0 t1 6 2 0.3451662 0.0860158 0.2000077 0.5263594
## p gn a1 t2 7 2 0.3451662 0.0860158 0.2000077 0.5263594
## p gn a2 t3 8 2 0.3451662 0.0860158 0.2000077 0.5263594
## p gn a3 t4 9 2 0.3451662 0.0860158 0.2000077 0.5263594
## p gn a4 t5 10 2 0.3451662 0.0860158 0.2000077 0.5263594
## fixed note group age time Age Time BrowCat
## p gb a0 t1 b 0 1 0 0 b
## p gb a1 t2 b 1 2 1 1 b
## p gb a2 t3 b 2 3 2 2 b
## p gb a3 t4 b 3 4 3 3 b
## p gb a4 t5 b 4 5 4 4 b
## p gn a0 t1 n 0 1 0 0 n
## p gn a1 t2 n 1 2 1 1 n
## p gn a2 t3 n 2 3 2 2 n
## p gn a3 t4 n 3 4 3 3 n
## p gn a4 t5 n 4 5 4 4 n
# some simple models
# Fit a model
# Note that formula have an equal sign
mod.fit4 <- RMark::mark(weta.data, ddl=weta.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~BrowCat),
p =list(formula=~Observer2+Observer3) # need to specify rest of obsevers after the intercept
)
)
##
## Output summary for Occupancy model
## Name : p(~Observer2 + Observer3)Psi(~BrowCat)
##
## Npar : 5
## -2lnL: 252.0421
## AICc : 262.9512
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.2452637 0.3614719 -1.9537487 -0.5367787
## p:Observer2 0.7498089 0.4427021 -0.1178871 1.6175050
## p:Observer3 1.0294214 0.4324361 0.1818466 1.8769963
## Psi:(Intercept) 1.1178712 0.6269319 -0.1109154 2.3466577
## Psi:BrowCatn -1.1792046 0.7017345 -2.5546041 0.1961950
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3027271 0.3205963 0.3205963 0.2943 0.2921418
## Group:BrowCatn 0.3027271 0.3205963 0.3205963 0.2943 0.2921418
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7535936
## Group:BrowCatn 0.4846715
summary(mod.fit4)
## Output summary for Occupancy model
## Name : p(~Observer2 + Observer3)Psi(~BrowCat)
##
## Npar : 5
## -2lnL: 252.0421
## AICc : 262.9512
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.2452637 0.3614719 -1.9537487 -0.5367787
## p:Observer2 0.7498089 0.4427021 -0.1178871 1.6175050
## p:Observer3 1.0294214 0.4324361 0.1818466 1.8769963
## Psi:(Intercept) 1.1178712 0.6269319 -0.1109154 2.3466577
## Psi:BrowCatn -1.1792046 0.7017345 -2.5546041 0.1961950
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3027271 0.3205963 0.3205963 0.2943 0.2921418
## Group:BrowCatn 0.3027271 0.3205963 0.3205963 0.2943 0.2921418
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7535936
## Group:BrowCatn 0.4846715
# Look the objects returned in more details
names(mod.fit4)
## [1] "data" "model" "title" "model.name"
## [5] "links" "mixtures" "call" "parameters"
## [9] "time.intervals" "number.of.groups" "group.labels" "nocc"
## [13] "begin.time" "covariates" "fixed" "design.matrix"
## [17] "pims" "design.data" "strata.labels" "mlogit.list"
## [21] "profile.int" "simplify" "model.parameters" "results"
## [25] "output"
names(mod.fit4$results)
## [1] "lnl" "deviance" "deviance.df" "npar"
## [5] "n" "AICc" "beta" "real"
## [9] "beta.vcv" "derived" "derived.vcv" "covariate.values"
## [13] "singular" "real.vcv"
# look at estimates on beta and original scale
mod.fit4$results$beta # on the logit scale
## estimate se lcl ucl
## p:(Intercept) -1.2452637 0.3614719 -1.9537487 -0.5367787
## p:Observer2 0.7498089 0.4427021 -0.1178871 1.6175050
## p:Observer3 1.0294214 0.4324361 0.1818466 1.8769963
## Psi:(Intercept) 1.1178712 0.6269319 -0.1109154 2.3466577
## Psi:BrowCatn -1.1792046 0.7017345 -2.5546041 0.1961950
mod.fit4$results$real# on the regular 0-1 scale for each site
## estimate se lcl ucl fixed note
## p gb a0 t1 0.3027271 0.0546476 0.2072189 0.4189906
## p gb a1 t2 0.3205963 0.0533763 0.2259472 0.4327286
## p gb a2 t3 0.3205963 0.0533763 0.2259472 0.4327286
## p gb a3 t4 0.2943000 0.0545457 0.1995112 0.4110010
## p gb a4 t5 0.2921418 0.0547807 0.1971480 0.4095578
## Psi gb a0 t1 0.7535936 0.1164152 0.4722995 0.9126682
## Psi gn a0 t1 0.4846715 0.1084301 0.2865447 0.6877365
# derived variabldes is the occupancy probability
names(mod.fit4$results$derived)
## [1] "Occupancy"
mod.fit4$results$derived$Occupancy
## estimate lcl ucl
## 1 0.7535936 0.4722995 0.9126682
## 2 0.4846715 0.2865447 0.6877365
# get the two psi values and their covariance
get.real(mod.fit4, "Psi", se=TRUE, vcv=TRUE)
## $estimates
## all.diff.index par.index estimate se lcl ucl
## Psi gb a0 t1 11 6 0.7535936 0.1164152 0.4722995 0.9126682
## Psi gn a0 t1 12 7 0.4846715 0.1084301 0.2865447 0.6877365
## group age time Age Time BrowCat
## Psi gb a0 t1 b 0 1 0 0 b
## Psi gn a0 t1 n 0 1 0 0 n
##
## $vcv.real
## 6 7
## 6 0.013552489 0.002065716
## 7 0.002065716 0.011757088
get.real(mod.fit4, "Psi", pim=TRUE)
## 1
## Group:BrowCatb 0.7535936
## Group:BrowCatn 0.4846715
# Estimate of p are much ore complicated when you have site x visit covariates
# These don't make much sense as they are computed at the average
# value of the observer for time 1:5
get.real(mod.fit4, "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p gb a0 t1 1 1 0.3027271 0.0546476 0.2072189 0.4189906
## p gb a1 t2 2 2 0.3205963 0.0533763 0.2259472 0.4327286
## p gb a2 t3 3 3 0.3205963 0.0533763 0.2259472 0.4327286
## p gb a3 t4 4 4 0.2943000 0.0545457 0.1995112 0.4110010
## p gb a4 t5 5 5 0.2921418 0.0547807 0.1971480 0.4095578
## p gn a0 t1 6 1 0.3027271 0.0546476 0.2072189 0.4189906
## p gn a1 t2 7 2 0.3205963 0.0533763 0.2259472 0.4327286
## p gn a2 t3 8 3 0.3205963 0.0533763 0.2259472 0.4327286
## p gn a3 t4 9 4 0.2943000 0.0545457 0.1995112 0.4110010
## p gn a4 t5 10 5 0.2921418 0.0547807 0.1971480 0.4095578
## fixed note group age time Age Time BrowCat
## p gb a0 t1 b 0 1 0 0 b
## p gb a1 t2 b 1 2 1 1 b
## p gb a2 t3 b 2 3 2 2 b
## p gb a3 t4 b 3 4 3 3 b
## p gb a4 t5 b 4 5 4 4 b
## p gn a0 t1 n 0 1 0 0 n
## p gn a1 t2 n 1 2 1 1 n
## p gn a2 t3 n 2 3 2 2 n
## p gn a3 t4 n 3 4 3 3 n
## p gn a4 t5 n 4 5 4 4 n
# get the estimated directly from the beta matrix
mod.fit4$results$beta # on the logit scale
## estimate se lcl ucl
## p:(Intercept) -1.2452637 0.3614719 -1.9537487 -0.5367787
## p:Observer2 0.7498089 0.4427021 -0.1178871 1.6175050
## p:Observer3 1.0294214 0.4324361 0.1818466 1.8769963
## Psi:(Intercept) 1.1178712 0.6269319 -0.1109154 2.3466577
## Psi:BrowCatn -1.1792046 0.7017345 -2.5546041 0.1961950
obs.logit <- matrix(c(1,0,0,0,0,
1,1,0,0,0,
1,0,1,0,0), ncol=5, byrow=TRUE) %*% mod.fit4$results$beta$estimate
expit <- function(x){1/(1+exp(-x))}
expit(obs.logit)
## [,1]
## [1,] 0.2235211
## [2,] 0.3786094
## [3,] 0.4462479
# Or we can get covariate predictions
weta.ddl$p # see the index numbers
## par.index model.index group age time Age Time BrowCat
## 1 1 1 b 0 1 0 0 b
## 2 2 2 b 1 2 1 1 b
## 3 3 3 b 2 3 2 2 b
## 4 4 4 b 3 4 3 3 b
## 5 5 5 b 4 5 4 4 b
## 6 6 6 n 0 1 0 0 n
## 7 7 7 n 1 2 1 1 n
## 8 8 8 n 2 3 2 2 n
## 9 9 9 n 3 4 3 3 n
## 10 10 10 n 4 5 4 4 n
# use the model.index index nubers
# create a special matrix, or use the data matrix for each site x visit
obs.df <-data.frame(Observer21=c(0,1,0),
Observer22=c(0,1,0),
Observer23=c(0,1,0),
Observer24=c(0,1,0),
Observer25=c(0,1,0),
Observer31=c(0,0,1),
Observer32=c(0,0,1),
Observer33=c(0,0,1),
Observer34=c(0,0,1),
Observer35=c(0,0,1))
obs.p <- covariate.predictions(mod.fit4, indices=1:10, data=obs.df)
obs.p$estimates
## vcv.index model.index par.index Observer21 Observer22 Observer23 Observer24
## 1 1 1 1 0 0 0 0
## 2 2 2 2 0 0 0 0
## 3 3 3 3 0 0 0 0
## 4 4 4 4 0 0 0 0
## 5 5 5 5 0 0 0 0
## 6 1 1 6 0 0 0 0
## 7 2 2 7 0 0 0 0
## 8 3 3 8 0 0 0 0
## 9 4 4 9 0 0 0 0
## 10 5 5 10 0 0 0 0
## 11 6 1 1 1 1 1 1
## 12 7 2 2 1 1 1 1
## 13 8 3 3 1 1 1 1
## 14 9 4 4 1 1 1 1
## 15 10 5 5 1 1 1 1
## 16 6 1 6 1 1 1 1
## 17 7 2 7 1 1 1 1
## 18 8 3 8 1 1 1 1
## 19 9 4 9 1 1 1 1
## 20 10 5 10 1 1 1 1
## 21 11 1 1 0 0 0 0
## 22 12 2 2 0 0 0 0
## 23 13 3 3 0 0 0 0
## 24 14 4 4 0 0 0 0
## 25 15 5 5 0 0 0 0
## 26 11 1 6 0 0 0 0
## 27 12 2 7 0 0 0 0
## 28 13 3 8 0 0 0 0
## 29 14 4 9 0 0 0 0
## 30 15 5 10 0 0 0 0
## Observer25 Observer31 Observer32 Observer33 Observer34 Observer35 estimate
## 1 0 0 0 0 0 0 0.2235211
## 2 0 0 0 0 0 0 0.2235211
## 3 0 0 0 0 0 0 0.2235211
## 4 0 0 0 0 0 0 0.2235211
## 5 0 0 0 0 0 0 0.2235211
## 6 0 0 0 0 0 0 0.2235211
## 7 0 0 0 0 0 0 0.2235211
## 8 0 0 0 0 0 0 0.2235211
## 9 0 0 0 0 0 0 0.2235211
## 10 0 0 0 0 0 0 0.2235211
## 11 1 0 0 0 0 0 0.3786094
## 12 1 0 0 0 0 0 0.3786094
## 13 1 0 0 0 0 0 0.3786094
## 14 1 0 0 0 0 0 0.3786094
## 15 1 0 0 0 0 0 0.3786094
## 16 1 0 0 0 0 0 0.3786094
## 17 1 0 0 0 0 0 0.3786094
## 18 1 0 0 0 0 0 0.3786094
## 19 1 0 0 0 0 0 0.3786094
## 20 1 0 0 0 0 0 0.3786094
## 21 0 1 1 1 1 1 0.4462479
## 22 0 1 1 1 1 1 0.4462479
## 23 0 1 1 1 1 1 0.4462479
## 24 0 1 1 1 1 1 0.4462479
## 25 0 1 1 1 1 1 0.4462479
## 26 0 1 1 1 1 1 0.4462479
## 27 0 1 1 1 1 1 0.4462479
## 28 0 1 1 1 1 1 0.4462479
## 29 0 1 1 1 1 1 0.4462479
## 30 0 1 1 1 1 1 0.4462479
## se lcl ucl fixed
## 1 0.06273686 0.1241466 0.3689342
## 2 0.06273686 0.1241466 0.3689342
## 3 0.06273686 0.1241466 0.3689342
## 4 0.06273686 0.1241466 0.3689342
## 5 0.06273686 0.1241466 0.3689342
## 6 0.06273686 0.1241466 0.3689342
## 7 0.06273686 0.1241466 0.3689342
## 8 0.06273686 0.1241466 0.3689342
## 9 0.06273686 0.1241466 0.3689342
## 10 0.06273686 0.1241466 0.3689342
## 11 0.07998876 0.2383343 0.5426272
## 12 0.07998876 0.2383343 0.5426272
## 13 0.07998876 0.2383343 0.5426272
## 14 0.07998876 0.2383343 0.5426272
## 15 0.07998876 0.2383343 0.5426272
## 16 0.07998876 0.2383343 0.5426272
## 17 0.07998876 0.2383343 0.5426272
## 18 0.07998876 0.2383343 0.5426272
## 19 0.07998876 0.2383343 0.5426272
## 20 0.07998876 0.2383343 0.5426272
## 21 0.08081049 0.2980099 0.6047049
## 22 0.08081049 0.2980099 0.6047049
## 23 0.08081049 0.2980099 0.6047049
## 24 0.08081049 0.2980099 0.6047049
## 25 0.08081049 0.2980099 0.6047049
## 26 0.08081049 0.2980099 0.6047049
## 27 0.08081049 0.2980099 0.6047049
## 28 0.08081049 0.2980099 0.6047049
## 29 0.08081049 0.2980099 0.6047049
## 30 0.08081049 0.2980099 0.6047049
mod.fit5 <- RMark::mark(weta.data, ddl = weta.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~BrowCat),
p =list(formula=~time+Observer2+Observer3) # need to specify rest of obsevers after the intercept
)
)
##
## Output summary for Occupancy model
## Name : p(~time + Observer2 + Observer3)Psi(~BrowCat)
##
## Npar : 9
## -2lnL: 239.5994
## AICc : 260.5027
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.2997191 0.5246876 -2.3281067 -0.2713315
## p:time2 -0.1548920 0.5504436 -1.2337616 0.9239776
## p:time3 -0.9431914 0.5997280 -2.1186584 0.2322756
## p:time4 -0.0706548 0.5629214 -1.1739809 1.0326712
## p:time5 1.0354546 0.5639405 -0.0698689 2.1407781
## p:Observer2 0.7264616 0.4679005 -0.1906235 1.6435467
## p:Observer3 1.0700176 0.4601062 0.1682095 1.9718258
## Psi:(Intercept) 1.1933139 0.6757955 -0.1312453 2.5178732
## Psi:BrowCatn -1.1693205 0.7412936 -2.6222560 0.2836150
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2928846 0.2777622 0.1488197 0.2697074 0.5249573
## Group:BrowCatn 0.2928846 0.2777622 0.1488197 0.2697074 0.5249573
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7673332
## Group:BrowCatn 0.5059981
summary(mod.fit5)
## Output summary for Occupancy model
## Name : p(~time + Observer2 + Observer3)Psi(~BrowCat)
##
## Npar : 9
## -2lnL: 239.5994
## AICc : 260.5027
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.2997191 0.5246876 -2.3281067 -0.2713315
## p:time2 -0.1548920 0.5504436 -1.2337616 0.9239776
## p:time3 -0.9431914 0.5997280 -2.1186584 0.2322756
## p:time4 -0.0706548 0.5629214 -1.1739809 1.0326712
## p:time5 1.0354546 0.5639405 -0.0698689 2.1407781
## p:Observer2 0.7264616 0.4679005 -0.1906235 1.6435467
## p:Observer3 1.0700176 0.4601062 0.1682095 1.9718258
## Psi:(Intercept) 1.1933139 0.6757955 -0.1312453 2.5178732
## Psi:BrowCatn -1.1693205 0.7412936 -2.6222560 0.2836150
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2928846 0.2777622 0.1488197 0.2697074 0.5249573
## Group:BrowCatn 0.2928846 0.2777622 0.1488197 0.2697074 0.5249573
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7673332
## Group:BrowCatn 0.5059981
mod.fit6 <- RMark::mark(weta.data, ddl = weta.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~BrowCat),
p =list(formula=~time+BrowCat+Observer2+Observer3) # need to specify rest of obsevers after the intercept
)
)
##
## Output summary for Occupancy model
## Name : p(~time + BrowCat + Observer2 + Observer3)Psi(~BrowCat)
##
## Npar : 10
## -2lnL: 239.5993
## AICc : 263.2059
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.3015575 0.5570566 -2.3933884 -0.2097266
## p:time2 -0.1547381 0.5506907 -1.2340920 0.9246157
## p:time3 -0.9437797 0.6027612 -2.1251916 0.2376321
## p:time4 -0.0708168 0.5631904 -1.1746700 1.0330364
## p:time5 1.0348548 0.5672484 -0.0769521 2.1466617
## p:BrowCatn 0.0052093 0.5293594 -1.0323352 1.0427538
## p:Observer2 0.7265862 0.4680886 -0.1908674 1.6440398
## p:Observer3 1.0703415 0.4612908 0.1662114 1.9744715
## Psi:(Intercept) 1.1957600 0.7216658 -0.2187050 2.6102249
## Psi:BrowCatn -1.1740835 0.8861534 -2.9109442 0.5627772
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2925269 0.2774494 0.1485283 0.2693320 0.5243721
## Group:BrowCatn 0.2936062 0.2784949 0.1491883 0.2703584 0.5256712
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7677696
## Group:BrowCatn 0.5054189
summary(mod.fit6)
## Output summary for Occupancy model
## Name : p(~time + BrowCat + Observer2 + Observer3)Psi(~BrowCat)
##
## Npar : 10
## -2lnL: 239.5993
## AICc : 263.2059
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.3015575 0.5570566 -2.3933884 -0.2097266
## p:time2 -0.1547381 0.5506907 -1.2340920 0.9246157
## p:time3 -0.9437797 0.6027612 -2.1251916 0.2376321
## p:time4 -0.0708168 0.5631904 -1.1746700 1.0330364
## p:time5 1.0348548 0.5672484 -0.0769521 2.1466617
## p:BrowCatn 0.0052093 0.5293594 -1.0323352 1.0427538
## p:Observer2 0.7265862 0.4680886 -0.1908674 1.6440398
## p:Observer3 1.0703415 0.4612908 0.1662114 1.9744715
## Psi:(Intercept) 1.1957600 0.7216658 -0.2187050 2.6102249
## Psi:BrowCatn -1.1740835 0.8861534 -2.9109442 0.5627772
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2925269 0.2774494 0.1485283 0.2693320 0.5243721
## Group:BrowCatn 0.2936062 0.2784949 0.1491883 0.2703584 0.5256712
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7677696
## Group:BrowCatn 0.5054189
mod.fit7 <- RMark::mark(weta.data, ddl = weta.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~time+BrowCat+Observer2+Observer3) # need to specify rest of obsevers after the intercept
)
)
##
## Output summary for Occupancy model
## Name : p(~time + BrowCat + Observer2 + Observer3)Psi(~1)
##
## Npar : 9
## -2lnL: 241.4147
## AICc : 262.318
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.1981930 0.5332003 -2.2432655 -0.1531205
## p:time2 -0.1651115 0.5413288 -1.2261160 0.8958929
## p:time3 -0.9041531 0.5926705 -2.0657874 0.2574811
## p:time4 -0.0471995 0.5539723 -1.1329851 1.0385862
## p:time5 1.0834171 0.5521178 0.0012663 2.1655679
## p:BrowCatn -0.4698721 0.4303053 -1.3132705 0.3735264
## p:Observer2 0.7341510 0.4618493 -0.1710737 1.6393758
## p:Observer3 1.0500605 0.4556464 0.1569936 1.9431275
## Psi:(Intercept) 0.7375142 0.4801906 -0.2036595 1.6786878
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3134259 0.2957352 0.1670429 0.2944888 0.5613735
## Group:BrowCatn 0.2220045 0.2079112 0.1113917 0.2069262 0.4444469
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.676452
## Group:BrowCatn 0.676452
summary(mod.fit7)
## Output summary for Occupancy model
## Name : p(~time + BrowCat + Observer2 + Observer3)Psi(~1)
##
## Npar : 9
## -2lnL: 241.4147
## AICc : 262.318
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.1981930 0.5332003 -2.2432655 -0.1531205
## p:time2 -0.1651115 0.5413288 -1.2261160 0.8958929
## p:time3 -0.9041531 0.5926705 -2.0657874 0.2574811
## p:time4 -0.0471995 0.5539723 -1.1329851 1.0385862
## p:time5 1.0834171 0.5521178 0.0012663 2.1655679
## p:BrowCatn -0.4698721 0.4303053 -1.3132705 0.3735264
## p:Observer2 0.7341510 0.4618493 -0.1710737 1.6393758
## p:Observer3 1.0500605 0.4556464 0.1569936 1.9431275
## Psi:(Intercept) 0.7375142 0.4801906 -0.2036595 1.6786878
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3134259 0.2957352 0.1670429 0.2944888 0.5613735
## Group:BrowCatn 0.2220045 0.2079112 0.1113917 0.2069262 0.4444469
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.676452
## Group:BrowCatn 0.676452
mod.fit8 <- RMark::mark(weta.data, ddl = weta.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~time+Observer2+Observer3) # need to specify rest of obsevers after the intercept
)
)
##
## Output summary for Occupancy model
## Name : p(~time + Observer2 + Observer3)Psi(~1)
##
## Npar : 8
## -2lnL: 242.5466
## AICc : 260.8323
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.3427886 0.5254504 -2.3726715 -0.3129058
## p:time2 -0.1376620 0.5473992 -1.2105645 0.9352406
## p:time3 -0.9640979 0.5967633 -2.1337540 0.2055583
## p:time4 -0.0470485 0.5606903 -1.1460015 1.0519046
## p:time5 1.0484585 0.5618030 -0.0526754 2.1495924
## p:Observer2 0.7480961 0.4660299 -0.1653226 1.6615148
## p:Observer3 1.1043723 0.4588471 0.2050319 2.0037127
## Psi:(Intercept) 0.5804494 0.4093418 -0.2218605 1.3827594
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2867266 0.2757033 0.1427895 0.2681741 0.5202915
## Group:BrowCatn 0.2867266 0.2757033 0.1427895 0.2681741 0.5202915
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.6411708
## Group:BrowCatn 0.6411708
summary(mod.fit8)
## Output summary for Occupancy model
## Name : p(~time + Observer2 + Observer3)Psi(~1)
##
## Npar : 8
## -2lnL: 242.5466
## AICc : 260.8323
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.3427886 0.5254504 -2.3726715 -0.3129058
## p:time2 -0.1376620 0.5473992 -1.2105645 0.9352406
## p:time3 -0.9640979 0.5967633 -2.1337540 0.2055583
## p:time4 -0.0470485 0.5606903 -1.1460015 1.0519046
## p:time5 1.0484585 0.5618030 -0.0526754 2.1495924
## p:Observer2 0.7480961 0.4660299 -0.1653226 1.6615148
## p:Observer3 1.1043723 0.4588471 0.2050319 2.0037127
## Psi:(Intercept) 0.5804494 0.4093418 -0.2218605 1.3827594
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2867266 0.2757033 0.1427895 0.2681741 0.5202915
## Group:BrowCatn 0.2867266 0.2757033 0.1427895 0.2681741 0.5202915
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.6411708
## Group:BrowCatn 0.6411708
######################################################333
######################################################333
######################################################333
model.set <- collect.models(type="Occupancy")
model.set
## model npar AICc
## 5 p(~time + Observer2 + Observer3)Psi(~BrowCat) 9 260.5027
## 8 p(~time + Observer2 + Observer3)Psi(~1) 8 260.8323
## 7 p(~time + BrowCat + Observer2 + Observer3)Psi(~1) 9 262.3180
## 4 p(~Observer2 + Observer3)Psi(~BrowCat) 5 262.9512
## 6 p(~time + BrowCat + Observer2 + Observer3)Psi(~BrowCat) 10 263.2059
## 2 p(~1)Psi(~BrowCat) 3 264.6172
## 1 p(~1)Psi(~1) 2 265.9611
## 3 p(~BrowCat)Psi(~BrowCat) 4 266.8576
## DeltaAICc weight Deviance
## 5 0.0000000 0.32904560 239.5994
## 8 0.3296085 0.27905031 242.5466
## 7 1.8152900 0.13276111 241.4147
## 4 2.4485251 0.09673099 252.0421
## 6 2.7032416 0.08516372 239.5993
## 2 4.1145554 0.04205249 -150.1160
## 1 5.4583972 0.02147732 -146.5931
## 3 6.3549091 0.01371845 -150.1197
psi.ma <- RMark::model.average(model.set, parameter="Psi", vcv=TRUE)
names(psi.ma)
## [1] "estimates" "vcv.real"
head(psi.ma$estimates)
## par.index estimate se lcl ucl fixed note
## Psi gb a0 t1 11 0.7150495 0.1255525 0.4285700 0.8935722
## Psi gn a0 t1 12 0.5652565 0.1321877 0.3117893 0.7886512
## group age time Age Time BrowCat
## Psi gb a0 t1 b 0 1 0 0 b
## Psi gn a0 t1 n 0 1 0 0 n
ggplot(data=psi.ma$estimates, aes(x=BrowCat, y=estimate))+
geom_point()+
geom_errorbar(aes(ymin=lcl,
ymax=ucl), width=0.2)+
ylim(0,1)+
ylab("Occupancy (with 95% ci)")+xlab("Browse Category")
# or use covariate predictions to estimate the occupancy at different
# browse categories
# use the get.real to figure out the index numbers
# by looking at the alt.diff.index numbers
get.real(mod.fit8, parameter="Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gb a0 t1 11 6 0.6411708 0.0941776 0.4447613 0.7994338
## Psi gn a0 t1 12 6 0.6411708 0.0941776 0.4447613 0.7994338
## fixed note group age time Age Time BrowCat
## Psi gb a0 t1 b 0 1 0 0 b
## Psi gn a0 t1 n 0 1 0 0 n
psi.ma2 <- covariate.predictions(model.set, indices=11:12)
head(psi.ma2$estimates)
## vcv.index model.index par.index index estimate se lcl ucl
## 1 2 2 11 11 0.7150495 0.1255525 0.4285700 0.8935722
## 2 2 2 12 12 0.5652565 0.1321877 0.3117893 0.7886512
## fixed
## 1
## 2
psi.ma2$estimates <- merge(psi.ma2$estimates,
weta.ddl$Psi[,c("model.index","BrowCat")],
by.x="par.index", by.y="model.index")
psi.ma2$estimates
## par.index vcv.index model.index index estimate se lcl ucl
## 1 11 2 2 11 0.7150495 0.1255525 0.4285700 0.8935722
## 2 12 2 2 12 0.5652565 0.1321877 0.3117893 0.7886512
## fixed BrowCat
## 1 b
## 2 n
ggplot(data=psi.ma2$estimates, aes(x=BrowCat, y=estimate))+
geom_point()+
geom_errorbar(aes(ymin=lcl,
ymax=ucl), width=0.2)+
ylim(0,1)+
ylab("Occupancy (with 95% ci)")+xlab("Browse Category")
# Much more difficult to compute model averaged values of p
# cleanup
cleanup(ask=FALSE)