# Multi state single season
# Nichols et al (2007) Ecology 88:1395-1400 looked at breeding and
# non-breeding California Spotted Owls in s = 54 sites over k = 5
# surveys.
# Hand-fed mice con rmed occupancy; breeding status only
# confrmed when owl took mouse to nest.
# . = site not visited.
# 0 = owl not detected.
# 1 = owl detected, but no detection of young.
# 2 = owl detected, along with evidence of young.
# Territories that were expected to be occupied were selectively monitored
# and so psi is expected to be close to 1 and not really of interest.
# See MacKenzie et al, Section 5.7.1 for more details
# 2018-10-01 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
# General model averaging framework
library(ggplot2)
library(plyr)
library(readxl)
library(reshape2)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
# get the data
occ.data <- readxl::read_excel(file.path("..","CalSpottedOwl.xls"),
sheet="RawData",
skip=10, na='.')
head(occ.data)
## # A tibble: 6 x 6
## Site `1` `2` `3` `4` `5`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 NA 2 2 1
## 2 2 0 1 1 NA NA
## 3 3 1 NA 1 NA NA
## 4 4 1 NA 1 NA NA
## 5 5 NA 1 2 1 2
## 6 6 1 1 1 1 NA
# create the capture history
input.history <- data.frame(ch =apply(occ.data[,-1],1,paste,sep="", collapse=""),
freq=1)
input.history[1:5,]
## ch freq
## 1 1NA221 1
## 2 011NANA 1
## 3 1NA1NANA 1
## 4 1NA1NANA 1
## 5 NA1212 1
# need to deal with missing values
input.history$ch <- gsub("NA",".",input.history$ch)
head(input.history)
## ch freq
## 1 1.221 1
## 2 011.. 1
## 3 1.1.. 1
## 4 1.1.. 1
## 5 .1212 1
## 6 1111. 1
# For biological reasons, we believe that errors in identifying states
# varies between the first 2 visits and the last 3 visit (early bs late
# breeding season). We create a site-time covariates.
# and append to the input.history data.frame
# The covariate matrix would be arranged originally as
# 1 2 3 4 5
# site1 [E E L L L]
# site2 [E E L L L]
# : [: : : : :]
# siteN [E E L L L]
#
# and this is stacked by columns
# create the Rmark data file
owl.data <- process.data(data=input.history,
model="MSOccupancy") # this parameterization is more stable
summary(owl.data)
## Length Class Mode
## data 2 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 54 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 4 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 1 -none- numeric
## group.covariates 0 -none- NULL
## nstrata 1 -none- numeric
## strata.labels 0 -none- NULL
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
# Time-specific covariates are added to the ddl's.
# We need to modify the design data to create this design matrix
# See the dipper example in Section C.6 of the RMark manual
owl.ddl=make.design.data(owl.data)
names(owl.ddl)
## [1] "Psi1" "Psi2" "p1" "p2" "Delta" "pimtypes"
owl.ddl$Delta
## par.index model.index group age time Age Time
## 1 1 13 1 0 1 0 0
## 2 2 14 1 1 2 1 1
## 3 3 15 1 2 3 2 2
## 4 4 16 1 3 4 3 3
## 5 5 17 1 4 5 4 4
owl.ddl$Delta$Per <- factor(c("E","E","L","L","L"))
owl.ddl$Delta
## par.index model.index group age time Age Time Per
## 1 1 13 1 0 1 0 0 E
## 2 2 14 1 1 2 1 1 E
## 3 3 15 1 2 3 2 2 L
## 4 4 16 1 3 4 3 3 L
## 5 5 17 1 4 5 4 4 L
str(owl.ddl$Delta)
## 'data.frame': 5 obs. of 8 variables:
## $ par.index : int 1 2 3 4 5
## $ model.index: num 13 14 15 16 17
## $ group : Factor w/ 1 level "1": 1 1 1 1 1
## $ age : Factor w/ 5 levels "0","1","2","3",..: 1 2 3 4 5
## $ time : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5
## $ Age : num 0 1 2 3 4
## $ Time : num 0 1 2 3 4
## $ Per : Factor w/ 2 levels "E","L": 1 1 2 2 2
setup.parameters("MSOccupancy", check=TRUE)
## [1] "Psi1" "Psi2" "p1" "p2" "Delta"
# Use the multi-state model, but only for a single season
# The real parameters are
# Psi1(i) = probability that a site i is occupied regardless
# of reproductive state, Pr(true state = 1 or 2);
# Psi2(i) = probability that young occurred, given that the site i is occupied,
# Pr(true state = 2 | true state = 1 or 2);
#
#
# p1(i, t) = probability that occupancy is detected for site i,
# period t, given that true state = 1, Pr(detection | true state = 1);
# p2(i, t) = probability that occupancy is detected for site i,
# period t, given that true state = 2, Pr(detection | true state = 2);
# Delta(i, t) = probability that evidence of successful reproduction is found,
# given detection of occupancy at site i, period t, with successful
# reproduction, Pr(classified state 2|true state = 2).
# Thus, psi1 and psi2 are single parameters relating to a site,
# whereas p1, p2, and delta are time-specific parameters,
# and thus the number of values in the PIM for each of these parameters equals the number of occasions.
#
# The product psi1*psi2 is the unconditional probability that a site is
# occupied by successful breeders, and is provided as a derived parameter.
# define the list of models to fit
# Notice the commas between the column and the placement of the quotes
# Notice the difference between Time and time in the model terms
model.list.csv <- textConnection("
Psi1, Psi2, p1, p2, Delta
~1, ~1, ~1, ~SHARE, ~1
~1, ~1, ~time, ~1, ~1
~1, ~1, ~1, ~1, ~1
~1, ~1, ~1, ~SHARE, ~Per
~1, ~1, ~time, ~SHARE, ~Per
~1, ~1, ~1, ~1, ~Per")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## Psi1 Psi2 p1 p2 Delta
## 1 ~1 ~1 ~1 ~SHARE ~1
## 2 ~1 ~1 ~time ~1 ~1
## 3 ~1 ~1 ~1 ~1 ~1
## 4 ~1 ~1 ~1 ~SHARE ~Per
## 5 ~1 ~1 ~time ~SHARE ~Per
## 6 ~1 ~1 ~1 ~1 ~Per
model.list$model.number <- 1:nrow(model.list)
model.list
## Psi1 Psi2 p1 p2 Delta model.number
## 1 ~1 ~1 ~1 ~SHARE ~1 1
## 2 ~1 ~1 ~time ~1 ~1 2
## 3 ~1 ~1 ~1 ~1 ~1 3
## 4 ~1 ~1 ~1 ~SHARE ~Per 4
## 5 ~1 ~1 ~time ~SHARE ~Per 5
## 6 ~1 ~1 ~1 ~1 ~Per 6
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing
rm(list=myobs)
# fit the models
model.fits <- plyr::dlply(model.list, "model.number", function(x,input.data, input.ddl){
cat("\n\n***** Starting ", unlist(x), "\n")
#browser()
model.parameters=list(
Psi1 =list(formula=as.formula(eval(x$Psi1 ))),
Psi2 =list(formula=as.formula(eval(x$Psi2))),
p1 =list(formula=as.formula(eval(x$p1 ))),
p2 =list(formula=as.formula(eval(x$p2 ))),
Delta =list(formula=as.formula(eval(x$Delta)))
)
if(x$p2 == '~SHARE'){
model.parameters$p1 =list(formula=as.formula(eval(x$p1 )),share=TRUE)
model.parameters$p2 = NULL
}
fit <- RMark::mark(input.data, ddl=input.ddl,
model="MSOccupancy",
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.data=owl.data, input.ddl=owl.ddl)
##
##
## ***** Starting ~1 ~1 ~1 ~SHARE ~1 1
##
## Output summary for MSOccupancy model
## Name : Psi1(~1)Psi2(~1)p1(~1)p2()Delta(~1)
##
## Npar : 4
## -2lnL: 327.3784
## AICc : 336.1948
##
## Beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.6729404 1.5637454 0.6079994 6.7378813
## Psi2:(Intercept) 0.2180147 0.4551343 -0.6740485 1.1100779
## p1:(Intercept) 1.0710730 0.1910325 0.6966493 1.4454966
## Delta:(Intercept) -0.3105595 0.3455643 -0.9878655 0.3667466
##
##
## Real Parameter Psi1
## 1
## 0.9752276
##
##
## Real Parameter Psi2
## 1
## 0.5542888
##
##
## Real Parameter p1
## 1 2 3 4 5
## 0.7448009 0.7448009 0.7448009 0.7448009 0.7448009
##
##
## Real Parameter p2
## 1 2 3 4 5
## 0.7448009 0.7448009 0.7448009 0.7448009 0.7448009
##
##
## Real Parameter Delta
## 1 2 3 4 5
## 0.4229782 0.4229782 0.4229782 0.4229782 0.4229782
##
##
## ***** Starting ~1 ~1 ~time ~1 ~1 2
##
## Output summary for MSOccupancy model
## Name : Psi1(~1)Psi2(~1)p1(~time)p2(~1)Delta(~1)
##
## Npar : 9
## -2lnL: 320.6246
## AICc : 342.7155
##
## Beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.7912334 1.5661059 0.7216658 6.8608011
## Psi2:(Intercept) 0.1487211 0.4854325 -0.8027267 1.1001689
## p1:(Intercept) 0.0831658 0.5817291 -1.0570233 1.2233549
## p1:time2 1.2954023 0.9417281 -0.5503848 3.1411893
## p1:time3 2.6843873 2.1259017 -1.4823801 6.8511548
## p1:time4 0.8812162 0.8464726 -0.7778702 2.5403026
## p1:time5 0.6262029 0.9067362 -1.1510002 2.4034059
## p2:(Intercept) 1.0822448 0.3410182 0.4138492 1.7506404
## Delta:(Intercept) -0.2806986 0.3398185 -0.9467429 0.3853456
##
##
## Real Parameter Psi1
## 1
## 0.9779303
##
##
## Real Parameter Psi2
## 1
## 0.5371119
##
##
## Real Parameter p1
## 1 2 3 4 5
## 0.5207795 0.7987609 0.9408971 0.7239983 0.6702616
##
##
## Real Parameter p2
## 1 2 3 4 5
## 0.7469186 0.7469186 0.7469186 0.7469186 0.7469186
##
##
## Real Parameter Delta
## 1 2 3 4 5
## 0.4302825 0.4302825 0.4302825 0.4302825 0.4302825
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 3
##
## Output summary for MSOccupancy model
## Name : Psi1(~1)Psi2(~1)p1(~1)p2(~1)Delta(~1)
##
## Npar : 5
## -2lnL: 327.3627
## AICc : 338.6127
##
## Beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.6773359 1.5713052 0.5975777 6.7570941
## Psi2:(Intercept) 0.2503774 0.5367397 -0.8016323 1.3023872
## p1:(Intercept) 1.1154524 0.4126504 0.3066577 1.9242472
## p2:(Intercept) 1.0353686 0.3408050 0.3673909 1.7033463
## Delta:(Intercept) -0.3215614 0.3588996 -1.0250046 0.3818818
##
##
## Real Parameter Psi1
## 1
## 0.9753336
##
##
## Real Parameter Psi2
## 1
## 0.5622694
##
##
## Real Parameter p1
## 1 2 3 4 5
## 0.7531442 0.7531442 0.7531442 0.7531442 0.7531442
##
##
## Real Parameter p2
## 1 2 3 4 5
## 0.7379554 0.7379554 0.7379554 0.7379554 0.7379554
##
##
## Real Parameter Delta
## 1 2 3 4 5
## 0.4202953 0.4202953 0.4202953 0.4202953 0.4202953
##
##
## ***** Starting ~1 ~1 ~1 ~SHARE ~Per 4
##
## Note: only 4 parameters counted of 5 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for MSOccupancy model
## Name : Psi1(~1)Psi2(~1)p1(~1)p2()Delta(~Per)
##
## Npar : 5 (unadjusted=4)
## -2lnL: 278.3492
## AICc : 289.5992 (unadjusted=287.16553)
##
## Beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.6729407 1.5637457 0.6079990 6.7378823
## Psi2:(Intercept) -0.2025604 0.3239793 -0.8375598 0.4324391
## p1:(Intercept) 1.0710730 0.1910325 0.6966494 1.4454967
## Delta:(Intercept) -19.6484380 0.0000000 -19.6484380 -19.6484380
## Delta:PerL 21.5347190 0.0000000 21.5347190 21.5347190
##
##
## Real Parameter Psi1
## 1
## 0.9752276
##
##
## Real Parameter Psi2
## 1
## 0.4495323
##
##
## Real Parameter p1
## 1 2 3 4 5
## 0.7448009 0.7448009 0.7448009 0.7448009 0.7448009
##
##
## Real Parameter p2
## 1 2 3 4 5
## 0.7448009 0.7448009 0.7448009 0.7448009 0.7448009
##
##
## Real Parameter Delta
## 1 2 3 4 5
## 2.929489e-09 2.929489e-09 0.868331 0.868331 0.868331
##
##
## ***** Starting ~1 ~1 ~time ~SHARE ~Per 5
##
## Note: only 8 parameters counted of 9 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for MSOccupancy model
## Name : Psi1(~1)Psi2(~1)p1(~time)p2()Delta(~Per)
##
## Npar : 9 (unadjusted=8)
## -2lnL: 261.355
## AICc : 283.4459 (unadjusted=280.55501)
##
## Beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.8395132 1.5703482 0.7616306 6.9173958
## Psi2:(Intercept) -0.2025609 0.3239784 -0.8375584 0.4324367
## p1:(Intercept) 0.3187289 0.3369260 -0.3416460 0.9791039
## p1:time2 1.2235204 0.5673946 0.1114269 2.3356138
## p1:time3 2.1563260 0.6702794 0.8425784 3.4700736
## p1:time4 0.5143787 0.4948296 -0.4554873 1.4842447
## p1:time5 0.1692565 0.5741595 -0.9560962 1.2946092
## Delta:(Intercept) -22.5559030 0.0000000 -22.5559030 -22.5559030
## Delta:PerL 24.4422220 0.0000000 24.4422220 24.4422220
##
##
## Real Parameter Psi1
## 1
## 0.9789486
##
##
## Real Parameter Psi2
## 1
## 0.4495322
##
##
## Real Parameter p1
## 1 2 3 4 5
## 0.5790144 0.8237915 0.9223745 0.6970116 0.6196317
##
##
## Real Parameter p2
## 1 2 3 4 5
## 0.5790144 0.8237915 0.9223745 0.6970116 0.6196317
##
##
## Real Parameter Delta
## 1 2 3 4 5
## 1.599911e-10 1.599911e-10 0.8683353 0.8683353 0.8683353
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~Per 6
##
## Note: only 5 parameters counted of 6 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for MSOccupancy model
## Name : Psi1(~1)Psi2(~1)p1(~1)p2(~1)Delta(~Per)
##
## Npar : 6 (unadjusted=5)
## -2lnL: 278.24
## AICc : 292.0273 (unadjusted=289.49005)
##
## Beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.6819710 1.5789064 0.5873145 6.7766276
## Psi2:(Intercept) -0.2327187 0.3334419 -0.8862648 0.4208273
## p1:(Intercept) 1.0060685 0.2704961 0.4758961 1.5362410
## p2:(Intercept) 1.1533787 0.3171119 0.5318394 1.7749180
## Delta:(Intercept) -21.3583780 0.0000000 -21.3583780 -21.3583780
## Delta:PerL 23.2665180 0.0000000 23.2665180 23.2665180
##
##
## Real Parameter Psi1
## 1
## 0.9754448
##
##
## Real Parameter Psi2
## 1
## 0.4420815
##
##
## Real Parameter p1
## 1 2 3 4 5
## 0.7322501 0.7322501 0.7322501 0.7322501 0.7322501
##
##
## Real Parameter p2
## 1 2 3 4 5
## 0.7601275 0.7601275 0.7601275 0.7601275 0.7601275
##
##
## Real Parameter Delta
## 1 2 3 4 5
## 5.298759e-10 5.298759e-10 0.87081 0.87081 0.87081
# examine individual model results
model.number <-2
summary(model.fits[[model.number]])
## Output summary for MSOccupancy model
## Name : Psi1(~1)Psi2(~1)p1(~time)p2(~1)Delta(~1)
##
## Npar : 9
## -2lnL: 320.6246
## AICc : 342.7155
##
## Beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.7912334 1.5661059 0.7216658 6.8608011
## Psi2:(Intercept) 0.1487211 0.4854325 -0.8027267 1.1001689
## p1:(Intercept) 0.0831658 0.5817291 -1.0570233 1.2233549
## p1:time2 1.2954023 0.9417281 -0.5503848 3.1411893
## p1:time3 2.6843873 2.1259017 -1.4823801 6.8511548
## p1:time4 0.8812162 0.8464726 -0.7778702 2.5403026
## p1:time5 0.6262029 0.9067362 -1.1510002 2.4034059
## p2:(Intercept) 1.0822448 0.3410182 0.4138492 1.7506404
## Delta:(Intercept) -0.2806986 0.3398185 -0.9467429 0.3853456
##
##
## Real Parameter Psi1
## 1
## 0.9779303
##
##
## Real Parameter Psi2
## 1
## 0.5371119
##
##
## Real Parameter p1
## 1 2 3 4 5
## 0.5207795 0.7987609 0.9408971 0.7239983 0.6702616
##
##
## Real Parameter p2
## 1 2 3 4 5
## 0.7469186 0.7469186 0.7469186 0.7469186 0.7469186
##
##
## Real Parameter Delta
## 1 2 3 4 5
## 0.4302825 0.4302825 0.4302825 0.4302825 0.4302825
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## Psi1 g1 a0 t1 0.9779303 0.0338007 0.6729737 0.9989530
## Psi2 g1 a0 t1 0.5371119 0.1206896 0.3094426 0.7502918
## p1 g1 a0 t1 0.5207795 0.1451811 0.2578787 0.7726534
## p1 g1 a1 t2 0.7987609 0.1207545 0.4765555 0.9453695
## p1 g1 a2 t3 0.9408971 0.1200944 0.1876662 0.9990893
## p1 g1 a3 t4 0.7239983 0.1327817 0.4162905 0.9060885
## p1 g1 a4 t5 0.6702616 0.1574566 0.3346987 0.8914605
## p2 g1 a0 t1 0.7469186 0.0644631 0.6020105 0.8520336
## Delta g1 a0 t1 0.4302825 0.0833029 0.2795403 0.5951618
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## Psi1:(Intercept) 3.7912334 1.5661059 0.7216658 6.8608011
## Psi2:(Intercept) 0.1487211 0.4854325 -0.8027267 1.1001689
## p1:(Intercept) 0.0831658 0.5817291 -1.0570233 1.2233549
## p1:time2 1.2954023 0.9417281 -0.5503848 3.1411893
## p1:time3 2.6843873 2.1259017 -1.4823801 6.8511548
## p1:time4 0.8812162 0.8464726 -0.7778702 2.5403026
## p1:time5 0.6262029 0.9067362 -1.1510002 2.4034059
## p2:(Intercept) 1.0822448 0.3410182 0.4138492 1.7506404
## Delta:(Intercept) -0.2806986 0.3398185 -0.9467429 0.3853456
model.fits[[model.number]]$results$derived
## $`psi1*psi2`
## estimate se lcl ucl
## 1 0.525258 0.1183269 0.3038697 0.7371435
get.real(model.fits[[model.number]], "Psi1", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi1 g1 a0 t1 1 1 0.9779303 0.0338007 0.6729737
## ucl fixed note group age time Age Time
## Psi1 g1 a0 t1 0.998953 1 0 1 0 0
get.real(model.fits[[model.number]], "Psi2", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi2 g1 a0 t1 2 2 0.5371119 0.1206896 0.3094426
## ucl fixed note group age time Age Time
## Psi2 g1 a0 t1 0.7502918 1 0 1 0 0
get.real(model.fits[[model.number]], "p1", se=TRUE)
## all.diff.index par.index estimate se lcl
## p1 g1 a0 t1 3 3 0.5207795 0.1451811 0.2578787
## p1 g1 a1 t2 4 4 0.7987609 0.1207545 0.4765555
## p1 g1 a2 t3 5 5 0.9408971 0.1200944 0.1876662
## p1 g1 a3 t4 6 6 0.7239983 0.1327817 0.4162905
## p1 g1 a4 t5 7 7 0.6702616 0.1574566 0.3346987
## ucl fixed note group age time Age Time
## p1 g1 a0 t1 0.7726534 1 0 1 0 0
## p1 g1 a1 t2 0.9453695 1 1 2 1 1
## p1 g1 a2 t3 0.9990893 1 2 3 2 2
## p1 g1 a3 t4 0.9060885 1 3 4 3 3
## p1 g1 a4 t5 0.8914605 1 4 5 4 4
get.real(model.fits[[model.number]], "p2", se=TRUE)
## all.diff.index par.index estimate se lcl
## p2 g1 a0 t1 8 8 0.7469186 0.0644631 0.6020105
## p2 g1 a1 t2 9 8 0.7469186 0.0644631 0.6020105
## p2 g1 a2 t3 10 8 0.7469186 0.0644631 0.6020105
## p2 g1 a3 t4 11 8 0.7469186 0.0644631 0.6020105
## p2 g1 a4 t5 12 8 0.7469186 0.0644631 0.6020105
## ucl fixed note group age time Age Time
## p2 g1 a0 t1 0.8520336 1 0 1 0 0
## p2 g1 a1 t2 0.8520336 1 1 2 1 1
## p2 g1 a2 t3 0.8520336 1 2 3 2 2
## p2 g1 a3 t4 0.8520336 1 3 4 3 3
## p2 g1 a4 t5 0.8520336 1 4 5 4 4
get.real(model.fits[[model.number]], "Delta", se=TRUE)
## all.diff.index par.index estimate se lcl
## Delta g1 a0 t1 13 9 0.4302825 0.0833029 0.2795403
## Delta g1 a1 t2 14 9 0.4302825 0.0833029 0.2795403
## Delta g1 a2 t3 15 9 0.4302825 0.0833029 0.2795403
## Delta g1 a3 t4 16 9 0.4302825 0.0833029 0.2795403
## Delta g1 a4 t5 17 9 0.4302825 0.0833029 0.2795403
## ucl fixed note group age time Age Time
## Delta g1 a0 t1 0.5951618 1 0 1 0 0
## Delta g1 a1 t2 0.5951618 1 1 2 1 1
## Delta g1 a2 t3 0.5951618 1 2 3 2 2
## Delta g1 a3 t4 0.5951618 1 3 4 3 3
## Delta g1 a4 t5 0.5951618 1 4 5 4 4
# collect models and make AICc table
model.set <- RMark::collect.models( type="MSOccupancy")
model.set
## model npar AICc DeltaAICc
## 5 Psi1(~1)Psi2(~1)p1(~time)p2()Delta(~Per) 9 283.4459 0.000000
## 4 Psi1(~1)Psi2(~1)p1(~1)p2()Delta(~Per) 5 289.5992 6.153281
## 6 Psi1(~1)Psi2(~1)p1(~1)p2(~1)Delta(~Per) 6 292.0273 8.581365
## 1 Psi1(~1)Psi2(~1)p1(~1)p2()Delta(~1) 4 336.1948 52.748857
## 3 Psi1(~1)Psi2(~1)p1(~1)p2(~1)Delta(~1) 5 338.6127 55.166751
## 2 Psi1(~1)Psi2(~1)p1(~time)p2(~1)Delta(~1) 9 342.7155 59.269570
## weight Deviance
## 5 9.435658e-01 -126.45304
## 4 4.351152e-02 -109.45885
## 6 1.292268e-02 -109.56800
## 1 3.315151e-12 -60.42960
## 3 9.896108e-13 -60.44538
## 2 1.272180e-13 -67.18347
names(model.set)
## [1] "m...001" "m...002" "m...003" "m...004" "m...005"
## [6] "m...006" "model.table"
model.set$model.table
## Psi1 Psi2 p1 p2 Delta model npar
## 5 ~1 ~1 ~time ~Per Psi1(~1)Psi2(~1)p1(~time)p2()Delta(~Per) 9
## 4 ~1 ~1 ~1 ~Per Psi1(~1)Psi2(~1)p1(~1)p2()Delta(~Per) 5
## 6 ~1 ~1 ~1 ~1 ~Per Psi1(~1)Psi2(~1)p1(~1)p2(~1)Delta(~Per) 6
## 1 ~1 ~1 ~1 ~1 Psi1(~1)Psi2(~1)p1(~1)p2()Delta(~1) 4
## 3 ~1 ~1 ~1 ~1 ~1 Psi1(~1)Psi2(~1)p1(~1)p2(~1)Delta(~1) 5
## 2 ~1 ~1 ~time ~1 ~1 Psi1(~1)Psi2(~1)p1(~time)p2(~1)Delta(~1) 9
## AICc DeltaAICc weight Deviance
## 5 283.4459 0.000000 9.435658e-01 -126.45304
## 4 289.5992 6.153281 4.351152e-02 -109.45885
## 6 292.0273 8.581365 1.292268e-02 -109.56800
## 1 336.1948 52.748857 3.315151e-12 -60.42960
## 3 338.6127 55.166751 9.896108e-13 -60.44538
## 2 342.7155 59.269570 1.272180e-13 -67.18347
# model averaged values at average covariate value
Psi1.ma <- RMark::model.average(model.set, param="Psi1")
## 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.
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
Psi1.ma
## par.index estimate se fixed note group age time
## Psi1 g1 a0 t1 1 0.9787414 0.03270321 1 0 1
## Age Time
## Psi1 g1 a0 t1 0 0
Psi2.ma <- RMark::model.average(model.set, param="Psi2")
## 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.
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
Psi2.ma
## par.index estimate se fixed note group age time
## Psi2 g1 a0 t1 2 0.4494359 0.08020097 1 0 1
## Age Time
## Psi2 g1 a0 t1 0 0
# model average derived parameters
source(file.path("..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
RMark.model.average.derived(model.set, "psi1*psi2") # Notice that lowercase Psi1 and lowercase Psi2
## estimate se lcl ucl
## 1 0.4398819 0.07985599 0.2833671 0.5963968
cleanup(ask=FALSE)