# multi season multi-state model
# Robust design, Multi State, Reproduction parameterization
# Using RMark
# Visits to hawk territories to establish occupancy and breeding
# 6 years x up to 9 visits/year
# 2018-09-27 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library(car)
## Loading required package: carData
library(ggplot2)
library(readxl)
library(reshape2)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(RMark)
# Robust design, Multi State, Reproduction parameterization
setup.parameters("RDMSOccRepro", check=TRUE) # returns a vector of parameter names (case sensitive)
## [1] "Phi0" "Psi" "R" "p" "Delta"
# phi0 = initial occupancy probabilities in each state parameterized as phi0[1]*(1-phi0[1]), phi0[1](phi0[2]])
# 1 - Phi0[1]
# Phi0[1](1 - Phi0[2])
# Phi0[1]*Phi0[2]
# psi = p(occupied in year t+1 | not occupied, occupied but not breeding, occupied and breeding in year t)
# r = p(reproduction in year t+1 | not occupied, occupied but not breeding, occupied and breeding in year t+1)
# p = p(detection in year t | state 1, state 2 in yar t
# delta=p(detect in state 2 | in state 2)
# The parameters Psi[0], Psi[1], Psi[2], R[0,2], R[1,2], and R[2,2] are used to construct
# the Phi_t matrix for transitions each interval (midway down second column page 826 of MacKenzie et al. 2009):
#
# 1 - Psi[0](t) Psi[0](t)*{1 - R[0,2](t)} Psi[0](t)*R[0,2](t)
# 1 - Psi[1](t) Psi[1](t)*{1 - R[1.2](t)} Psi[1](t)*R[1,2](t)
# 1 - Psi[2](t) Psi[2](t)*{1 - R[2,2](t)} Psi[2](t)*R[2,2](t)
# These 3 PIMs for each primary interval are used to construct the detection matrix
# (top of first column page 826 MacKenzie et al. 2009):
# Observed State
# 0 1 2
#True State 0| 1 0 0
# 1| 1 - p[1] p[1] 0
# 2| 1 - p[2] p[2]*(1 - Delta[2,2]) p[2]*Delta[2,2]
input.history <- read.csv(file.path("..","hawk.csv"),
as.is=TRUE, strip.white=TRUE, header=TRUE)
head(input.history)
## TerritoryN tsh V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16
## 1 1 9 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 2 2 3 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 3 10 NA NA NA NA NA NA NA NA NA 0 2 0 NA NA NA NA
## 4 4 2 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 5 7 NA NA NA NA NA NA NA NA NA 0 NA NA NA NA NA NA
## 6 6 7 0 0 NA NA NA NA NA NA NA 0 0 NA NA NA NA NA
## V17 V18 V19 V20 V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34
## 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 2 NA NA 0 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 NA NA 0 NA NA NA NA NA NA NA NA 0 1 NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## V35 V36 V37 V38 V39 V40 V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52
## 1 NA NA NA NA NA NA NA NA NA NA NA 0 0 NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## V53 V54
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
hawk.history <- data.frame(ch =apply(input.history[,-(1:2)],1,paste,sep="", collapse=""),
freq =1,
tsh =input.history$tsh,
Territory = input.history$TerritoryN,
stringsAsFactors=FALSE)
head(hawk.history)
## ch
## 1 NANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANA00NANANANANANANA
## 2 NANANANANANANANANANANANANANANANANANA01NANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANA
## 3 NANANANANANANANANA020NANANANANANA0NANANANANANANANA01NANANANANANANANANANANANANANANANANANANANANANANANANA
## 4 1NANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANA
## 5 NANANANANANANANANA0NANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANA
## 6 00NANANANANANANA00NANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANANA
## freq tsh Territory
## 1 1 9 1
## 2 1 3 2
## 3 1 10 3
## 4 1 2 4
## 5 1 7 5
## 6 1 7 6
# fix up the NA in the history
hawk.history$ch <- gsub("NA",".", hawk.history$ch)
head(hawk.history)
## ch freq tsh
## 1 .............................................00....... 1 9
## 2 ..................01.................................. 1 3
## 3 .........020......0........01......................... 1 10
## 4 1..................................................... 1 2
## 5 .........0............................................ 1 7
## 6 00.......00........................................... 1 7
## Territory
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
# make a plot of the imputed occupancy
plotdata <- reshape2::melt(input.history,
id.vars=c("TerritoryN","tsh"),
value.name="ObsOcc",
variable.name="Visit")
plotdata$Visit <- as.numeric(substring(plotdata$Visit,2))
head(plotdata)
## TerritoryN tsh Visit ObsOcc
## 1 1 9 1 NA
## 2 2 3 1 NA
## 3 3 10 1 NA
## 4 4 2 1 1
## 5 5 7 1 NA
## 6 6 7 1 0
ggplot(data=plotdata, aes(x=Visit, y=TerritoryN, color=as.factor(ObsOcc), shape=as.factor(ObsOcc)))+
ggtitle("Observed occupancy state", subtitle="Not corrected for false negatives")+
geom_point()+
scale_shape_manual(name="Observed\nOccupancy", values=c(1, 16, 19))+
scale_color_manual(name="Observed\nOccupancy", values=c("black","green","red"))+
theme_bw()
## Warning: Removed 7275 rows containing missing values (geom_point).
hawk.data <- process.data(data=hawk.history,
model="RDMSOccRepro",
time.intervals=c( rep( c(rep(0,8),1),5),rep(0,8))) # 6 years x 9 visit occasion/year
summary(hawk.data)
## Length Class Mode
## data 4 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 148 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 6 -none- numeric
## time.intervals 53 -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 2 -none- character
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
# Robust design, Multi State, Reproduction parameterization
setup.parameters("RDMSOccRepro", check=TRUE) # returns a vector of parameter names (case sensitive)
## [1] "Phi0" "Psi" "R" "p" "Delta"
model.list.csv <- textConnection("
Phi0, Psi, R, p, Delta
~stratum, ~stratum, ~stratum, ~stratum, ~1
~stratum*tsh, ~stratum, ~stratum, ~stratum, ~1
~stratum*tsh, ~stratum*tsh, ~stratum*tsh, ~stratum, ~1
~stratum*tsh, ~stratum, ~1, ~stratum, ~1
~stratum, ~stratum*tsh, ~stratum, ~stratum, ~1
~stratum*tsh, ~stratum*tsh, ~stratum, ~stratum, ~1")
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
## Phi0 Psi R p Delta model.number
## 1 ~stratum ~stratum ~stratum ~stratum ~1 1
## 2 ~stratum*tsh ~stratum ~stratum ~stratum ~1 2
## 3 ~stratum*tsh ~stratum*tsh ~stratum*tsh ~stratum ~1 3
## 4 ~stratum*tsh ~stratum ~1 ~stratum ~1 4
## 5 ~stratum ~stratum*tsh ~stratum ~stratum ~1 5
## 6 ~stratum*tsh ~stratum*tsh ~stratum ~stratum ~1 6
# 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,hawk.data){
cat("\n\n***** Starting ", unlist(x), "\n")
#browser()
#fit <- myoccMod(model=list(as.formula(paste("psi",x$psi)),
fit <- RMark::mark(hawk.data,
model="RDMSOccRepro",
model.parameters=list(
Phi0 =list(formula=as.formula(eval(x$Phi0))),
Psi =list(formula=as.formula(eval(x$Psi))),
R =list(formula=as.formula(eval(x$R))),
p =list(formula=as.formula(eval(x$p))),
Delta =list(formula=as.formula(eval(x$Delta))))
,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
},hawk.data=hawk.data)
##
##
## ***** Starting ~stratum ~stratum ~stratum ~stratum ~1 1
##
## Note: only 10 parameters counted of 11 specified parameters
## AICc and parameter count have been adjusted upward
##
##
## ***** Starting ~stratum*tsh ~stratum ~stratum ~stratum ~1 2
##
## Note: only 12 parameters counted of 13 specified parameters
##
## AICc and parameter count have been adjusted upward
##
##
## ***** Starting ~stratum*tsh ~stratum*tsh ~stratum*tsh ~stratum ~1 3
##
## Note: only 14 parameters counted of 19 specified parameters
##
## AICc and parameter count have been adjusted upward
##
##
## ***** Starting ~stratum*tsh ~stratum ~1 ~stratum ~1 4
##
##
## ***** Starting ~stratum ~stratum*tsh ~stratum ~stratum ~1 5
##
## Note: only 13 parameters counted of 14 specified parameters
##
## AICc and parameter count have been adjusted upward
##
##
## ***** Starting ~stratum*tsh ~stratum*tsh ~stratum ~stratum ~1 6
##
## Note: only 15 parameters counted of 16 specified parameters
##
## AICc and parameter count have been adjusted upward
model.number <-2
summary(model.fits[[model.number]])
## Output summary for RDMSOccRepro model
## Name : Phi0(~stratum * tsh)Psi(~stratum)R(~stratum)p(~stratum)Delta(~1)
##
## Npar : 13 (unadjusted=12)
## -2lnL: 1051.256
## AICc : 1078.336 (unadjusted=1076.1794)
##
## Beta
## estimate se lcl ucl
## Phi0:(Intercept) -0.3587343 0.7833774 -1.8941540 1.1766853
## Phi0:stratum2 -0.5475147 1.1934986 -2.8867720 1.7917426
## Phi0:tsh 0.1699270 0.1126001 -0.0507691 0.3906231
## Phi0:stratum2:tsh -0.1427356 0.1272311 -0.3921087 0.1066374
## Psi:(Intercept) -1.8552009 1.2655619 -4.3357023 0.6253004
## Psi:stratum1 5.0219247 3.4909944 -1.8204244 11.8642740
## Psi:stratum2 4.6902478 1.6887155 1.3803654 8.0001302
## R:(Intercept) -13.9643060 178.4505700 -363.7274300 335.7988200
## R:stratum1 13.1130640 178.4503000 -336.6495400 362.8756700
## R:stratum2 16.3052170 178.4529200 -333.4625100 366.0729500
## p:(Intercept) -1.0543689 0.2653280 -1.5744118 -0.5343260
## p:stratum2 0.9272896 0.3319148 0.2767365 1.5778427
## Delta:(Intercept) 2.0181457 0.3859670 1.2616503 2.7746410
##
##
## Real Parameter Phi0
##
##
## 0.8191298 0.352711
##
##
## Real Parameter Psi
## Stratum:0
## 2 3 4 5 6
## 0.1352634 0.1352634 0.1352634 0.1352634 0.1352634
##
## Stratum:1
## 2 3 4 5 6
## 0.9595627 0.9595627 0.9595627 0.9595627 0.9595627
##
## Stratum:2
## 2 3 4 5 6
## 0.9445406 0.9445406 0.9445406 0.9445406 0.9445406
##
##
## Real Parameter R
## Stratum:0 To:2
## 2 3 4 5 6
## 8.617446e-07 8.617446e-07 8.617446e-07 8.617446e-07 8.617446e-07
##
## Stratum:1 To:2
## 2 3 4 5 6
## 0.2991725 0.2991725 0.2991725 0.2991725 0.2991725
##
## Stratum:2 To:2
## 2 3 4 5 6
## 0.912209 0.912209 0.912209 0.912209 0.912209
##
##
## Real Parameter p
## Session:1 Stratum:1
## 1 2 3 4 5 6 7 8
## 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387
## 9
## 0.258387
##
## Session:2 Stratum:1
## 1 2 3 4 5 6 7 8
## 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387
## 9
## 0.258387
##
## Session:3 Stratum:1
## 1 2 3 4 5 6 7 8
## 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387
## 9
## 0.258387
##
## Session:4 Stratum:1
## 1 2 3 4 5 6 7 8
## 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387
## 9
## 0.258387
##
## Session:5 Stratum:1
## 1 2 3 4 5 6 7 8
## 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387
## 9
## 0.258387
##
## Session:6 Stratum:1
## 1 2 3 4 5 6 7 8
## 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387 0.258387
## 9
## 0.258387
##
## Session:1 Stratum:2
## 1 2 3 4 5 6 7
## 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729
## 8 9
## 0.4682729 0.4682729
##
## Session:2 Stratum:2
## 1 2 3 4 5 6 7
## 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729
## 8 9
## 0.4682729 0.4682729
##
## Session:3 Stratum:2
## 1 2 3 4 5 6 7
## 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729
## 8 9
## 0.4682729 0.4682729
##
## Session:4 Stratum:2
## 1 2 3 4 5 6 7
## 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729
## 8 9
## 0.4682729 0.4682729
##
## Session:5 Stratum:2
## 1 2 3 4 5 6 7
## 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729
## 8 9
## 0.4682729 0.4682729
##
## Session:6 Stratum:2
## 1 2 3 4 5 6 7
## 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729 0.4682729
## 8 9
## 0.4682729 0.4682729
##
##
## Real Parameter Delta
## Session:1 Stratum:2 To:2
## 2 3 4 5 6 7 8
## 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891
## 9 10
## 0.8826891 0.8826891
##
## Session:2 Stratum:2 To:2
## 2 3 4 5 6 7 8
## 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891
## 9 10
## 0.8826891 0.8826891
##
## Session:3 Stratum:2 To:2
## 2 3 4 5 6 7 8
## 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891
## 9 10
## 0.8826891 0.8826891
##
## Session:4 Stratum:2 To:2
## 2 3 4 5 6 7 8
## 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891
## 9 10
## 0.8826891 0.8826891
##
## Session:5 Stratum:2 To:2
## 2 3 4 5 6 7 8
## 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891
## 9 10
## 0.8826891 0.8826891
##
## Session:6 Stratum:2 To:2
## 2 3 4 5 6 7 8
## 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891 0.8826891
## 9 10
## 0.8826891 0.8826891
model.fits[[model.number]]$results$real
## estimate se lcl ucl
## Phi0 s1 g1 8.191298e-01 0.1378001000 4.224964e-01 0.9655590
## Phi0 s2 g1 3.527110e-01 0.1048636000 1.813252e-01 0.5727565
## Psi s0 g1 a1 t2 1.352634e-01 0.1480292000 1.292350e-02 0.6514231
## Psi s1 g1 a1 t2 9.595627e-01 0.0942040000 1.691351e-01 0.9996386
## Psi s2 g1 a1 t2 9.445406e-01 0.0430566000 7.727724e-01 0.9884112
## R s0 to2 g1 a1 t2 8.617446e-07 0.0001537787 1.084387e-158 1.0000000
## R s1 to2 g1 a1 t2 2.991725e-01 0.1010079000 1.424030e-01 0.5232308
## R s2 to2 g1 a1 t2 9.122090e-01 0.0581681000 7.144877e-01 0.9773468
## p s1 g1 s1 t1 2.583870e-01 0.0508430000 1.715884e-01 0.3695085
## p s2 g1 s1 t1 4.682729e-01 0.0442617000 3.833167e-01 0.5551091
## Delta s2 to2 g1 a0 s1 t2 8.826891e-01 0.0399665000 7.793101e-01 0.9412900
## fixed note
## Phi0 s1 g1
## Phi0 s2 g1
## Psi s0 g1 a1 t2
## Psi s1 g1 a1 t2
## Psi s2 g1 a1 t2
## R s0 to2 g1 a1 t2
## R s1 to2 g1 a1 t2
## R s2 to2 g1 a1 t2
## p s1 g1 s1 t1
## p s2 g1 s1 t1
## Delta s2 to2 g1 a0 s1 t2
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## Phi0:(Intercept) -0.3587343 0.7833774 -1.8941540 1.1766853
## Phi0:stratum2 -0.5475147 1.1934986 -2.8867720 1.7917426
## Phi0:tsh 0.1699270 0.1126001 -0.0507691 0.3906231
## Phi0:stratum2:tsh -0.1427356 0.1272311 -0.3921087 0.1066374
## Psi:(Intercept) -1.8552009 1.2655619 -4.3357023 0.6253004
## Psi:stratum1 5.0219247 3.4909944 -1.8204244 11.8642740
## Psi:stratum2 4.6902478 1.6887155 1.3803654 8.0001302
## R:(Intercept) -13.9643060 178.4505700 -363.7274300 335.7988200
## R:stratum1 13.1130640 178.4503000 -336.6495400 362.8756700
## R:stratum2 16.3052170 178.4529200 -333.4625100 366.0729500
## p:(Intercept) -1.0543689 0.2653280 -1.5744118 -0.5343260
## p:stratum2 0.9272896 0.3319148 0.2767365 1.5778427
## Delta:(Intercept) 2.0181457 0.3859670 1.2616503 2.7746410
model.fits[[model.number]]$results$derived
## $`Phi_t transition matrix`
## estimate se lcl ucl
## 1 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 2 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 3 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 4 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 5 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 6 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 7 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 8 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 9 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 10 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 11 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 12 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 13 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 14 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 15 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 16 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 17 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 18 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 19 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 20 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 21 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 22 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 23 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 24 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 25 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 26 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 27 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 28 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 29 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 30 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 31 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 32 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 33 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 34 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 35 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 36 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 37 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 38 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 39 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 40 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 41 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 42 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 43 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 44 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 45 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
##
## $`1 - psi`
## estimate se lcl ucl
## 1 0.1808702 0.13780011 0.03444101 0.5775036
## 2 0.1938687 0.12206621 0.04945676 0.5264262
## 3 0.2062692 0.10909243 0.06576981 0.4896109
## 4 0.2174035 0.09823965 0.08221840 0.4627831
## 5 0.2270524 0.09079958 0.09629465 0.4474522
## 6 0.2352285 0.08755032 0.10594033 0.4439508
##
## $psi
## estimate se lcl ucl
## 1 0.8191298 0.13780011 0.4224964 0.9655590
## 2 0.8061313 0.12206621 0.4735738 0.9505432
## 3 0.7937308 0.10909243 0.5103891 0.9342302
## 4 0.7825965 0.09823965 0.5372169 0.9177816
## 5 0.7729476 0.09079958 0.5525478 0.9037053
## 6 0.7647715 0.08755032 0.5560492 0.8940597
##
## $psi_1
## estimate se lcl ucl
## 1 0.5302137 0.14253911 0.2688235 0.7760170
## 2 0.4049850 0.11184065 0.2151061 0.6283031
## 3 0.3318347 0.09504404 0.1765230 0.5350144
## 4 0.2893569 0.08237797 0.1566046 0.4717036
## 5 0.2648962 0.07435934 0.1456504 0.4323632
## 6 0.2509801 0.07050788 0.1384150 0.4113791
##
## $psi_2
## estimate se lcl ucl
## 1 0.2889161 0.08381683 0.1544268 0.4747680
## 2 0.4011464 0.06564123 0.2816580 0.5336666
## 3 0.4618961 0.06501905 0.3395189 0.5890441
## 4 0.4932396 0.06758508 0.3642402 0.6231454
## 5 0.5080514 0.07142366 0.3710072 0.6438961
## 6 0.5137915 0.07583258 0.3682341 0.6570465
get.real(model.fits[[model.number]], "Phi0", se=TRUE)
## all.diff.index par.index estimate se lcl
## Phi0 s1 g1 1 1 0.8191298 0.1378001 0.4224964
## Phi0 s2 g1 2 2 0.3527110 0.1048636 0.1813252
## ucl fixed note stratum group
## Phi0 s1 g1 0.9655590 1 1
## Phi0 s2 g1 0.5727565 2 1
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi s0 g1 a1 t2 3 3 0.1352634 0.1480292 0.0129235
## Psi s0 g1 a2 t3 4 3 0.1352634 0.1480292 0.0129235
## Psi s0 g1 a3 t4 5 3 0.1352634 0.1480292 0.0129235
## Psi s0 g1 a4 t5 6 3 0.1352634 0.1480292 0.0129235
## Psi s0 g1 a5 t6 7 3 0.1352634 0.1480292 0.0129235
## Psi s1 g1 a1 t2 8 4 0.9595627 0.0942040 0.1691351
## Psi s1 g1 a2 t3 9 4 0.9595627 0.0942040 0.1691351
## Psi s1 g1 a3 t4 10 4 0.9595627 0.0942040 0.1691351
## Psi s1 g1 a4 t5 11 4 0.9595627 0.0942040 0.1691351
## Psi s1 g1 a5 t6 12 4 0.9595627 0.0942040 0.1691351
## Psi s2 g1 a1 t2 13 5 0.9445406 0.0430566 0.7727724
## Psi s2 g1 a2 t3 14 5 0.9445406 0.0430566 0.7727724
## Psi s2 g1 a3 t4 15 5 0.9445406 0.0430566 0.7727724
## Psi s2 g1 a4 t5 16 5 0.9445406 0.0430566 0.7727724
## Psi s2 g1 a5 t6 17 5 0.9445406 0.0430566 0.7727724
## ucl fixed note group age time stratum Age Time 0
## Psi s0 g1 a1 t2 0.6514231 1 1 2 0 1 0 1
## Psi s0 g1 a2 t3 0.6514231 1 2 3 0 2 1 1
## Psi s0 g1 a3 t4 0.6514231 1 3 4 0 3 2 1
## Psi s0 g1 a4 t5 0.6514231 1 4 5 0 4 3 1
## Psi s0 g1 a5 t6 0.6514231 1 5 6 0 5 4 1
## Psi s1 g1 a1 t2 0.9996386 1 1 2 1 1 0 0
## Psi s1 g1 a2 t3 0.9996386 1 2 3 1 2 1 0
## Psi s1 g1 a3 t4 0.9996386 1 3 4 1 3 2 0
## Psi s1 g1 a4 t5 0.9996386 1 4 5 1 4 3 0
## Psi s1 g1 a5 t6 0.9996386 1 5 6 1 5 4 0
## Psi s2 g1 a1 t2 0.9884112 1 1 2 2 1 0 0
## Psi s2 g1 a2 t3 0.9884112 1 2 3 2 2 1 0
## Psi s2 g1 a3 t4 0.9884112 1 3 4 2 3 2 0
## Psi s2 g1 a4 t5 0.9884112 1 4 5 2 4 3 0
## Psi s2 g1 a5 t6 0.9884112 1 5 6 2 5 4 0
## 1 2
## Psi s0 g1 a1 t2 0 0
## Psi s0 g1 a2 t3 0 0
## Psi s0 g1 a3 t4 0 0
## Psi s0 g1 a4 t5 0 0
## Psi s0 g1 a5 t6 0 0
## Psi s1 g1 a1 t2 1 0
## Psi s1 g1 a2 t3 1 0
## Psi s1 g1 a3 t4 1 0
## Psi s1 g1 a4 t5 1 0
## Psi s1 g1 a5 t6 1 0
## Psi s2 g1 a1 t2 0 1
## Psi s2 g1 a2 t3 0 1
## Psi s2 g1 a3 t4 0 1
## Psi s2 g1 a4 t5 0 1
## Psi s2 g1 a5 t6 0 1
get.real(model.fits[[model.number]], "R", se=TRUE)
## all.diff.index par.index estimate se
## R s0 to2 g1 a1 t2 18 6 8.617446e-07 0.0001537787
## R s0 to2 g1 a2 t3 19 6 8.617446e-07 0.0001537787
## R s0 to2 g1 a3 t4 20 6 8.617446e-07 0.0001537787
## R s0 to2 g1 a4 t5 21 6 8.617446e-07 0.0001537787
## R s0 to2 g1 a5 t6 22 6 8.617446e-07 0.0001537787
## R s1 to2 g1 a1 t2 23 7 2.991725e-01 0.1010079000
## R s1 to2 g1 a2 t3 24 7 2.991725e-01 0.1010079000
## R s1 to2 g1 a3 t4 25 7 2.991725e-01 0.1010079000
## R s1 to2 g1 a4 t5 26 7 2.991725e-01 0.1010079000
## R s1 to2 g1 a5 t6 27 7 2.991725e-01 0.1010079000
## R s2 to2 g1 a1 t2 28 8 9.122090e-01 0.0581681000
## R s2 to2 g1 a2 t3 29 8 9.122090e-01 0.0581681000
## R s2 to2 g1 a3 t4 30 8 9.122090e-01 0.0581681000
## R s2 to2 g1 a4 t5 31 8 9.122090e-01 0.0581681000
## R s2 to2 g1 a5 t6 32 8 9.122090e-01 0.0581681000
## lcl ucl fixed note group age time
## R s0 to2 g1 a1 t2 1.084387e-158 1.0000000 1 1 2
## R s0 to2 g1 a2 t3 1.084387e-158 1.0000000 1 2 3
## R s0 to2 g1 a3 t4 1.084387e-158 1.0000000 1 3 4
## R s0 to2 g1 a4 t5 1.084387e-158 1.0000000 1 4 5
## R s0 to2 g1 a5 t6 1.084387e-158 1.0000000 1 5 6
## R s1 to2 g1 a1 t2 1.424030e-01 0.5232308 1 1 2
## R s1 to2 g1 a2 t3 1.424030e-01 0.5232308 1 2 3
## R s1 to2 g1 a3 t4 1.424030e-01 0.5232308 1 3 4
## R s1 to2 g1 a4 t5 1.424030e-01 0.5232308 1 4 5
## R s1 to2 g1 a5 t6 1.424030e-01 0.5232308 1 5 6
## R s2 to2 g1 a1 t2 7.144877e-01 0.9773468 1 1 2
## R s2 to2 g1 a2 t3 7.144877e-01 0.9773468 1 2 3
## R s2 to2 g1 a3 t4 7.144877e-01 0.9773468 1 3 4
## R s2 to2 g1 a4 t5 7.144877e-01 0.9773468 1 4 5
## R s2 to2 g1 a5 t6 7.144877e-01 0.9773468 1 5 6
## stratum tostratum Age Time 0 1 2 to0 to1 to2
## R s0 to2 g1 a1 t2 0 2 1 0 1 0 0 0 0 1
## R s0 to2 g1 a2 t3 0 2 2 1 1 0 0 0 0 1
## R s0 to2 g1 a3 t4 0 2 3 2 1 0 0 0 0 1
## R s0 to2 g1 a4 t5 0 2 4 3 1 0 0 0 0 1
## R s0 to2 g1 a5 t6 0 2 5 4 1 0 0 0 0 1
## R s1 to2 g1 a1 t2 1 2 1 0 0 1 0 0 0 1
## R s1 to2 g1 a2 t3 1 2 2 1 0 1 0 0 0 1
## R s1 to2 g1 a3 t4 1 2 3 2 0 1 0 0 0 1
## R s1 to2 g1 a4 t5 1 2 4 3 0 1 0 0 0 1
## R s1 to2 g1 a5 t6 1 2 5 4 0 1 0 0 0 1
## R s2 to2 g1 a1 t2 2 2 1 0 0 0 1 0 0 1
## R s2 to2 g1 a2 t3 2 2 2 1 0 0 1 0 0 1
## R s2 to2 g1 a3 t4 2 2 3 2 0 0 1 0 0 1
## R s2 to2 g1 a4 t5 2 2 4 3 0 0 1 0 0 1
## R s2 to2 g1 a5 t6 2 2 5 4 0 0 1 0 0 1
temp <- get.real(model.fits[[model.number]], "p", se=TRUE)
temp[ temp$time==1 & temp$session==1,]
## all.diff.index par.index estimate se lcl
## p s1 g1 s1 t1 33 9 0.2583870 0.0508430 0.1715884
## p s2 g1 s1 t1 87 10 0.4682729 0.0442617 0.3833167
## ucl fixed note group time stratum session Time 1 2
## p s1 g1 s1 t1 0.3695085 1 1 1 1 0 1 0
## p s2 g1 s1 t1 0.5551091 1 1 2 1 0 0 1
# collect models and make AICc table
model.set <- RMark::collect.models( type="RDMSOccRepro")
model.set
## model
## 2 Phi0(~stratum * tsh)Psi(~stratum)R(~stratum)p(~stratum)Delta(~1)
## 1 Phi0(~stratum)Psi(~stratum)R(~stratum)p(~stratum)Delta(~1)
## 5 Phi0(~stratum)Psi(~stratum * tsh)R(~stratum)p(~stratum)Delta(~1)
## 6 Phi0(~stratum * tsh)Psi(~stratum * tsh)R(~stratum)p(~stratum)Delta(~1)
## 3 Phi0(~stratum * tsh)Psi(~stratum * tsh)R(~stratum * tsh)p(~stratum)Delta(~1)
## 4 Phi0(~stratum * tsh)Psi(~stratum)R(~1)p(~stratum)Delta(~1)
## npar AICc DeltaAICc weight Deviance
## 2 13 1078.336 0.000000 0.4803850272 1051.256
## 1 11 1079.522 1.185242 0.2655934453 1056.743
## 5 14 1080.577 2.240881 0.1566708753 1051.327
## 6 16 1081.720 3.383224 0.0884976512 1048.091
## 3 19 1086.460 8.123954 0.0082698057 1046.164
## 4 11 1091.764 13.427642 0.0005831952 1068.985
# get model averaged values for the initial occupancy probabilities at different values of the covariates
get.real(model.fits[[model.number]], "Phi0", se=TRUE)
## all.diff.index par.index estimate se lcl
## Phi0 s1 g1 1 1 0.8191298 0.1378001 0.4224964
## Phi0 s2 g1 2 2 0.3527110 0.1048636 0.1813252
## ucl fixed note stratum group
## Phi0 s1 g1 0.9655590 1 1
## Phi0 s2 g1 0.5727565 2 1
Phi0.ma <- RMark::model.average(model.set, param="Phi0")
Phi0.ma
## par.index estimate se fixed note stratum group
## Phi0 s1 g1 1 0.7753864 0.1397825 1 1
## Phi0 s2 g1 2 0.3671753 0.1081650 2 1
range(hawk.history$tsh, na.rm=TRUE)
## [1] 1 36
Phi0.pred <- covariate.predictions(model.set, data=hawk.history, indices=1:4)
Phi0.pred$estimates[1:10,]
## vcv.index model.index par.index
## 1 1 1 1
## 2 2 2 2
## 3 3 3 3
## 4 3 3 4
## 5 4 1 1
## 6 5 2 2
## 7 6 3 3
## 8 6 3 4
## 9 7 1 1
## 10 8 2 2
## ch freq tsh
## 1 .............................................00....... 1 9
## 2 .............................................00....... 1 9
## 3 .............................................00....... 1 9
## 4 .............................................00....... 1 9
## 5 ..................01.................................. 1 3
## 6 ..................01.................................. 1 3
## 7 ..................01.................................. 1 3
## 8 ..................01.................................. 1 3
## 9 .........020......0........01......................... 1 10
## 10 .........020......0........01......................... 1 10
## Territory estimate se lcl ucl fixed
## 1 1 0.7440811 0.1331486 0.42477990 0.9196618
## 2 1 0.3588518 0.1117965 0.17759431 0.5919493
## 3 1 0.2072014 0.1791005 0.02992026 0.6889224
## 4 1 0.2072014 0.1791005 0.02992026 0.6889224
## 5 2 0.6194517 0.1662538 0.29005767 0.8664059
## 6 2 0.3350442 0.1345756 0.13361208 0.6221011
## 7 2 0.1926305 0.1749167 0.02564555 0.6838214
## 8 2 0.1926305 0.1749167 0.02564555 0.6838214
## 9 3 0.7604688 0.1363573 0.42263789 0.9322930
## 10 3 0.3629915 0.1096319 0.18366763 0.5907064
# Extract the estimated psi values and combine back with standardized values
# This is tricky because the region creates different group
plotdata1 <- Phi0.pred$estimates[ Phi0.pred$estimates$par.index==1, ]
plotdata1$Source <- "Initial occupancy (breeding or not breeding)"
plotdata2 <- Phi0.pred$estimates[ Phi0.pred$estimates$par.index==2, ]
plotdata2$Source <- "Conditional p(breeding | occupied)"
plotdata <- rbind(plotdata1, plotdata2)
plotdata.ci <- plyr::ddply(plotdata, c("Source"), function(x){
# select 5 points at random for ci
x <- x[ sample(nrow(x),5),]
x
})
head(plotdata)
## vcv.index model.index par.index
## 1 1 1 1
## 5 4 1 1
## 9 7 1 1
## 13 10 1 1
## 17 13 1 1
## 21 16 1 1
## ch freq tsh
## 1 .............................................00....... 1 9
## 5 ..................01.................................. 1 3
## 9 .........020......0........01......................... 1 10
## 13 1..................................................... 1 2
## 17 .........0............................................ 1 7
## 21 00.......00........................................... 1 7
## Territory estimate se lcl ucl fixed
## 1 1 0.7440811 0.1331486 0.4247799 0.9196618
## 5 2 0.6194517 0.1662538 0.2900577 0.8664059
## 9 3 0.7604688 0.1363573 0.4226379 0.9322930
## 13 4 0.5960933 0.1839560 0.2481854 0.8683834
## 17 5 0.7070213 0.1306724 0.4120556 0.8925825
## 21 6 0.7070213 0.1306724 0.4120556 0.8925825
## Source
## 1 Initial occupancy (breeding or not breeding)
## 5 Initial occupancy (breeding or not breeding)
## 9 Initial occupancy (breeding or not breeding)
## 13 Initial occupancy (breeding or not breeding)
## 17 Initial occupancy (breeding or not breeding)
## 21 Initial occupancy (breeding or not breeding)
i.occ.plot <- ggplot(data=plotdata, aes(x=tsh, y=estimate))+
ggtitle(paste('Model averaged initial occupancy probability from multi-state occupancy model',
'\nVariation in trend represents effects from other models',sep=""))+
geom_point()+
ylim(0,1)+
ylab("Initial occupancy probability\nSelected 95% ci")+
xlab("Area of suitable habitat ")+
theme(legend.position=c(1,0), legend.justification=c(1,0))+
facet_wrap(~Source, ncol=1)+
geom_errorbar(data=plotdata.ci, aes(ymin=lcl, ymax=ucl), width=0.1)
i.occ.plot
# Model averaged estimates of occupancy over time given previous state
Psi.ma <- RMark::model.average(model.set, param="Psi")
Psi.ma
## par.index estimate se fixed note group age time
## Psi s0 g1 a1 t2 3 0.2129985 0.1817714 1 1 2
## Psi s0 g1 a2 t3 4 0.2129985 0.1817714 1 2 3
## Psi s0 g1 a3 t4 5 0.2129985 0.1817714 1 3 4
## Psi s0 g1 a4 t5 6 0.2129985 0.1817714 1 4 5
## Psi s0 g1 a5 t6 7 0.2129985 0.1817714 1 5 6
## Psi s1 g1 a1 t2 8 0.9054559 0.1459711 1 1 2
## Psi s1 g1 a2 t3 9 0.9054559 0.1459711 1 2 3
## Psi s1 g1 a3 t4 10 0.9054559 0.1459711 1 3 4
## Psi s1 g1 a4 t5 11 0.9054559 0.1459711 1 4 5
## Psi s1 g1 a5 t6 12 0.9054559 0.1459711 1 5 6
## Psi s2 g1 a1 t2 13 0.9422494 0.0438346 1 1 2
## Psi s2 g1 a2 t3 14 0.9422494 0.0438346 1 2 3
## Psi s2 g1 a3 t4 15 0.9422494 0.0438346 1 3 4
## Psi s2 g1 a4 t5 16 0.9422494 0.0438346 1 4 5
## Psi s2 g1 a5 t6 17 0.9422494 0.0438346 1 5 6
## stratum Age Time 0 1 2
## Psi s0 g1 a1 t2 0 1 0 1 0 0
## Psi s0 g1 a2 t3 0 2 1 1 0 0
## Psi s0 g1 a3 t4 0 3 2 1 0 0
## Psi s0 g1 a4 t5 0 4 3 1 0 0
## Psi s0 g1 a5 t6 0 5 4 1 0 0
## Psi s1 g1 a1 t2 1 1 0 0 1 0
## Psi s1 g1 a2 t3 1 2 1 0 1 0
## Psi s1 g1 a3 t4 1 3 2 0 1 0
## Psi s1 g1 a4 t5 1 4 3 0 1 0
## Psi s1 g1 a5 t6 1 5 4 0 1 0
## Psi s2 g1 a1 t2 2 1 0 0 0 1
## Psi s2 g1 a2 t3 2 2 1 0 0 1
## Psi s2 g1 a3 t4 2 3 2 0 0 1
## Psi s2 g1 a4 t5 2 4 3 0 0 1
## Psi s2 g1 a5 t6 2 5 4 0 0 1
R.ma <- RMark::model.average(model.set, param="R")
R.ma
## par.index estimate se fixed note group
## R s0 to2 g1 a1 t2 18 0.0004069587 0.01692707 1
## R s0 to2 g1 a2 t3 19 0.0004069587 0.01692707 1
## R s0 to2 g1 a3 t4 20 0.0004069587 0.01692707 1
## R s0 to2 g1 a4 t5 21 0.0004069587 0.01692707 1
## R s0 to2 g1 a5 t6 22 0.0004069587 0.01692707 1
## R s1 to2 g1 a1 t2 23 0.3448114871 0.12678938 1
## R s1 to2 g1 a2 t3 24 0.3448114871 0.12678938 1
## R s1 to2 g1 a3 t4 25 0.3448114871 0.12678938 1
## R s1 to2 g1 a4 t5 26 0.3448114871 0.12678938 1
## R s1 to2 g1 a5 t6 27 0.3448114871 0.12678938 1
## R s2 to2 g1 a1 t2 28 0.9216555001 0.05655124 1
## R s2 to2 g1 a2 t3 29 0.9216555001 0.05655124 1
## R s2 to2 g1 a3 t4 30 0.9216555001 0.05655124 1
## R s2 to2 g1 a4 t5 31 0.9216555001 0.05655124 1
## R s2 to2 g1 a5 t6 32 0.9216555001 0.05655124 1
## age time stratum tostratum Age Time 0 1 2 to0 to1 to2
## R s0 to2 g1 a1 t2 1 2 0 2 1 0 1 0 0 0 0 1
## R s0 to2 g1 a2 t3 2 3 0 2 2 1 1 0 0 0 0 1
## R s0 to2 g1 a3 t4 3 4 0 2 3 2 1 0 0 0 0 1
## R s0 to2 g1 a4 t5 4 5 0 2 4 3 1 0 0 0 0 1
## R s0 to2 g1 a5 t6 5 6 0 2 5 4 1 0 0 0 0 1
## R s1 to2 g1 a1 t2 1 2 1 2 1 0 0 1 0 0 0 1
## R s1 to2 g1 a2 t3 2 3 1 2 2 1 0 1 0 0 0 1
## R s1 to2 g1 a3 t4 3 4 1 2 3 2 0 1 0 0 0 1
## R s1 to2 g1 a4 t5 4 5 1 2 4 3 0 1 0 0 0 1
## R s1 to2 g1 a5 t6 5 6 1 2 5 4 0 1 0 0 0 1
## R s2 to2 g1 a1 t2 1 2 2 2 1 0 0 0 1 0 0 1
## R s2 to2 g1 a2 t3 2 3 2 2 2 1 0 0 1 0 0 1
## R s2 to2 g1 a3 t4 3 4 2 2 3 2 0 0 1 0 0 1
## R s2 to2 g1 a4 t5 4 5 2 2 4 3 0 0 1 0 0 1
## R s2 to2 g1 a5 t6 5 6 2 2 5 4 0 0 1 0 0 1
# Get the transition matrix
model.number <- 2
names(model.fits[[model.number]]$results$derived)
## [1] "Phi_t transition matrix" "1 - psi"
## [3] "psi" "psi_1"
## [5] "psi_2"
model.fits[[model.number]]$results$real
## estimate se lcl ucl
## Phi0 s1 g1 8.191298e-01 0.1378001000 4.224964e-01 0.9655590
## Phi0 s2 g1 3.527110e-01 0.1048636000 1.813252e-01 0.5727565
## Psi s0 g1 a1 t2 1.352634e-01 0.1480292000 1.292350e-02 0.6514231
## Psi s1 g1 a1 t2 9.595627e-01 0.0942040000 1.691351e-01 0.9996386
## Psi s2 g1 a1 t2 9.445406e-01 0.0430566000 7.727724e-01 0.9884112
## R s0 to2 g1 a1 t2 8.617446e-07 0.0001537787 1.084387e-158 1.0000000
## R s1 to2 g1 a1 t2 2.991725e-01 0.1010079000 1.424030e-01 0.5232308
## R s2 to2 g1 a1 t2 9.122090e-01 0.0581681000 7.144877e-01 0.9773468
## p s1 g1 s1 t1 2.583870e-01 0.0508430000 1.715884e-01 0.3695085
## p s2 g1 s1 t1 4.682729e-01 0.0442617000 3.833167e-01 0.5551091
## Delta s2 to2 g1 a0 s1 t2 8.826891e-01 0.0399665000 7.793101e-01 0.9412900
## fixed note
## Phi0 s1 g1
## Phi0 s2 g1
## Psi s0 g1 a1 t2
## Psi s1 g1 a1 t2
## Psi s2 g1 a1 t2
## R s0 to2 g1 a1 t2
## R s1 to2 g1 a1 t2
## R s2 to2 g1 a1 t2
## p s1 g1 s1 t1
## p s2 g1 s1 t1
## Delta s2 to2 g1 a0 s1 t2
model.fits[[model.number]]$results$derived$'Phi_t transition matrix'
## estimate se lcl ucl
## 1 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 2 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 3 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 4 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 5 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 6 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 7 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 8 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 9 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 10 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 11 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 12 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 13 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 14 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 15 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 16 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 17 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 18 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 19 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 20 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 21 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 22 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 23 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 24 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 25 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 26 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 27 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 28 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 29 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 30 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 31 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 32 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 33 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 34 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 35 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 36 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
## 37 8.647366e-01 1.480293e-01 3.485768e-01 0.9870765
## 38 1.352633e-01 1.480291e-01 1.292346e-02 0.6514228
## 39 1.165625e-07 2.080159e-05 1.443614e-159 1.0000000
## 40 4.043735e-02 9.420433e-02 3.613749e-04 0.8308672
## 41 6.724879e-01 1.478236e-01 3.552520e-01 0.8844181
## 42 2.870747e-01 8.204277e-02 1.550755e-01 0.4690564
## 43 5.545943e-02 4.305668e-02 1.158873e-02 0.2272284
## 44 8.292213e-02 5.592352e-02 2.094567e-02 0.2764928
## 45 8.616184e-01 5.983614e-02 6.995550e-01 0.9433431
# Model average the derived parameter
source(file.path("..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
trans.ma <- RMark.model.average.derived(model.set, "Phi_t transition matrix")
## Loading required package: plyr
head(trans.ma)
## estimate se lcl ucl
## 1 7.870016e-01 0.181771423 0.430736111 1.143266998
## 2 2.129552e-01 0.181785910 -0.143338621 0.569249052
## 3 4.322983e-05 0.002666943 -0.005183882 0.005270342
## 4 9.454409e-02 0.145971330 -0.191554463 0.380642637
## 5 5.964856e-01 0.189775431 0.224532621 0.968438640
## 6 3.089703e-01 0.088309138 0.135887553 0.482053011
temp <- matrix(trans.ma[1:9,"estimate"], 3,3, byrow=TRUE)
temp
## [,1] [,2] [,3]
## [1,] 0.78700155 0.21295522 4.322983e-05
## [2,] 0.09454409 0.59648563 3.089703e-01
## [3,] 0.05775064 0.07381982 8.684295e-01
# Estimate the average level 1 occupancy over time
# This is a derived parameter and so we need to model average as before.
psi1.ma <- RMark.model.average.derived(model.set, parameter="psi")
psi1.ma$Year <- 1:6
psi1.ma$Source="Model"
psi1.ma$Level="Level 1 occupancy"
psi1.ma
## estimate se lcl ucl Year Source Level
## 1 0.7753864 0.13978254 0.5014176 1.0493551 1 Model Level 1 occupancy
## 2 0.7656593 0.11563622 0.5390164 0.9923021 2 Model Level 1 occupancy
## 3 0.7637376 0.10110703 0.5655714 0.9619037 3 Model Level 1 occupancy
## 4 0.7631890 0.09056797 0.5856790 0.9406990 4 Model Level 1 occupancy
## 5 0.7625875 0.08472481 0.5965299 0.9286451 5 Model Level 1 occupancy
## 6 0.7617404 0.08312087 0.5988265 0.9246543 6 Model Level 1 occupancy
psi2.ma <- RMark.model.average.derived(model.set, parameter="psi_2")
psi2.ma$Year <- 1:6
psi2.ma$Source="Model"
psi2.ma$Level="Level 2 occupancy (Breeding)"
psi2.ma
## estimate se lcl ucl Year Source
## 1 0.2839436 0.08057250 0.1260244 0.4418628 1 Model
## 2 0.3975798 0.06570988 0.2687908 0.5263688 2 Model
## 3 0.4582409 0.06602333 0.3288375 0.5876442 3 Model
## 4 0.4918840 0.07011454 0.3544620 0.6293060 4 Model
## 5 0.5107890 0.07598203 0.3618670 0.6597110 5 Model
## 6 0.5213612 0.08255414 0.3595580 0.6831643 6 Model
## Level
## 1 Level 2 occupancy (Breeding)
## 2 Level 2 occupancy (Breeding)
## 3 Level 2 occupancy (Breeding)
## 4 Level 2 occupancy (Breeding)
## 5 Level 2 occupancy (Breeding)
## 6 Level 2 occupancy (Breeding)
plotdata<- plyr::rbind.fill(psi1.ma, psi2.ma )
plotdata
## estimate se lcl ucl Year Source
## 1 0.7753864 0.13978254 0.5014176 1.0493551 1 Model
## 2 0.7656593 0.11563622 0.5390164 0.9923021 2 Model
## 3 0.7637376 0.10110703 0.5655714 0.9619037 3 Model
## 4 0.7631890 0.09056797 0.5856790 0.9406990 4 Model
## 5 0.7625875 0.08472481 0.5965299 0.9286451 5 Model
## 6 0.7617404 0.08312087 0.5988265 0.9246543 6 Model
## 7 0.2839436 0.08057250 0.1260244 0.4418628 1 Model
## 8 0.3975798 0.06570988 0.2687908 0.5263688 2 Model
## 9 0.4582409 0.06602333 0.3288375 0.5876442 3 Model
## 10 0.4918840 0.07011454 0.3544620 0.6293060 4 Model
## 11 0.5107890 0.07598203 0.3618670 0.6597110 5 Model
## 12 0.5213612 0.08255414 0.3595580 0.6831643 6 Model
## Level
## 1 Level 1 occupancy
## 2 Level 1 occupancy
## 3 Level 1 occupancy
## 4 Level 1 occupancy
## 5 Level 1 occupancy
## 6 Level 1 occupancy
## 7 Level 2 occupancy (Breeding)
## 8 Level 2 occupancy (Breeding)
## 9 Level 2 occupancy (Breeding)
## 10 Level 2 occupancy (Breeding)
## 11 Level 2 occupancy (Breeding)
## 12 Level 2 occupancy (Breeding)
plot.occ1 <- ggplot(data=plotdata, aes(x=Year, y=estimate, linetype=Source))+
ggtitle("Model averaged level 1 and 2 occupancy from multi-state model", subtitle="At a global average tsh")+
geom_line()+
ylim(0,1)+ylab("Average level 1 occupancy over all territories \n(95% ci)")+
theme(legend.position=c(0,.5), legend.justification=c(0,1), legend.box="horizontal")+
scale_linetype_discrete(name="Estimate")+
geom_errorbar(aes(ymin=pmax(0,lcl), ymax=pmin(1,ucl),width=.1))+
facet_wrap(~Level, ncol=1)
plot.occ1
remove.mark(collect.models(), 1:nrow(collect.models()$model.table))
## NULL
cleanup(ask=FALSE)