# Multi species single season
# Co-occurence of Bull and Brook trout
# Data provided by Parks Canada
# 2018-12-27 code submitted by Neil Faught
library(car)
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
# get the data
occ.data <- readxl::read_excel("../Trout2.xlsx",
sheet="fish_sp_stacked_juv_delta")
head(occ.data)
## # A tibble: 6 x 7
## site watershed species presence rep temp discharge
## <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 1.00 W3 BKTR 1 1.00 10.2 8.93
## 2 1.00 W3 BKTR 1 2.00 10.2 8.93
## 3 1.00 W3 BLTR 1 1.00 10.2 8.93
## 4 1.00 W3 BLTR 1 2.00 10.2 8.93
## 5 2.00 W3 BKTR 1 1.00 9.41 10.4
## 6 2.00 W3 BKTR 1 2.00 9.41 10.4
# check site code
xtabs(~site+rep, data=occ.data, exclude=NULL, na.action=na.pass)
## rep
## site 1 2
## 1 2 2
## 2 2 2
## 3 4 4
## 4 2 2
## 5 2 2
## 6 2 2
## 7 6 6
## 8 2 2
## 9 4 4
## 10 2 2
## 11 4 4
## 12 2 2
## 13 2 2
## 14 2 2
## 15 2 2
## 16 4 4
## 19 4 4
## 21 2 2
## 22 2 2
## 24 2 2
## 25 2 2
## 27 2 2
## 35 2 2
## 39 2 2
## 40 2 2
## 41 2 2
## 42 2 2
## 47 2 2
## 51 4 4
## 54 4 4
## 55 4 4
## 58 2 2
## 60 2 2
## 63 2 2
## 64 2 2
## 66 2 2
## 67 4 4
## 70 2 2
## 72 2 2
## 73 2 2
## 74 2 2
## 77 2 2
## 82 2 2
## 83 2 2
## 84 2 2
## 89 2 2
## 94 2 2
## 100 2 2
## 103 2 2
## 107 2 2
## 108 2 2
## 109 2 2
## 112 2 2
## 113 2 2
## 114 2 2
## 115 2 2
## 116 2 2
## 117 4 4
## 118 2 2
## 119 2 2
## 120 4 4
## 121 2 2
## 122 2 2
## 123 2 2
## 124 2 2
## 125 2 2
## 127 4 4
## 128 2 2
## 129 2 2
## 130 2 2
## 131 2 2
## 133 4 4
## 136 2 2
## 137 2 2
## 139 2 2
## 141 2 2
## 143 4 4
## 144 2 2
## 145 2 2
## 147 2 2
## 148 4 4
## 149 2 2
## 150 2 2
## 152 2 2
## 154 4 4
## 155 4 4
## 157 2 2
## 158 2 2
## 161 2 2
## 163 2 2
## 165 2 2
## 166 2 2
## 168 2 2
## 169 2 2
## 170 4 4
## 173 2 2
## 175 2 2
## 176 2 2
## 188 2 2
## 189 2 2
## 191 2 2
## 198 2 2
## 199 2 2
## 200 2 2
## 201 2 2
## 202 2 2
## 203 2 2
## 204 2 2
## 205 2 2
## 206 2 2
## 207 2 2
## 209 2 2
## 210 2 2
## 211 2 2
## 212 2 2
## 213 2 2
## 215 2 2
## 216 2 2
## 217 2 2
## 218 2 2
## 220 2 2
## 221 2 2
## 222 2 2
## 224 2 2
## 225 2 2
## 226 2 2
## 227 2 2
## 229 2 2
## 230 2 2
## 231 2 2
## 234 2 2
## 236 2 2
## 237 2 2
## 238 2 2
## 240 2 2
## 242 2 2
## 243 2 2
## 245 2 2
## 246 2 2
## 249 2 2
## 253 2 2
## 422 2 2
## 453 2 2
## 470 2 2
## 534 2 2
## 598 2 2
## 666 2 2
## 789 2 2
## 805 2 2
## 821 2 2
## 825 2 2
## 837 2 2
## 853 2 2
## 869 2 2
## 873 2 2
## 885 2 2
## 889 2 2
## 901 2 2
## 917 2 2
## 933 2 2
## 937 2 2
## 949 2 2
## 985 2 2
# some site numbers are "reused", but a combination of watershed x site is unique
nmeasure <- plyr::ddply(occ.data, c("watershed","site"), plyr::summarize, count=length(rep))
nmeasure[ nmeasure$count != 4,]
## [1] watershed site count
## <0 rows> (or 0-length row.names)
xtabs(~presence, data=occ.data, exclude=NULL, na.action=na.pass)
## presence
## 0 1 NA
## 432 298 2
occ.data$presence <- as.numeric(occ.data$presence)
## Warning: NAs introduced by coercion
# create the detection records for each rep
# Species A = BKTR; Species B=BLTR
# it is believed that BKTR have a preference for warmer water (see later)
# so we code them as the "dominant" species (Species A) so that
# a linear logisitic regression makes sense for it.
dhist1 <- plyr::ddply(occ.data, c("watershed","site","rep"), plyr::summarize,
history=sum(presence*1*(species=="BKTR")+
presence*2*(species=="BLTR")),
temp = mean(temp),
discharge=mean(discharge))
history <- reshape2::dcast(dhist1, watershed+site+temp+discharge~rep,
variable.name="rep",
value.var="history")
head(history)
## watershed site temp discharge 1 2
## 1 W1 3 4.743874 0.27503 0 0
## 2 W1 7 4.743874 0.47006 0 0
## 3 W1 9 6.884988 1.02842 0 0
## 4 W1 19 4.743874 0.60002 0 0
## 5 W1 21 7.848305 0.25631 2 2
## 6 W1 24 5.364172 0.44846 0 0
history <- plyr::rename(history, c("1"="r1", "2"="r2"))
xtabs(~r1+r2, data=history, exclude=NULL, na.action=na.pass)
## r2
## r1 0 1 2 3 <NA>
## 0 62 1 4 0 0
## 1 0 14 1 4 0
## 2 6 0 53 1 0
## 3 1 4 4 27 1
input.data <- history[,c("r1","r2")]
input.data[1:5,]
## r1 r2
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 2 2
input.chistory <- data.frame(lapply(input.data, as.character), stringsAsFactors=FALSE)
input.chistory <- data.frame(lapply(input.chistory,
car::recode, " '0'='00'; '1'='10';'2'='01';'3'='11';", as.numeric=FALSE, as.factor=FALSE),
stringsAsFactors=FALSE)
input.history <- data.frame(ch=apply(input.chistory, 1, paste, sep="",collapse=""), freq=1)
# Change any NA to . in the chapter history
select <- grepl("NA", input.history$ch)
input.history[ select,]
## ch freq
## 163 11NA 1
input.history$ch <- gsub("NA","..", input.history$ch, fixed=TRUE)
input.history[ select,]
## ch freq
## 163 11.. 1
site.covar <- history[,c("watershed","temp","discharge")]
names(site.covar)
## [1] "watershed" "temp" "discharge"
head(site.covar)
## watershed temp discharge
## 1 W1 4.743874 0.27503
## 2 W1 4.743874 0.47006
## 3 W1 6.884988 1.02842
## 4 W1 4.743874 0.60002
## 5 W1 7.848305 0.25631
## 6 W1 5.364172 0.44846
temp.mean <- mean(site.covar$temp)
temp.sd <- sd (site.covar$temp)
site.covar$temp.std <- (site.covar$temp - temp.mean)/ temp.sd
discharge.mean <- mean(site.covar$discharge)
discharge.sd <- sd (site.covar$discharge)
site.covar$discharge.std <- (site.covar$discharge - discharge.mean)/ discharge.sd
# Combine site covariates with input history
input.history = cbind(input.history, site.covar)
head(input.history)
## ch freq watershed temp discharge temp.std discharge.std
## 1 0000 1 W1 4.743874 0.27503 -1.12739888 -0.8011166
## 2 0000 1 W1 4.743874 0.47006 -1.12739888 -0.7435560
## 3 0000 1 W1 6.884988 1.02842 -0.04090698 -0.5787631
## 4 0000 1 W1 4.743874 0.60002 -1.12739888 -0.7051999
## 5 0101 1 W1 7.848305 0.25631 0.44792090 -0.8066415
## 6 0000 1 W1 5.364172 0.44846 -0.81263325 -0.7499309
# Creat RMark object
trout.data <- process.data(data=input.history,
model="2SpecConOccup") # this parameterization is more stable
summary(trout.data)
## Length Class Mode
## data 7 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 183 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 2 -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
# if there are time-specific covariates add them to the ddl's
trout.ddl <- RMark::make.design.data(trout.data) # need for covariate predictions
# What are the parameter names for Single Season Single Species models
setup.parameters("2SpecConOccup", check=TRUE)
## [1] "PsiA" "PsiBA" "PsiBa" "pA" "pB" "rA" "rBA" "rBa"
# Use the psiAB parameterization
# Parameters of the model are
# psiA - occupancy probability of species A
# psiBA - occupancy probability of species B if species A is present
# psiBa - occupancy probability of species B if species A is absent
#
# If species are independent thatn psiBA = psiBa.
# Alternatively, nu = odds(B|A)/odds(B|a) = 1.
#
# Detection parameters
# pA - probability of detection if A is alone in the site
# pB - probability of detection if B is alone in the site
# rA - probability of detecting A only given both on the site
# rBA - probability of detecting B given that A was detected and both are on the site
# rBa - probability of detecting B given that A not detected and both are on the site
# Ifspecies do not interact, then
# rBA = rBa
# and
# rho = odds(rBA)/odds(rBa) = 1
setup.parameters("2SpecConOccup", check=TRUE)
## [1] "PsiA" "PsiBA" "PsiBa" "pA" "pB" "rA" "rBA" "rBa"
# Get the list of models. NOtice NO equal signs here
model.list.csv <- textConnection("
PsiA, PsiBA, PsiBa, pA, pB, rA, rBA, rBa
~1, ~1, ~1, ~1, ~1, ~1, ~1, ~1
~1, ~1, ~1, ~1, ~SHARE, ~1, ~1, ~SHARE
~1, ~1, ~SHARE, ~1, ~1, ~1, ~1, ~1
~temp.std, ~temp.std, ~temp.std, ~1, ~SHARE, ~1, ~1, ~SHARE")
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
## PsiA PsiBA PsiBa pA pB rA rBA rBa model.number
## 1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 1
## 2 ~1 ~1 ~1 ~1 ~SHARE ~1 ~1 ~SHARE 2
## 3 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 3
## 4 ~temp.std ~temp.std ~temp.std ~1 ~SHARE ~1 ~1 ~SHARE 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.data, input.ddl){
cat("\n\n***** Starting ", unlist(x), "\n")
#browser()
model.parameters=list(
PsiA =list(formula=as.formula(eval(x$PsiA ))),
PsiBA =list(formula=as.formula(eval(x$PsiBA))),
PsiBa =list(formula=as.formula(eval(x$PsiBa))),
pA =list(formula=as.formula(eval(x$pA ))),
pB =list(formula=as.formula(eval(x$pB ))),
rA =list(formula=as.formula(eval(x$rA ))),
rBA =list(formula=as.formula(eval(x$rBA))),
rBa =list(formula=as.formula(eval(x$rBa)))
)
if(x$rBa == '~SHARE'){
model.parameters$rBA =list(formula=as.formula(eval(x$rBA)), share=TRUE)
model.parameters$rBa = NULL
}
if(x$PsiBa == '~SHARE'){
model.parameters$PsiBA =list(formula=as.formula(eval(x$PsiBA)), share=TRUE)
model.parameters$PsiBa = NULL
}
if(x$pB == '~SHARE'){
model.parameters$pA =list(formula=as.formula(eval(x$pA)), share=TRUE)
model.parameters$pB = NULL
}
fit <- RMark::mark(input.data, ddl=input.ddl,
model="2SpecConOccup",
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=trout.data,input.ddl=trout.ddl)
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 1
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
##
## Npar : 8
## -2lnL: 653.798
## AICc : 670.6256
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) -0.7584873 0.1594990 -1.0711053 -0.4458693
## PsiBA:(Intercept) 1.1318242 0.3199195 0.5047821 1.7588664
## PsiBa:(Intercept) 0.0257210 0.1810942 -0.3292236 0.3806657
## pA:(Intercept) 3.4524326 1.2273624 1.0468022 5.8580630
## pB:(Intercept) 2.3658321 0.3480175 1.6837177 3.0479465
## rA:(Intercept) 2.3114218 0.4143267 1.4993416 3.1235021
## rBA:(Intercept) 1.8914804 0.3784509 1.1497166 2.6332441
## rBa:(Intercept) 1.6493564 1.0812927 -0.4699774 3.7686901
##
##
## Real Parameter PsiA
## 1
## 0.3189748
##
##
## Real Parameter PsiBA
## 1
## 0.7561754
##
##
## Real Parameter PsiBa
## 1
## 0.5064299
##
##
## Real Parameter pA
## 1 2
## 0.9693036 0.9693036
##
##
## Real Parameter pB
## 1 2
## 0.9141844 0.9141844
##
##
## Real Parameter rA
## 1 2
## 0.9098186 0.9098186
##
##
## Real Parameter rBA
## 1 2
## 0.8689242 0.8689242
##
##
## Real Parameter rBa
## 1 2
## 0.838804 0.838804
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~SHARE ~1 ~1 ~SHARE 2
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB()rA(~1)rBA(~1)rBa()
##
## Npar : 6
## -2lnL: 654.8385
## AICc : 667.3158
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) -0.7568657 0.1595781 -1.0696388 -0.4440926
## PsiBA:(Intercept) 1.1257237 0.3193645 0.4997692 1.7516782
## PsiBa:(Intercept) 0.0235540 0.1808124 -0.3308382 0.3779463
## pA:(Intercept) 2.4991782 0.3298836 1.8526063 3.1457501
## rA:(Intercept) 2.3254034 0.4141395 1.5136900 3.1371168
## rBA:(Intercept) 1.8711050 0.3615426 1.1624816 2.5797285
##
##
## Real Parameter PsiA
## 1
## 0.3193271
##
##
## Real Parameter PsiBA
## 1
## 0.7550489
##
##
## Real Parameter PsiBa
## 1
## 0.5058882
##
##
## Real Parameter pA
## 1 2
## 0.9240842 0.9240842
##
##
## Real Parameter pB
## 1 2
## 0.9240842 0.9240842
##
##
## Real Parameter rA
## 1 2
## 0.9109592 0.9109592
##
##
## Real Parameter rBA
## 1 2
## 0.8665861 0.8665861
##
##
## Real Parameter rBa
## 1 2
## 0.8665861 0.8665861
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 3
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
##
## Npar : 7
## -2lnL: 663.8794
## AICc : 678.5194
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) -0.7613256 0.1592714 -1.0734977 -0.4491536
## PsiBA:(Intercept) 0.3422651 0.1521158 0.0441182 0.6404120
## pA:(Intercept) 3.3385660 1.1198843 1.1435927 5.5335393
## pB:(Intercept) 2.3240805 0.3540616 1.6301197 3.0180414
## rA:(Intercept) 2.3477257 0.4065485 1.5508906 3.1445609
## rBA:(Intercept) 1.9753220 0.3632630 1.2633266 2.6873174
## rBa:(Intercept) 1.7333484 1.0822363 -0.3878349 3.8545317
##
##
## Real Parameter PsiA
## 1
## 0.3183585
##
##
## Real Parameter PsiBA
## 1
## 0.5847406
##
##
## Real Parameter PsiBa
## 1
## 0.5847406
##
##
## Real Parameter pA
## 1 2
## 0.9657284 0.9657284
##
##
## Real Parameter pB
## 1 2
## 0.9108518 0.9108518
##
##
## Real Parameter rA
## 1 2
## 0.9127533 0.9127533
##
##
## Real Parameter rBA
## 1 2
## 0.8781816 0.8781816
##
##
## Real Parameter rBa
## 1 2
## 0.8498402 0.8498402
##
##
## ***** Starting ~temp.std ~temp.std ~temp.std ~1 ~SHARE ~1 ~1 ~SHARE 4
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~temp.std)PsiBA(~temp.std)PsiBa(~temp.std)pA(~1)pB()rA(~1)rBA(~1)rBa()
##
## Npar : 9
## -2lnL: 521.3936
## AICc : 540.434
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) -1.1535590 0.2448024 -1.6333717 -0.6737463
## PsiA:temp.std 2.5263631 0.4057293 1.7311335 3.3215926
## PsiBA:(Intercept) 3.8249070 1.1487507 1.5733556 6.0764585
## PsiBA:temp.std -2.0174878 0.6754506 -3.3413709 -0.6936047
## PsiBa:(Intercept) 0.4976977 0.2753907 -0.0420682 1.0374635
## PsiBa:temp.std 0.9257765 0.3689921 0.2025519 1.6490011
## pA:(Intercept) 2.4436669 0.3379554 1.7812743 3.1060595
## rA:(Intercept) 2.3495746 0.4072935 1.5512793 3.1478699
## rBA:(Intercept) 1.9541310 0.3471786 1.2736609 2.6346011
##
##
## Real Parameter PsiA
## 1
## 0.2398396
##
##
## Real Parameter PsiBA
## 1
## 0.9786455
##
##
## Real Parameter PsiBa
## 1
## 0.6219181
##
##
## Real Parameter pA
## 1 2
## 0.9200971 0.9200971
##
##
## Real Parameter pB
## 1 2
## 0.9200971 0.9200971
##
##
## Real Parameter rA
## 1 2
## 0.9129004 0.9129004
##
##
## Real Parameter rBA
## 1 2
## 0.8758964 0.8758964
##
##
## Real Parameter rBa
## 1 2
## 0.8758964 0.8758964
# examine individual model results
model.number <-4
summary(model.fits[[model.number]])
## Output summary for 2SpecConOccup model
## Name : PsiA(~temp.std)PsiBA(~temp.std)PsiBa(~temp.std)pA(~1)pB()rA(~1)rBA(~1)rBa()
##
## Npar : 9
## -2lnL: 521.3936
## AICc : 540.434
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) -1.1535590 0.2448024 -1.6333717 -0.6737463
## PsiA:temp.std 2.5263631 0.4057293 1.7311335 3.3215926
## PsiBA:(Intercept) 3.8249070 1.1487507 1.5733556 6.0764585
## PsiBA:temp.std -2.0174878 0.6754506 -3.3413709 -0.6936047
## PsiBa:(Intercept) 0.4976977 0.2753907 -0.0420682 1.0374635
## PsiBa:temp.std 0.9257765 0.3689921 0.2025519 1.6490011
## pA:(Intercept) 2.4436669 0.3379554 1.7812743 3.1060595
## rA:(Intercept) 2.3495746 0.4072935 1.5512793 3.1478699
## rBA:(Intercept) 1.9541310 0.3471786 1.2736609 2.6346011
##
##
## Real Parameter PsiA
## 1
## 0.2398396
##
##
## Real Parameter PsiBA
## 1
## 0.9786455
##
##
## Real Parameter PsiBa
## 1
## 0.6219181
##
##
## Real Parameter pA
## 1 2
## 0.9200971 0.9200971
##
##
## Real Parameter pB
## 1 2
## 0.9200971 0.9200971
##
##
## Real Parameter rA
## 1 2
## 0.9129004 0.9129004
##
##
## Real Parameter rBA
## 1 2
## 0.8758964 0.8758964
##
##
## Real Parameter rBa
## 1 2
## 0.8758964 0.8758964
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.2398396 0.0446315 0.1633690 0.3376585
## PsiBA g1 a0 t1 0.9786455 0.0240071 0.8282615 0.9977090
## PsiBa g1 a0 t1 0.6219181 0.0647543 0.4894845 0.7383603
## pA g1 a0 t1 0.9200971 0.0248460 0.8558542 0.9571420
## rA g1 a0 t1 0.9129004 0.0323852 0.8250984 0.9588247
## rBA g1 a0 t1 0.8758964 0.0377390 0.7813688 0.9330555
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## PsiA:(Intercept) -1.1535590 0.2448024 -1.6333717 -0.6737463
## PsiA:temp.std 2.5263631 0.4057293 1.7311335 3.3215926
## PsiBA:(Intercept) 3.8249070 1.1487507 1.5733556 6.0764585
## PsiBA:temp.std -2.0174878 0.6754506 -3.3413709 -0.6936047
## PsiBa:(Intercept) 0.4976977 0.2753907 -0.0420682 1.0374635
## PsiBa:temp.std 0.9257765 0.3689921 0.2025519 1.6490011
## pA:(Intercept) 2.4436669 0.3379554 1.7812743 3.1060595
## rA:(Intercept) 2.3495746 0.4072935 1.5512793 3.1478699
## rBA:(Intercept) 1.9541310 0.3471786 1.2736609 2.6346011
model.fits[[model.number]]$results$derived
## $`Species Interaction Factor`
## estimate se lcl ucl
## 1 1.383292 0.1048245 1.177836 1.588749
##
## $`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.7074755 0.0521949 0.5960015 0.7985855
##
## $`Occupancy of Both Species`
## estimate se lcl ucl
## 1 0.234718 0.04403106 0.1594523 0.3314995
get.real(model.fits[[model.number]], "PsiA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiA g1 a0 t1 1 1 0.2398396 0.0446315 0.163369
## ucl fixed note group age time Age Time
## PsiA g1 a0 t1 0.3376585 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBA g1 a0 t1 2 2 0.9786455 0.0240071 0.8282615
## ucl fixed note group age time Age Time
## PsiBA g1 a0 t1 0.997709 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBa g1 a0 t1 3 3 0.6219181 0.0647543 0.4894845
## ucl fixed note group age time Age Time
## PsiBa g1 a0 t1 0.7383603 1 0 1 0 0
get.real(model.fits[[model.number]], "pA", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## pA g1 a0 t1 4 4 0.9200971 0.024846 0.8558542 0.957142
## pA g1 a1 t2 5 4 0.9200971 0.024846 0.8558542 0.957142
## fixed note group age time Age Time
## pA g1 a0 t1 1 0 1 0 0
## pA g1 a1 t2 1 1 2 1 1
get.real(model.fits[[model.number]], "pB", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## pB g1 a0 t1 6 4 0.9200971 0.024846 0.8558542 0.957142
## pB g1 a1 t2 7 4 0.9200971 0.024846 0.8558542 0.957142
## fixed note group age time Age Time
## pB g1 a0 t1 1 0 1 0 0
## pB g1 a1 t2 1 1 2 1 1
get.real(model.fits[[model.number]], "rA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rA g1 a0 t1 8 5 0.9129004 0.0323852 0.8250984
## rA g1 a1 t2 9 5 0.9129004 0.0323852 0.8250984
## ucl fixed note group age time Age Time
## rA g1 a0 t1 0.9588247 1 0 1 0 0
## rA g1 a1 t2 0.9588247 1 1 2 1 1
get.real(model.fits[[model.number]], "rBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBA g1 a0 t1 10 6 0.8758964 0.037739 0.7813688
## rBA g1 a1 t2 11 6 0.8758964 0.037739 0.7813688
## ucl fixed note group age time Age Time
## rBA g1 a0 t1 0.9330555 1 0 1 0 0
## rBA g1 a1 t2 0.9330555 1 1 2 1 1
get.real(model.fits[[model.number]], "rBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBa g1 a0 t1 12 6 0.8758964 0.037739 0.7813688
## rBa g1 a1 t2 13 6 0.8758964 0.037739 0.7813688
## ucl fixed note group age time Age Time
## rBa g1 a0 t1 0.9330555 1 0 1 0 0
## rBa g1 a1 t2 0.9330555 1 1 2 1 1
# collect models and make AICc table
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model
## 4 PsiA(~temp.std)PsiBA(~temp.std)PsiBa(~temp.std)pA(~1)pB()rA(~1)rBA(~1)rBa()
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB()rA(~1)rBA(~1)rBa()
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
## 3 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
## npar AICc DeltaAICc weight Deviance
## 4 9 540.4340 0.0000 1 521.3936
## 2 6 667.3158 126.8818 0 654.8385
## 1 8 670.6256 130.1915 0 653.7980
## 3 7 678.5194 138.0854 0 663.8794
# No need for model averaging as 1 model contains 100% of weight
# cleanup
cleanup(ask=FALSE)