# Single Species Single Season Occupancy
## Weta Example with model averaging with a list of models to fit
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)
xtabs(~BrowCat, data=input.history,exclude=NULL, na.action=na.pass)
## BrowCat
## b n
## 35 37
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"
# is survey covariates are present modify the ddl
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"
# define the list of models to fit
# Notice the commas between the column and the placement of the quotes
model.list.csv <- textConnection("
p, psi
~1, ~1
~1, ~BrowCat
~BrowCat, ~BrowCat
~time, ~BrowCat
~Observer2+Observer3, ~BrowCat
~Observer2+Observer3+time, ~BrowCat
~Observer2+Observer3+time+BrowCat, ~BrowCat
~Observer2+Observer3+time+BrowCat, ~1
~Observer2+Observer3+time, ~1")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p psi
## 1 ~1 ~1
## 2 ~1 ~BrowCat
## 3 ~BrowCat ~BrowCat
## 4 ~time ~BrowCat
## 5 ~Observer2+Observer3 ~BrowCat
## 6 ~Observer2+Observer3+time ~BrowCat
## 7 ~Observer2+Observer3+time+BrowCat ~BrowCat
## 8 ~Observer2+Observer3+time+BrowCat ~1
## 9 ~Observer2+Observer3+time ~1
# fit the models
model.fits <- plyr::alply(model.list, 1, 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)))
))
fit
},input.data=weta.data, input.ddl=weta.ddl)
##
##
## ***** Starting ~1 ~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.0851236 -0.1585387
## Psi:(Intercept) 0.4751241 0.3745433 -0.2589809 1.2092290
##
##
## 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
##
##
## ***** Starting ~1 ~BrowCat
##
## 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.0860904 -0.1584188
## Psi:(Intercept) 1.1492333 0.6557671 -0.1360704 2.4345369
## Psi:BrowCatn -1.2252776 0.7205601 -2.6375755 0.1870203
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3492689 0.3492689 0.3492689 0.3492689 0.3492689
## Group:BrowCatn 0.3492689 0.3492689 0.3492689 0.3492689 0.3492689
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7593708
## Group:BrowCatn 0.4809981
##
##
## ***** Starting ~BrowCat ~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.6107333 0.3022129 -1.2030705 -0.0183960
## p:BrowCatn -0.0296219 0.4859586 -0.9821007 0.9228569
## Psi:(Intercept) 1.1336538 0.6935116 -0.2256290 2.4929366
## Psi:BrowCatn -1.1980367 0.8416146 -2.8476013 0.4515279
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3518920 0.3518920 0.3518920 0.3518920 0.3518920
## Group:BrowCatn 0.3451663 0.3451663 0.3451663 0.3451663 0.3451663
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7565126
## Group:BrowCatn 0.4839098
##
##
## ***** Starting ~time ~BrowCat
##
## Output summary for Occupancy model
## Name : p(~time)Psi(~BrowCat)
##
## Npar : 7
## -2lnL: 245.4368
## AICc : 261.1868
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.6100128 0.4312909 -1.4553430 0.2353174
## p:time2 -0.1551359 0.5283110 -1.1906256 0.8803538
## p:time3 -0.9792834 0.5881659 -2.1320887 0.1735218
## p:time4 -0.1827226 0.5424148 -1.2458556 0.8804103
## p:time5 0.9811151 0.5456454 -0.0883498 2.0505801
## Psi:(Intercept) 1.2078313 0.6959862 -0.1563016 2.5719642
## Psi:BrowCatn -1.2351490 0.7508901 -2.7068936 0.2365956
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3520563 0.3175295 0.1694829 0.3115816 0.5917253
## Group:BrowCatn 0.3520563 0.3175295 0.1694829 0.3115816 0.5917253
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.769915
## Group:BrowCatn 0.493171
##
##
## ***** Starting ~Observer2+Observer3 ~BrowCat
##
## 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.2452636 0.3614719 -1.9537485 -0.5367787
## p:Observer2 0.7498088 0.4427020 -0.1178871 1.6175047
## p:Observer3 1.0294213 0.4324361 0.1818466 1.8769961
## Psi:(Intercept) 1.1178714 0.6269319 -0.1109151 2.3466579
## Psi:BrowCatn -1.1792049 0.7017344 -2.5546043 0.1961945
##
##
## 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.7535937
## Group:BrowCatn 0.4846714
##
##
## ***** Starting ~Observer2+Observer3+time ~BrowCat
##
## Output summary for Occupancy model
## Name : p(~Observer2 + Observer3 + time)Psi(~BrowCat)
##
## Npar : 9
## -2lnL: 239.5994
## AICc : 260.5027
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.2997188 0.5246875 -2.3281064 -0.2713312
## p:Observer2 0.7264614 0.4679006 -0.1906238 1.6435466
## p:Observer3 1.0700174 0.4601062 0.1682092 1.9718256
## p:time2 -0.1548921 0.5504438 -1.2337620 0.9239778
## p:time3 -0.9431914 0.5997280 -2.1186584 0.2322755
## p:time4 -0.0706551 0.5629214 -1.1739809 1.0326708
## p:time5 1.0354544 0.5639406 -0.0698692 2.1407780
## Psi:(Intercept) 1.1933137 0.6757954 -0.1312453 2.5178727
## Psi:BrowCatn -1.1693203 0.7412935 -2.6222556 0.2836150
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2928847 0.2777622 0.1488197 0.2697074 0.5249573
## Group:BrowCatn 0.2928847 0.2777622 0.1488197 0.2697074 0.5249573
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7673332
## Group:BrowCatn 0.5059981
##
##
## ***** Starting ~Observer2+Observer3+time+BrowCat ~BrowCat
##
## Output summary for Occupancy model
## Name : p(~Observer2 + Observer3 + time + BrowCat)Psi(~BrowCat)
##
## Npar : 10
## -2lnL: 239.5993
## AICc : 263.2059
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.3015574 0.5570564 -2.3933879 -0.2097269
## p:Observer2 0.7265861 0.4680886 -0.1908676 1.6440397
## p:Observer3 1.0703413 0.4612908 0.1662113 1.9744713
## p:time2 -0.1547380 0.5506905 -1.2340915 0.9246155
## p:time3 -0.9437799 0.6027611 -2.1251916 0.2376319
## p:time4 -0.0708170 0.5631902 -1.1746698 1.0330358
## p:time5 1.0348549 0.5672482 -0.0769517 2.1466615
## p:BrowCatn 0.0052092 0.5293593 -1.0323350 1.0427535
## Psi:(Intercept) 1.1957600 0.7216659 -0.2187051 2.6102251
## Psi:BrowCatn -1.1740836 0.8861535 -2.9109444 0.5627773
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2925270 0.2774494 0.1485283 0.2693320 0.5243721
## Group:BrowCatn 0.2936062 0.2784949 0.1491883 0.2703583 0.5256712
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.7677697
## Group:BrowCatn 0.5054189
##
##
## ***** Starting ~Observer2+Observer3+time+BrowCat ~1
##
## Output summary for Occupancy model
## Name : p(~Observer2 + Observer3 + time + BrowCat)Psi(~1)
##
## Npar : 9
## -2lnL: 241.4147
## AICc : 262.318
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.1981930 0.5332006 -2.2432662 -0.1531197
## p:Observer2 0.7341509 0.4618495 -0.1710740 1.6393759
## p:Observer3 1.0500603 0.4556465 0.1569932 1.9431274
## p:time2 -0.1651113 0.5413295 -1.2261171 0.8958945
## p:time3 -0.9041529 0.5926709 -2.0657878 0.2574820
## p:time4 -0.0471992 0.5539729 -1.1329862 1.0385878
## p:time5 1.0834173 0.5521181 0.0012657 2.1655688
## p:BrowCatn -0.4698723 0.4303054 -1.3132710 0.3735263
## Psi:(Intercept) 0.7375142 0.4801906 -0.2036593 1.6786878
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.3134259 0.2957352 0.1670429 0.2944888 0.5613736
## Group:BrowCatn 0.2220044 0.2079112 0.1113917 0.2069262 0.4444468
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.676452
## Group:BrowCatn 0.676452
##
##
## ***** Starting ~Observer2+Observer3+time ~1
##
## Output summary for Occupancy model
## Name : p(~Observer2 + Observer3 + time)Psi(~1)
##
## Npar : 8
## -2lnL: 242.5466
## AICc : 260.8323
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.3427885 0.5254503 -2.3726710 -0.3129060
## p:Observer2 0.7480961 0.4660298 -0.1653224 1.6615146
## p:Observer3 1.1043723 0.4588471 0.2050320 2.0037126
## p:time2 -0.1376622 0.5473991 -1.2105644 0.9352400
## p:time3 -0.9640982 0.5967632 -2.1337541 0.2055577
## p:time4 -0.0470487 0.5606901 -1.1460012 1.0519039
## p:time5 1.0484583 0.5618028 -0.0526753 2.1495918
## Psi:(Intercept) 0.5804493 0.4093417 -0.2218605 1.3827591
##
##
## Real Parameter p
## 1 2 3 4 5
## Group:BrowCatb 0.2867266 0.2757032 0.1427895 0.2681741 0.5202915
## Group:BrowCatn 0.2867266 0.2757032 0.1427895 0.2681741 0.5202915
##
##
## Real Parameter Psi
## 1
## Group:BrowCatb 0.6411708
## Group:BrowCatn 0.6411708
model.1 = model.fits[[1]]
model.2 = model.fits[[2]]
model.3 = model.fits[[3]]
model.4 = model.fits[[4]]
model.5 = model.fits[[5]]
model.6 = model.fits[[6]]
model.7 = model.fits[[7]]
model.8 = model.fits[[8]]
model.9 = model.fits[[9]]
#AICc model averaging
model.set <- collect.models(type="Occupancy")
model.set
## model npar AICc
## 6 p(~Observer2 + Observer3 + time)Psi(~BrowCat) 9 260.5027
## 9 p(~Observer2 + Observer3 + time)Psi(~1) 8 260.8323
## 4 p(~time)Psi(~BrowCat) 7 261.1868
## 8 p(~Observer2 + Observer3 + time + BrowCat)Psi(~1) 9 262.3180
## 5 p(~Observer2 + Observer3)Psi(~BrowCat) 5 262.9512
## 7 p(~Observer2 + Observer3 + time + BrowCat)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
## 6 0.0000000 0.26671085 239.5994
## 9 0.3296085 0.22618672 242.5466
## 4 0.6841742 0.18944107 -162.9434
## 8 1.8152900 0.10761070 241.4147
## 5 2.4485251 0.07840616 252.0421
## 7 2.7032416 0.06903022 239.5993
## 2 4.1145554 0.03408602 -150.1160
## 1 5.4583972 0.01740864 -146.5931
## 3 6.3549091 0.01111961 -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.7254433 0.1269608 0.4310130 0.9021168
## Psi gn a0 t1 12 0.5516005 0.1316221 0.3024059 0.7773251
## 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(model.fits[[8]], parameter="Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gb a0 t1 11 11 0.676452 0.1050968 0.4492604 0.8427307
## Psi gn a0 t1 12 11 0.676452 0.1050968 0.4492604 0.8427307
## 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.7254433 0.1269608 0.4310130 0.9021168
## 2 2 2 12 12 0.5516005 0.1316221 0.3024059 0.7773251
## 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.7254433 0.1269608 0.4310130 0.9021168
## 2 12 2 2 12 0.5516005 0.1316221 0.3024059 0.7773251
## 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")
p.ma <- RMark::model.average(model.set, parameter="p")
p.ma
## par.index estimate se fixed note group age time Age Time
## p gb a0 t1 1 0.3092201 0.09405035 b 0 1 0 0
## p gb a1 t2 2 0.2946092 0.08736754 b 1 2 1 1
## p gb a2 t3 3 0.1793613 0.08986047 b 2 3 2 2
## p gb a3 t4 4 0.2868748 0.08818353 b 3 4 3 3
## p gb a4 t5 5 0.5112050 0.13504099 b 4 5 4 4
## p gn a0 t1 6 0.2993819 0.09867159 n 0 1 0 0
## p gn a1 t2 7 0.2851558 0.09288497 n 1 2 1 1
## p gn a2 t3 8 0.1733434 0.09120876 n 2 3 2 2
## p gn a3 t4 9 0.2774482 0.09185280 n 3 4 3 3
## p gn a4 t5 10 0.4986373 0.13684062 n 4 5 4 4
## BrowCat
## p gb a0 t1 b
## p gb a1 t2 b
## p gb a2 t3 b
## p gb a3 t4 b
## p gb a4 t5 b
## p gn a0 t1 n
## p gn a1 t2 n
## p gn a2 t3 n
## p gn a3 t4 n
## p gn a4 t5 n
# cleanup
cleanup(ask=FALSE)