# Single Species, Multi Season Occupancy analyais
# Salamanders 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 <- readxl::read_excel(file.path("..","SalamanderMultiSeason.xls"),
sheet="Detections",
skip=2, na="-",
col_names=FALSE)
head(input.data)
## # A tibble: 6 x 22
## X__1 X__2 X__3 X__4 X__5 X__6 X__7 X__8 X__9 X__10 X__11 X__12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00 1.00 0 1.00 0 1.00 1.00 0 1.00 1.00 1.00 0
## 2 2.00 0 0 0 0 0 0 0 0 0 0 0
## 3 3.00 0 1.00 0 1.00 1.00 1.00 0 0 0 1.00 1.00
## 4 4.00 1.00 1.00 1.00 1.00 0 0 0 1.00 1.00 1.00 1.00
## 5 5.00 1.00 1.00 1.00 1.00 0 1.00 0 1.00 0 0 0
## 6 6.00 0 0 0 0 0 0 0 0 1.00 0 0
## # ... with 10 more variables: X__13 <dbl>, X__14 <dbl>, X__15 <dbl>, X__16
## # <dbl>, X__17 <dbl>, X__18 <dbl>, X__19 <dbl>, X__20 <lgl>, X__21
## # <dbl>, X__22 <dbl>
input.history <- input.data[, 10:19]
head(input.history)
## # A tibble: 6 x 10
## X__10 X__11 X__12 X__13 X__14 X__15 X__16 X__17 X__18 X__19
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00 1.00 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 1.00 1.00 0 0 0 0 1.00 0 0
## 4 1.00 1.00 1.00 1.00 0 1.00 1.00 1.00 1.00 0
## 5 0 0 0 0 1.00 0 0 0 0 0
## 6 1.00 0 0 0 0 0 0 1.00 0 0
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 39
ncol(input.history)
## [1] 10
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 0
site.covar <- data.frame(Site=1:nrow(input.data),
Elevation=unlist(input.data[,21]),
Prox=car::recode(unlist(input.data[,22]),
"1='Yes'; 0='No';"))
row.names(site.covar) <- NULL
head(site.covar)
## Site Elevation Prox
## 1 1 841.248 Yes
## 2 2 853.440 No
## 3 3 780.288 Yes
## 4 4 707.136 Yes
## 5 5 670.560 No
## 6 6 670.560 Yes
#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 ch
## 1 1 1100000000
## 2 1 0000000000
## 3 1 0110000100
## 4 1 1111011110
## 5 1 0000100000
## 6 1 1000000100
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## freq ch Site Elevation Prox
## 1 1 1100000000 1 841.248 Yes
## 2 1 0000000000 2 853.440 No
## 3 1 0110000100 3 780.288 Yes
## 4 1 1111011110 4 707.136 Yes
## 5 1 0000100000 5 670.560 No
## 6 1 1000000100 6 670.560 Yes
#Create the RMark data structure
#Two Seasons, with 5 visits per season
max.visit.per.year <- 5
n.season <- 2
salamander.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)),
group="Prox")
summary(salamander.data)
## Length Class Mode
## data 6 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 2 data.frame list
## nocc 1 -none- numeric
## nocc.secondary 2 -none- numeric
## time.intervals 9 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 2 -none- numeric
## group.covariates 1 data.frame list
## nstrata 1 -none- numeric
## strata.labels 0 -none- NULL
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
# add visit level covariates in the loop below because of different data types
# 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.
model.list.csv <- textConnection("
p, Psi, Gamma, Epsilon, model.type
~1, ~1, ~1, ~1, RDOccupEG
~1, ~Elevation, ~1, ~1, RDOccupEG
~1, ~Prox, ~Prox, ~1, 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 ~1 ~Elevation ~1 ~1 RDOccupEG 2
## 3 ~1 ~Prox ~Prox ~1 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")
max.visit.per.year <- 5
n.season <- 2
# 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)),
group = "Prox")
# 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: 316.773
## AICc : 325.321
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.1254141 0.4123724 -0.6828358 0.9336639
## Epsilon:(Intercept) -1.8185760 1.2131228 -4.1962968 0.5591448
## Gamma:(Intercept) -1.2370966 0.8684571 -2.9392725 0.4650793
## p:(Intercept) -0.8921720 0.2030460 -1.2901421 -0.4942018
##
##
## Real Parameter Psi
## Group:ProxNo
## 1
## 0.5313125
##
## Group:ProxYes
## 1
## 0.5313125
##
##
## Real Parameter Epsilon
## Group:ProxNo
## 1
## 0.1396048
##
## Group:ProxYes
## 1
## 0.1396048
##
##
## Real Parameter Gamma
## Group:ProxNo
## 1
## 0.2249418
##
## Group:ProxYes
## 1
## 0.2249418
##
##
## Real Parameter p
## Session:1Group:ProxNo
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
##
## Session:2Group:ProxNo
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
##
## Session:1Group:ProxYes
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
##
## Session:2Group:ProxYes
## 1 2 3 4 5
## 0.2906618 0.2906618 0.2906618 0.2906618 0.2906618
##
##
## ***** Starting ~1 ~Elevation ~1 ~1 RDOccupEG 2
##
## Output summary for RDOccupEG model
## Name : Psi(~Elevation)Epsilon(~1)Gamma(~1)p(~1)
##
## Npar : 5
## -2lnL: 315.6731
## AICc : 326.5065
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 3.2246213 3.0532662 -2.759781 9.2090232
## Psi:Elevation -0.0036358 0.0035267 -0.010548 0.0032765
## Epsilon:(Intercept) -1.8757640 1.2673522 -4.359774 0.6082464
## Gamma:(Intercept) -1.2994714 0.9262164 -3.114855 0.5159127
## p:(Intercept) -0.8989659 0.2047022 -1.300182 -0.4977495
##
##
## Real Parameter Psi
## Group:ProxNo
## 1
## 0.5364704
##
## Group:ProxYes
## 1
## 0.5364704
##
##
## Real Parameter Epsilon
## Group:ProxNo
## 1
## 0.1328762
##
## Group:ProxYes
## 1
## 0.1328762
##
##
## Real Parameter Gamma
## Group:ProxNo
## 1
## 0.214254
##
## Group:ProxYes
## 1
## 0.214254
##
##
## Real Parameter p
## Session:1Group:ProxNo
## 1 2 3 4 5
## 0.2892631 0.2892631 0.2892631 0.2892631 0.2892631
##
## Session:2Group:ProxNo
## 1 2 3 4 5
## 0.2892631 0.2892631 0.2892631 0.2892631 0.2892631
##
## Session:1Group:ProxYes
## 1 2 3 4 5
## 0.2892631 0.2892631 0.2892631 0.2892631 0.2892631
##
## Session:2Group:ProxYes
## 1 2 3 4 5
## 0.2892631 0.2892631 0.2892631 0.2892631 0.2892631
##
##
## ***** Starting ~1 ~Prox ~Prox ~1 RDOccupEG 3
##
## Note: only 5 parameters counted of 6 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for RDOccupEG model
## Name : Psi(~Prox)Epsilon(~1)Gamma(~Prox)p(~1)
##
## Npar : 6 (unadjusted=5)
## -2lnL: 296.6992
## AICc : 309.8823 (unadjusted=307.53254)
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -1.3866309 0.5747103 -2.513063 -0.2601987
## Psi:ProxYes 5.4785383 4.5030334 -3.347407 14.3044840
## Epsilon:(Intercept) -1.5971488 0.9920702 -3.541606 0.3473088
## Gamma:(Intercept) -0.9553679 0.6083314 -2.147698 0.2369617
## Gamma:ProxYes -16.4490720 0.0000000 -16.449072 -16.4490720
## p:(Intercept) -0.8408081 0.1845680 -1.202561 -0.4790549
##
##
## Real Parameter Psi
## Group:ProxNo
## 1
## 0.1999462
##
## Group:ProxYes
## 1
## 0.9835672
##
##
## Real Parameter Epsilon
## Group:ProxNo
## 1
## 0.1683805
##
## Group:ProxYes
## 1
## 0.1683805
##
##
## Real Parameter Gamma
## Group:ProxNo
## 1
## 0.2778066
##
## Group:ProxYes
## 1
## 2.76279e-08
##
##
## Real Parameter p
## Session:1Group:ProxNo
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
##
## Session:2Group:ProxNo
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
##
## Session:1Group:ProxYes
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
##
## Session:2Group:ProxYes
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
# examine individula model results
model.number <-3
summary(model.fits[[model.number]])
## Output summary for RDOccupEG model
## Name : Psi(~Prox)Epsilon(~1)Gamma(~Prox)p(~1)
##
## Npar : 6 (unadjusted=5)
## -2lnL: 296.6992
## AICc : 309.8823 (unadjusted=307.53254)
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -1.3866309 0.5747103 -2.513063 -0.2601987
## Psi:ProxYes 5.4785383 4.5030334 -3.347407 14.3044840
## Epsilon:(Intercept) -1.5971488 0.9920702 -3.541606 0.3473088
## Gamma:(Intercept) -0.9553679 0.6083314 -2.147698 0.2369617
## Gamma:ProxYes -16.4490720 0.0000000 -16.449072 -16.4490720
## p:(Intercept) -0.8408081 0.1845680 -1.202561 -0.4790549
##
##
## Real Parameter Psi
## Group:ProxNo
## 1
## 0.1999462
##
## Group:ProxYes
## 1
## 0.9835672
##
##
## Real Parameter Epsilon
## Group:ProxNo
## 1
## 0.1683805
##
## Group:ProxYes
## 1
## 0.1683805
##
##
## Real Parameter Gamma
## Group:ProxNo
## 1
## 0.2778066
##
## Group:ProxYes
## 1
## 2.76279e-08
##
##
## Real Parameter p
## Session:1Group:ProxNo
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
##
## Session:2Group:ProxNo
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
##
## Session:1Group:ProxYes
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
##
## Session:2Group:ProxYes
## 1 2 3 4 5
## 0.3013646 0.3013646 0.3013646 0.3013646 0.3013646
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed
## Psi gNo a0 t1 1.999462e-01 0.0919351 7.494750e-02 4.353149e-01
## Psi gYes a0 t1 9.835672e-01 0.0724389 9.080500e-03 9.999974e-01
## Epsilon gNo a0 t1 1.683805e-01 0.1389181 2.815130e-02 5.859648e-01
## Gamma gNo a0 t1 2.778066e-01 0.1220496 1.045466e-01 5.589648e-01
## Gamma gYes a0 t1 2.762790e-08 0.0000000 2.762790e-08 2.762790e-08
## p gNo s1 t1 3.013646e-01 0.0388597 2.310199e-01 3.824753e-01
## note
## Psi gNo a0 t1
## Psi gYes a0 t1
## Epsilon gNo a0 t1
## Gamma gNo a0 t1
## Gamma gYes a0 t1
## p gNo s1 t1
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## Psi:(Intercept) -1.3866309 0.5747103 -2.513063 -0.2601987
## Psi:ProxYes 5.4785383 4.5030334 -3.347407 14.3044840
## Epsilon:(Intercept) -1.5971488 0.9920702 -3.541606 0.3473088
## Gamma:(Intercept) -0.9553679 0.6083314 -2.147698 0.2369617
## Gamma:ProxYes -16.4490720 0.0000000 -16.449072 -16.4490720
## p:(Intercept) -0.8408081 0.1845680 -1.202561 -0.4790549
#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.fits[[model.number]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.1999462 0.09193508 0.01975339 0.3801389
## 2 0.3885393 0.10907887 0.17474475 0.6023339
## 3 0.9835672 0.07243916 0.84158646 1.1255480
## 4 0.8179537 0.13780894 0.54784815 1.0880592
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.9432198 0.8460224 0.2850159 3.601424
## 2 0.8316195 0.1389181 0.5593400 1.103899
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 2.54256831 1.4722414 -0.3430249 5.4281615
## 2 0.07506785 0.3269494 -0.5657531 0.7158888
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.1999462 0.09193508 0.01975339 0.3801389
## 2 0.3885393 0.10907887 0.17474475 0.6023339
## 3 0.9835672 0.07243916 0.84158646 1.1255480
## 4 0.8179537 0.13780894 0.54784815 1.0880592
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 1.9432198 0.8460224 0.2850159 3.601424
## 2 0.8316195 0.1389181 0.5593400 1.103899
model.fits[[model.number]]$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 2.54256831 1.4722414 -0.3430249 5.4281615
## 2 0.07506785 0.3269494 -0.5657531 0.7158888
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## Warning in get.real(model.fits[[model.number]], "Psi", se = TRUE):
## Improper V-C matrix for beta estimates. Some variances non-positive.
## all.diff.index par.index estimate se lcl
## Psi gNo a0 t1 1 1 0.1999462 0.0919351 0.0749475
## Psi gYes a0 t1 2 2 0.9835672 0.0724389 0.0090805
## ucl fixed note group age time Age Time Prox
## Psi gNo a0 t1 0.4353149 No 0 1 0 0 No
## Psi gYes a0 t1 0.9999974 Yes 0 1 0 0 Yes
get.real(model.fits[[model.number]], "p", se=TRUE)
## Warning in get.real(model.fits[[model.number]], "p", se = TRUE):
## Improper V-C matrix for beta estimates. Some variances non-positive.
## all.diff.index par.index estimate se lcl
## p gNo s1 t1 7 6 0.3013646 0.0388597 0.2310199
## p gNo s1 t2 8 6 0.3013646 0.0388597 0.2310199
## p gNo s1 t3 9 6 0.3013646 0.0388597 0.2310199
## p gNo s1 t4 10 6 0.3013646 0.0388597 0.2310199
## p gNo s1 t5 11 6 0.3013646 0.0388597 0.2310199
## p gNo s2 t1 12 6 0.3013646 0.0388597 0.2310199
## p gNo s2 t2 13 6 0.3013646 0.0388597 0.2310199
## p gNo s2 t3 14 6 0.3013646 0.0388597 0.2310199
## p gNo s2 t4 15 6 0.3013646 0.0388597 0.2310199
## p gNo s2 t5 16 6 0.3013646 0.0388597 0.2310199
## p gYes s1 t1 17 6 0.3013646 0.0388597 0.2310199
## p gYes s1 t2 18 6 0.3013646 0.0388597 0.2310199
## p gYes s1 t3 19 6 0.3013646 0.0388597 0.2310199
## p gYes s1 t4 20 6 0.3013646 0.0388597 0.2310199
## p gYes s1 t5 21 6 0.3013646 0.0388597 0.2310199
## p gYes s2 t1 22 6 0.3013646 0.0388597 0.2310199
## p gYes s2 t2 23 6 0.3013646 0.0388597 0.2310199
## p gYes s2 t3 24 6 0.3013646 0.0388597 0.2310199
## p gYes s2 t4 25 6 0.3013646 0.0388597 0.2310199
## p gYes s2 t5 26 6 0.3013646 0.0388597 0.2310199
## ucl fixed note group time session Time Prox
## p gNo s1 t1 0.3824753 No 1 1 0 No
## p gNo s1 t2 0.3824753 No 2 1 1 No
## p gNo s1 t3 0.3824753 No 3 1 2 No
## p gNo s1 t4 0.3824753 No 4 1 3 No
## p gNo s1 t5 0.3824753 No 5 1 4 No
## p gNo s2 t1 0.3824753 No 1 2 0 No
## p gNo s2 t2 0.3824753 No 2 2 1 No
## p gNo s2 t3 0.3824753 No 3 2 2 No
## p gNo s2 t4 0.3824753 No 4 2 3 No
## p gNo s2 t5 0.3824753 No 5 2 4 No
## p gYes s1 t1 0.3824753 Yes 1 1 0 Yes
## p gYes s1 t2 0.3824753 Yes 2 1 1 Yes
## p gYes s1 t3 0.3824753 Yes 3 1 2 Yes
## p gYes s1 t4 0.3824753 Yes 4 1 3 Yes
## p gYes s1 t5 0.3824753 Yes 5 1 4 Yes
## p gYes s2 t1 0.3824753 Yes 1 2 0 Yes
## p gYes s2 t2 0.3824753 Yes 2 2 1 Yes
## p gYes s2 t3 0.3824753 Yes 3 2 2 Yes
## p gYes s2 t4 0.3824753 Yes 4 2 3 Yes
## p gYes s2 t5 0.3824753 Yes 5 2 4 Yes
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
model.set
## model npar AICc DeltaAICc
## 3 Psi(~Prox)Epsilon(~1)Gamma(~Prox)p(~1) 6 309.8823 0.00000
## 1 Psi(~1)Epsilon(~1)Gamma(~1)p(~1) 4 325.3210 15.43866
## 2 Psi(~Elevation)Epsilon(~1)Gamma(~1)p(~1) 5 326.5065 16.62414
## weight Deviance
## 3 0.9993107818 152.4704
## 1 0.0004438527 172.5442
## 2 0.0002453654 315.6731
# model averaging in the usual way
RMark::model.average(model.set, "Psi")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note group age time
## Psi gNo a0 t1 1 0.2001758 0.09235896 No 0 1
## Psi gYes a0 t1 2 0.9832568 0.07342334 Yes 0 1
## Age Time Prox
## Psi gNo a0 t1 0 0 No
## Psi gYes a0 t1 0 0 Yes
RMark::model.average(model.set, "p")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note group time session
## p gNo s1 t1 7 0.3013569 0.038863 No 1 1
## p gNo s1 t2 8 0.3013569 0.038863 No 2 1
## p gNo s1 t3 9 0.3013569 0.038863 No 3 1
## p gNo s1 t4 10 0.3013569 0.038863 No 4 1
## p gNo s1 t5 11 0.3013569 0.038863 No 5 1
## p gNo s2 t1 12 0.3013569 0.038863 No 1 2
## p gNo s2 t2 13 0.3013569 0.038863 No 2 2
## p gNo s2 t3 14 0.3013569 0.038863 No 3 2
## p gNo s2 t4 15 0.3013569 0.038863 No 4 2
## p gNo s2 t5 16 0.3013569 0.038863 No 5 2
## p gYes s1 t1 17 0.3013569 0.038863 Yes 1 1
## p gYes s1 t2 18 0.3013569 0.038863 Yes 2 1
## p gYes s1 t3 19 0.3013569 0.038863 Yes 3 1
## p gYes s1 t4 20 0.3013569 0.038863 Yes 4 1
## p gYes s1 t5 21 0.3013569 0.038863 Yes 5 1
## p gYes s2 t1 22 0.3013569 0.038863 Yes 1 2
## p gYes s2 t2 23 0.3013569 0.038863 Yes 2 2
## p gYes s2 t3 24 0.3013569 0.038863 Yes 3 2
## p gYes s2 t4 25 0.3013569 0.038863 Yes 4 2
## p gYes s2 t5 26 0.3013569 0.038863 Yes 5 2
## Time Prox
## p gNo s1 t1 0 No
## p gNo s1 t2 1 No
## p gNo s1 t3 2 No
## p gNo s1 t4 3 No
## p gNo s1 t5 4 No
## p gNo s2 t1 0 No
## p gNo s2 t2 1 No
## p gNo s2 t3 2 No
## p gNo s2 t4 3 No
## p gNo s2 t5 4 No
## p gYes s1 t1 0 Yes
## p gYes s1 t2 1 Yes
## p gYes s1 t3 2 Yes
## p gYes s1 t4 3 Yes
## p gYes s1 t5 4 Yes
## p gYes s2 t1 0 Yes
## p gYes s2 t2 1 Yes
## p gYes s2 t3 2 Yes
## p gYes s2 t4 3 Yes
## p gYes s2 t5 4 Yes
# 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.5313125 0.1026888 0.3300425 0.7325825
## 2 0.5625661 0.1038644 0.3589919 0.7661403
## 3 0.5313125 0.1026888 0.3300425 0.7325825
## 4 0.5625661 0.1038644 0.3589919 0.7661403
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.058823 0.2007517 0.6653501 1.452297
## 2 1.058823 0.2007517 0.6653501 1.452297
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 1.134474 0.4746753 0.2041102 2.064837
## 2 1.134474 0.4746753 0.2041102 2.064837
##
## $all_psi
## estimate se lcl ucl
## 1 0.5313125 0.1026888 0.3300425 0.7325825
## 2 0.5625661 0.1038644 0.3589919 0.7661403
## 3 0.5313125 0.1026888 0.3300425 0.7325825
## 4 0.5625661 0.1038644 0.3589919 0.7661403
RMark.model.average.derived(model.set, "all_psi")
## Loading required package: plyr
## Loading required package: boot
## estimate se lcl ucl
## 1 0.2001758 0.09235897 0.074756046 0.4367003
## 2 0.3886598 0.10917199 0.205299223 0.6100681
## 3 0.9832568 0.07342358 0.009299462 0.9999973
## 4 0.8177781 0.13795070 0.422371342 0.9649662
RMark.model.average.derived(model.set, "lambda Rate of Change")
## estimate se lcl ucl
## 1 1.9426087 0.8460673 0.2843473 3.600870
## 2 0.8317745 0.1390949 0.5591535 1.104396
RMark.model.average.derived(model.set, "log odds lambda")
## estimate se lcl ucl
## 1 2.54159427 1.472254 -0.3439699 5.4271584
## 2 0.07579446 0.328241 -0.5675461 0.7191351
# model averaging of derived parameters such as the occupancy at each time step
psi.ma <- RMark.model.average.derived(model.set, "all_psi")
# Note that due to the categorical covariate, there are 4 estimates, but only 2 years
# worth of data. Need to account for this when making Year column
psi.ma$Year <- rep(1:(nrow(psi.ma)/2),2)
psi.ma$parameter <- 'psi'
psi.ma$Prox = c(rep("No",2),rep("Yes",2))
psi.ma
## estimate se lcl ucl Year parameter Prox
## 1 0.2001758 0.09235897 0.074756046 0.4367003 1 psi No
## 2 0.3886598 0.10917199 0.205299223 0.6100681 2 psi No
## 3 0.9832568 0.07342358 0.009299462 0.9999973 1 psi Yes
## 4 0.8177781 0.13795070 0.422371342 0.9649662 2 psi Yes
# likely more interested in colonization and extinction probabilities
epsilon.ma <- RMark::model.average(model.set, "Epsilon", vcv=TRUE)$estimates
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
epsilon.ma$Year <- rep(1:(nrow(epsilon.ma)/2),2)
epsilon.ma$parameter <- 'epsilon'
epsilon.ma
## par.index estimate se lcl ucl fixed
## Epsilon gNo a0 t1 3 0.168359 0.1389254 0.02813988 0.5859917
## Epsilon gYes a0 t1 4 0.168359 0.1389254 0.02813988 0.5859917
## note group age time Age Time Prox Year parameter
## Epsilon gNo a0 t1 No 0 1 0 0 No 1 epsilon
## Epsilon gYes a0 t1 Yes 0 1 0 0 Yes 1 epsilon
gamma.ma <- RMark::model.average(model.set, "Gamma", vcv=TRUE)$estimates
##
## Model 3dropped from the model averaging because one or more beta variances are not positive
gamma.ma$Year <- rep(1:(nrow(gamma.ma)/2),2)
gamma.ma$parameter <- 'gamma'
gamma.ma
## par.index estimate se lcl ucl fixed
## Gamma gNo a0 t1 5 0.2211369 0.1531188 0.0473591 0.6185439
## Gamma gYes a0 t1 6 0.2211369 0.1531188 0.0473591 0.6185439
## note group age time Age Time Prox Year parameter
## Gamma gNo a0 t1 No 0 1 0 0 No 1 gamma
## Gamma gYes a0 t1 Yes 0 1 0 0 Yes 1 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 and 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)+
facet_wrap(~Prox)
cleanup(ask=FALSE)