# Single Species, Multi Season Occupancy analyais
# Northern Spotted Owl (Strix occidentalis caurina) in California.
# s=55 sites visited up to K=5 times per season between 1997 and 2001 (Y=5).
# Detection probabilities relatively constant within years, but likely different among years.
# Using RMark and several models including different parameterizations of the multi-season model
# RMark package
library(plyr)
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 RMark additional functions
source(file.path("..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
# 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.csv(file.path("..","NSO.csv"), header=FALSE, skip=2, na.strings="-")
input.data$V1 <- NULL # drop the site number
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.data)
## [1] 55
ncol(input.data)
## [1] 40
range(input.data, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.data))
## [1] 752
# Extract the history records and create a capture history
input.history <- data.frame(freq=1,
ch=apply(input.data,1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 0111NANANANA010NANANANANA11NANANANANANA011NANANANANA0101NANANANA
## 2 1 00NANANANANANA000NANANANANA000NANANANANA000110NANA0011NANANANA
## 3 1 1000111NA0100NANANANA0001NANANANA000000NANA00000NANANA
## 4 1 111NANANANANA10NANANANANANA011NANANANANA00000100000011NANA
## 5 1 000000NANA0000000NA000000NANA000000NANA0111NANANANA
## 6 1 1NANANANANANANA000111NANA0111NANANANA01111NANANA011111NANA
# 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 0111....010.....11......011.....0101....
## 2 1 00......000.....000.....000110..0011....
## 3 1 1000111.0100....0001....000000..00000...
## 4 1 111.....10......011.....00000100000011..
## 5 1 000000..0000000.000000..000000..0111....
## 6 1 1.......000111..0111....01111...011111..
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 55
# Create the data structure
max.visit.per.year <- 8
n.season <- 5
nso.data <- process.data(data=input.history, model="RDOccupEG",
time.intervals=c( rep( c(rep(0,max.visit.per.year-1),1),n.season-1),
rep(0,max.visit.per.year-1)))
summary(nso.data)
## Length Class Mode
## data 2 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 55 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 5 -none- numeric
## time.intervals 39 -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("RDOccupEG", check=TRUE)
## [1] "Psi" "Epsilon" "Gamma" "p"
#there are other parameterizations available
setup.parameters("RDOccupPE", check=TRUE) # psi, epsilon, p
## [1] "Psi" "Epsilon" "p"
setup.parameters("RDOccupPG", check=TRUE) # psi, gamma, p
## [1] "Psi" "Gamma" "p"
# time-specific covariates must be added to the ddl's
# look at builtin variables for p
nso.ddl <- make.design.data(nso.data)
nso.ddl$p
## par.index model.index group time session Time
## 1 1 10 1 1 1 0
## 2 2 11 1 2 1 1
## 3 3 12 1 3 1 2
## 4 4 13 1 4 1 3
## 5 5 14 1 5 1 4
## 6 6 15 1 6 1 5
## 7 7 16 1 7 1 6
## 8 8 17 1 8 1 7
## 9 9 18 1 1 2 0
## 10 10 19 1 2 2 1
## 11 11 20 1 3 2 2
## 12 12 21 1 4 2 3
## 13 13 22 1 5 2 4
## 14 14 23 1 6 2 5
## 15 15 24 1 7 2 6
## 16 16 25 1 8 2 7
## 17 17 26 1 1 3 0
## 18 18 27 1 2 3 1
## 19 19 28 1 3 3 2
## 20 20 29 1 4 3 3
## 21 21 30 1 5 3 4
## 22 22 31 1 6 3 5
## 23 23 32 1 7 3 6
## 24 24 33 1 8 3 7
## 25 25 34 1 1 4 0
## 26 26 35 1 2 4 1
## 27 27 36 1 3 4 2
## 28 28 37 1 4 4 3
## 29 29 38 1 5 4 4
## 30 30 39 1 6 4 5
## 31 31 40 1 7 4 6
## 32 32 41 1 8 4 7
## 33 33 42 1 1 5 0
## 34 34 43 1 2 5 1
## 35 35 44 1 3 5 2
## 36 36 45 1 4 5 3
## 37 37 46 1 5 5 4
## 38 38 47 1 6 5 5
## 39 39 48 1 7 5 6
## 40 40 49 1 8 5 7
str(nso.ddl$p)
## 'data.frame': 40 obs. of 6 variables:
## $ par.index : int 1 2 3 4 5 6 7 8 9 10 ...
## $ model.index: num 10 11 12 13 14 15 16 17 18 19 ...
## $ group : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## $ time : Factor w/ 8 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 1 2 ...
## $ session : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ Time : num 0 1 2 3 4 5 6 7 0 1 ...
nso.ddl$Gamma
## par.index model.index group age time Age Time
## 1 1 6 1 0 1 0 0
## 2 2 7 1 1 2 1 1
## 3 3 8 1 2 3 2 2
## 4 4 9 1 3 4 3 3
# define the list of models to fit
# Notice the commas between the column and the placement of the quotes
# Define the models.
# model.type can be RDOccupPE, RDOccupPG, RDOccupEG
# The random occupancy model cannot be fit here. See the other code in this directory.
#
# Notice that the last 3 models are all identical.
model.list.csv <- textConnection("
p, Psi, Gamma, Epsilon, model.type
~1, ~1, ~1, ~1, RDOccupEG
~session, ~1, ~1, ~1, RDOccupEG
~session, ~1, ~time, ~time, RDOccupEG
~session, ~time, ~time, ~1, RDOccupPG
~session, ~time, ~1, ~time, RDOccupPE")
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 Gamma Epsilon model.type model.number
## 1 ~1 ~1 ~1 ~1 RDOccupEG 1
## 2 ~session ~1 ~1 ~1 RDOccupEG 2
## 3 ~session ~1 ~time ~time RDOccupEG 3
## 4 ~session ~time ~time ~1 RDOccupPG 4
## 5 ~session ~time ~1 ~time RDOccupPE 5
# 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.history){
cat("\n\n***** Starting ", unlist(x), "\n")
max.visit.per.year <- 8
n.season <- 5
# we need to process the data in the loop to allow for different data types
input.data <- process.data(data=input.history, model=x$model.type,
time.intervals=c( rep( c(rep(0,max.visit.per.year-1),1),n.season-1),
rep(0,max.visit.per.year-1)))
# set up the ddls as needed for time-varying covariates.
# you need to do this in the loop because different paraemterizations have different ddl structures
input.ddl <- make.design.data(input.data)
model.parameters=list(
Psi =list(formula=as.formula(eval(x$Psi))),
p =list(formula=as.formula(eval(x$p))),
Epsilon=list(formula=as.formula(eval(x$Epsilon))),
Gamma =list(formula=as.formula(eval(x$Gamma)))
)
if(x$model.type == "RDOccupPG"){ # psi, gamma, p formulation
model.parameters$Epsilon= NULL
}
if(x$model.type == "RDOccupPE"){ # psi, epsilon, p formulation
model.parameters$Gamma = NULL
}
fit <- RMark::mark(input.data, ddl=input.ddl,
model=x$model.type,
model.parameters=model.parameters
#,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.history=input.history)
##
##
## ***** Starting ~1 ~1 ~1 ~1 RDOccupEG 1
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~1)
##
## Npar : 4
## -2lnL: 1355.315
## AICc : 1363.464
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5372049 0.2889444 -0.0291261 1.1035359
## Epsilon:(Intercept) -1.7288379 0.2595729 -2.2376007 -1.2200751
## Gamma:(Intercept) -1.4883083 0.2843058 -2.0455478 -0.9310689
## p:(Intercept) -0.0210406 0.0743374 -0.1667418 0.1246607
##
##
## Real Parameter Psi
##
## 1
## 0.631162
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.1507363 0.1507363 0.1507363 0.1507363
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1841758 0.1841758 0.1841758 0.1841758
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:2
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:4
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
## Session:5
## 1 2 3 4 5 6 7
## 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401 0.4947401
## 8
## 0.4947401
##
##
## ***** Starting ~session ~1 ~1 ~1 RDOccupEG 2
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~session)
##
## Npar : 8
## -2lnL: 1337.523
## AICc : 1354.064
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5096825 0.2853276 -0.0495596 1.0689246
## Epsilon:(Intercept) -1.7963671 0.2687696 -2.3231556 -1.2695787
## Gamma:(Intercept) -1.5242751 0.2929070 -2.0983728 -0.9501774
## p:(Intercept) 0.3663506 0.1630641 0.0467449 0.6859562
## p:session2 -0.2764264 0.2294433 -0.7261352 0.1732824
## p:session3 -0.7406950 0.2449777 -1.2208513 -0.2605387
## p:session4 -0.8353622 0.2386106 -1.3030390 -0.3676854
## p:session5 -0.2197454 0.2267268 -0.6641300 0.2246391
##
##
## Real Parameter Psi
##
## 1
## 0.624732
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.1422939 0.1422939 0.1422939 0.1422939
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1788329 0.1788329 0.1788329 0.1788329
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5905769 0.5905769 0.5905769 0.5905769 0.5905769 0.5905769 0.5905769
## 8
## 0.5905769
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659
## 8
## 0.5224659
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4074917 0.4074917 0.4074917 0.4074917 0.4074917 0.4074917 0.4074917
## 8
## 0.4074917
##
## Session:4
## 1 2 3 4 5 6 7
## 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502
## 8
## 0.3848502
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858
## 8
## 0.5365858
##
##
## ***** Starting ~session ~1 ~time ~time RDOccupEG 3
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~time)Gamma(~time)p(~session)
##
## Npar : 14
## -2lnL: 1327.639
## AICc : 1357.255
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5301137 0.2856343 -0.0297295 1.0899569
## Epsilon:(Intercept) -2.3314311 0.6232547 -3.5530103 -1.1098519
## Epsilon:time2 0.4654276 0.8563187 -1.2129571 2.1438124
## Epsilon:time3 1.1739899 0.7855113 -0.3656123 2.7135920
## Epsilon:time4 0.3276693 0.8600529 -1.3580344 2.0133729
## Gamma:(Intercept) -2.1158505 0.7697003 -3.6244632 -0.6072378
## Gamma:time2 -0.4832565 1.2828779 -2.9976973 2.0311843
## Gamma:time3 1.6524519 0.8902763 -0.0924897 3.3973934
## Gamma:time4 0.0876487 1.1419229 -2.1505203 2.3258178
## p:(Intercept) 0.3612754 0.1634009 0.0410095 0.6815412
## p:session2 -0.2839547 0.2291275 -0.7330447 0.1651352
## p:session3 -0.7063170 0.2460574 -1.1885895 -0.2240446
## p:session4 -0.8327931 0.2453479 -1.3136751 -0.3519112
## p:session5 -0.2168013 0.2274616 -0.6626260 0.2290234
##
##
## Real Parameter Psi
##
## 1
## 0.6295096
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.0885531 0.1340048 0.2391325 0.1188085
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1075658 0.0691959 0.3861799 0.1162736
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5893491 0.5893491 0.5893491 0.5893491 0.5893491 0.5893491 0.5893491
## 8
## 0.5893491
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205
## 8
## 0.5193205
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4145853 0.4145853 0.4145853 0.4145853 0.4145853 0.4145853 0.4145853
## 8
## 0.4145853
##
## Session:4
## 1 2 3 4 5 6 7
## 0.3842571 0.3842571 0.3842571 0.3842571 0.3842571 0.3842571 0.3842571
## 8
## 0.3842571
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558
## 8
## 0.5360558
##
##
## ***** Starting ~session ~time ~time ~1 RDOccupPG 4
##
## Output summary for RDOccupPG model
## Name : Psi(~time)Gamma(~time)p(~session)
##
## Npar : 14
## -2lnL: 1327.639
## AICc : 1357.255
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5301144 0.2856351 -0.0297304 1.0899593
## Psi:time2 -0.0675737 0.1873771 -0.4348329 0.2996855
## Psi:time3 -0.2965581 0.2608942 -0.8079107 0.2147945
## Psi:time4 -0.1441849 0.3654834 -0.8605323 0.5721625
## Psi:time5 -0.2416053 0.3622116 -0.9515401 0.4683295
## Gamma:(Intercept) -2.1158488 0.7697151 -3.6244904 -0.6072073
## Gamma:time2 -0.4832636 1.2829016 -2.9977508 2.0312237
## Gamma:time3 1.6524497 0.8902884 -0.0925155 3.3974150
## Gamma:time4 0.0876476 1.1419584 -2.1505908 2.3258861
## p:(Intercept) 0.3612749 0.1634013 0.0410083 0.6815415
## p:session2 -0.2839545 0.2291280 -0.7330454 0.1651364
## p:session3 -0.7063163 0.2460577 -1.1885894 -0.2240432
## p:session4 -0.8327929 0.2453485 -1.3136759 -0.3519098
## p:session5 -0.2168008 0.2274620 -0.6626264 0.2290248
##
##
## Real Parameter Psi
##
## 1 2 3 4 5
## 0.6295098 0.6136167 0.5581251 0.5953024 0.5716311
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1075659 0.0691956 0.3861798 0.1162736
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7 8
## 0.589349 0.589349 0.589349 0.589349 0.589349 0.589349 0.589349 0.589349
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205 0.5193205
## 8
## 0.5193205
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4145854 0.4145854 0.4145854 0.4145854 0.4145854 0.4145854 0.4145854
## 8
## 0.4145854
##
## Session:4
## 1 2 3 4 5 6 7 8
## 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558
## 8
## 0.5360558
##
##
## ***** Starting ~session ~time ~1 ~time RDOccupPE 5
##
## Output summary for RDOccupPE model
## Name : Psi(~time)Epsilon(~time)p(~session)
##
## Npar : 14
## -2lnL: 1327.639
## AICc : 1357.255
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5301141 0.2856355 -0.0297316 1.0899598
## Psi:time2 -0.0675735 0.1873787 -0.4348357 0.2996887
## Psi:time3 -0.2965574 0.2608972 -0.8079159 0.2148011
## Psi:time4 -0.1441839 0.3654839 -0.8605323 0.5721646
## Psi:time5 -0.2416053 0.3622126 -0.9515420 0.4683315
## Epsilon:(Intercept) -2.3314297 0.6232653 -3.5530297 -1.1098297
## Epsilon:time2 0.4654277 0.8563252 -1.2129697 2.1438251
## Epsilon:time3 1.1739874 0.7855217 -0.3656352 2.7136101
## Epsilon:time4 0.3276682 0.8600714 -1.3580718 2.0134082
## p:(Intercept) 0.3612748 0.1634012 0.0410084 0.6815413
## p:session2 -0.2839539 0.2291279 -0.7330447 0.1651368
## p:session3 -0.7063166 0.2460577 -1.1885898 -0.2240434
## p:session4 -0.8327929 0.2453484 -1.3136758 -0.3519100
## p:session5 -0.2168007 0.2274620 -0.6626261 0.2290248
##
##
## Real Parameter Psi
##
## 1 2 3 4 5
## 0.6295097 0.6136167 0.5581252 0.5953026 0.571631
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.0885532 0.134005 0.2391323 0.1188086
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7 8
## 0.589349 0.589349 0.589349 0.589349 0.589349 0.589349 0.589349 0.589349
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5193206 0.5193206 0.5193206 0.5193206 0.5193206 0.5193206 0.5193206
## 8
## 0.5193206
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4145853 0.4145853 0.4145853 0.4145853 0.4145853 0.4145853 0.4145853
## 8
## 0.4145853
##
## Session:4
## 1 2 3 4 5 6 7 8
## 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257 0.384257
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558 0.5360558
## 8
## 0.5360558
# examine individula model results
model.number <-2
summary(model.fits[[model.number]])
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~session)
##
## Npar : 8
## -2lnL: 1337.523
## AICc : 1354.064
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5096825 0.2853276 -0.0495596 1.0689246
## Epsilon:(Intercept) -1.7963671 0.2687696 -2.3231556 -1.2695787
## Gamma:(Intercept) -1.5242751 0.2929070 -2.0983728 -0.9501774
## p:(Intercept) 0.3663506 0.1630641 0.0467449 0.6859562
## p:session2 -0.2764264 0.2294433 -0.7261352 0.1732824
## p:session3 -0.7406950 0.2449777 -1.2208513 -0.2605387
## p:session4 -0.8353622 0.2386106 -1.3030390 -0.3676854
## p:session5 -0.2197454 0.2267268 -0.6641300 0.2246391
##
##
## Real Parameter Psi
##
## 1
## 0.624732
##
##
## Real Parameter Epsilon
##
## 1 2 3 4
## 0.1422939 0.1422939 0.1422939 0.1422939
##
##
## Real Parameter Gamma
##
## 1 2 3 4
## 0.1788329 0.1788329 0.1788329 0.1788329
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5905769 0.5905769 0.5905769 0.5905769 0.5905769 0.5905769 0.5905769
## 8
## 0.5905769
##
## Session:2
## 1 2 3 4 5 6 7
## 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659 0.5224659
## 8
## 0.5224659
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4074917 0.4074917 0.4074917 0.4074917 0.4074917 0.4074917 0.4074917
## 8
## 0.4074917
##
## Session:4
## 1 2 3 4 5 6 7
## 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502 0.3848502
## 8
## 0.3848502
##
## Session:5
## 1 2 3 4 5 6 7
## 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858 0.5365858
## 8
## 0.5365858
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## Psi g1 a0 t1 0.6247320 0.0668928 0.4876126 0.7443924
## Epsilon g1 a0 t1 0.1422939 0.0328023 0.0892233 0.2193294
## Gamma g1 a0 t1 0.1788329 0.0430139 0.1092551 0.2788491
## p g1 s1 t1 0.5905769 0.0394282 0.5116841 0.6650668
## p g1 s2 t1 0.5224659 0.0404720 0.4432412 0.6005761
## p g1 s3 t1 0.4074917 0.0441805 0.3245375 0.4960771
## p g1 s4 t1 0.3848502 0.0412514 0.3077762 0.4681714
## p g1 s5 t1 0.5365858 0.0391612 0.4595707 0.6118942
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## Psi:(Intercept) 0.5096825 0.2853276 -0.0495596 1.0689246
## Epsilon:(Intercept) -1.7963671 0.2687696 -2.3231556 -1.2695787
## Gamma:(Intercept) -1.5242751 0.2929070 -2.0983728 -0.9501774
## p:(Intercept) 0.3663506 0.1630641 0.0467449 0.6859562
## p:session2 -0.2764264 0.2294433 -0.7261352 0.1732824
## p:session3 -0.7406950 0.2449777 -1.2208513 -0.2605387
## p:session4 -0.8353622 0.2386106 -1.3030390 -0.3676854
## p:session5 -0.2197454 0.2267268 -0.6641300 0.2246391
model.fits[[model.number]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.6247320 0.06689275 0.4936223 0.7558418
## 2 0.6029467 0.05143305 0.5021380 0.7037555
## 3 0.5881573 0.05153012 0.4871582 0.6891563
## 4 0.5781171 0.05696019 0.4664751 0.6897591
## 5 0.5713011 0.06259527 0.4486144 0.6939879
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 0.9651286 0.05067869 0.8657983 1.064459
## 2 0.9754714 0.03700607 0.9029395 1.048003
## 3 0.9829295 0.02652773 0.9309351 1.034924
## 4 0.9882100 0.01877819 0.9514048 1.025015
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 0.9121744 0.12994057 0.6574909 1.166858
## 2 0.9404418 0.08808530 0.7677946 1.113089
## 3 0.9595372 0.06022556 0.8414951 1.077579
## 4 0.9724982 0.04135490 0.8914426 1.053554
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.6247320 0.06689275 0.4936223 0.7558418
## 2 0.6029467 0.05143305 0.5021380 0.7037555
## 3 0.5881573 0.05153012 0.4871582 0.6891563
## 4 0.5781171 0.05696019 0.4664751 0.6897591
## 5 0.5713011 0.06259527 0.4486144 0.6939879
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 0.9651286 0.05067869 0.8657983 1.064459
## 2 0.9754714 0.03700607 0.9029395 1.048003
## 3 0.9829295 0.02652773 0.9309351 1.034924
## 4 0.9882100 0.01877819 0.9514048 1.025015
model.fits[[model.number]]$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 0.9121744 0.12994057 0.6574909 1.166858
## 2 0.9404418 0.08808530 0.7677946 1.113089
## 3 0.9595372 0.06022556 0.8414951 1.077579
## 4 0.9724982 0.04135490 0.8914426 1.053554
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 1 1 0.624732 0.0668928 0.4876126
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.7443924 1 0 1 0 0
get.real(model.fits[[model.number]], "p", se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 s1 t1 10 4 0.5905769 0.0394282 0.5116841
## p g1 s1 t2 11 4 0.5905769 0.0394282 0.5116841
## p g1 s1 t3 12 4 0.5905769 0.0394282 0.5116841
## p g1 s1 t4 13 4 0.5905769 0.0394282 0.5116841
## p g1 s1 t5 14 4 0.5905769 0.0394282 0.5116841
## p g1 s1 t6 15 4 0.5905769 0.0394282 0.5116841
## p g1 s1 t7 16 4 0.5905769 0.0394282 0.5116841
## p g1 s1 t8 17 4 0.5905769 0.0394282 0.5116841
## p g1 s2 t1 18 5 0.5224659 0.0404720 0.4432412
## p g1 s2 t2 19 5 0.5224659 0.0404720 0.4432412
## p g1 s2 t3 20 5 0.5224659 0.0404720 0.4432412
## p g1 s2 t4 21 5 0.5224659 0.0404720 0.4432412
## p g1 s2 t5 22 5 0.5224659 0.0404720 0.4432412
## p g1 s2 t6 23 5 0.5224659 0.0404720 0.4432412
## p g1 s2 t7 24 5 0.5224659 0.0404720 0.4432412
## p g1 s2 t8 25 5 0.5224659 0.0404720 0.4432412
## p g1 s3 t1 26 6 0.4074917 0.0441805 0.3245375
## p g1 s3 t2 27 6 0.4074917 0.0441805 0.3245375
## p g1 s3 t3 28 6 0.4074917 0.0441805 0.3245375
## p g1 s3 t4 29 6 0.4074917 0.0441805 0.3245375
## p g1 s3 t5 30 6 0.4074917 0.0441805 0.3245375
## p g1 s3 t6 31 6 0.4074917 0.0441805 0.3245375
## p g1 s3 t7 32 6 0.4074917 0.0441805 0.3245375
## p g1 s3 t8 33 6 0.4074917 0.0441805 0.3245375
## p g1 s4 t1 34 7 0.3848502 0.0412514 0.3077762
## p g1 s4 t2 35 7 0.3848502 0.0412514 0.3077762
## p g1 s4 t3 36 7 0.3848502 0.0412514 0.3077762
## p g1 s4 t4 37 7 0.3848502 0.0412514 0.3077762
## p g1 s4 t5 38 7 0.3848502 0.0412514 0.3077762
## p g1 s4 t6 39 7 0.3848502 0.0412514 0.3077762
## p g1 s4 t7 40 7 0.3848502 0.0412514 0.3077762
## p g1 s4 t8 41 7 0.3848502 0.0412514 0.3077762
## p g1 s5 t1 42 8 0.5365858 0.0391612 0.4595707
## p g1 s5 t2 43 8 0.5365858 0.0391612 0.4595707
## p g1 s5 t3 44 8 0.5365858 0.0391612 0.4595707
## p g1 s5 t4 45 8 0.5365858 0.0391612 0.4595707
## p g1 s5 t5 46 8 0.5365858 0.0391612 0.4595707
## p g1 s5 t6 47 8 0.5365858 0.0391612 0.4595707
## p g1 s5 t7 48 8 0.5365858 0.0391612 0.4595707
## p g1 s5 t8 49 8 0.5365858 0.0391612 0.4595707
## ucl fixed note group time session Time
## p g1 s1 t1 0.6650668 1 1 1 0
## p g1 s1 t2 0.6650668 1 2 1 1
## p g1 s1 t3 0.6650668 1 3 1 2
## p g1 s1 t4 0.6650668 1 4 1 3
## p g1 s1 t5 0.6650668 1 5 1 4
## p g1 s1 t6 0.6650668 1 6 1 5
## p g1 s1 t7 0.6650668 1 7 1 6
## p g1 s1 t8 0.6650668 1 8 1 7
## p g1 s2 t1 0.6005761 1 1 2 0
## p g1 s2 t2 0.6005761 1 2 2 1
## p g1 s2 t3 0.6005761 1 3 2 2
## p g1 s2 t4 0.6005761 1 4 2 3
## p g1 s2 t5 0.6005761 1 5 2 4
## p g1 s2 t6 0.6005761 1 6 2 5
## p g1 s2 t7 0.6005761 1 7 2 6
## p g1 s2 t8 0.6005761 1 8 2 7
## p g1 s3 t1 0.4960771 1 1 3 0
## p g1 s3 t2 0.4960771 1 2 3 1
## p g1 s3 t3 0.4960771 1 3 3 2
## p g1 s3 t4 0.4960771 1 4 3 3
## p g1 s3 t5 0.4960771 1 5 3 4
## p g1 s3 t6 0.4960771 1 6 3 5
## p g1 s3 t7 0.4960771 1 7 3 6
## p g1 s3 t8 0.4960771 1 8 3 7
## p g1 s4 t1 0.4681714 1 1 4 0
## p g1 s4 t2 0.4681714 1 2 4 1
## p g1 s4 t3 0.4681714 1 3 4 2
## p g1 s4 t4 0.4681714 1 4 4 3
## p g1 s4 t5 0.4681714 1 5 4 4
## p g1 s4 t6 0.4681714 1 6 4 5
## p g1 s4 t7 0.4681714 1 7 4 6
## p g1 s4 t8 0.4681714 1 8 4 7
## p g1 s5 t1 0.6118942 1 1 5 0
## p g1 s5 t2 0.6118942 1 2 5 1
## p g1 s5 t3 0.6118942 1 3 5 2
## p g1 s5 t4 0.6118942 1 4 5 3
## p g1 s5 t5 0.6118942 1 5 5 4
## p g1 s5 t6 0.6118942 1 6 5 5
## p g1 s5 t7 0.6118942 1 7 5 6
## p g1 s5 t8 0.6118942 1 8 5 7
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
## Warning in model.table(x, type, pf = 2, adjust = adjust): Model list contains models of differing types
model.set
## model npar AICc DeltaAICc
## 2 Psi(~1)Epsilon(~1)Gamma(~1)p(~session) 8 1354.064 0.000000
## 3 Psi(~1)Epsilon(~time)Gamma(~time)p(~session) 14 1357.254 3.190531
## 4 Psi(~time)Gamma(~time)p(~session) 14 1357.254 3.190531
## 5 Psi(~time)Epsilon(~time)p(~session) 14 1357.254 3.190531
## 1 Psi(~1)Epsilon(~1)Gamma(~1)p(~1) 4 1363.463 9.399495
## weight Deviance
## 2 0.618176302 896.7160
## 3 0.125399931 886.8325
## 4 0.125399931 886.8325
## 5 0.125399931 886.8325
## 1 0.005623905 914.5087
names(model.set)
## [1] "m...001" "m...002" "m...003" "m...004" "m...005"
## [6] "model.table"
model.set$model.table
## model npar AICc DeltaAICc
## 2 Psi(~1)Epsilon(~1)Gamma(~1)p(~session) 8 1354.064 0.000000
## 3 Psi(~1)Epsilon(~time)Gamma(~time)p(~session) 14 1357.254 3.190531
## 4 Psi(~time)Gamma(~time)p(~session) 14 1357.254 3.190531
## 5 Psi(~time)Epsilon(~time)p(~session) 14 1357.254 3.190531
## 1 Psi(~1)Epsilon(~1)Gamma(~1)p(~1) 4 1363.463 9.399495
## weight Deviance
## 2 0.618176302 896.7160
## 3 0.125399931 886.8325
## 4 0.125399931 886.8325
## 5 0.125399931 886.8325
## 1 0.005623905 914.5087
# The usual model averaging takes place
# RMark::model.average(model.set, "Psi") # oops a problem with different model structures
RMark::model.average(model.set, "p") # we can model the average p's
## par.index estimate se fixed note group time session
## p g1 s1 t1 10 0.5895760 0.04003120 1 1 1
## p g1 s1 t2 11 0.5895760 0.04003120 1 2 1
## p g1 s1 t3 12 0.5895760 0.04003120 1 3 1
## p g1 s1 t4 13 0.5895760 0.04003120 1 4 1
## p g1 s1 t5 14 0.5895760 0.04003120 1 5 1
## p g1 s1 t6 15 0.5895760 0.04003120 1 6 1
## p g1 s1 t7 16 0.5895760 0.04003120 1 7 1
## p g1 s1 t8 17 0.5895760 0.04003120 1 8 1
## p g1 s2 t1 18 0.5211267 0.04038531 1 1 2
## p g1 s2 t2 19 0.5211267 0.04038531 1 2 2
## p g1 s2 t3 20 0.5211267 0.04038531 1 3 2
## p g1 s2 t4 21 0.5211267 0.04038531 1 4 2
## p g1 s2 t5 22 0.5211267 0.04038531 1 5 2
## p g1 s2 t6 23 0.5211267 0.04038531 1 6 2
## p g1 s2 t7 24 0.5211267 0.04038531 1 7 2
## p g1 s2 t8 25 0.5211267 0.04038531 1 8 2
## p g1 s3 t1 26 0.4106510 0.04483895 1 1 3
## p g1 s3 t2 27 0.4106510 0.04483895 1 2 3
## p g1 s3 t3 28 0.4106510 0.04483895 1 3 3
## p g1 s3 t4 29 0.4106510 0.04483895 1 4 3
## p g1 s3 t5 30 0.4106510 0.04483895 1 5 3
## p g1 s3 t6 31 0.4106510 0.04483895 1 6 3
## p g1 s3 t7 32 0.4106510 0.04483895 1 7 3
## p g1 s3 t8 33 0.4106510 0.04483895 1 8 3
## p g1 s4 t1 34 0.3852451 0.04274576 1 1 4
## p g1 s4 t2 35 0.3852451 0.04274576 1 2 4
## p g1 s4 t3 36 0.3852451 0.04274576 1 3 4
## p g1 s4 t4 37 0.3852451 0.04274576 1 4 4
## p g1 s4 t5 38 0.3852451 0.04274576 1 5 4
## p g1 s4 t6 39 0.3852451 0.04274576 1 6 4
## p g1 s4 t7 40 0.3852451 0.04274576 1 7 4
## p g1 s4 t8 41 0.3852451 0.04274576 1 8 4
## p g1 s5 t1 42 0.5361511 0.03927290 1 1 5
## p g1 s5 t2 43 0.5361511 0.03927290 1 2 5
## p g1 s5 t3 44 0.5361511 0.03927290 1 3 5
## p g1 s5 t4 45 0.5361511 0.03927290 1 4 5
## p g1 s5 t5 46 0.5361511 0.03927290 1 5 5
## p g1 s5 t6 47 0.5361511 0.03927290 1 6 5
## p g1 s5 t7 48 0.5361511 0.03927290 1 7 5
## p g1 s5 t8 49 0.5361511 0.03927290 1 8 5
## Time
## p g1 s1 t1 0
## p g1 s1 t2 1
## p g1 s1 t3 2
## p g1 s1 t4 3
## p g1 s1 t5 4
## p g1 s1 t6 5
## p g1 s1 t7 6
## p g1 s1 t8 7
## p g1 s2 t1 0
## p g1 s2 t2 1
## p g1 s2 t3 2
## p g1 s2 t4 3
## p g1 s2 t5 4
## p g1 s2 t6 5
## p g1 s2 t7 6
## p g1 s2 t8 7
## p g1 s3 t1 0
## p g1 s3 t2 1
## p g1 s3 t3 2
## p g1 s3 t4 3
## p g1 s3 t5 4
## p g1 s3 t6 5
## p g1 s3 t7 6
## p g1 s3 t8 7
## p g1 s4 t1 0
## p g1 s4 t2 1
## p g1 s4 t3 2
## p g1 s4 t4 3
## p g1 s4 t5 4
## p g1 s4 t6 5
## p g1 s4 t7 6
## p g1 s4 t8 7
## p g1 s5 t1 0
## p g1 s5 t2 1
## p g1 s5 t3 2
## p g1 s5 t4 3
## p g1 s5 t5 4
## p g1 s5 t6 5
## p g1 s5 t7 6
## p g1 s5 t8 7
# Model average the derived parameters
# Because RMarks stores psi in different places, we standarize the dervied parameters.
# We need to do this here because collect.models re-extracts the output from MARK and wipes anything else
model.set[-length(model.set)] <- plyr::llply(model.set[-length(model.set)], function(x){RMark.add.derived(x)})
model.set[[1]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.6311620 0.06726525 0.4993221 0.7630018
## 2 0.6039540 0.05095844 0.5040754 0.7038325
## 3 0.5858583 0.05112760 0.4856482 0.6860684
## 4 0.5738230 0.05661538 0.4628569 0.6847892
## 5 0.5658186 0.06213794 0.4440282 0.6876089
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 0.9568922 0.05137574 0.8561958 1.057589
## 2 0.9700379 0.03736524 0.8968021 1.043274
## 3 0.9794571 0.02655675 0.9274059 1.031508
## 4 0.9860506 0.01858500 0.9496240 1.022477
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 0.8911547 0.13254333 0.6313697 1.150940
## 2 0.9276527 0.08800695 0.7551590 1.100146
## 3 0.9517972 0.05914705 0.8358690 1.067725
## 4 0.9678720 0.03998058 0.8895100 1.046234
##
## $all_psi
## estimate se lcl ucl
## 1 0.6311620 0.06726525 0.4993221 0.7630018
## 2 0.6039540 0.05095844 0.5040754 0.7038325
## 3 0.5858583 0.05112760 0.4856482 0.6860684
## 4 0.5738230 0.05661538 0.4628569 0.6847892
## 5 0.5658186 0.06213794 0.4440282 0.6876089
RMark.model.average.derived(model.set, "all_psi")
## estimate se lcl ucl
## 1 0.6265656 0.06683242 0.4955764 0.7575547
## 2 0.6069665 0.05796474 0.4933577 0.7205753
## 3 0.5768462 0.06098531 0.4573172 0.6963752
## 4 0.5845582 0.06340090 0.4602947 0.7088217
## 5 0.5713944 0.06463779 0.4447067 0.6980822
RMark.model.average.derived(model.set, "lambda Rate of Change")
## estimate se lcl ucl
## 1 0.9687031 0.05851943 0.8540071 1.083399
## 2 0.9506472 0.06513859 0.8229779 1.078316
## 3 1.0143912 0.10329114 0.8119443 1.216838
## 4 0.9776741 0.05898438 0.8620668 1.093281
RMark.model.average.derived(model.set, "log odds lambda")
## estimate se lcl ucl
## 1 0.9205150 0.1489909 0.6284982 1.212532
## 2 0.8857825 0.1414350 0.6085750 1.162990
## 3 1.0366371 0.2633759 0.5204298 1.552844
## 4 0.9478970 0.1342364 0.6847986 1.210995
cleanup(ask=FALSE)