# Single Species, Multi Season Occupancy analyais
# nightingales over four seasons with covariates on elevation and stream type.
# We will use only the last 2 years (unsure why)
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
# RMark package
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 RPResence 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("..","nightingales.csv"), header=TRUE, as.is=TRUE, strip.white=TRUE)
head(input.data)
## Site V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19
## 1 1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1
## 2 2 0 1 1 1 1 1 1 0 0 0 0 0 0 1 0 1 0 0 0
## 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 4 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
## 5 5 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0
## 6 6 1 1 1 1 1 1 1 1 0 0 1 0 0 1 0 1 0 0 1
## V20 V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37
## 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 0
## 2 0 1 1 1 0 0 0 0 1 1 1 1 0 0 0 1 1 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0
## 6 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1
## V38 V39 V40 V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52 V53 V54 V55
## 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 1 0 1 0 1 1 1 0 0 0 0 1 0 1 0 0 1
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 1 1 0 0 0 0 1 0 1 0 0 0 1 1 1 0 1 0
## V56 V57 V58 V59 V60 V61 V62 V63 V64 V65 V66 V67 V68 V69 V70 V71 V72 V73
## 1 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 0
## 2 0 0 0 1 0 1 0 0 0 0 0 1 1 1 1 1 1 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 1 0 0 1 1 0 1 0 0 0 0 1 1 1 1 1 1 0
## V74 V75 V76 V77 V78 V79 V80
## 1 0 1 0 0 0 0 0
## 2 0 0 0 0 0 0 0
## 3 0 1 0 0 1 0 0
## 4 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0
## 6 0 1 1 1 1 0 0
input.history <- input.data[, -1]
head(input.history)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20
## 1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 1
## 2 0 1 1 1 1 1 1 0 0 0 0 0 0 1 0 1 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0
## 6 1 1 1 1 1 1 1 1 0 0 1 0 0 1 0 1 0 0 1 1
## V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37 V38
## 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 0 1
## 2 1 1 1 0 0 0 0 1 1 1 1 0 0 0 1 1 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1
## 6 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1
## V39 V40 V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52 V53 V54 V55 V56
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 1 0 1 0 1 1 1 0 0 0 0 1 0 1 0 0 1 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 1 0 0 0 0 1 0 1 0 0 0 1 1 1 0 1 0 1
## V57 V58 V59 V60 V61 V62 V63 V64 V65 V66 V67 V68 V69 V70 V71 V72 V73 V74
## 1 0 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 0 0
## 2 0 0 1 0 1 0 0 0 0 0 1 1 1 1 1 1 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 1 1 0 1 0 0 0 0 1 1 1 1 1 1 0 0
## V75 V76 V77 V78 V79 V80
## 1 1 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 1 0 0 1 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 1 1 1 1 0 0
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 55
ncol(input.history)
## [1] 80
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 0
# Site covariate. Here there are none
site.covar <- NULL
row.names(site.covar) <- NULL
head(site.covar)
## NULL
#Format the capture history to be used by RMark
input.history <- data.frame(freq=1,
ch=apply(input.history,1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history)
## freq
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
## ch
## 1 01010000000100000111111111101111111001000000000000000000000100000011001000100000
## 2 01111110000001010000111000011110001100101011100001010010001010000011111100000000
## 3 00000000000000000000000000000000000000000000000000000000000000000000000000100100
## 4 00111110000000000000000000000000000000000000000000000000000000000000000000000000
## 5 00001000000001010000000000000000001101100000000000000000000000000000000000000000
## 6 11111111001001010011110000000000001111100001010001110101001101000011111100111100
#Add site covariates to input history
#input.history = cbind(input.history,site.covar)
#head(input.history)
#Create the RMark data structure
# 10 years with 8 visits per year
max.visit.per.year <- 8
n.season <- 10
nightingale.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(nightingale.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 10 -none- numeric
## time.intervals 79 -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
# all time-specific covariates must be added to the ddl in the loop below
# Get the parameter names
# What are the parameter names for Single Season Single Species models
setup.parameters("RDOccupEG", check=TRUE)
## [1] "Psi" "Epsilon" "Gamma" "p"
# 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~time
#The random occupancy model cannot be fit here. See the other code in this directory.
# The reserved keywords can be found by looking at the ddl structures earlier
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")
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
# 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")
# we need to process the data in the loop to allow for different data types
# notice that max.visits.per.year and n.season are defined outside the function (bad programming practise)
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: 3810.899
## AICc : 3818.972
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 1.4269742 0.3482761 0.7443531 2.1095953
## Epsilon:(Intercept) -0.8329071 0.1350572 -1.0976192 -0.5681950
## Gamma:(Intercept) -1.0798478 0.1529212 -1.3795734 -0.7801223
## p:(Intercept) -0.0524793 0.0425940 -0.1359635 0.0310049
##
##
## Real Parameter Psi
##
## 1
## 0.8064294
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.3030307 0.3030307 0.3030307 0.3030307 0.3030307 0.3030307 0.3030307
## 8 9
## 0.3030307 0.3030307
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7
## 0.2535348 0.2535348 0.2535348 0.2535348 0.2535348 0.2535348 0.2535348
## 8 9
## 0.2535348 0.2535348
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:2
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:3
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:4
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:5
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:6
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:7
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:8
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:9
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
## Session:10
## 1 2 3 4 5 6 7
## 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832 0.4868832
## 8
## 0.4868832
##
##
## ***** Starting ~session ~1 ~1 ~1 RDOccupEG 2
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~session)
##
## Npar : 13
## -2lnL: 3566.397
## AICc : 3593.077
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 1.4171968 0.3457931 0.7394424 2.0949512
## Epsilon:(Intercept) -0.9217753 0.1417699 -1.1996443 -0.6439063
## Gamma:(Intercept) -1.1269845 0.1601751 -1.4409276 -0.8130413
## p:(Intercept) 0.0330090 0.1090620 -0.1807525 0.2467706
## p:session2 -1.1916082 0.1728642 -1.5304220 -0.8527944
## p:session3 1.0092192 0.1744336 0.6673293 1.3511090
## p:session4 1.1243693 0.2043931 0.7237588 1.5249798
## p:session5 -0.2825578 0.1725432 -0.6207426 0.0556269
## p:session6 -0.3061072 0.1701695 -0.6396395 0.0274251
## p:session7 -0.0019924 0.1893119 -0.3730439 0.3690590
## p:session8 -1.1188374 0.2294534 -1.5685660 -0.6691089
## p:session9 -0.0942170 0.1919391 -0.4704176 0.2819835
## p:session10 -0.1154180 0.1891406 -0.4861336 0.2552977
##
##
## Real Parameter Psi
##
## 1
## 0.8048986
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.2845963 0.2845963 0.2845963 0.2845963 0.2845963 0.2845963 0.2845963
## 8 9
## 0.2845963 0.2845963
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7 8
## 0.244718 0.244718 0.244718 0.244718 0.244718 0.244718 0.244718 0.244718
## 9
## 0.244718
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5082515 0.5082515 0.5082515 0.5082515 0.5082515 0.5082515 0.5082515
## 8
## 0.5082515
##
## Session:2
## 1 2 3 4 5 6 7
## 0.2389219 0.2389219 0.2389219 0.2389219 0.2389219 0.2389219 0.2389219
## 8
## 0.2389219
##
## Session:3
## 1 2 3 4 5 6 7
## 0.7392797 0.7392797 0.7392797 0.7392797 0.7392797 0.7392797 0.7392797
## 8
## 0.7392797
##
## Session:4
## 1 2 3 4 5 6 7 8
## 0.760856 0.760856 0.760856 0.760856 0.760856 0.760856 0.760856 0.760856
##
## Session:5
## 1 2 3 4 5 6 7
## 0.4379346 0.4379346 0.4379346 0.4379346 0.4379346 0.4379346 0.4379346
## 8
## 0.4379346
##
## Session:6
## 1 2 3 4 5 6 7
## 0.4321466 0.4321466 0.4321466 0.4321466 0.4321466 0.4321466 0.4321466
## 8
## 0.4321466
##
## Session:7
## 1 2 3 4 5 6 7
## 0.5077535 0.5077535 0.5077535 0.5077535 0.5077535 0.5077535 0.5077535
## 8
## 0.5077535
##
## Session:8
## 1 2 3 4 5 6 7
## 0.2524046 0.2524046 0.2524046 0.2524046 0.2524046 0.2524046 0.2524046
## 8
## 0.2524046
##
## Session:9
## 1 2 3 4 5 6 7
## 0.4847028 0.4847028 0.4847028 0.4847028 0.4847028 0.4847028 0.4847028
## 8
## 0.4847028
##
## Session:10
## 1 2 3 4 5 6 7
## 0.4794094 0.4794094 0.4794094 0.4794094 0.4794094 0.4794094 0.4794094
## 8
## 0.4794094
##
##
## ***** Starting ~session ~1 ~time ~time RDOccupEG 3
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~time)Gamma(~time)p(~session)
##
## Npar : 29
## -2lnL: 3533.994
## AICc : 3595.34
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 1.4031461000 0.3417650 0.7332868 2.0730055
## Epsilon:(Intercept) -2.3102572000 0.8125807 -3.9029154 -0.7175990
## Epsilon:time2 1.3255168000 0.9147279 -0.4673499 3.1183835
## Epsilon:time3 1.7841484000 0.8846823 0.0501712 3.5181257
## Epsilon:time4 -0.0769521000 1.1160191 -2.2643495 2.1104453
## Epsilon:time5 0.9076214000 0.9464242 -0.9473700 2.7626129
## Epsilon:time6 2.1019332000 0.8898330 0.3578605 3.8460059
## Epsilon:time7 1.7793104000 0.9777432 -0.1370662 3.6956871
## Epsilon:time8 1.4738256000 0.9531283 -0.3943060 3.3419571
## Epsilon:time9 2.0096698000 0.9260584 0.1945953 3.8247442
## Gamma:(Intercept) 0.4967845000 0.7463589 -0.9660790 1.9596481
## Gamma:time2 -2.6410196000 3.0152805 -8.5509695 3.2689304
## Gamma:time3 -3.4414334000 1.2689191 -5.9285148 -0.9543520
## Gamma:time4 -1.5811830000 0.8514516 -3.2500282 0.0876621
## Gamma:time5 -1.2970472000 0.8614353 -2.9854605 0.3913661
## Gamma:time6 -1.8751978000 0.9133852 -3.6654329 -0.0849627
## Gamma:time7 -1.7071197000 0.8672590 -3.4069473 -0.0072922
## Gamma:time8 -1.9943784000 0.9116515 -3.7812154 -0.2075414
## Gamma:time9 -1.3691970000 0.8371373 -3.0099862 0.2715921
## p:(Intercept) 0.0386232000 0.1079577 -0.1729740 0.2502203
## p:session2 -1.3130052000 0.1825759 -1.6708540 -0.9551564
## p:session3 1.0035789000 0.1737524 0.6630242 1.3441337
## p:session4 1.1187878000 0.2037933 0.7193529 1.5182227
## p:session5 -0.2903862000 0.1722299 -0.6279569 0.0471844
## p:session6 -0.3212176000 0.1711305 -0.6566334 0.0141981
## p:session7 -0.0007891196 0.1871691 -0.3676406 0.3660623
## p:session8 -1.1042485000 0.2414087 -1.5774095 -0.6310875
## p:session9 -0.0960527000 0.1905299 -0.4694914 0.2773860
## p:session10 -0.1171112000 0.1876904 -0.4849844 0.2507621
##
##
## Real Parameter Psi
##
## 1
## 0.8026827
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.090277 0.2719522 0.3714249 0.0841533 0.1973982 0.4481065 0.3702961
## 8 9
## 0.3022869 0.4254139
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7
## 0.6217034 0.1048712 0.04999 0.2526745 0.3099693 0.201264 0.2296417
## 8 9
## 0.1827847 0.2947526
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5096546 0.5096546 0.5096546 0.5096546 0.5096546 0.5096546 0.5096546
## 8
## 0.5096546
##
## Session:2
## 1 2 3 4 5 6 7 8
## 0.218508 0.218508 0.218508 0.218508 0.218508 0.218508 0.218508 0.218508
##
## Session:3
## 1 2 3 4 5 6 7
## 0.7392747 0.7392747 0.7392747 0.7392747 0.7392747 0.7392747 0.7392747
## 8
## 0.7392747
##
## Session:4
## 1 2 3 4 5 6 7 8
## 0.760862 0.760862 0.760862 0.760862 0.760862 0.760862 0.760862 0.760862
##
## Session:5
## 1 2 3 4 5 6 7
## 0.4373896 0.4373896 0.4373896 0.4373896 0.4373896 0.4373896 0.4373896
## 8
## 0.4373896
##
## Session:6
## 1 2 3 4 5 6 7
## 0.4298178 0.4298178 0.4298178 0.4298178 0.4298178 0.4298178 0.4298178
## 8
## 0.4298178
##
## Session:7
## 1 2 3 4 5 6 7
## 0.5094574 0.5094574 0.5094574 0.5094574 0.5094574 0.5094574 0.5094574
## 8
## 0.5094574
##
## Session:8
## 1 2 3 4 5 6 7
## 0.2562359 0.2562359 0.2562359 0.2562359 0.2562359 0.2562359 0.2562359
## 8
## 0.2562359
##
## Session:9
## 1 2 3 4 5 6 7
## 0.4856466 0.4856466 0.4856466 0.4856466 0.4856466 0.4856466 0.4856466
## 8
## 0.4856466
##
## Session:10
## 1 2 3 4 5 6 7
## 0.4803881 0.4803881 0.4803881 0.4803881 0.4803881 0.4803881 0.4803881
## 8
## 0.4803881
# 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 : 13
## -2lnL: 3566.397
## AICc : 3593.077
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 1.4171968 0.3457931 0.7394424 2.0949512
## Epsilon:(Intercept) -0.9217753 0.1417699 -1.1996443 -0.6439063
## Gamma:(Intercept) -1.1269845 0.1601751 -1.4409276 -0.8130413
## p:(Intercept) 0.0330090 0.1090620 -0.1807525 0.2467706
## p:session2 -1.1916082 0.1728642 -1.5304220 -0.8527944
## p:session3 1.0092192 0.1744336 0.6673293 1.3511090
## p:session4 1.1243693 0.2043931 0.7237588 1.5249798
## p:session5 -0.2825578 0.1725432 -0.6207426 0.0556269
## p:session6 -0.3061072 0.1701695 -0.6396395 0.0274251
## p:session7 -0.0019924 0.1893119 -0.3730439 0.3690590
## p:session8 -1.1188374 0.2294534 -1.5685660 -0.6691089
## p:session9 -0.0942170 0.1919391 -0.4704176 0.2819835
## p:session10 -0.1154180 0.1891406 -0.4861336 0.2552977
##
##
## Real Parameter Psi
##
## 1
## 0.8048986
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.2845963 0.2845963 0.2845963 0.2845963 0.2845963 0.2845963 0.2845963
## 8 9
## 0.2845963 0.2845963
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7 8
## 0.244718 0.244718 0.244718 0.244718 0.244718 0.244718 0.244718 0.244718
## 9
## 0.244718
##
##
## Real Parameter p
## Session:1
## 1 2 3 4 5 6 7
## 0.5082515 0.5082515 0.5082515 0.5082515 0.5082515 0.5082515 0.5082515
## 8
## 0.5082515
##
## Session:2
## 1 2 3 4 5 6 7
## 0.2389219 0.2389219 0.2389219 0.2389219 0.2389219 0.2389219 0.2389219
## 8
## 0.2389219
##
## Session:3
## 1 2 3 4 5 6 7
## 0.7392797 0.7392797 0.7392797 0.7392797 0.7392797 0.7392797 0.7392797
## 8
## 0.7392797
##
## Session:4
## 1 2 3 4 5 6 7 8
## 0.760856 0.760856 0.760856 0.760856 0.760856 0.760856 0.760856 0.760856
##
## Session:5
## 1 2 3 4 5 6 7
## 0.4379346 0.4379346 0.4379346 0.4379346 0.4379346 0.4379346 0.4379346
## 8
## 0.4379346
##
## Session:6
## 1 2 3 4 5 6 7
## 0.4321466 0.4321466 0.4321466 0.4321466 0.4321466 0.4321466 0.4321466
## 8
## 0.4321466
##
## Session:7
## 1 2 3 4 5 6 7
## 0.5077535 0.5077535 0.5077535 0.5077535 0.5077535 0.5077535 0.5077535
## 8
## 0.5077535
##
## Session:8
## 1 2 3 4 5 6 7
## 0.2524046 0.2524046 0.2524046 0.2524046 0.2524046 0.2524046 0.2524046
## 8
## 0.2524046
##
## Session:9
## 1 2 3 4 5 6 7
## 0.4847028 0.4847028 0.4847028 0.4847028 0.4847028 0.4847028 0.4847028
## 8
## 0.4847028
##
## Session:10
## 1 2 3 4 5 6 7
## 0.4794094 0.4794094 0.4794094 0.4794094 0.4794094 0.4794094 0.4794094
## 8
## 0.4794094
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## Psi g1 a0 t1 0.8048986 0.0543023 0.6768739 0.8904115
## Epsilon g1 a0 t1 0.2845963 0.0288645 0.2315385 0.3443641
## Gamma g1 a0 t1 0.2447180 0.0296053 0.1914017 0.3072428
## p g1 s1 t1 0.5082515 0.0272581 0.4549345 0.5613815
## p g1 s2 t1 0.2389219 0.0244116 0.1943910 0.2899820
## p g1 s3 t1 0.7392797 0.0262392 0.6846889 0.7873544
## p g1 s4 t1 0.7608560 0.0314533 0.6939316 0.8170059
## p g1 s5 t1 0.4379346 0.0329098 0.3748187 0.5031253
## p g1 s6 t1 0.4321466 0.0320578 0.3707159 0.4957377
## p g1 s7 t1 0.5077535 0.0386790 0.4323425 0.5828133
## p g1 s8 t1 0.2524046 0.0381034 0.1851843 0.3340227
## p g1 s9 t1 0.4847028 0.0394506 0.4083501 0.5617763
## p g1 s10 t1 0.4794094 0.0385668 0.4048526 0.5548949
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## Psi:(Intercept) 1.4171968 0.3457931 0.7394424 2.0949512
## Epsilon:(Intercept) -0.9217753 0.1417699 -1.1996443 -0.6439063
## Gamma:(Intercept) -1.1269845 0.1601751 -1.4409276 -0.8130413
## p:(Intercept) 0.0330090 0.1090620 -0.1807525 0.2467706
## p:session2 -1.1916082 0.1728642 -1.5304220 -0.8527944
## p:session3 1.0092192 0.1744336 0.6673293 1.3511090
## p:session4 1.1243693 0.2043931 0.7237588 1.5249798
## p:session5 -0.2825578 0.1725432 -0.6207426 0.0556269
## p:session6 -0.3061072 0.1701695 -0.6396395 0.0274251
## p:session7 -0.0019924 0.1893119 -0.3730439 0.3690590
## p:session8 -1.1188374 0.2294534 -1.5685660 -0.6691089
## p:session9 -0.0942170 0.1919391 -0.4704176 0.2819835
## p:session10 -0.1154180 0.1891406 -0.4861336 0.2552977
model.fits[[model.number]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.8048986 0.05430225 0.6984662 0.9113310
## 2 0.6235723 0.03478308 0.5553974 0.6917471
## 3 0.5382245 0.03348381 0.4725963 0.6038528
## 4 0.4980526 0.03482287 0.4297998 0.5663054
## 5 0.4791443 0.03607055 0.4084460 0.5498425
## 6 0.4702444 0.03690755 0.3979056 0.5425832
## 7 0.4660553 0.03740399 0.3927435 0.5393671
## 8 0.4640836 0.03767857 0.3902336 0.5379336
## 9 0.4631555 0.03782397 0.3890205 0.5372905
## 10 0.4627187 0.03789885 0.3884369 0.5370004
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 0.7747215 0.0357469083 0.7046576 0.8447855
## 2 0.8631310 0.0253220195 0.8134999 0.9127622
## 3 0.9253621 0.0174611120 0.8911383 0.9595859
## 4 0.9620354 0.0114210969 0.9396501 0.9844208
## 5 0.9814254 0.0070037914 0.9676980 0.9951529
## 6 0.9910918 0.0040766966 0.9831014 0.9990821
## 7 0.9957693 0.0022864650 0.9912879 1.0002508
## 8 0.9980002 0.0012492620 0.9955517 1.0004488
## 9 0.9990568 0.0006696603 0.9977443 1.0003694
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 0.4015359 0.103297251 0.1990733 0.6039985
## 2 0.7036027 0.049225447 0.6071208 0.8000846
## 3 0.8513034 0.030703752 0.7911240 0.9114827
## 4 0.9271112 0.019809109 0.8882853 0.9659370
## 5 0.9649375 0.012223197 0.9409800 0.9888950
## 6 0.9833162 0.007183656 0.9692362 0.9973961
## 7 0.9921057 0.004065852 0.9841367 1.0000748
## 8 0.9962749 0.002238512 0.9918875 1.0006624
## 9 0.9982446 0.001207345 0.9958782 1.0006110
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.8048986 0.05430225 0.6984662 0.9113310
## 2 0.6235723 0.03478308 0.5553974 0.6917471
## 3 0.5382245 0.03348381 0.4725963 0.6038528
## 4 0.4980526 0.03482287 0.4297998 0.5663054
## 5 0.4791443 0.03607055 0.4084460 0.5498425
## 6 0.4702444 0.03690755 0.3979056 0.5425832
## 7 0.4660553 0.03740399 0.3927435 0.5393671
## 8 0.4640836 0.03767857 0.3902336 0.5379336
## 9 0.4631555 0.03782397 0.3890205 0.5372905
## 10 0.4627187 0.03789885 0.3884369 0.5370004
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 0.7747215 0.0357469083 0.7046576 0.8447855
## 2 0.8631310 0.0253220195 0.8134999 0.9127622
## 3 0.9253621 0.0174611120 0.8911383 0.9595859
## 4 0.9620354 0.0114210969 0.9396501 0.9844208
## 5 0.9814254 0.0070037914 0.9676980 0.9951529
## 6 0.9910918 0.0040766966 0.9831014 0.9990821
## 7 0.9957693 0.0022864650 0.9912879 1.0002508
## 8 0.9980002 0.0012492620 0.9955517 1.0004488
## 9 0.9990568 0.0006696603 0.9977443 1.0003694
model.fits[[model.number]]$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 0.4015359 0.103297251 0.1990733 0.6039985
## 2 0.7036027 0.049225447 0.6071208 0.8000846
## 3 0.8513034 0.030703752 0.7911240 0.9114827
## 4 0.9271112 0.019809109 0.8882853 0.9659370
## 5 0.9649375 0.012223197 0.9409800 0.9888950
## 6 0.9833162 0.007183656 0.9692362 0.9973961
## 7 0.9921057 0.004065852 0.9841367 1.0000748
## 8 0.9962749 0.002238512 0.9918875 1.0006624
## 9 0.9982446 0.001207345 0.9958782 1.0006110
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 1 1 0.8048986 0.0543023 0.6768739
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.8904115 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 20 4 0.5082515 0.0272581 0.4549345
## p g1 s1 t2 21 4 0.5082515 0.0272581 0.4549345
## p g1 s1 t3 22 4 0.5082515 0.0272581 0.4549345
## p g1 s1 t4 23 4 0.5082515 0.0272581 0.4549345
## p g1 s1 t5 24 4 0.5082515 0.0272581 0.4549345
## p g1 s1 t6 25 4 0.5082515 0.0272581 0.4549345
## p g1 s1 t7 26 4 0.5082515 0.0272581 0.4549345
## p g1 s1 t8 27 4 0.5082515 0.0272581 0.4549345
## p g1 s2 t1 28 5 0.2389219 0.0244116 0.1943910
## p g1 s2 t2 29 5 0.2389219 0.0244116 0.1943910
## p g1 s2 t3 30 5 0.2389219 0.0244116 0.1943910
## p g1 s2 t4 31 5 0.2389219 0.0244116 0.1943910
## p g1 s2 t5 32 5 0.2389219 0.0244116 0.1943910
## p g1 s2 t6 33 5 0.2389219 0.0244116 0.1943910
## p g1 s2 t7 34 5 0.2389219 0.0244116 0.1943910
## p g1 s2 t8 35 5 0.2389219 0.0244116 0.1943910
## p g1 s3 t1 36 6 0.7392797 0.0262392 0.6846889
## p g1 s3 t2 37 6 0.7392797 0.0262392 0.6846889
## p g1 s3 t3 38 6 0.7392797 0.0262392 0.6846889
## p g1 s3 t4 39 6 0.7392797 0.0262392 0.6846889
## p g1 s3 t5 40 6 0.7392797 0.0262392 0.6846889
## p g1 s3 t6 41 6 0.7392797 0.0262392 0.6846889
## p g1 s3 t7 42 6 0.7392797 0.0262392 0.6846889
## p g1 s3 t8 43 6 0.7392797 0.0262392 0.6846889
## p g1 s4 t1 44 7 0.7608560 0.0314533 0.6939316
## p g1 s4 t2 45 7 0.7608560 0.0314533 0.6939316
## p g1 s4 t3 46 7 0.7608560 0.0314533 0.6939316
## p g1 s4 t4 47 7 0.7608560 0.0314533 0.6939316
## p g1 s4 t5 48 7 0.7608560 0.0314533 0.6939316
## p g1 s4 t6 49 7 0.7608560 0.0314533 0.6939316
## p g1 s4 t7 50 7 0.7608560 0.0314533 0.6939316
## p g1 s4 t8 51 7 0.7608560 0.0314533 0.6939316
## p g1 s5 t1 52 8 0.4379346 0.0329098 0.3748187
## p g1 s5 t2 53 8 0.4379346 0.0329098 0.3748187
## p g1 s5 t3 54 8 0.4379346 0.0329098 0.3748187
## p g1 s5 t4 55 8 0.4379346 0.0329098 0.3748187
## p g1 s5 t5 56 8 0.4379346 0.0329098 0.3748187
## p g1 s5 t6 57 8 0.4379346 0.0329098 0.3748187
## p g1 s5 t7 58 8 0.4379346 0.0329098 0.3748187
## p g1 s5 t8 59 8 0.4379346 0.0329098 0.3748187
## p g1 s6 t1 60 9 0.4321466 0.0320578 0.3707159
## p g1 s6 t2 61 9 0.4321466 0.0320578 0.3707159
## p g1 s6 t3 62 9 0.4321466 0.0320578 0.3707159
## p g1 s6 t4 63 9 0.4321466 0.0320578 0.3707159
## p g1 s6 t5 64 9 0.4321466 0.0320578 0.3707159
## p g1 s6 t6 65 9 0.4321466 0.0320578 0.3707159
## p g1 s6 t7 66 9 0.4321466 0.0320578 0.3707159
## p g1 s6 t8 67 9 0.4321466 0.0320578 0.3707159
## p g1 s7 t1 68 10 0.5077535 0.0386790 0.4323425
## p g1 s7 t2 69 10 0.5077535 0.0386790 0.4323425
## p g1 s7 t3 70 10 0.5077535 0.0386790 0.4323425
## p g1 s7 t4 71 10 0.5077535 0.0386790 0.4323425
## p g1 s7 t5 72 10 0.5077535 0.0386790 0.4323425
## p g1 s7 t6 73 10 0.5077535 0.0386790 0.4323425
## p g1 s7 t7 74 10 0.5077535 0.0386790 0.4323425
## p g1 s7 t8 75 10 0.5077535 0.0386790 0.4323425
## p g1 s8 t1 76 11 0.2524046 0.0381034 0.1851843
## p g1 s8 t2 77 11 0.2524046 0.0381034 0.1851843
## p g1 s8 t3 78 11 0.2524046 0.0381034 0.1851843
## p g1 s8 t4 79 11 0.2524046 0.0381034 0.1851843
## p g1 s8 t5 80 11 0.2524046 0.0381034 0.1851843
## p g1 s8 t6 81 11 0.2524046 0.0381034 0.1851843
## p g1 s8 t7 82 11 0.2524046 0.0381034 0.1851843
## p g1 s8 t8 83 11 0.2524046 0.0381034 0.1851843
## p g1 s9 t1 84 12 0.4847028 0.0394506 0.4083501
## p g1 s9 t2 85 12 0.4847028 0.0394506 0.4083501
## p g1 s9 t3 86 12 0.4847028 0.0394506 0.4083501
## p g1 s9 t4 87 12 0.4847028 0.0394506 0.4083501
## p g1 s9 t5 88 12 0.4847028 0.0394506 0.4083501
## p g1 s9 t6 89 12 0.4847028 0.0394506 0.4083501
## p g1 s9 t7 90 12 0.4847028 0.0394506 0.4083501
## p g1 s9 t8 91 12 0.4847028 0.0394506 0.4083501
## p g1 s10 t1 92 13 0.4794094 0.0385668 0.4048526
## p g1 s10 t2 93 13 0.4794094 0.0385668 0.4048526
## p g1 s10 t3 94 13 0.4794094 0.0385668 0.4048526
## p g1 s10 t4 95 13 0.4794094 0.0385668 0.4048526
## p g1 s10 t5 96 13 0.4794094 0.0385668 0.4048526
## p g1 s10 t6 97 13 0.4794094 0.0385668 0.4048526
## p g1 s10 t7 98 13 0.4794094 0.0385668 0.4048526
## p g1 s10 t8 99 13 0.4794094 0.0385668 0.4048526
## ucl fixed note group time session Time
## p g1 s1 t1 0.5613815 1 1 1 0
## p g1 s1 t2 0.5613815 1 2 1 1
## p g1 s1 t3 0.5613815 1 3 1 2
## p g1 s1 t4 0.5613815 1 4 1 3
## p g1 s1 t5 0.5613815 1 5 1 4
## p g1 s1 t6 0.5613815 1 6 1 5
## p g1 s1 t7 0.5613815 1 7 1 6
## p g1 s1 t8 0.5613815 1 8 1 7
## p g1 s2 t1 0.2899820 1 1 2 0
## p g1 s2 t2 0.2899820 1 2 2 1
## p g1 s2 t3 0.2899820 1 3 2 2
## p g1 s2 t4 0.2899820 1 4 2 3
## p g1 s2 t5 0.2899820 1 5 2 4
## p g1 s2 t6 0.2899820 1 6 2 5
## p g1 s2 t7 0.2899820 1 7 2 6
## p g1 s2 t8 0.2899820 1 8 2 7
## p g1 s3 t1 0.7873544 1 1 3 0
## p g1 s3 t2 0.7873544 1 2 3 1
## p g1 s3 t3 0.7873544 1 3 3 2
## p g1 s3 t4 0.7873544 1 4 3 3
## p g1 s3 t5 0.7873544 1 5 3 4
## p g1 s3 t6 0.7873544 1 6 3 5
## p g1 s3 t7 0.7873544 1 7 3 6
## p g1 s3 t8 0.7873544 1 8 3 7
## p g1 s4 t1 0.8170059 1 1 4 0
## p g1 s4 t2 0.8170059 1 2 4 1
## p g1 s4 t3 0.8170059 1 3 4 2
## p g1 s4 t4 0.8170059 1 4 4 3
## p g1 s4 t5 0.8170059 1 5 4 4
## p g1 s4 t6 0.8170059 1 6 4 5
## p g1 s4 t7 0.8170059 1 7 4 6
## p g1 s4 t8 0.8170059 1 8 4 7
## p g1 s5 t1 0.5031253 1 1 5 0
## p g1 s5 t2 0.5031253 1 2 5 1
## p g1 s5 t3 0.5031253 1 3 5 2
## p g1 s5 t4 0.5031253 1 4 5 3
## p g1 s5 t5 0.5031253 1 5 5 4
## p g1 s5 t6 0.5031253 1 6 5 5
## p g1 s5 t7 0.5031253 1 7 5 6
## p g1 s5 t8 0.5031253 1 8 5 7
## p g1 s6 t1 0.4957377 1 1 6 0
## p g1 s6 t2 0.4957377 1 2 6 1
## p g1 s6 t3 0.4957377 1 3 6 2
## p g1 s6 t4 0.4957377 1 4 6 3
## p g1 s6 t5 0.4957377 1 5 6 4
## p g1 s6 t6 0.4957377 1 6 6 5
## p g1 s6 t7 0.4957377 1 7 6 6
## p g1 s6 t8 0.4957377 1 8 6 7
## p g1 s7 t1 0.5828133 1 1 7 0
## p g1 s7 t2 0.5828133 1 2 7 1
## p g1 s7 t3 0.5828133 1 3 7 2
## p g1 s7 t4 0.5828133 1 4 7 3
## p g1 s7 t5 0.5828133 1 5 7 4
## p g1 s7 t6 0.5828133 1 6 7 5
## p g1 s7 t7 0.5828133 1 7 7 6
## p g1 s7 t8 0.5828133 1 8 7 7
## p g1 s8 t1 0.3340227 1 1 8 0
## p g1 s8 t2 0.3340227 1 2 8 1
## p g1 s8 t3 0.3340227 1 3 8 2
## p g1 s8 t4 0.3340227 1 4 8 3
## p g1 s8 t5 0.3340227 1 5 8 4
## p g1 s8 t6 0.3340227 1 6 8 5
## p g1 s8 t7 0.3340227 1 7 8 6
## p g1 s8 t8 0.3340227 1 8 8 7
## p g1 s9 t1 0.5617763 1 1 9 0
## p g1 s9 t2 0.5617763 1 2 9 1
## p g1 s9 t3 0.5617763 1 3 9 2
## p g1 s9 t4 0.5617763 1 4 9 3
## p g1 s9 t5 0.5617763 1 5 9 4
## p g1 s9 t6 0.5617763 1 6 9 5
## p g1 s9 t7 0.5617763 1 7 9 6
## p g1 s9 t8 0.5617763 1 8 9 7
## p g1 s10 t1 0.5548949 1 1 10 0
## p g1 s10 t2 0.5548949 1 2 10 1
## p g1 s10 t3 0.5548949 1 3 10 2
## p g1 s10 t4 0.5548949 1 4 10 3
## p g1 s10 t5 0.5548949 1 5 10 4
## p g1 s10 t6 0.5548949 1 6 10 5
## p g1 s10 t7 0.5548949 1 7 10 6
## p g1 s10 t8 0.5548949 1 8 10 7
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
model.set
## model npar AICc DeltaAICc
## 2 Psi(~1)Epsilon(~1)Gamma(~1)p(~session) 13 3593.077 0.000000
## 3 Psi(~1)Epsilon(~time)Gamma(~time)p(~session) 29 3595.340 2.263749
## 1 Psi(~1)Epsilon(~1)Gamma(~1)p(~1) 4 3818.972 225.895790
## weight Deviance
## 2 0.7561847 3125.591
## 3 0.2438153 3093.187
## 1 0.0000000 3370.092
# model averaging in the usual way
RMark::model.average(model.set, "Psi")
## par.index estimate se fixed note group age time
## Psi g1 a0 t1 1 0.8043583 0.05426861 1 0 1
## Age Time
## Psi g1 a0 t1 0 0
RMark::model.average(model.set, "p")
## par.index estimate se fixed note group time
## p g1 s1 t1 20 0.5085936 0.02719706 1 1
## p g1 s1 t2 21 0.5085936 0.02719706 1 2
## p g1 s1 t3 22 0.5085936 0.02719706 1 3
## p g1 s1 t4 23 0.5085936 0.02719706 1 4
## p g1 s1 t5 24 0.5085936 0.02719706 1 5
## p g1 s1 t6 25 0.5085936 0.02719706 1 6
## p g1 s1 t7 26 0.5085936 0.02719706 1 7
## p g1 s1 t8 27 0.5085936 0.02719706 1 8
## p g1 s2 t1 28 0.2339447 0.02610728 1 1
## p g1 s2 t2 29 0.2339447 0.02610728 1 2
## p g1 s2 t3 30 0.2339447 0.02610728 1 3
## p g1 s2 t4 31 0.2339447 0.02610728 1 4
## p g1 s2 t5 32 0.2339447 0.02610728 1 5
## p g1 s2 t6 33 0.2339447 0.02610728 1 6
## p g1 s2 t7 34 0.2339447 0.02610728 1 7
## p g1 s2 t8 35 0.2339447 0.02610728 1 8
## p g1 s3 t1 36 0.7392785 0.02623974 1 1
## p g1 s3 t2 37 0.7392785 0.02623974 1 2
## p g1 s3 t3 38 0.7392785 0.02623974 1 3
## p g1 s3 t4 39 0.7392785 0.02623974 1 4
## p g1 s3 t5 40 0.7392785 0.02623974 1 5
## p g1 s3 t6 41 0.7392785 0.02623974 1 6
## p g1 s3 t7 42 0.7392785 0.02623974 1 7
## p g1 s3 t8 43 0.7392785 0.02623974 1 8
## p g1 s4 t1 44 0.7608575 0.03145255 1 1
## p g1 s4 t2 45 0.7608575 0.03145255 1 2
## p g1 s4 t3 46 0.7608575 0.03145255 1 3
## p g1 s4 t4 47 0.7608575 0.03145255 1 4
## p g1 s4 t5 48 0.7608575 0.03145255 1 5
## p g1 s4 t6 49 0.7608575 0.03145255 1 6
## p g1 s4 t7 50 0.7608575 0.03145255 1 7
## p g1 s4 t8 51 0.7608575 0.03145255 1 8
## p g1 s5 t1 52 0.4378017 0.03293816 1 1
## p g1 s5 t2 53 0.4378017 0.03293816 1 2
## p g1 s5 t3 54 0.4378017 0.03293816 1 3
## p g1 s5 t4 55 0.4378017 0.03293816 1 4
## p g1 s5 t5 56 0.4378017 0.03293816 1 5
## p g1 s5 t6 57 0.4378017 0.03293816 1 6
## p g1 s5 t7 58 0.4378017 0.03293816 1 7
## p g1 s5 t8 59 0.4378017 0.03293816 1 8
## p g1 s6 t1 60 0.4315789 0.03219185 1 1
## p g1 s6 t2 61 0.4315789 0.03219185 1 2
## p g1 s6 t3 62 0.4315789 0.03219185 1 3
## p g1 s6 t4 63 0.4315789 0.03219185 1 4
## p g1 s6 t5 64 0.4315789 0.03219185 1 5
## p g1 s6 t6 65 0.4315789 0.03219185 1 6
## p g1 s6 t7 66 0.4315789 0.03219185 1 7
## p g1 s6 t8 67 0.4315789 0.03219185 1 8
## p g1 s7 t1 68 0.5081690 0.03857226 1 1
## p g1 s7 t2 69 0.5081690 0.03857226 1 2
## p g1 s7 t3 70 0.5081690 0.03857226 1 3
## p g1 s7 t4 71 0.5081690 0.03857226 1 4
## p g1 s7 t5 72 0.5081690 0.03857226 1 5
## p g1 s7 t6 73 0.5081690 0.03857226 1 6
## p g1 s7 t7 74 0.5081690 0.03857226 1 7
## p g1 s7 t8 75 0.5081690 0.03857226 1 8
## p g1 s8 t1 76 0.2533388 0.03890318 1 1
## p g1 s8 t2 77 0.2533388 0.03890318 1 2
## p g1 s8 t3 78 0.2533388 0.03890318 1 3
## p g1 s8 t4 79 0.2533388 0.03890318 1 4
## p g1 s8 t5 80 0.2533388 0.03890318 1 5
## p g1 s8 t6 81 0.2533388 0.03890318 1 6
## p g1 s8 t7 82 0.2533388 0.03890318 1 7
## p g1 s8 t8 83 0.2533388 0.03890318 1 8
## p g1 s9 t1 84 0.4849329 0.03939556 1 1
## p g1 s9 t2 85 0.4849329 0.03939556 1 2
## p g1 s9 t3 86 0.4849329 0.03939556 1 3
## p g1 s9 t4 87 0.4849329 0.03939556 1 4
## p g1 s9 t5 88 0.4849329 0.03939556 1 5
## p g1 s9 t6 89 0.4849329 0.03939556 1 6
## p g1 s9 t7 90 0.4849329 0.03939556 1 7
## p g1 s9 t8 91 0.4849329 0.03939556 1 8
## p g1 s10 t1 92 0.4796480 0.03851016 1 1
## p g1 s10 t2 93 0.4796480 0.03851016 1 2
## p g1 s10 t3 94 0.4796480 0.03851016 1 3
## p g1 s10 t4 95 0.4796480 0.03851016 1 4
## p g1 s10 t5 96 0.4796480 0.03851016 1 5
## p g1 s10 t6 97 0.4796480 0.03851016 1 6
## p g1 s10 t7 98 0.4796480 0.03851016 1 7
## p g1 s10 t8 99 0.4796480 0.03851016 1 8
## session Time
## p g1 s1 t1 1 0
## p g1 s1 t2 1 1
## p g1 s1 t3 1 2
## p g1 s1 t4 1 3
## p g1 s1 t5 1 4
## p g1 s1 t6 1 5
## p g1 s1 t7 1 6
## p g1 s1 t8 1 7
## p g1 s2 t1 2 0
## p g1 s2 t2 2 1
## p g1 s2 t3 2 2
## p g1 s2 t4 2 3
## p g1 s2 t5 2 4
## p g1 s2 t6 2 5
## p g1 s2 t7 2 6
## p g1 s2 t8 2 7
## p g1 s3 t1 3 0
## p g1 s3 t2 3 1
## p g1 s3 t3 3 2
## p g1 s3 t4 3 3
## p g1 s3 t5 3 4
## p g1 s3 t6 3 5
## p g1 s3 t7 3 6
## p g1 s3 t8 3 7
## p g1 s4 t1 4 0
## p g1 s4 t2 4 1
## p g1 s4 t3 4 2
## p g1 s4 t4 4 3
## p g1 s4 t5 4 4
## p g1 s4 t6 4 5
## p g1 s4 t7 4 6
## p g1 s4 t8 4 7
## p g1 s5 t1 5 0
## p g1 s5 t2 5 1
## p g1 s5 t3 5 2
## p g1 s5 t4 5 3
## p g1 s5 t5 5 4
## p g1 s5 t6 5 5
## p g1 s5 t7 5 6
## p g1 s5 t8 5 7
## p g1 s6 t1 6 0
## p g1 s6 t2 6 1
## p g1 s6 t3 6 2
## p g1 s6 t4 6 3
## p g1 s6 t5 6 4
## p g1 s6 t6 6 5
## p g1 s6 t7 6 6
## p g1 s6 t8 6 7
## p g1 s7 t1 7 0
## p g1 s7 t2 7 1
## p g1 s7 t3 7 2
## p g1 s7 t4 7 3
## p g1 s7 t5 7 4
## p g1 s7 t6 7 5
## p g1 s7 t7 7 6
## p g1 s7 t8 7 7
## p g1 s8 t1 8 0
## p g1 s8 t2 8 1
## p g1 s8 t3 8 2
## p g1 s8 t4 8 3
## p g1 s8 t5 8 4
## p g1 s8 t6 8 5
## p g1 s8 t7 8 6
## p g1 s8 t8 8 7
## p g1 s9 t1 9 0
## p g1 s9 t2 9 1
## p g1 s9 t3 9 2
## p g1 s9 t4 9 3
## p g1 s9 t5 9 4
## p g1 s9 t6 9 5
## p g1 s9 t7 9 6
## p g1 s9 t8 9 7
## p g1 s10 t1 10 0
## p g1 s10 t2 10 1
## p g1 s10 t3 10 2
## p g1 s10 t4 10 3
## p g1 s10 t5 10 4
## p g1 s10 t6 10 5
## p g1 s10 t7 10 6
## p g1 s10 t8 10 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)})
# NOTE: Because a group was used in the process.data step above, the first half
# of each derived parameter table will be for Prox=No, and the second
# will be for Prox = Yes
model.set[[1]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.8064294 0.05436625 0.6998716 0.9129873
## 2 0.6111334 0.03374706 0.5449892 0.6772776
## 3 0.5245324 0.03247428 0.4608828 0.5881820
## 4 0.4861306 0.03377543 0.4199307 0.5523304
## 5 0.4691018 0.03492109 0.4006565 0.5375472
## 6 0.4615507 0.03564158 0.3916932 0.5314082
## 7 0.4582023 0.03604050 0.3875629 0.5288417
## 8 0.4567175 0.03624646 0.3856744 0.5277606
## 9 0.4560591 0.03634844 0.3848161 0.5273020
## 10 0.4557671 0.03639766 0.3844277 0.5271065
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 0.7578263 0.0361205511 0.6870300 0.8286226
## 2 0.8582945 0.0254302233 0.8084512 0.9081377
## 3 0.9267884 0.0170986289 0.8932751 0.9603017
## 4 0.9649709 0.0106762572 0.9440455 0.9858964
## 5 0.9839030 0.0061833101 0.9717837 0.9960223
## 6 0.9927453 0.0033886004 0.9861036 0.9993869
## 7 0.9967595 0.0017887345 0.9932536 1.0002654
## 8 0.9985584 0.0009200834 0.9967550 1.0003617
## 9 0.9993598 0.0004644587 0.9984495 1.0002702
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 0.3772319 0.1003416557 0.1805622 0.5739015
## 2 0.7019659 0.0478415780 0.6081964 0.7957354
## 3 0.8575288 0.0295032549 0.7997024 0.9153552
## 4 0.9340192 0.0182909121 0.8981690 0.9698694
## 5 0.9701050 0.0106865935 0.9491592 0.9910507
## 6 0.9866099 0.0059184067 0.9750098 0.9982100
## 7 0.9940353 0.0031530440 0.9878554 1.0002153
## 8 0.9973497 0.0016339140 0.9941472 1.0005521
## 9 0.9988237 0.0008295864 0.9971977 1.0004497
##
## $all_psi
## estimate se lcl ucl
## 1 0.8064294 0.05436625 0.6998716 0.9129873
## 2 0.6111334 0.03374706 0.5449892 0.6772776
## 3 0.5245324 0.03247428 0.4608828 0.5881820
## 4 0.4861306 0.03377543 0.4199307 0.5523304
## 5 0.4691018 0.03492109 0.4006565 0.5375472
## 6 0.4615507 0.03564158 0.3916932 0.5314082
## 7 0.4582023 0.03604050 0.3875629 0.5288417
## 8 0.4567175 0.03624646 0.3856744 0.5277606
## 9 0.4560591 0.03634844 0.3848161 0.5273020
## 10 0.4557671 0.03639766 0.3844277 0.5271065
RMark.model.average.derived(model.set, "all_psi")
## Loading required package: plyr
## Loading required package: boot
## estimate se lcl ucl
## 1 0.8043583 0.05426861 0.6765265 0.8898946
## 2 0.6794870 0.10874275 0.4434377 0.8494191
## 3 0.5621561 0.06041388 0.4424676 0.6750219
## 4 0.4785788 0.05631776 0.3709666 0.5882166
## 5 0.4915453 0.05077127 0.3936548 0.5900885
## 6 0.4948277 0.06345185 0.3732741 0.6169959
## 7 0.4503140 0.05383814 0.3484874 0.5564806
## 8 0.4460846 0.05884892 0.3355256 0.5622482
## 9 0.4437925 0.05746169 0.3358059 0.5573646
## 10 0.4479471 0.05327823 0.3472296 0.5531241
RMark.model.average.derived(model.set, "lambda Rate of Change")
## estimate se lcl ucl
## 1 0.8449029 0.13678674 0.5768058 1.1129999
## 2 0.8346043 0.06809584 0.7011389 0.9680697
## 3 0.8599617 0.12379990 0.6173183 1.1026050
## 4 1.0364892 0.15236857 0.7378523 1.3351261
## 5 1.0048453 0.07857538 0.8508404 1.1588502
## 6 0.9208624 0.13658238 0.6531659 1.1885589
## 7 0.9899801 0.09920378 0.7955443 1.1844159
## 8 0.9944131 0.09360450 0.8109516 1.1778745
## 9 1.0109739 0.10702913 0.8012007 1.2207472
RMark.model.average.derived(model.set, "log odds lambda")
## estimate se lcl ucl
## 1 0.6511369 0.6119301 -0.5482242 1.8504979
## 2 0.6056451 0.1954195 0.2226299 0.9886603
## 3 0.7438725 0.1981265 0.3555517 1.1321933
## 4 1.0836015 0.3275662 0.4415836 1.7256195
## 5 1.0175283 0.1870411 0.6509344 1.3841222
## 6 0.8664133 0.2204494 0.4343405 1.2984862
## 7 0.9828441 0.1608442 0.6675953 1.2980929
## 8 0.9905711 0.1512557 0.6941153 1.2870269
## 9 1.0182207 0.1815792 0.6623321 1.3741094
# model averaging of derived parameters such as the occupancy at each time step
psi.ma <- RMark.model.average.derived(model.set, "all_psi")
psi.ma$Year <- 1:nrow(psi.ma)
psi.ma$parameter <- 'psi'
psi.ma
## estimate se lcl ucl Year parameter
## 1 0.8043583 0.05426861 0.6765265 0.8898946 1 psi
## 2 0.6794870 0.10874275 0.4434377 0.8494191 2 psi
## 3 0.5621561 0.06041388 0.4424676 0.6750219 3 psi
## 4 0.4785788 0.05631776 0.3709666 0.5882166 4 psi
## 5 0.4915453 0.05077127 0.3936548 0.5900885 5 psi
## 6 0.4948277 0.06345185 0.3732741 0.6169959 6 psi
## 7 0.4503140 0.05383814 0.3484874 0.5564806 7 psi
## 8 0.4460846 0.05884892 0.3355256 0.5622482 8 psi
## 9 0.4437925 0.05746169 0.3358059 0.5573646 9 psi
## 10 0.4479471 0.05327823 0.3472296 0.5531241 10 psi
# likely more interested in colonization and extinction probabilities
epsilon.ma <- RMark::model.average(model.set, "Epsilon", vcv=TRUE)$estimates
epsilon.ma$Year <- 1:nrow(epsilon.ma)
epsilon.ma$parameter <- 'epsilon'
epsilon.ma
## par.index estimate se lcl ucl fixed
## Epsilon g1 a0 t1 2 0.2372183 0.09315391 0.10183381 0.4603413
## Epsilon g1 a1 t2 3 0.2815135 0.04278359 0.20561773 0.3722941
## Epsilon g1 a2 t3 4 0.3057664 0.06038582 0.20140406 0.4347660
## Epsilon g1 a3 t4 5 0.2357252 0.09426047 0.09959533 0.4623731
## Epsilon g1 a4 t5 6 0.2633361 0.05892969 0.16464131 0.3933364
## Epsilon g1 a5 t6 7 0.3244626 0.08672015 0.18111833 0.5105265
## Epsilon g1 a6 t7 8 0.3054912 0.07683739 0.17783324 0.4721612
## Epsilon g1 a7 t8 9 0.2889095 0.05813098 0.18919153 0.4143285
## Epsilon g1 a8 t9 10 0.3189298 0.08461770 0.17913514 0.5012072
## note group age time Age Time Year parameter
## Epsilon g1 a0 t1 1 0 1 0 0 1 epsilon
## Epsilon g1 a1 t2 1 1 2 1 1 2 epsilon
## Epsilon g1 a2 t3 1 2 3 2 2 3 epsilon
## Epsilon g1 a3 t4 1 3 4 3 3 4 epsilon
## Epsilon g1 a4 t5 1 4 5 4 4 5 epsilon
## Epsilon g1 a5 t6 1 5 6 5 5 6 epsilon
## Epsilon g1 a6 t7 1 6 7 6 6 7 epsilon
## Epsilon g1 a7 t8 1 7 8 7 7 8 epsilon
## Epsilon g1 a8 t9 1 8 9 8 8 9 epsilon
gamma.ma <- RMark::model.average(model.set, "Gamma", vcv=TRUE)$estimates
gamma.ma$Year <- 1:nrow(gamma.ma)
gamma.ma$parameter <- 'gamma'
gamma.ma
## par.index estimate se lcl ucl fixed
## Gamma g1 a0 t1 11 0.3366328 0.18541180 0.09065496 0.7209124
## Gamma g1 a1 t2 12 0.2106212 0.14013238 0.04865404 0.5819475
## Gamma g1 a2 t3 13 0.1972404 0.09073580 0.07400049 0.4303404
## Gamma g1 a3 t4 14 0.2466579 0.04619805 0.16745076 0.3476851
## Gamma g1 a4 t5 15 0.2606273 0.05925787 0.16172837 0.3917416
## Gamma g1 a5 t6 16 0.2341233 0.05251326 0.14689357 0.3517921
## Gamma g1 a6 t7 17 0.2410422 0.04683228 0.16128216 0.3440653
## Gamma g1 a7 t8 18 0.2296177 0.05348732 0.14147339 0.3502738
## Gamma g1 a8 t9 19 0.2569172 0.05136910 0.16946425 0.3694264
## note group age time Age Time Year parameter
## Gamma g1 a0 t1 1 0 1 0 0 1 gamma
## Gamma g1 a1 t2 1 1 2 1 1 2 gamma
## Gamma g1 a2 t3 1 2 3 2 2 3 gamma
## Gamma g1 a3 t4 1 3 4 3 3 4 gamma
## Gamma g1 a4 t5 1 4 5 4 4 5 gamma
## Gamma g1 a5 t6 1 5 6 5 5 6 gamma
## Gamma g1 a6 t7 1 6 7 6 6 7 gamma
## Gamma g1 a7 t8 1 7 8 7 7 8 gamma
## Gamma g1 a8 t9 1 8 9 8 8 9 gamma
all.est <- plyr::rbind.fill(psi.ma, epsilon.ma, gamma.ma)
ggplot(data=all.est, aes(x=Year,y=estimate, color=parameter))+
ggtitle("Estimated occupancy, extinction, colonization, over time")+
geom_point(position=position_dodge(w=0.2))+
geom_line(position=position_dodge(w=0.2))+
ylim(0,1)+
geom_errorbar(aes(ymin=lcl, ymax=ucl), width=.1,position=position_dodge(w=0.2))+
scale_x_continuous(breaks=1:10)
cleanup(ask=FALSE)