# Single Species, MultiSeason Occupancy analyais - multi-model averaging
# Grand Skinks
# Data has been collected on 352 tors over a 5 year (the ‘seasons’) period,^
# although not all tors (rock piles) were surveyed each year, with up to 3 surveys^
# of each tor per year.
# The 15 columns are in 5 blocks of 3.
# There is also a site-specific covariate Pasture indicating whether the surrounding^
# matrix is either predominately the modified habitat (farm pasture, Pasture =1) or
# "native" grassland (tussock, Pasture = 0).
# 2018-08-18 Code submitted by Carl James Schwarz (cschwarz.stat.sfu.ca@gmail.com)
# Using the RMark package
library(car)
library(readxl)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(ggplot2)
# Get the RMark additional functions
source(file.path("..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
# get the data read in
# Data for detections should be a data frame with rows corresponding to sites
# and columns to visits.
# The usual 1=detected; 0=not detected; NA=not visited is used.
input.data <- readxl::read_excel(file.path("..","GrandSkinks.xls"),
sheet="DetectionHistory",
skip=3, na="-",
col_names=FALSE)
head(input.data)
## # A tibble: 6 x 18
## 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 0 0 0 0 0 0 0 0 NA 0 0
## 2 2.00 0 0 0 0 0 NA 0 0 0 0 0
## 3 3.00 0 0 0 0 0 NA 0 0 0 0 0
## 4 4.00 0 NA NA 0 NA NA 0 NA NA 0 0
## 5 5.00 0 NA NA NA NA NA NA NA NA NA NA
## 6 6.00 0 NA NA NA NA NA NA NA NA NA NA
## # ... with 6 more variables: X__13 <dbl>, X__14 <dbl>, X__15 <dbl>, X__16
## # <dbl>, X__17 <lgl>, X__18 <dbl>
input.history <- input.data[, 2:16]
head(input.history)
## # A tibble: 6 x 15
## X__2 X__3 X__4 X__5 X__6 X__7 X__8 X__9 X__10 X__11 X__12 X__13
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 0 0 0 0 NA 0 0 0
## 2 0 0 0 0 0 NA 0 0 0 0 0 0
## 3 0 0 0 0 0 NA 0 0 0 0 0 0
## 4 0 NA NA 0 NA NA 0 NA NA 0 0 NA
## 5 0 NA NA NA NA NA NA NA NA NA NA NA
## 6 0 NA NA NA NA NA NA NA NA NA NA NA
## # ... with 3 more variables: X__14 <dbl>, X__15 <dbl>, X__16 <dbl>
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 352
ncol(input.history)
## [1] 15
range(input.history, na.rm=TRUE)
## [1] 0 1
#Looks like there are some missing values
sum(is.na(input.history))
## [1] 2969
site.covar <- data.frame(Site=1:nrow(input.data),
Pasture=car::recode(input.data$X__18,
"1='Modified'; 0='Native';"))
head(site.covar)
## Site Pasture
## 1 1 Modified
## 2 2 Modified
## 3 3 Modified
## 4 4 Modified
## 5 5 Modified
## 6 6 Modified
#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 00000000NA000000
## 2 1 00000NA000000000
## 3 1 00000NA000000000
## 4 1 0NANA0NANA0NANA00NA0NANA
## 5 1 0NANANANANANANANANANANA0NANA
## 6 1 0NANANANANANANANANANANA0NANA
# Change any NA to . in the chapter history
input.history$ch <- gsub("NA",".", input.history$ch, fixed=TRUE)
head(input.history)
## freq ch
## 1 1 00000000.000000
## 2 1 00000.000000000
## 3 1 00000.000000000
## 4 1 0..0..0..00.0..
## 5 1 0...........0..
## 6 1 0...........0..
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## freq ch Site Pasture
## 1 1 00000000.000000 1 Modified
## 2 1 00000.000000000 2 Modified
## 3 1 00000.000000000 3 Modified
## 4 1 0..0..0..00.0.. 4 Modified
## 5 1 0...........0.. 5 Modified
## 6 1 0...........0.. 6 Modified
#Create the RMark data structure
#Five Seasons, with 3 visits per season
max.visit.per.year <- 3
n.season <- 5
skink.data <- process.data(data=input.history,
model="RDOccupEG",
groups = "Pasture",
time.intervals=c( rep( c(rep(0,max.visit.per.year-1),1),n.season-1),
rep(0,max.visit.per.year-1)))
summary(skink.data)
## Length Class Mode
## data 5 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 2 data.frame list
## nocc 1 -none- numeric
## nocc.secondary 5 -none- numeric
## time.intervals 14 -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
# any time specific covariates need to be added to the ddl's in the loop below
# 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
# What are the parameter names for Multi Season Single Species models
setup.parameters("RDOccupEG", check=TRUE)
## [1] "Psi" "Epsilon" "Gamma" "p"
#there are other parameterizations available
setup.parameters("RDOccupPE", check=TRUE) # psi, epsilon, p
## [1] "Psi" "Epsilon" "p"
setup.parameters("RDOccupPG", check=TRUE) # psi, gamma, p
## [1] "Psi" "Gamma" "p"
# 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
~session, ~1, ~time, ~time, RDOccupEG
~session, ~Pasture, ~time, ~time, RDOccupEG
~session, ~Pasture, ~time*Pasture, ~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 ~time ~time RDOccupEG 2
## 3 ~session ~Pasture ~time ~time RDOccupEG 3
## 4 ~session ~Pasture ~time*Pasture ~time RDOccupEG 4
# 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 <- 3
n.season <- 5
# we need to process the data in the loop to allow for different data types
input.data <- process.data(data=input.history,
model=x$model.type,
time.intervals=c( rep( c(rep(0,max.visit.per.year-1),1),n.season-1),
rep(0,max.visit.per.year-1)),
group = "Pasture")
# 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: 1767.015
## AICc : 1775.043
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.4284162 0.1318158 -0.6867753 -0.1700572
## Epsilon:(Intercept) -2.1970555 0.2028550 -2.5946513 -1.7994598
## Gamma:(Intercept) -2.6132999 0.1864896 -2.9788196 -2.2477802
## p:(Intercept) 0.7961677 0.1077368 0.5850035 1.0073319
##
##
## Real Parameter Psi
## Group:PastureModified
## 1
## 0.3945046
##
## Group:PastureNative
## 1
## 0.3945046
##
##
## Real Parameter Epsilon
## Group:PastureModified
## 1 2 3 4
## 0.1000152 0.1000152 0.1000152 0.1000152
##
## Group:PastureNative
## 1 2 3 4
## 0.1000152 0.1000152 0.1000152 0.1000152
##
##
## Real Parameter Gamma
## Group:PastureModified
## 1 2 3 4
## 0.0682873 0.0682873 0.0682873 0.0682873
##
## Group:PastureNative
## 1 2 3 4
## 0.0682873 0.0682873 0.0682873 0.0682873
##
##
## Real Parameter p
## Session:1Group:PastureModified
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:2Group:PastureModified
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:3Group:PastureModified
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:4Group:PastureModified
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:5Group:PastureModified
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:1Group:PastureNative
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:2Group:PastureNative
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:3Group:PastureNative
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:4Group:PastureNative
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
## Session:5Group:PastureNative
## 1 2 3
## 0.6891541 0.6891541 0.6891541
##
##
## ***** Starting ~session ~1 ~time ~time RDOccupEG 2
##
## Note: only 13 parameters counted of 14 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~time)Gamma(~time)p(~session)
##
## Npar : 14 (unadjusted=13)
## -2lnL: 1746.601
## AICc : 1774.894 (unadjusted=1772.8549)
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.4759619 0.1447491 -0.7596701 -0.1922536
## Epsilon:(Intercept) -2.6094002 0.4617956 -3.5145196 -1.7042808
## Epsilon:time2 0.0640641 0.6757473 -1.2604007 1.3885289
## Epsilon:time3 0.7863089 0.5766713 -0.3439668 1.9165847
## Epsilon:time4 1.0252830 0.6207833 -0.1914523 2.2420182
## Gamma:(Intercept) -1.9164899 0.3300351 -2.5633587 -1.2696211
## Gamma:time2 -9.4583808 0.0000000 -9.4583808 -9.4583808
## Gamma:time3 -0.6370415 0.4996370 -1.6163300 0.3422470
## Gamma:time4 -0.3289008 0.4648862 -1.2400777 0.5822762
## p:(Intercept) 0.8248251 0.2720551 0.2915971 1.3580532
## p:session2 -0.2676570 0.3230589 -0.9008524 0.3655383
## p:session3 -0.0250358 0.3313501 -0.6744821 0.6244104
## p:session4 0.8511617 0.4194278 0.0290832 1.6732403
## p:session5 -0.1425150 0.3797182 -0.8867627 0.6017327
##
##
## Real Parameter Psi
## Group:PastureModified
## 1
## 0.3832061
##
## Group:PastureNative
## 1
## 0.3832061
##
##
## Real Parameter Epsilon
## Group:PastureModified
## 1 2 3 4
## 0.0685359 0.0727404 0.1390634 0.1702132
##
## Group:PastureNative
## 1 2 3 4
## 0.0685359 0.0727404 0.1390634 0.1702132
##
##
## Real Parameter Gamma
## Group:PastureModified
## 1 2 3 4
## 0.1282535 1.148025e-05 0.0721896 0.0957478
##
## Group:PastureNative
## 1 2 3 4
## 0.1282535 1.148025e-05 0.0721896 0.0957478
##
##
## Real Parameter p
## Session:1Group:PastureModified
## 1 2 3
## 0.6952596 0.6952596 0.6952596
##
## Session:2Group:PastureModified
## 1 2 3
## 0.635797 0.635797 0.635797
##
## Session:3Group:PastureModified
## 1 2 3
## 0.6899294 0.6899294 0.6899294
##
## Session:4Group:PastureModified
## 1 2 3
## 0.8423724 0.8423724 0.8423724
##
## Session:5Group:PastureModified
## 1 2 3
## 0.6642541 0.6642541 0.6642541
##
## Session:1Group:PastureNative
## 1 2 3
## 0.6952596 0.6952596 0.6952596
##
## Session:2Group:PastureNative
## 1 2 3
## 0.635797 0.635797 0.635797
##
## Session:3Group:PastureNative
## 1 2 3
## 0.6899294 0.6899294 0.6899294
##
## Session:4Group:PastureNative
## 1 2 3
## 0.8423724 0.8423724 0.8423724
##
## Session:5Group:PastureNative
## 1 2 3
## 0.6642541 0.6642541 0.6642541
##
##
## ***** Starting ~session ~Pasture ~time ~time RDOccupEG 3
##
## Output summary for RDOccupEG model
## Name : Psi(~Pasture)Epsilon(~time)Gamma(~time)p(~session)
##
## Npar : 15
## -2lnL: 1721.523
## AICc : 1751.858
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -1.2415714 0.2392899 -1.7105796 -0.7725631
## Psi:PastureNative 1.3614364 0.2890357 0.7949264 1.9279464
## Epsilon:(Intercept) -2.5512963 0.4529245 -3.4390283 -1.6635643
## Epsilon:time2 0.0478822 0.6600699 -1.2458549 1.3416193
## Epsilon:time3 0.7447768 0.5692779 -0.3710080 1.8605616
## Epsilon:time4 0.9678044 0.6139786 -0.2355937 2.1712025
## Gamma:(Intercept) -2.3394815 0.4974169 -3.3144185 -1.3645444
## Gamma:time2 -1.7282194 1.3249985 -4.3252166 0.8687778
## Gamma:time3 -0.2326664 0.6245554 -1.4567950 0.9914621
## Gamma:time4 0.0997778 0.5940123 -1.0644862 1.2640419
## p:(Intercept) 0.6045438 0.2844120 0.0470964 1.1619913
## p:session2 0.0317154 0.3400432 -0.6347693 0.6982002
## p:session3 0.1796942 0.3390874 -0.4849171 0.8443055
## p:session4 1.0606995 0.4257265 0.2262756 1.8951233
## p:session5 0.0786483 0.3889799 -0.6837524 0.8410490
##
##
## Real Parameter Psi
## Group:PastureModified
## 1
## 0.2241626
##
## Group:PastureNative
## 1
## 0.5299304
##
##
## Real Parameter Epsilon
## Group:PastureModified
## 1 2 3 4
## 0.0723394 0.0756192 0.1410593 0.1703015
##
## Group:PastureNative
## 1 2 3 4
## 0.0723394 0.0756192 0.1410593 0.1703015
##
##
## Real Parameter Gamma
## Group:PastureModified
## 1 2 3 4
## 0.0879055 0.0168286 0.0709526 0.0962413
##
## Group:PastureNative
## 1 2 3 4
## 0.0879055 0.0168286 0.0709526 0.0962413
##
##
## Real Parameter p
## Session:1Group:PastureModified
## 1 2 3
## 0.6466952 0.6466952 0.6466952
##
## Session:2Group:PastureModified
## 1 2 3
## 0.6539074 0.6539074 0.6539074
##
## Session:3Group:PastureModified
## 1 2 3
## 0.6865928 0.6865928 0.6865928
##
## Session:4Group:PastureModified
## 1 2 3
## 0.8409406 0.8409406 0.8409406
##
## Session:5Group:PastureModified
## 1 2 3
## 0.6644508 0.6644508 0.6644508
##
## Session:1Group:PastureNative
## 1 2 3
## 0.6466952 0.6466952 0.6466952
##
## Session:2Group:PastureNative
## 1 2 3
## 0.6539074 0.6539074 0.6539074
##
## Session:3Group:PastureNative
## 1 2 3
## 0.6865928 0.6865928 0.6865928
##
## Session:4Group:PastureNative
## 1 2 3
## 0.8409406 0.8409406 0.8409406
##
## Session:5Group:PastureNative
## 1 2 3
## 0.6644508 0.6644508 0.6644508
##
##
## ***** Starting ~session ~Pasture ~time*Pasture ~time RDOccupEG 4
##
## Note: only 18 parameters counted of 19 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for RDOccupEG model
## Name : Psi(~Pasture)Epsilon(~time)Gamma(~time * Pasture)p(~session)
##
## Npar : 19 (unadjusted=18)
## -2lnL: 1712.075
## AICc : 1750.607 (unadjusted=1748.5538)
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -1.2116258 0.2359949 -1.6741759 -0.7490757
## Psi:PastureNative 1.1807539 0.2909556 0.6104809 1.7510268
## Epsilon:(Intercept) -2.5239920 0.4497066 -3.4054169 -1.6425670
## Epsilon:time2 0.0298064 0.6562341 -1.2564125 1.3160254
## Epsilon:time3 0.7337478 0.5682024 -0.3799289 1.8474245
## Epsilon:time4 0.9431004 0.6117750 -0.2559786 2.1421794
## Gamma:(Intercept) -2.7293984 0.5774503 -3.8612011 -1.5975957
## Gamma:time2 -1.1243139 1.2922289 -3.6570825 1.4084548
## Gamma:time3 -1.0093895 1.2034924 -3.3682347 1.3494557
## Gamma:time4 0.2316498 0.7452488 -1.2290380 1.6923375
## Gamma:PastureNative 1.2643446 0.6733714 -0.0554635 2.5841526
## Gamma:time2:PastureNative -12.1146750 0.0000000 -12.1146750 -12.1146750
## Gamma:time3:PastureNative 0.5325512 1.3144784 -2.0438266 3.1089289
## Gamma:time4:PastureNative -0.7361701 0.9171421 -2.5337687 1.0614284
## p:(Intercept) 0.8072217 0.2711390 0.2757893 1.3386542
## p:session2 -0.1804033 0.3269619 -0.8212487 0.4604421
## p:session3 -0.0194137 0.3314947 -0.6691433 0.6303160
## p:session4 0.8518957 0.4236257 0.0215894 1.6822021
## p:session5 -0.1268057 0.3793323 -0.8702970 0.6166855
##
##
## Real Parameter Psi
## Group:PastureModified
## 1
## 0.2294135
##
## Group:PastureNative
## 1
## 0.4922826
##
##
## Real Parameter Epsilon
## Group:PastureModified
## 1 2 3 4
## 0.0741933 0.0762668 0.1430428 0.1706693
##
## Group:PastureNative
## 1 2 3 4
## 0.0741933 0.0762668 0.1430428 0.1706693
##
##
## Real Parameter Gamma
## Group:PastureModified
## 1 2 3 4
## 0.0612608 0.0207607 0.0232304 0.0760162
##
## Group:PastureNative
## 1 2 3 4
## 0.1876956 4.11259e-07 0.1254401 0.1224346
##
##
## Real Parameter p
## Session:1Group:PastureModified
## 1 2 3
## 0.6915172 0.6915172 0.6915172
##
## Session:2Group:PastureModified
## 1 2 3
## 0.6517677 0.6517677 0.6517677
##
## Session:3Group:PastureModified
## 1 2 3
## 0.6873605 0.6873605 0.6873605
##
## Session:4Group:PastureModified
## 1 2 3
## 0.8401195 0.8401195 0.8401195
##
## Session:5Group:PastureModified
## 1 2 3
## 0.6638315 0.6638315 0.6638315
##
## Session:1Group:PastureNative
## 1 2 3
## 0.6915172 0.6915172 0.6915172
##
## Session:2Group:PastureNative
## 1 2 3
## 0.6517677 0.6517677 0.6517677
##
## Session:3Group:PastureNative
## 1 2 3
## 0.6873605 0.6873605 0.6873605
##
## Session:4Group:PastureNative
## 1 2 3
## 0.8401195 0.8401195 0.8401195
##
## Session:5Group:PastureNative
## 1 2 3
## 0.6638315 0.6638315 0.6638315
# examine individula model results
model.number <-4
summary(model.fits[[model.number]])
## Output summary for RDOccupEG model
## Name : Psi(~Pasture)Epsilon(~time)Gamma(~time * Pasture)p(~session)
##
## Npar : 19 (unadjusted=18)
## -2lnL: 1712.075
## AICc : 1750.607 (unadjusted=1748.5538)
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -1.2116258 0.2359949 -1.6741759 -0.7490757
## Psi:PastureNative 1.1807539 0.2909556 0.6104809 1.7510268
## Epsilon:(Intercept) -2.5239920 0.4497066 -3.4054169 -1.6425670
## Epsilon:time2 0.0298064 0.6562341 -1.2564125 1.3160254
## Epsilon:time3 0.7337478 0.5682024 -0.3799289 1.8474245
## Epsilon:time4 0.9431004 0.6117750 -0.2559786 2.1421794
## Gamma:(Intercept) -2.7293984 0.5774503 -3.8612011 -1.5975957
## Gamma:time2 -1.1243139 1.2922289 -3.6570825 1.4084548
## Gamma:time3 -1.0093895 1.2034924 -3.3682347 1.3494557
## Gamma:time4 0.2316498 0.7452488 -1.2290380 1.6923375
## Gamma:PastureNative 1.2643446 0.6733714 -0.0554635 2.5841526
## Gamma:time2:PastureNative -12.1146750 0.0000000 -12.1146750 -12.1146750
## Gamma:time3:PastureNative 0.5325512 1.3144784 -2.0438266 3.1089289
## Gamma:time4:PastureNative -0.7361701 0.9171421 -2.5337687 1.0614284
## p:(Intercept) 0.8072217 0.2711390 0.2757893 1.3386542
## p:session2 -0.1804033 0.3269619 -0.8212487 0.4604421
## p:session3 -0.0194137 0.3314947 -0.6691433 0.6303160
## p:session4 0.8518957 0.4236257 0.0215894 1.6822021
## p:session5 -0.1268057 0.3793323 -0.8702970 0.6166855
##
##
## Real Parameter Psi
## Group:PastureModified
## 1
## 0.2294135
##
## Group:PastureNative
## 1
## 0.4922826
##
##
## Real Parameter Epsilon
## Group:PastureModified
## 1 2 3 4
## 0.0741933 0.0762668 0.1430428 0.1706693
##
## Group:PastureNative
## 1 2 3 4
## 0.0741933 0.0762668 0.1430428 0.1706693
##
##
## Real Parameter Gamma
## Group:PastureModified
## 1 2 3 4
## 0.0612608 0.0207607 0.0232304 0.0760162
##
## Group:PastureNative
## 1 2 3 4
## 0.1876956 4.11259e-07 0.1254401 0.1224346
##
##
## Real Parameter p
## Session:1Group:PastureModified
## 1 2 3
## 0.6915172 0.6915172 0.6915172
##
## Session:2Group:PastureModified
## 1 2 3
## 0.6517677 0.6517677 0.6517677
##
## Session:3Group:PastureModified
## 1 2 3
## 0.6873605 0.6873605 0.6873605
##
## Session:4Group:PastureModified
## 1 2 3
## 0.8401195 0.8401195 0.8401195
##
## Session:5Group:PastureModified
## 1 2 3
## 0.6638315 0.6638315 0.6638315
##
## Session:1Group:PastureNative
## 1 2 3
## 0.6915172 0.6915172 0.6915172
##
## Session:2Group:PastureNative
## 1 2 3
## 0.6517677 0.6517677 0.6517677
##
## Session:3Group:PastureNative
## 1 2 3
## 0.6873605 0.6873605 0.6873605
##
## Session:4Group:PastureNative
## 1 2 3
## 0.8401195 0.8401195 0.8401195
##
## Session:5Group:PastureNative
## 1 2 3
## 0.6638315 0.6638315 0.6638315
model.fits[[model.number]]$results$real
## estimate se lcl ucl
## Psi gModified a0 t1 2.294135e-01 0.0417199 1.578682e-01 3.210227e-01
## Psi gNative a0 t1 4.922826e-01 0.0473034 4.008754e-01 5.842087e-01
## Epsilon gModified a0 t1 7.419330e-02 0.0308897 3.212660e-02 1.621161e-01
## Epsilon gModified a1 t2 7.626680e-02 0.0314465 3.327650e-02 1.653000e-01
## Epsilon gModified a2 t3 1.430428e-01 0.0419580 7.862880e-02 2.461294e-01
## Epsilon gModified a3 t4 1.706693e-01 0.0585849 8.377340e-02 3.165582e-01
## Gamma gModified a0 t1 6.126080e-02 0.0332079 2.060900e-02 1.683179e-01
## Gamma gModified a1 t2 2.076070e-02 0.0225531 2.404400e-03 1.571807e-01
## Gamma gModified a2 t3 2.323040e-02 0.0229242 3.272400e-03 1.469623e-01
## Gamma gModified a3 t4 7.601620e-02 0.0329348 3.177450e-02 1.709802e-01
## Gamma gNative a0 t1 1.876956e-01 0.0614644 9.490190e-02 3.373987e-01
## Gamma gNative a1 t2 4.112590e-07 0.0000000 4.112590e-07 4.112590e-07
## Gamma gNative a2 t3 1.254401e-01 0.0426906 6.270200e-02 2.352003e-01
## Gamma gNative a3 t4 1.224346e-01 0.0466704 5.620330e-02 2.463430e-01
## p gModified s1 t1 6.915172e-01 0.0578397 5.685136e-01 7.922685e-01
## p gModified s2 t1 6.517677e-01 0.0423119 5.649875e-01 7.295239e-01
## p gModified s3 t1 6.873605e-01 0.0414540 6.010204e-01 7.624027e-01
## p gModified s4 t1 8.401195e-01 0.0439474 7.345498e-01 9.089105e-01
## p gModified s5 t1 6.638315e-01 0.0593554 5.396912e-01 7.688328e-01
## fixed note
## Psi gModified a0 t1
## Psi gNative a0 t1
## Epsilon gModified a0 t1
## Epsilon gModified a1 t2
## Epsilon gModified a2 t3
## Epsilon gModified a3 t4
## Gamma gModified a0 t1
## Gamma gModified a1 t2
## Gamma gModified a2 t3
## Gamma gModified a3 t4
## Gamma gNative a0 t1
## Gamma gNative a1 t2
## Gamma gNative a2 t3
## Gamma gNative a3 t4
## p gModified s1 t1
## p gModified s2 t1
## p gModified s3 t1
## p gModified s4 t1
## p gModified s5 t1
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## Psi:(Intercept) -1.2116258 0.2359949 -1.6741759 -0.7490757
## Psi:PastureNative 1.1807539 0.2909556 0.6104809 1.7510268
## Epsilon:(Intercept) -2.5239920 0.4497066 -3.4054169 -1.6425670
## Epsilon:time2 0.0298064 0.6562341 -1.2564125 1.3160254
## Epsilon:time3 0.7337478 0.5682024 -0.3799289 1.8474245
## Epsilon:time4 0.9431004 0.6117750 -0.2559786 2.1421794
## Gamma:(Intercept) -2.7293984 0.5774503 -3.8612011 -1.5975957
## Gamma:time2 -1.1243139 1.2922289 -3.6570825 1.4084548
## Gamma:time3 -1.0093895 1.2034924 -3.3682347 1.3494557
## Gamma:time4 0.2316498 0.7452488 -1.2290380 1.6923375
## Gamma:PastureNative 1.2643446 0.6733714 -0.0554635 2.5841526
## Gamma:time2:PastureNative -12.1146750 0.0000000 -12.1146750 -12.1146750
## Gamma:time3:PastureNative 0.5325512 1.3144784 -2.0438266 3.1089289
## Gamma:time4:PastureNative -0.7361701 0.9171421 -2.5337687 1.0614284
## p:(Intercept) 0.8072217 0.2711390 0.2757893 1.3386542
## p:session2 -0.1804033 0.3269619 -0.8212487 0.4604421
## p:session3 -0.0194137 0.3314947 -0.6691433 0.6303160
## p:session4 0.8518957 0.4236257 0.0215894 1.6822021
## p:session5 -0.1268057 0.3793323 -0.8702970 0.6166855
#NOTE: Because a group was used in the process.data step above, the first half
#of each derived parameter table will be for Modified Pastures, and the second
#will be for Native pastures
model.fits[[model.number]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.2294135 0.04171988 0.1476425 0.3111845
## 2 0.2595993 0.04023803 0.1807327 0.3384658
## 3 0.2551717 0.03906549 0.1786034 0.3317401
## 4 0.2359739 0.03486851 0.1676317 0.3043162
## 5 0.2537788 0.03845114 0.1784145 0.3291430
## 6 0.4922826 0.04730337 0.3995680 0.5849972
## 7 0.5510549 0.03938525 0.4738598 0.6282500
## 8 0.5090279 0.03888970 0.4328041 0.5852517
## 9 0.4978027 0.03859611 0.4221543 0.5734511
## 10 0.4743294 0.04562270 0.3849089 0.5637499
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.1315780 0.13398752 0.8689624 1.394194
## 2 0.9829447 0.07403097 0.8378440 1.128045
## 3 0.9247652 0.08195122 0.7641408 1.085390
## 4 1.0754525 0.13639162 0.8081250 1.342780
## 5 1.1193872 0.08918858 0.9445776 1.294197
## 6 0.9237335 0.03144617 0.8620990 0.985368
## 7 0.9779479 0.06348516 0.8535169 1.102379
## 8 0.9528462 0.08161287 0.7928850 1.112807
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 1.1777119 0.18252613 0.8199606 1.5354631
## 2 0.9771017 0.09913543 0.7827963 1.1714072
## 3 0.9015285 0.10627959 0.6932205 1.1098365
## 4 1.1011128 0.18516768 0.7381842 1.4640415
## 5 1.2659283 0.20450086 0.8651066 1.6667500
## 6 0.8446623 0.06067709 0.7257352 0.9635894
## 7 0.9560887 0.12498040 0.7111271 1.2010503
## 8 0.9102979 0.15037607 0.6155608 1.2050350
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.2294135 0.04171988 0.1476425 0.3111845
## 2 0.2595993 0.04023803 0.1807327 0.3384658
## 3 0.2551717 0.03906549 0.1786034 0.3317401
## 4 0.2359739 0.03486851 0.1676317 0.3043162
## 5 0.2537788 0.03845114 0.1784145 0.3291430
## 6 0.4922826 0.04730337 0.3995680 0.5849972
## 7 0.5510549 0.03938525 0.4738598 0.6282500
## 8 0.5090279 0.03888970 0.4328041 0.5852517
## 9 0.4978027 0.03859611 0.4221543 0.5734511
## 10 0.4743294 0.04562270 0.3849089 0.5637499
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 1.1315780 0.13398752 0.8689624 1.394194
## 2 0.9829447 0.07403097 0.8378440 1.128045
## 3 0.9247652 0.08195122 0.7641408 1.085390
## 4 1.0754525 0.13639162 0.8081250 1.342780
## 5 1.1193872 0.08918858 0.9445776 1.294197
## 6 0.9237335 0.03144617 0.8620990 0.985368
## 7 0.9779479 0.06348516 0.8535169 1.102379
## 8 0.9528462 0.08161287 0.7928850 1.112807
model.fits[[model.number]]$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 1.1777119 0.18252613 0.8199606 1.5354631
## 2 0.9771017 0.09913543 0.7827963 1.1714072
## 3 0.9015285 0.10627959 0.6932205 1.1098365
## 4 1.1011128 0.18516768 0.7381842 1.4640415
## 5 1.2659283 0.20450086 0.8651066 1.6667500
## 6 0.8446623 0.06067709 0.7257352 0.9635894
## 7 0.9560887 0.12498040 0.7111271 1.2010503
## 8 0.9102979 0.15037607 0.6155608 1.2050350
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 gModified a0 t1 1 1 0.2294135 0.0417199 0.1578682
## Psi gNative a0 t1 2 2 0.4922826 0.0473034 0.4008754
## ucl fixed note group age time Age Time
## Psi gModified a0 t1 0.3210227 Modified 0 1 0 0
## Psi gNative a0 t1 0.5842087 Native 0 1 0 0
## Pasture
## Psi gModified a0 t1 Modified
## Psi gNative a0 t1 Native
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 gModified s1 t1 19 15 0.6915172 0.0578397 0.5685136
## p gModified s1 t2 20 15 0.6915172 0.0578397 0.5685136
## p gModified s1 t3 21 15 0.6915172 0.0578397 0.5685136
## p gModified s2 t1 22 16 0.6517677 0.0423119 0.5649875
## p gModified s2 t2 23 16 0.6517677 0.0423119 0.5649875
## p gModified s2 t3 24 16 0.6517677 0.0423119 0.5649875
## p gModified s3 t1 25 17 0.6873605 0.0414540 0.6010204
## p gModified s3 t2 26 17 0.6873605 0.0414540 0.6010204
## p gModified s3 t3 27 17 0.6873605 0.0414540 0.6010204
## p gModified s4 t1 28 18 0.8401195 0.0439474 0.7345498
## p gModified s4 t2 29 18 0.8401195 0.0439474 0.7345498
## p gModified s4 t3 30 18 0.8401195 0.0439474 0.7345498
## p gModified s5 t1 31 19 0.6638315 0.0593554 0.5396912
## p gModified s5 t2 32 19 0.6638315 0.0593554 0.5396912
## p gModified s5 t3 33 19 0.6638315 0.0593554 0.5396912
## p gNative s1 t1 34 15 0.6915172 0.0578397 0.5685136
## p gNative s1 t2 35 15 0.6915172 0.0578397 0.5685136
## p gNative s1 t3 36 15 0.6915172 0.0578397 0.5685136
## p gNative s2 t1 37 16 0.6517677 0.0423119 0.5649875
## p gNative s2 t2 38 16 0.6517677 0.0423119 0.5649875
## p gNative s2 t3 39 16 0.6517677 0.0423119 0.5649875
## p gNative s3 t1 40 17 0.6873605 0.0414540 0.6010204
## p gNative s3 t2 41 17 0.6873605 0.0414540 0.6010204
## p gNative s3 t3 42 17 0.6873605 0.0414540 0.6010204
## p gNative s4 t1 43 18 0.8401195 0.0439474 0.7345498
## p gNative s4 t2 44 18 0.8401195 0.0439474 0.7345498
## p gNative s4 t3 45 18 0.8401195 0.0439474 0.7345498
## p gNative s5 t1 46 19 0.6638315 0.0593554 0.5396912
## p gNative s5 t2 47 19 0.6638315 0.0593554 0.5396912
## p gNative s5 t3 48 19 0.6638315 0.0593554 0.5396912
## ucl fixed note group time session Time
## p gModified s1 t1 0.7922685 Modified 1 1 0
## p gModified s1 t2 0.7922685 Modified 2 1 1
## p gModified s1 t3 0.7922685 Modified 3 1 2
## p gModified s2 t1 0.7295239 Modified 1 2 0
## p gModified s2 t2 0.7295239 Modified 2 2 1
## p gModified s2 t3 0.7295239 Modified 3 2 2
## p gModified s3 t1 0.7624027 Modified 1 3 0
## p gModified s3 t2 0.7624027 Modified 2 3 1
## p gModified s3 t3 0.7624027 Modified 3 3 2
## p gModified s4 t1 0.9089105 Modified 1 4 0
## p gModified s4 t2 0.9089105 Modified 2 4 1
## p gModified s4 t3 0.9089105 Modified 3 4 2
## p gModified s5 t1 0.7688328 Modified 1 5 0
## p gModified s5 t2 0.7688328 Modified 2 5 1
## p gModified s5 t3 0.7688328 Modified 3 5 2
## p gNative s1 t1 0.7922685 Native 1 1 0
## p gNative s1 t2 0.7922685 Native 2 1 1
## p gNative s1 t3 0.7922685 Native 3 1 2
## p gNative s2 t1 0.7295239 Native 1 2 0
## p gNative s2 t2 0.7295239 Native 2 2 1
## p gNative s2 t3 0.7295239 Native 3 2 2
## p gNative s3 t1 0.7624027 Native 1 3 0
## p gNative s3 t2 0.7624027 Native 2 3 1
## p gNative s3 t3 0.7624027 Native 3 3 2
## p gNative s4 t1 0.9089105 Native 1 4 0
## p gNative s4 t2 0.9089105 Native 2 4 1
## p gNative s4 t3 0.9089105 Native 3 4 2
## p gNative s5 t1 0.7688328 Native 1 5 0
## p gNative s5 t2 0.7688328 Native 2 5 1
## p gNative s5 t3 0.7688328 Native 3 5 2
## Pasture
## p gModified s1 t1 Modified
## p gModified s1 t2 Modified
## p gModified s1 t3 Modified
## p gModified s2 t1 Modified
## p gModified s2 t2 Modified
## p gModified s2 t3 Modified
## p gModified s3 t1 Modified
## p gModified s3 t2 Modified
## p gModified s3 t3 Modified
## p gModified s4 t1 Modified
## p gModified s4 t2 Modified
## p gModified s4 t3 Modified
## p gModified s5 t1 Modified
## p gModified s5 t2 Modified
## p gModified s5 t3 Modified
## p gNative s1 t1 Native
## p gNative s1 t2 Native
## p gNative s1 t3 Native
## p gNative s2 t1 Native
## p gNative s2 t2 Native
## p gNative s2 t3 Native
## p gNative s3 t1 Native
## p gNative s3 t2 Native
## p gNative s3 t3 Native
## p gNative s4 t1 Native
## p gNative s4 t2 Native
## p gNative s4 t3 Native
## p gNative s5 t1 Native
## p gNative s5 t2 Native
## p gNative s5 t3 Native
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
model.set
## model npar
## 4 Psi(~Pasture)Epsilon(~time)Gamma(~time * Pasture)p(~session) 19
## 3 Psi(~Pasture)Epsilon(~time)Gamma(~time)p(~session) 15
## 2 Psi(~1)Epsilon(~time)Gamma(~time)p(~session) 14
## 1 Psi(~1)Epsilon(~1)Gamma(~1)p(~1) 4
## AICc DeltaAICc weight Deviance
## 4 1750.607 0.000000 6.514053e-01 -1130.993
## 3 1751.858 1.250483 3.485880e-01 -1121.545
## 2 1774.894 24.286679 3.467896e-06 -1096.467
## 1 1775.043 24.435807 3.218720e-06 -1076.053
# 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.
## 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
## Psi gModified a0 t1 1 0.2275842 0.04176069 Modified
## Psi gNative a0 t1 2 0.5054055 0.05157205 Native
## age time Age Time Pasture
## Psi gModified a0 t1 0 1 0 0 Modified
## Psi gNative a0 t1 0 1 0 0 Native
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.
## 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
## p gModified s1 t1 19 0.6758927 0.06408927 Modified
## p gModified s1 t2 20 0.6758927 0.06408927 Modified
## p gModified s1 t3 21 0.6758927 0.06408927 Modified
## p gModified s2 t1 22 0.6525136 0.04278745 Modified
## p gModified s2 t2 23 0.6525136 0.04278745 Modified
## p gModified s2 t3 24 0.6525136 0.04278745 Modified
## p gModified s3 t1 25 0.6870929 0.04138822 Modified
## p gModified s3 t2 26 0.6870929 0.04138822 Modified
## p gModified s3 t3 27 0.6870929 0.04138822 Modified
## p gModified s4 t1 28 0.8404052 0.04358461 Modified
## p gModified s4 t2 29 0.8404052 0.04358461 Modified
## p gModified s4 t3 30 0.8404052 0.04358461 Modified
## p gModified s5 t1 31 0.6640475 0.05932827 Modified
## p gModified s5 t2 32 0.6640475 0.05932827 Modified
## p gModified s5 t3 33 0.6640475 0.05932827 Modified
## p gNative s1 t1 34 0.6758927 0.06408927 Native
## p gNative s1 t2 35 0.6758927 0.06408927 Native
## p gNative s1 t3 36 0.6758927 0.06408927 Native
## p gNative s2 t1 37 0.6525136 0.04278745 Native
## p gNative s2 t2 38 0.6525136 0.04278745 Native
## p gNative s2 t3 39 0.6525136 0.04278745 Native
## p gNative s3 t1 40 0.6870929 0.04138822 Native
## p gNative s3 t2 41 0.6870929 0.04138822 Native
## p gNative s3 t3 42 0.6870929 0.04138822 Native
## p gNative s4 t1 43 0.8404052 0.04358461 Native
## p gNative s4 t2 44 0.8404052 0.04358461 Native
## p gNative s4 t3 45 0.8404052 0.04358461 Native
## p gNative s5 t1 46 0.6640475 0.05932827 Native
## p gNative s5 t2 47 0.6640475 0.05932827 Native
## p gNative s5 t3 48 0.6640475 0.05932827 Native
## time session Time Pasture
## p gModified s1 t1 1 1 0 Modified
## p gModified s1 t2 2 1 1 Modified
## p gModified s1 t3 3 1 2 Modified
## p gModified s2 t1 1 2 0 Modified
## p gModified s2 t2 2 2 1 Modified
## p gModified s2 t3 3 2 2 Modified
## p gModified s3 t1 1 3 0 Modified
## p gModified s3 t2 2 3 1 Modified
## p gModified s3 t3 3 3 2 Modified
## p gModified s4 t1 1 4 0 Modified
## p gModified s4 t2 2 4 1 Modified
## p gModified s4 t3 3 4 2 Modified
## p gModified s5 t1 1 5 0 Modified
## p gModified s5 t2 2 5 1 Modified
## p gModified s5 t3 3 5 2 Modified
## p gNative s1 t1 1 1 0 Native
## p gNative s1 t2 2 1 1 Native
## p gNative s1 t3 3 1 2 Native
## p gNative s2 t1 1 2 0 Native
## p gNative s2 t2 2 2 1 Native
## p gNative s2 t3 3 2 2 Native
## p gNative s3 t1 1 3 0 Native
## p gNative s3 t2 2 3 1 Native
## p gNative s3 t3 3 3 2 Native
## p gNative s4 t1 1 4 0 Native
## p gNative s4 t2 2 4 1 Native
## p gNative s4 t3 3 4 2 Native
## p gNative s5 t1 1 5 0 Native
## p gNative s5 t2 2 5 1 Native
## p gNative s5 t3 3 5 2 Native
# 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 Modified Pastures, and the second
#will be for Native pastures
RMark.model.average.derived(model.set, "all_psi")
## Loading required package: plyr
## Loading required package: boot
##
## Attaching package: 'boot'
## The following object is masked from 'package:car':
##
## logit
## estimate se lcl ucl
## 1 0.2275842 0.04176069 0.1560900 0.3194300
## 2 0.2653687 0.04148887 0.1922637 0.3540854
## 3 0.2594516 0.03933514 0.1899901 0.3435379
## 4 0.2519136 0.04110273 0.1800643 0.3405256
## 5 0.2708868 0.04437625 0.1930214 0.3659201
## 6 0.5054055 0.05157204 0.4054658 0.6049151
## 7 0.5447314 0.04042789 0.4650320 0.6222035
## 8 0.5060471 0.03892168 0.4302097 0.5816072
## 9 0.4864803 0.04078150 0.4075635 0.5660770
## 10 0.4615033 0.04725620 0.3712208 0.5543838
RMark.model.average.derived(model.set, "lambda Rate of Change")
## estimate se lcl ucl
## 1 1.1665506 0.15610787 0.8605848 1.4725164
## 2 0.9779069 0.06986316 0.8409776 1.1148362
## 3 0.9695658 0.10433957 0.7650640 1.1740676
## 4 1.0753309 0.12732884 0.8257709 1.3248908
## 5 1.0797347 0.09571588 0.8921351 1.2673344
## 6 0.9291010 0.03368414 0.8630813 0.9951207
## 7 0.9611494 0.06330744 0.8370691 1.0852297
## 8 0.9484673 0.07878500 0.7940516 1.1028831
RMark.model.average.derived(model.set, "log odds lambda")
## estimate se lcl ucl
## 1 1.2274427 0.21643542 0.8032371 1.6516484
## 2 0.9700909 0.09399868 0.7858568 1.1543249
## 3 0.9617135 0.13966705 0.6879711 1.2354559
## 4 1.1034195 0.17610649 0.7582571 1.4485819
## 5 1.1774345 0.21597697 0.7541274 1.6007416
## 6 0.8563344 0.06552356 0.7279106 0.9847582
## 7 0.9256004 0.12134654 0.6877656 1.1634353
## 8 0.9045581 0.14211298 0.6260218 1.1830944
# 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 10 estimates, but only 5 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$Pasture = c(rep("Modified",5),rep("Native",5))
psi.ma
## estimate se lcl ucl Year parameter Pasture
## 1 0.2275842 0.04176069 0.1560900 0.3194300 1 psi Modified
## 2 0.2653687 0.04148887 0.1922637 0.3540854 2 psi Modified
## 3 0.2594516 0.03933514 0.1899901 0.3435379 3 psi Modified
## 4 0.2519136 0.04110273 0.1800643 0.3405256 4 psi Modified
## 5 0.2708868 0.04437625 0.1930214 0.3659201 5 psi Modified
## 6 0.5054055 0.05157204 0.4054658 0.6049151 1 psi Native
## 7 0.5447314 0.04042789 0.4650320 0.6222035 2 psi Native
## 8 0.5060471 0.03892168 0.4302097 0.5816072 3 psi Native
## 9 0.4864803 0.04078150 0.4075635 0.5660770 4 psi Native
## 10 0.4615033 0.04725620 0.3712208 0.5543838 5 psi Native
# 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.
## 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
## Epsilon gModified a0 t1 3 0.07354712 0.03073059 0.03175632
## Epsilon gModified a1 t2 4 0.07604111 0.03132233 0.03320869
## Epsilon gModified a2 t3 5 0.14235122 0.04165533 0.07837575
## Epsilon gModified a3 t4 6 0.17054083 0.05852733 0.08372883
## Epsilon gNative a0 t1 7 0.07354712 0.03073059 0.03175632
## Epsilon gNative a1 t2 8 0.07604111 0.03132233 0.03320869
## Epsilon gNative a2 t3 9 0.14235122 0.04165533 0.07837575
## Epsilon gNative a3 t4 10 0.17054083 0.05852733 0.08372883
## ucl fixed note group age time Age Time
## Epsilon gModified a0 t1 0.1611790 Modified 0 1 0 0
## Epsilon gModified a1 t2 0.1647068 Modified 1 2 1 1
## Epsilon gModified a2 t3 0.2446834 Modified 2 3 2 2
## Epsilon gModified a3 t4 0.3162912 Modified 3 4 3 3
## Epsilon gNative a0 t1 0.1611790 Native 0 1 0 0
## Epsilon gNative a1 t2 0.1647068 Native 1 2 1 1
## Epsilon gNative a2 t3 0.2446834 Native 2 3 2 2
## Epsilon gNative a3 t4 0.3162912 Native 3 4 3 3
## Pasture Year parameter
## Epsilon gModified a0 t1 Modified 1 epsilon
## Epsilon gModified a1 t2 Modified 2 epsilon
## Epsilon gModified a2 t3 Modified 3 epsilon
## Epsilon gModified a3 t4 Modified 4 epsilon
## Epsilon gNative a0 t1 Native 1 epsilon
## Epsilon gNative a1 t2 Native 2 epsilon
## Epsilon gNative a2 t3 Native 3 epsilon
## Epsilon gNative a3 t4 Native 4 epsilon
gamma.ma <- RMark::model.average(model.set, "Gamma", vcv=TRUE)$estimates
##
## Model 2dropped from the model averaging because one or more beta variances are not positive
##
## Model 4dropped 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
## Gamma gModified a0 t1 11 0.08790530 0.03988182 0.035080407
## Gamma gModified a1 t2 12 0.01682912 0.01884814 0.001832299
## Gamma gModified a2 t3 13 0.07095256 0.02383591 0.036233428
## Gamma gModified a3 t4 14 0.09624105 0.02828349 0.053300198
## Gamma gNative a0 t1 15 0.08790530 0.03988182 0.035080407
## Gamma gNative a1 t2 16 0.01682912 0.01884814 0.001832299
## Gamma gNative a2 t3 17 0.07095256 0.02383591 0.036233428
## Gamma gNative a3 t4 18 0.09624105 0.02828349 0.053300198
## ucl fixed note group age time Age Time
## Gamma gModified a0 t1 0.2034992 Modified 0 1 0 0
## Gamma gModified a1 t2 0.1376444 Modified 1 2 1 1
## Gamma gModified a2 t3 0.1343038 Modified 2 3 2 2
## Gamma gModified a3 t4 0.1676505 Modified 3 4 3 3
## Gamma gNative a0 t1 0.2034992 Native 0 1 0 0
## Gamma gNative a1 t2 0.1376444 Native 1 2 1 1
## Gamma gNative a2 t3 0.1343038 Native 2 3 2 2
## Gamma gNative a3 t4 0.1676505 Native 3 4 3 3
## Pasture Year parameter
## Gamma gModified a0 t1 Modified 1 gamma
## Gamma gModified a1 t2 Modified 2 gamma
## Gamma gModified a2 t3 Modified 3 gamma
## Gamma gModified a3 t4 Modified 4 gamma
## Gamma gNative a0 t1 Native 1 gamma
## Gamma gNative a1 t2 Native 2 gamma
## Gamma gNative a2 t3 Native 3 gamma
## Gamma gNative a3 t4 Native 4 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)+
facet_wrap(~Pasture)
cleanup(ask=FALSE)