# Multi species single season with a list of models
# Using RMark
# Co-occurence of Jordan's salamander (Plethodon jordani) (PJ) and
# members of Plethodon glutinosus (PG) in Great Smokey
# Mountains National Park (MacKenzie et al. 2004).
# 2018-09-27 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library(car)
## Loading required package: carData
library(ggplot2)
library(readxl)
library(reshape2)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
# get the data
PG.data <- readxl::read_excel(file.path("..","Salamander2.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "B3:F90")
PJ.data <- readxl::read_excel(file.path("..","Salamander2.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "H3:L90")
# MARK wants a 2 "digit" code xy for each visit where
# x = 0/1 if species A is not detected/detected
# and y = 0/1 if species B is not detected/detected
input.data <- PG.data + 2*PJ.data
input.data
## X__1 X__2 X__3 X__4 X__5
## 1 1 0 0 0 0
## 2 0 1 0 0 0
## 3 0 0 0 0 0
## 4 1 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 1 1 1
## 7 0 0 0 0 0
## 8 0 0 0 0 0
## 9 0 0 0 0 0
## 10 1 1 0 0 0
## 11 0 0 0 0 0
## 12 0 0 1 1 0
## 13 0 0 1 0 1
## 14 0 0 1 0 1
## 15 1 1 0 1 0
## 16 1 1 1 0 1
## 17 0 1 0 1 0
## 18 1 1 1 0 1
## 19 0 0 0 0 0
## 20 0 1 0 0 0
## 21 NA NA NA 1 0
## 22 0 0 0 1 1
## 23 0 1 1 1 1
## 24 0 0 0 0 0
## 25 1 1 1 3 3
## 26 3 1 1 1 1
## 27 1 1 1 1 1
## 28 0 0 0 0 0
## 29 0 0 0 0 0
## 30 0 0 0 0 0
## 31 0 1 1 0 1
## 32 2 0 2 0 0
## 33 1 1 0 0 0
## 34 0 1 0 0 0
## 35 1 1 1 1 1
## 36 NA 1 1 1 3
## 37 NA 2 2 2 2
## 38 0 0 0 0 0
## 39 1 1 1 3 1
## 40 1 1 1 1 1
## 41 1 1 1 1 1
## 42 0 1 0 0 1
## 43 1 0 1 0 0
## 44 1 1 0 1 0
## 45 1 1 1 1 1
## 46 1 0 1 0 2
## 47 0 1 0 0 0
## 48 NA 2 2 2 2
## 49 2 2 3 2 0
## 50 1 0 0 0 0
## 51 NA 1 1 1 1
## 52 NA 1 1 0 1
## 53 2 2 2 2 2
## 54 1 3 3 3 3
## 55 0 0 1 1 0
## 56 2 0 2 2 2
## 57 NA 2 0 2 2
## 58 2 2 2 2 2
## 59 2 2 2 2 2
## 60 NA 0 0 2 0
## 61 3 2 2 2 2
## 62 NA 2 0 0 2
## 63 NA 2 2 0 2
## 64 0 0 1 2 3
## 65 0 0 2 3 2
## 66 2 2 2 3 2
## 67 2 2 2 0 2
## 68 2 2 2 2 2
## 69 0 1 1 1 1
## 70 0 0 2 2 3
## 71 2 2 2 2 2
## 72 2 2 2 2 2
## 73 2 2 2 2 2
## 74 2 0 2 2 2
## 75 NA 2 2 2 2
## 76 0 3 3 2 0
## 77 2 2 2 2 2
## 78 NA 2 2 2 2
## 79 NA 2 2 2 2
## 80 2 3 3 2 2
## 81 0 2 2 2 2
## 82 2 2 2 2 2
## 83 2 2 2 2 2
## 84 2 2 2 2 2
## 85 0 1 1 1 0
## 86 2 2 2 2 2
## 87 0 0 2 2 2
## 88 2 2 2 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
## 21 NANANA1000 1
## 36 NA10101011 1
## 37 NA01010101 1
## 48 NA01010101 1
## 51 NA10101010 1
## 52 NA10100010 1
## 57 NA01000101 1
## 60 NA00000100 1
## 62 NA01000001 1
## 63 NA01010001 1
## 75 NA01010101 1
## 78 NA01010101 1
## 79 NA01010101 1
input.history$ch <- gsub("NA","..", input.history$ch, fixed=TRUE)
input.history[ select,]
## ch freq
## 21 ......1000 1
## 36 ..10101011 1
## 37 ..01010101 1
## 48 ..01010101 1
## 51 ..10101010 1
## 52 ..10100010 1
## 57 ..01000101 1
## 60 ..00000100 1
## 62 ..01000001 1
## 63 ..01010001 1
## 75 ..01010101 1
## 78 ..01010101 1
## 79 ..01010101 1
# Any site covariates?
site.covar <- readxl::read_excel(file.path("..","Salamander2.xls"),
sheet="RawData", col_names=TRUE,
range = "N2:O90")
names(site.covar)
## [1] "Elevation (m)" "Std. Elevation"
names(site.covar) <- make.names(names(site.covar))
names(site.covar)
## [1] "Elevation..m." "Std..Elevation"
# Standardize Elevation covariates
elevation.mean <- mean(site.covar$Elevation..m.)
elevation.std <- sd (site.covar$Elevation..m.)
site.covar$El <- (site.covar$Elevation..m. - elevation.mean)/elevation.std
# add the site covariate to the chapture hisotry
input.history <- cbind(input.history, site.covar)
head(input.history)
## ch freq Elevation..m. Std..Elevation El
## 1 1000000000 1 502.92 -0.34208 -2.148793
## 2 0010000000 1 518.16 -0.32684 -2.052179
## 3 0000000000 1 518.16 -0.32684 -2.052179
## 4 1000000000 1 518.16 -0.32684 -2.052179
## 5 0000000000 1 533.40 -0.31160 -1.955564
## 6 0000101010 1 548.64 -0.29636 -1.858950
Nsites <- nrow(input.history)
Nvisits<- ncol(input.history)
sal.data <- process.data(data=input.history,
model="2SpecConOccup") # this parameterization is more stable
summary(sal.data)
## Length Class Mode
## data 5 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 88 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 5 -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
sal.ddl <- RMark::make.design.data(sal.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, ~1, ~1, ~1, ~SHARE
~1, ~1, ~SHARE, ~1, ~1, ~1, ~1, ~1
~1, ~1, ~SHARE, ~1, ~1, ~1, ~1, ~SHARE
~El, ~El, ~El, ~1, ~1, ~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 ~1 ~1 ~1 ~SHARE 2
## 3 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 3
## 4 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 4
## 5 ~El ~El ~El ~1 ~1 ~1 ~1 ~SHARE 5
# 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
}
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=sal.data,input.ddl=sal.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: 735.8665
## AICc : 753.6893
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.2896860 0.2337644 -0.1684922 0.7478642
## PsiBA:(Intercept) -0.7644291 0.3305658 -1.4123382 -0.1165200
## PsiBa:(Intercept) 0.8494996 0.3707647 0.1228009 1.5761984
## pA:(Intercept) 0.1586284 0.1676121 -0.1698914 0.4871481
## pB:(Intercept) 2.2273183 0.4451948 1.3547365 3.0999001
## rA:(Intercept) -0.0473871 0.3067518 -0.6486207 0.5538464
## rBA:(Intercept) 0.0054705 0.3389494 -0.6588705 0.6698114
## rBa:(Intercept) 0.4332825 0.3854088 -0.3221187 1.1886836
##
##
## Real Parameter PsiA
## 1
## 0.5719193
##
##
## Real Parameter PsiBA
## 1
## 0.3176854
##
##
## Real Parameter PsiBa
## 1
## 0.7004622
##
##
## Real Parameter pA
## 1 2 3 4 5
## 0.5395741 0.5395741 0.5395741 0.5395741 0.5395741
##
##
## Real Parameter pB
## 1 2 3 4 5
## 0.902676 0.902676 0.902676 0.902676 0.902676
##
##
## Real Parameter rA
## 1 2 3 4 5
## 0.4881554 0.4881554 0.4881554 0.4881554 0.4881554
##
##
## Real Parameter rBA
## 1 2 3 4 5
## 0.5013676 0.5013676 0.5013676 0.5013676 0.5013676
##
##
## Real Parameter rBa
## 1 2 3 4 5
## 0.6066572 0.6066572 0.6066572 0.6066572 0.6066572
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~SHARE 2
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 7
## -2lnL: 736.624
## AICc : 752.024
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.3039969 0.2325338 -0.1517693 0.7597632
## PsiBA:(Intercept) -0.7428440 0.3261672 -1.3821317 -0.1035563
## PsiBa:(Intercept) 0.8374771 0.3707038 0.1108976 1.5640566
## pA:(Intercept) 0.1628002 0.1675176 -0.1655343 0.4911346
## pB:(Intercept) 2.3033987 0.4512448 1.4189589 3.1878385
## rA:(Intercept) -0.0895065 0.2925274 -0.6628602 0.4838471
## rBA:(Intercept) 0.2021934 0.2504401 -0.2886691 0.6930560
##
##
## Real Parameter PsiA
## 1
## 0.5754193
##
##
## Real Parameter PsiBA
## 1
## 0.3223826
##
##
## Real Parameter PsiBa
## 1
## 0.6979336
##
##
## Real Parameter pA
## 1 2 3 4 5
## 0.5406104 0.5406104 0.5406104 0.5406104 0.5406104
##
##
## Real Parameter pB
## 1 2 3 4 5
## 0.9091581 0.9091581 0.9091581 0.9091581 0.9091581
##
##
## Real Parameter rA
## 1 2 3 4 5
## 0.4776383 0.4776383 0.4776383 0.4776383 0.4776383
##
##
## Real Parameter rBA
## 1 2 3 4 5
## 0.5503769 0.5503769 0.5503769 0.5503769 0.5503769
##
##
## Real Parameter rBa
## 1 2 3 4 5
## 0.5503769 0.5503769 0.5503769 0.5503769 0.5503769
##
##
## ***** 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: 746.3286
## AICc : 761.7286
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.3842900 0.2367120 -0.0796656 0.8482456
## PsiBA:(Intercept) -0.0533568 0.2183830 -0.4813875 0.3746738
## pA:(Intercept) 0.1917279 0.1650127 -0.1316969 0.5151528
## pB:(Intercept) 2.6787106 0.5300501 1.6398123 3.7176089
## rA:(Intercept) -0.2928086 0.2786267 -0.8389169 0.2532997
## rBA:(Intercept) -0.0338729 0.3444471 -0.7089891 0.6412434
## rBa:(Intercept) 0.3306090 0.3358839 -0.3277235 0.9889415
##
##
## Real Parameter PsiA
## 1
## 0.5949074
##
##
## Real Parameter PsiBA
## 1
## 0.486664
##
##
## Real Parameter PsiBa
## 1
## 0.486664
##
##
## Real Parameter pA
## 1 2 3 4 5
## 0.5477857 0.5477857 0.5477857 0.5477857 0.5477857
##
##
## Real Parameter pB
## 1 2 3 4 5
## 0.9357587 0.9357587 0.9357587 0.9357587 0.9357587
##
##
## Real Parameter rA
## 1 2 3 4 5
## 0.4273164 0.4273164 0.4273164 0.4273164 0.4273164
##
##
## Real Parameter rBA
## 1 2 3 4 5
## 0.4915326 0.4915326 0.4915326 0.4915326 0.4915326
##
##
## Real Parameter rBa
## 1 2 3 4 5
## 0.5819075 0.5819075 0.5819075 0.5819075 0.5819075
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 4
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 6
## -2lnL: 746.9737
## AICc : 760.0108
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.3755761 0.2330026 -0.0811091 0.8322612
## PsiBA:(Intercept) -0.0509461 0.2187934 -0.4797811 0.3778890
## pA:(Intercept) 0.1974423 0.1652824 -0.1265111 0.5213957
## pB:(Intercept) 2.6712087 0.5066860 1.6781041 3.6643133
## rA:(Intercept) -0.2852436 0.2675671 -0.8096751 0.2391880
## rBA:(Intercept) 0.1542574 0.2546934 -0.3449418 0.6534565
##
##
## Real Parameter PsiA
## 1
## 0.5928057
##
##
## Real Parameter PsiBA
## 1
## 0.4872662
##
##
## Real Parameter PsiBa
## 1
## 0.4872662
##
##
## Real Parameter pA
## 1 2 3 4 5
## 0.5492008 0.5492008 0.5492008 0.5492008 0.5492008
##
##
## Real Parameter pB
## 1 2 3 4 5
## 0.9353062 0.9353062 0.9353062 0.9353062 0.9353062
##
##
## Real Parameter rA
## 1 2 3 4 5
## 0.4291687 0.4291687 0.4291687 0.4291687 0.4291687
##
##
## Real Parameter rBA
## 1 2 3 4 5
## 0.5384881 0.5384881 0.5384881 0.5384881 0.5384881
##
##
## Real Parameter rBa
## 1 2 3 4 5
## 0.5384881 0.5384881 0.5384881 0.5384881 0.5384881
##
##
## ***** Starting ~El ~El ~El ~1 ~1 ~1 ~1 ~SHARE 5
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~El)PsiBA(~El)PsiBa(~El)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 10
## -2lnL: 667.4714
## AICc : 690.3285
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.3724853 0.2443912 -0.1065215 0.8514920
## PsiA:El -0.6869444 0.2752634 -1.2264607 -0.1474280
## PsiBA:(Intercept) -0.6118099 0.3881654 -1.3726142 0.1489943
## PsiBA:El 2.4555417 0.7652059 0.9557382 3.9553453
## PsiBa:(Intercept) 3.8920068 4.3189268 -4.5730899 12.3571040
## PsiBa:El 18.4986560 18.0710010 -16.9205070 53.9178190
## pA:(Intercept) 0.1552531 0.1684578 -0.1749242 0.4854304
## pB:(Intercept) 2.4408254 0.4311967 1.5956798 3.2859710
## rA:(Intercept) -0.1421755 0.2623664 -0.6564136 0.3720626
## rBA:(Intercept) 0.1938594 0.2396024 -0.2757613 0.6634802
##
##
## Real Parameter PsiA
## 1
## 0.5920594
##
##
## Real Parameter PsiBA
## 1
## 0.3516464
##
##
## Real Parameter PsiBa
## 1
## 0.9800037
##
##
## Real Parameter pA
## 1 2 3 4 5
## 0.5387355 0.5387355 0.5387355 0.5387355 0.5387355
##
##
## Real Parameter pB
## 1 2 3 4 5
## 0.9198879 0.9198879 0.9198879 0.9198879 0.9198879
##
##
## Real Parameter rA
## 1 2 3 4 5
## 0.4645159 0.4645159 0.4645159 0.4645159 0.4645159
##
##
## Real Parameter rBA
## 1 2 3 4 5
## 0.5483136 0.5483136 0.5483136 0.5483136 0.5483136
##
##
## Real Parameter rBa
## 1 2 3 4 5
## 0.5483136 0.5483136 0.5483136 0.5483136 0.5483136
# examine individual model results
model.number <-2
summary(model.fits[[model.number]])
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 7
## -2lnL: 736.624
## AICc : 752.024
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.3039969 0.2325338 -0.1517693 0.7597632
## PsiBA:(Intercept) -0.7428440 0.3261672 -1.3821317 -0.1035563
## PsiBa:(Intercept) 0.8374771 0.3707038 0.1108976 1.5640566
## pA:(Intercept) 0.1628002 0.1675176 -0.1655343 0.4911346
## pB:(Intercept) 2.3033987 0.4512448 1.4189589 3.1878385
## rA:(Intercept) -0.0895065 0.2925274 -0.6628602 0.4838471
## rBA:(Intercept) 0.2021934 0.2504401 -0.2886691 0.6930560
##
##
## Real Parameter PsiA
## 1
## 0.5754193
##
##
## Real Parameter PsiBA
## 1
## 0.3223826
##
##
## Real Parameter PsiBa
## 1
## 0.6979336
##
##
## Real Parameter pA
## 1 2 3 4 5
## 0.5406104 0.5406104 0.5406104 0.5406104 0.5406104
##
##
## Real Parameter pB
## 1 2 3 4 5
## 0.9091581 0.9091581 0.9091581 0.9091581 0.9091581
##
##
## Real Parameter rA
## 1 2 3 4 5
## 0.4776383 0.4776383 0.4776383 0.4776383 0.4776383
##
##
## Real Parameter rBA
## 1 2 3 4 5
## 0.5503769 0.5503769 0.5503769 0.5503769 0.5503769
##
##
## Real Parameter rBa
## 1 2 3 4 5
## 0.5503769 0.5503769 0.5503769 0.5503769 0.5503769
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.5754193 0.0568108 0.4621303 0.6813023
## PsiBA g1 a0 t1 0.3223826 0.0712519 0.2006669 0.4741340
## PsiBa g1 a0 t1 0.6979336 0.0781526 0.5276960 0.8269347
## pA g1 a0 t1 0.5406104 0.0416031 0.4587107 0.6203737
## pB g1 a0 t1 0.9091581 0.0372681 0.8051752 0.9603740
## rA g1 a0 t1 0.4776383 0.0729856 0.3400974 0.6186559
## rBA g1 a0 t1 0.5503769 0.0619744 0.4283297 0.6666464
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.3039969 0.2325338 -0.1517693 0.7597632
## PsiBA:(Intercept) -0.7428440 0.3261672 -1.3821317 -0.1035563
## PsiBa:(Intercept) 0.8374771 0.3707038 0.1108976 1.5640566
## pA:(Intercept) 0.1628002 0.1675176 -0.1655343 0.4911346
## pB:(Intercept) 2.3033987 0.4512448 1.4189589 3.1878385
## rA:(Intercept) -0.0895065 0.2925274 -0.6628602 0.4838471
## rBA:(Intercept) 0.2021934 0.2504401 -0.2886691 0.6930560
model.fits[[model.number]]$results$derived
## $`Species Interaction Factor`
## estimate se lcl ucl
## 1 0.6690735 0.1114698 0.4505928 0.8875543
##
## $`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.4818343 0.05382968 0.3786511 0.5865903
##
## $`Occupancy of Both Species`
## estimate se lcl ucl
## 1 0.1855051 0.04641912 0.1108941 0.2937314
get.real(model.fits[[model.number]], "PsiA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiA g1 a0 t1 1 1 0.5754193 0.0568108 0.4621303
## ucl fixed note group age time Age Time
## PsiA g1 a0 t1 0.6813023 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.3223826 0.0712519 0.2006669
## ucl fixed note group age time Age Time
## PsiBA g1 a0 t1 0.474134 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.6979336 0.0781526 0.527696
## ucl fixed note group age time Age Time
## PsiBa g1 a0 t1 0.8269347 1 0 1 0 0
get.real(model.fits[[model.number]], "pA", se=TRUE)
## all.diff.index par.index estimate se lcl
## pA g1 a0 t1 4 4 0.5406104 0.0416031 0.4587107
## pA g1 a1 t2 5 4 0.5406104 0.0416031 0.4587107
## pA g1 a2 t3 6 4 0.5406104 0.0416031 0.4587107
## pA g1 a3 t4 7 4 0.5406104 0.0416031 0.4587107
## pA g1 a4 t5 8 4 0.5406104 0.0416031 0.4587107
## ucl fixed note group age time Age Time
## pA g1 a0 t1 0.6203737 1 0 1 0 0
## pA g1 a1 t2 0.6203737 1 1 2 1 1
## pA g1 a2 t3 0.6203737 1 2 3 2 2
## pA g1 a3 t4 0.6203737 1 3 4 3 3
## pA g1 a4 t5 0.6203737 1 4 5 4 4
get.real(model.fits[[model.number]], "pB", se=TRUE)
## all.diff.index par.index estimate se lcl
## pB g1 a0 t1 9 5 0.9091581 0.0372681 0.8051752
## pB g1 a1 t2 10 5 0.9091581 0.0372681 0.8051752
## pB g1 a2 t3 11 5 0.9091581 0.0372681 0.8051752
## pB g1 a3 t4 12 5 0.9091581 0.0372681 0.8051752
## pB g1 a4 t5 13 5 0.9091581 0.0372681 0.8051752
## ucl fixed note group age time Age Time
## pB g1 a0 t1 0.960374 1 0 1 0 0
## pB g1 a1 t2 0.960374 1 1 2 1 1
## pB g1 a2 t3 0.960374 1 2 3 2 2
## pB g1 a3 t4 0.960374 1 3 4 3 3
## pB g1 a4 t5 0.960374 1 4 5 4 4
get.real(model.fits[[model.number]], "rA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rA g1 a0 t1 14 6 0.4776383 0.0729856 0.3400974
## rA g1 a1 t2 15 6 0.4776383 0.0729856 0.3400974
## rA g1 a2 t3 16 6 0.4776383 0.0729856 0.3400974
## rA g1 a3 t4 17 6 0.4776383 0.0729856 0.3400974
## rA g1 a4 t5 18 6 0.4776383 0.0729856 0.3400974
## ucl fixed note group age time Age Time
## rA g1 a0 t1 0.6186559 1 0 1 0 0
## rA g1 a1 t2 0.6186559 1 1 2 1 1
## rA g1 a2 t3 0.6186559 1 2 3 2 2
## rA g1 a3 t4 0.6186559 1 3 4 3 3
## rA g1 a4 t5 0.6186559 1 4 5 4 4
get.real(model.fits[[model.number]], "rBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBA g1 a0 t1 19 7 0.5503769 0.0619744 0.4283297
## rBA g1 a1 t2 20 7 0.5503769 0.0619744 0.4283297
## rBA g1 a2 t3 21 7 0.5503769 0.0619744 0.4283297
## rBA g1 a3 t4 22 7 0.5503769 0.0619744 0.4283297
## rBA g1 a4 t5 23 7 0.5503769 0.0619744 0.4283297
## ucl fixed note group age time Age Time
## rBA g1 a0 t1 0.6666464 1 0 1 0 0
## rBA g1 a1 t2 0.6666464 1 1 2 1 1
## rBA g1 a2 t3 0.6666464 1 2 3 2 2
## rBA g1 a3 t4 0.6666464 1 3 4 3 3
## rBA g1 a4 t5 0.6666464 1 4 5 4 4
get.real(model.fits[[model.number]], "rBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBa g1 a0 t1 24 7 0.5503769 0.0619744 0.4283297
## rBa g1 a1 t2 25 7 0.5503769 0.0619744 0.4283297
## rBa g1 a2 t3 26 7 0.5503769 0.0619744 0.4283297
## rBa g1 a3 t4 27 7 0.5503769 0.0619744 0.4283297
## rBa g1 a4 t5 28 7 0.5503769 0.0619744 0.4283297
## ucl fixed note group age time Age Time
## rBa g1 a0 t1 0.6666464 1 0 1 0 0
## rBa g1 a1 t2 0.6666464 1 1 2 1 1
## rBa g1 a2 t3 0.6666464 1 2 3 2 2
## rBa g1 a3 t4 0.6666464 1 3 4 3 3
## rBa g1 a4 t5 0.6666464 1 4 5 4 4
# collect models and make AICc table
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model npar
## 5 PsiA(~El)PsiBA(~El)PsiBa(~El)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 10
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 7
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 6
## 3 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 7
## AICc DeltaAICc weight Deviance
## 5 690.3285 0.00000 1.000000e+00 667.4714
## 2 752.0240 61.69544 4.008723e-14 736.6240
## 1 753.6893 63.36079 1.743327e-14 735.8665
## 4 760.0107 69.68221 0.000000e+00 746.9737
## 3 761.7286 71.40005 0.000000e+00 746.3286
names(model.set)
## [1] "m...001" "m...002" "m...003" "m...004" "m...005"
## [6] "model.table"
model.set$model.table
## PsiA PsiBA PsiBa pA pB rA rBA rBa
## 5 ~El ~El ~El ~1 ~1 ~1 ~1
## 2 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 4 ~1 ~1 ~1 ~1 ~1 ~1
## 3 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## model npar
## 5 PsiA(~El)PsiBA(~El)PsiBa(~El)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 10
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 7
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 6
## 3 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 7
## AICc DeltaAICc weight Deviance
## 5 690.3285 0.00000 1.000000e+00 667.4714
## 2 752.0240 61.69544 4.008723e-14 736.6240
## 1 753.6893 63.36079 1.743327e-14 735.8665
## 4 760.0107 69.68221 0.000000e+00 746.9737
## 3 761.7286 71.40005 0.000000e+00 746.3286
# model averaged values at average covariate value
PsiA.ma <- RMark::model.average(model.set, param="PsiA")
PsiA.ma
## par.index estimate se fixed note group age time
## PsiA g1 a0 t1 1 0.5920594 0.0590266 1 0 1
## Age Time
## PsiA g1 a0 t1 0 0
# make predictions as a function of elevation
# We need to know the indices to the parameters
# Look at the all.diff.index columns to get the index numbers
get.real(model.fits[[model.number]], "PsiA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiA g1 a0 t1 1 1 0.5754193 0.0568108 0.4621303
## ucl fixed note group age time Age Time
## PsiA g1 a0 t1 0.6813023 1 0 1 0 0
PsiA.ma.e <- RMark::covariate.predictions(model.set, indices=1, data=input.history)$estimates
PsiA.ma.e$parameter <- 'PsiA'
get.real(model.fits[[model.number]], "PsiBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBa g1 a0 t1 3 3 0.6979336 0.0781526 0.527696
## ucl fixed note group age time Age Time
## PsiBa g1 a0 t1 0.8269347 1 0 1 0 0
PsiBa.ma.e <- RMark::covariate.predictions(model.set, indices=3, data=input.history)$estimates
PsiBa.ma.e$parameter <- 'PsiBa'
get.real(model.fits[[model.number]], "PsiBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBA g1 a0 t1 2 2 0.3223826 0.0712519 0.2006669
## ucl fixed note group age time Age Time
## PsiBA g1 a0 t1 0.474134 1 0 1 0 0
PsiBA.ma.e <- RMark::covariate.predictions(model.set, indices=2, data=input.history)$estimates
PsiBA.ma.e$parameter <- 'PsiBA'
# Not easy to get the se of dervied parameters at the covariate values
PsiB.ma.e <- data.frame(estimate=PsiBA.ma.e$estimate*PsiA.ma.e$estimate +
PsiBa.ma.e$estimate*(1-PsiA.ma.e$estimate),
Elevation..m.=PsiA.ma.e$Elevation..m., stringAsFactors=FALSE)
PsiB.ma.e$parameter <- "PsiB"
plotdata <- plyr::rbind.fill(PsiA.ma.e,
PsiBA.ma.e,
PsiBa.ma.e,
PsiB.ma.e)
# It is clear that the fit is not satisfactory with a huge se for PsiBa and very steep decline.
# Try reversing the species and you get a much better fit
ggplot(data=plotdata, aes(x=Elevation..m., y=estimate))+
ggtitle("Effects of elevation")+
geom_point()+
geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=0.2)+
facet_wrap(~parameter, ncol=2)
# Why is psiBa so bad? This may be a case of complete separation.
# The detection probabilities are all fairly large so that over 5 visits, there is almost certain detection
# if the species is present. Consequently, the observed occupancy is likely a pretty good indication of the
# actual occupancy.
RMark::model.average(model.set, param="pA")
## par.index estimate se fixed note group age time
## pA g1 a0 t1 4 0.5387355 0.04186169 1 0 1
## pA g1 a1 t2 5 0.5387355 0.04186169 1 1 2
## pA g1 a2 t3 6 0.5387355 0.04186169 1 2 3
## pA g1 a3 t4 7 0.5387355 0.04186169 1 3 4
## pA g1 a4 t5 8 0.5387355 0.04186169 1 4 5
## Age Time
## pA g1 a0 t1 0 0
## pA g1 a1 t2 1 1
## pA g1 a2 t3 2 2
## pA g1 a3 t4 3 3
## pA g1 a4 t5 4 4
RMark::model.average(model.set, param="pB")
## par.index estimate se fixed note group age time
## pB g1 a0 t1 9 0.9198879 0.03177666 1 0 1
## pB g1 a1 t2 10 0.9198879 0.03177666 1 1 2
## pB g1 a2 t3 11 0.9198879 0.03177666 1 2 3
## pB g1 a3 t4 12 0.9198879 0.03177666 1 3 4
## pB g1 a4 t5 13 0.9198879 0.03177666 1 4 5
## Age Time
## pB g1 a0 t1 0 0
## pB g1 a1 t2 1 1
## pB g1 a2 t3 2 2
## pB g1 a3 t4 3 3
## pB g1 a4 t5 4 4
RMark::model.average(model.set, param="rA")
## par.index estimate se fixed note group age time
## rA g1 a0 t1 14 0.4645159 0.06526124 1 0 1
## rA g1 a1 t2 15 0.4645159 0.06526124 1 1 2
## rA g1 a2 t3 16 0.4645159 0.06526124 1 2 3
## rA g1 a3 t4 17 0.4645159 0.06526124 1 3 4
## rA g1 a4 t5 18 0.4645159 0.06526124 1 4 5
## Age Time
## rA g1 a0 t1 0 0
## rA g1 a1 t2 1 1
## rA g1 a2 t3 2 2
## rA g1 a3 t4 3 3
## rA g1 a4 t5 4 4
RMark::model.average(model.set, param="rBA")
## par.index estimate se fixed note group age time
## rBA g1 a0 t1 19 0.5483136 0.05934132 1 0 1
## rBA g1 a1 t2 20 0.5483136 0.05934132 1 1 2
## rBA g1 a2 t3 21 0.5483136 0.05934132 1 2 3
## rBA g1 a3 t4 22 0.5483136 0.05934132 1 3 4
## rBA g1 a4 t5 23 0.5483136 0.05934132 1 4 5
## Age Time
## rBA g1 a0 t1 0 0
## rBA g1 a1 t2 1 1
## rBA g1 a2 t3 2 2
## rBA g1 a3 t4 3 3
## rBA g1 a4 t5 4 4
RMark::model.average(model.set, param="rBa")
## par.index estimate se fixed note group age time
## rBa g1 a0 t1 24 0.5483136 0.05934132 1 0 1
## rBa g1 a1 t2 25 0.5483136 0.05934132 1 1 2
## rBa g1 a2 t3 26 0.5483136 0.05934132 1 2 3
## rBa g1 a3 t4 27 0.5483136 0.05934132 1 3 4
## rBa g1 a4 t5 28 0.5483136 0.05934132 1 4 5
## Age Time
## rBa g1 a0 t1 0 0
## rBa g1 a1 t2 1 1
## rBa g1 a2 t3 2 2
## rBa g1 a3 t4 3 3
## rBa g1 a4 t5 4 4
# Compute the observed occupance of each species
obs.occ <- data.frame(Elevation..m. = site.covar$Elevation..m.,
ooPG = apply(PG.data,1,max, na.rm=TRUE),
ooPJ = apply(PJ.data,1,max, na.rm=TRUE))
head(obs.occ)
## Elevation..m. ooPG ooPJ
## 1 502.92 1 0
## 2 518.16 1 0
## 3 518.16 0 0
## 4 518.16 1 0
## 5 533.40 0 0
## 6 548.64 1 0
ggplot(data=obs.occ, aes(x=Elevation..m., y=ooPJ))+
ggtitle("Observed occupancy of PJ given PG")+
geom_point()+
geom_smooth(method="glm", method.args = list(family = quasibinomial(link = 'logit')))+
facet_wrap(~ooPG, ncol=1)
# and look at when the species are reversed
ggplot(data=obs.occ, aes(x=Elevation..m., y=ooPG))+
ggtitle("Observed occupancy of PG given PJ")+
geom_point()+
geom_smooth(method="glm", method.args = list(family = quasibinomial(link = 'logit')))+
facet_wrap(~ooPJ, ncol=1)
# cleanup
cleanup(ask=FALSE)