# Multi species single season
# Co-occurence of Bull and Brook trout
# Data provided by Parks Canada
# 2018-12-13 code submitted by Neil Faught
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
occ.data <- readxl::read_excel("../Trout2.xlsx",
sheet="fish_sp_stacked_juv_delta")
head(occ.data)
## # A tibble: 6 x 32
## site E N patch watershed species presence count presence_juv
## <dbl> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1 6.06e5 5.65e6 spray Spray BKTR 1 2 1
## 2 1 6.06e5 5.65e6 spray Spray BKTR 1 3 1
## 3 1 6.06e5 5.65e6 spray Spray BLTR 1 1 1
## 4 1 6.06e5 5.65e6 spray Spray BLTR 1 1 1
## 5 2 6.04e5 5.66e6 spray Spray BKTR 1 2 1
## 6 2 6.04e5 5.66e6 spray Spray BKTR 1 7 1
## # ... with 23 more variables: count_juv <chr>, rep <dbl>, temp <dbl>,
## # discharge <dbl>, width <dbl>, bf <dbl>, depth <dbl>, vel <dbl>,
## # cond <chr>, DO <chr>, turb <chr>, pools <dbl>, under <dbl>,
## # willow <dbl>, cover <dbl>, ind <dbl>, agg <dbl>, wol <dbl>,
## # reachslope <dbl>, cellslope <dbl>, `cont area` <dbl>,
## # power_reachslope <dbl>, power_cellslope <dbl>
# 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 Cascade 3 4.743874 0.27503 0 0
## 2 Cascade 7 4.743874 0.47006 0 0
## 3 Cascade 9 6.884988 1.02842 0 0
## 4 Cascade 19 4.743874 0.60002 0 0
## 5 Cascade 21 7.848305 0.25631 2 2
## 6 Cascade 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 Cascade 4.743874 0.27503
## 2 Cascade 4.743874 0.47006
## 3 Cascade 6.884988 1.02842
## 4 Cascade 4.743874 0.60002
## 5 Cascade 7.848305 0.25631
## 6 Cascade 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 Cascade 4.743874 0.27503 -1.12739888 -0.8011166
## 2 0000 1 Cascade 4.743874 0.47006 -1.12739888 -0.7435560
## 3 0000 1 Cascade 6.884988 1.02842 -0.04090698 -0.5787631
## 4 0000 1 Cascade 4.743874 0.60002 -1.12739888 -0.7051999
## 5 0101 1 Cascade 7.848305 0.25631 0.44792090 -0.8066415
## 6 0000 1 Cascade 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
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 #
# Model where occupancy is not independent, and there is detection interaction
mod1 = RMark::mark(trout.data, ddl=trout.ddl,
model="2SpecConOccup",
model.parameters=list(
PsiA =list(formula= ~1),
PsiBA =list(formula= ~1),
PsiBa =list(formula= ~1),
pA =list(formula= ~1),
pB =list(formula= ~1),
rA =list(formula= ~1),
rBA =list(formula= ~1),
rBa =list(formula= ~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.1318244 0.3199195 0.5047821 1.7588666
## PsiBa:(Intercept) 0.0257211 0.1810942 -0.3292236 0.3806658
## pA:(Intercept) 3.4524334 1.2273634 1.0468012 5.8580656
## pB:(Intercept) 2.3658321 0.3480175 1.6837178 3.0479465
## rA:(Intercept) 2.3114220 0.4143267 1.4993417 3.1235023
## rBA:(Intercept) 1.8914803 0.3784509 1.1497165 2.6332441
## rBa:(Intercept) 1.6493568 1.0812930 -0.4699775 3.7686910
##
##
## 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.9141845 0.9141845
##
##
## Real Parameter rA
## 1 2
## 0.9098186 0.9098186
##
##
## Real Parameter rBA
## 1 2
## 0.8689242 0.8689242
##
##
## Real Parameter rBa
## 1 2
## 0.8388041 0.8388041
summary(mod1)
## 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.1318244 0.3199195 0.5047821 1.7588666
## PsiBa:(Intercept) 0.0257211 0.1810942 -0.3292236 0.3806658
## pA:(Intercept) 3.4524334 1.2273634 1.0468012 5.8580656
## pB:(Intercept) 2.3658321 0.3480175 1.6837178 3.0479465
## rA:(Intercept) 2.3114220 0.4143267 1.4993417 3.1235023
## rBA:(Intercept) 1.8914803 0.3784509 1.1497165 2.6332441
## rBa:(Intercept) 1.6493568 1.0812930 -0.4699775 3.7686910
##
##
## 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.9141845 0.9141845
##
##
## Real Parameter rA
## 1 2
## 0.9098186 0.9098186
##
##
## Real Parameter rBA
## 1 2
## 0.8689242 0.8689242
##
##
## Real Parameter rBa
## 1 2
## 0.8388041 0.8388041
names(mod1)
## [1] "data" "model" "title"
## [4] "model.name" "links" "mixtures"
## [7] "call" "parameters" "time.intervals"
## [10] "number.of.groups" "group.labels" "nocc"
## [13] "begin.time" "covariates" "fixed"
## [16] "design.matrix" "pims" "design.data"
## [19] "strata.labels" "mlogit.list" "profile.int"
## [22] "simplify" "model.parameters" "results"
## [25] "output"
# Probability of occupancy of A
psiA = mod1$results$real[1,]
psiA
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.3189748 0.0346479 0.2551929 0.3903433
# If species occupy sites independently, then
# psiBA = psiBa
# i.e. regardlesss of occupancy by A, occupancy of B is
# the same
# Pr(Occupancy of B | A present)
psiBA = mod1$results$real[2,]
psiBA
## estimate se lcl ucl fixed note
## PsiBA g1 a0 t1 0.7561754 0.0589849 0.6235825 0.8530677
# Pr(Occupancy of B | A absent)
psiBa = mod1$results$real[3,]
psiBa
## estimate se lcl ucl fixed note
## PsiBa g1 a0 t1 0.5064299 0.0452661 0.4184295 0.5940337
# We can compute the marginal estimate of occupancy of B
psiB = psiA[1,1]*psiBA[1,1] +
(1-psiA[1,1])*psiBa[1,1]
psiB
## [1] 0.5860924
# RMark can generate this too, with confindence limits
mod1$results$derived$`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.5860924 0.03706475 0.5120898 0.656402
# Look at Species Interaction Factor (SIF), which is
# similar to odds ratio of occupancy of B|A to B|a.
# SIF has the form:
# SIF = psi(A and B)/(psi(A)*psi(B))
# = psi(B|A)*psi(A)/(psi(A)*(psi(A)*psi(B|A) + (1-psi(A))*psi(B|a)))
# If there is no interaction between species then this is 1
# If less than 1, then occur less often than expected if independent
# If greater than 1, then occur more often than expected if independent
mod1$results$derived$`Species Interaction Factor`
## estimate se lcl ucl
## 1 1.290198 0.0920699 1.109741 1.470655
# look at odds ratio of occupancy of B|A to B|a.
# If there is no interaction between species then this is 1
# Note, no confidence intervals (But RPrence analysis has them)
(psiBA[1,1]/(1-psiBA[1,1])) /
(psiBa[1,1]/(1-psiBa[1,1]))
## [1] 3.022557
# For example, the marginal estimates are and prob if indendent are
c(psiA[1,1],psiB, psiA[1,1]*psiB)
## [1] 0.3189748 0.5860924 0.1869487
# Actual estimation of occupancy of both species
psiBA[1,1]*psiA[1,1]
## [1] 0.2412009
# Probability of detection of species A if alone
pA = mod1$results$real[4,]
pA
## estimate se lcl ucl fixed note
## pA g1 a0 t1 0.9693036 0.0365191 0.7401602 0.9971514
# Probability of detection of species B if alone
pB = mod1$results$real[5,]
pB
## estimate se lcl ucl fixed note
## pB g1 a0 t1 0.9141845 0.0273024 0.8433962 0.9546938
# Probability of detection when both species present
# rA - just A; rBA B given A also detected; rBa B given A not detected
# Again, if species do not interact, then
# rBA = rBa
rA = mod1$results$real[6,]
rA
## estimate se lcl ucl fixed note
## rA g1 a0 t1 0.9098186 0.033995 0.8174763 0.9578518
rBA = mod1$results$real[7,]
rBA
## estimate se lcl ucl fixed note
## rBA g1 a0 t1 0.8689242 0.0431036 0.7594591 0.9329707
rBa = mod1$results$real[8,]
rBa
## estimate se lcl ucl fixed note
## rBa g1 a0 t1 0.8388041 0.1462036 0.3846216 0.9774385
# similarly, the odds ratio of detection of B|A vs B|a
# rho < 1, B is harder to detect if A detected compared to detecting B if A not detected
# Note: RMark doesn't automatically produce confidence interaval for rho
# RPresence does though
(rBA[1,1]/(1-rBA[1,1])) /
(rBa[1,1]/(1-rBa[1,1]))
## [1] 1.273951
###########################################################################
# The detection rates are all high and almost all the same.
# Fit a simpler model where pA = pB, rBA = rBa
# Model where pA = pB = rA = rBA = rBa cannot be fit in RMark. It can be
# fit in RPresence. See http://www.phidot.org/forum/viewtopic.php?f=21&t=2974 for details
mod1.psimp = RMark::mark(trout.data, ddl=trout.ddl,
model="2SpecConOccup",
model.parameters=list(
PsiA =list(formula= ~1),
PsiBA =list(formula= ~1),
PsiBa =list(formula= ~1),
pA =list(formula= ~1, share=TRUE),
rA =list(formula= ~1),
rBA =list(formula= ~1, share=TRUE)
))
##
## 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.7568656 0.1595781 -1.0696387 -0.4440925
## PsiBA:(Intercept) 1.1257238 0.3193645 0.4997693 1.7516783
## PsiBa:(Intercept) 0.0235540 0.1808124 -0.3308382 0.3779463
## pA:(Intercept) 2.4991784 0.3298836 1.8526065 3.1457503
## rA:(Intercept) 2.3254031 0.4141395 1.5136897 3.1371165
## rBA:(Intercept) 1.8711051 0.3615426 1.1624817 2.5797285
##
##
## Real Parameter PsiA
## 1
## 0.3193272
##
##
## 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
summary(mod1.psimp)
## 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.7568656 0.1595781 -1.0696387 -0.4440925
## PsiBA:(Intercept) 1.1257238 0.3193645 0.4997693 1.7516783
## PsiBa:(Intercept) 0.0235540 0.1808124 -0.3308382 0.3779463
## pA:(Intercept) 2.4991784 0.3298836 1.8526065 3.1457503
## rA:(Intercept) 2.3254031 0.4141395 1.5136897 3.1371165
## rBA:(Intercept) 1.8711051 0.3615426 1.1624817 2.5797285
##
##
## Real Parameter PsiA
## 1
## 0.3193272
##
##
## 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
names(mod1.psimp)
## [1] "data" "model" "title"
## [4] "model.name" "links" "mixtures"
## [7] "call" "parameters" "time.intervals"
## [10] "number.of.groups" "group.labels" "nocc"
## [13] "begin.time" "covariates" "fixed"
## [16] "design.matrix" "pims" "design.data"
## [19] "strata.labels" "mlogit.list" "profile.int"
## [22] "simplify" "model.parameters" "results"
## [25] "output"
# Probability of occupancy of A
psiA = mod1.psimp$results$real[1,]
psiA
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.3193272 0.0346855 0.2554718 0.3907662
# If species occupy sites independently, then
# psiBA = psiBa
# i.e. regardlesss of occupancy by A, occupancy of B is
# the same
# Pr(Occupancy of B | A present)
psiBA = mod1.psimp$results$real[2,]
psiBA
## estimate se lcl ucl fixed note
## PsiBA g1 a0 t1 0.7550489 0.0590665 0.6224051 0.8521644
# Pr(Occupancy of B | A absent)
psiBa = mod1.psimp$results$real[3,]
psiBa
## estimate se lcl ucl fixed note
## PsiBa g1 a0 t1 0.5058882 0.0451968 0.4180367 0.5933777
# We can compute the marginal estimate of occupancy of B
psiB = psiA[1,1]*psiBA[1,1] +
(1-psiA[1,1])*psiBa[1,1]
psiB
## [1] 0.585452
# RMark can generate this too, with confindence limits
mod1.psimp$results$derived$`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.585452 0.03700965 0.5115754 0.6556754
# Look at Species Interaction Factor (SIF), which is
# similar to odds ratio of occupancy of B|A to B|a.
# SIF has the form:
# SIF = psi(A and B)/(psi(A)*psi(B))
# = psi(B|A)*psi(A)/(psi(A)*(psi(A)*psi(B|A) + (1-psi(A))*psi(B|a)))
# If there is no interaction between species then this is 1
# If less than 1, then occur less often than expected if independent
# If greater than 1, then occur more often than expected if independent
mod1.psimp$results$derived$`Species Interaction Factor`
## estimate se lcl ucl
## 1 1.289685 0.09218856 1.108996 1.470375
# look at odds ratio of occupancy of B|A to B|a.
# If there is no interaction between species then this is 1
# Note, no confidence intervals (But RPrence analysis has them)
(psiBA[1,1]/(1-psiBA[1,1])) /
(psiBa[1,1]/(1-psiBa[1,1]))
## [1] 3.010692
# For example, the marginal estimates are and prob if indendent are
c(psiA[1,1],psiB, psiA[1,1]*psiB)
## [1] 0.3193272 0.5854520 0.1869507
# Actual estimation of occupancy of both species
psiBA[1,1]*psiA[1,1]
## [1] 0.2411077
# Probability of detection of species A if alone
pA = mod1.psimp$results$real[4,]
pA
## estimate se lcl ucl fixed note
## pA g1 a0 t1 0.9240842 0.0231422 0.8644328 0.9587409
# Probability of detection of species B if alone
pB = mod1.psimp$results$real[4,]
pB
## estimate se lcl ucl fixed note
## pA g1 a0 t1 0.9240842 0.0231422 0.8644328 0.9587409
# Probability of detection when both species present
# rA - just A; rBA B given A also detected; rBa B given A not detected
# Again, if species do not interact, then
# rBA = rBa
rA = mod1.psimp$results$real[5,]
rA
## estimate se lcl ucl fixed note
## rA g1 a0 t1 0.9109592 0.0335919 0.8196074 0.9583981
rBA = mod1.psimp$results$real[6,]
rBA
## estimate se lcl ucl fixed note
## rBA g1 a0 t1 0.8665861 0.0417996 0.7617834 0.9295455
rBa = mod1.psimp$results$real[6,]
rBa
## estimate se lcl ucl fixed note
## rBA g1 a0 t1 0.8665861 0.0417996 0.7617834 0.9295455
# similarly, the odds ratio of detection of B|A vs B|a
# rho < 1, B is harder to detect if A detected compared to detecting B if A not detected
# Note: RMark doesn't automatically produce confidence interaval for rho
# RPresence does though
(rBA[1,1]/(1-rBA[1,1])) /
(rBa[1,1]/(1-rBa[1,1]))
## [1] 1
# Creat AICc table for mod1 and mod1.psimp
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model npar AICc
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB()rA(~1)rBA(~1)rBa() 6 667.3158
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8 670.6256
## DeltaAICc weight Deviance
## 2 0.000000 0.8395518 654.8385
## 1 3.309793 0.1604482 653.7980
rm(model.set)
###########################################################################
# mod1.occind is same as model 1, but we assume that species occupancy is
# independent. This implies that a psiBA=psiBa.
mod1.occind = RMark::mark(trout.data, ddl=trout.ddl,
model="2SpecConOccup",
model.parameters=list(
PsiA =list(formula= ~1),
PsiBA =list(formula= ~1, share = TRUE),
pA =list(formula= ~1),
pB =list(formula= ~1),
rA =list(formula= ~1),
rBA =list(formula= ~1),
rBa =list(formula= ~1)
))
##
## 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.0734976 -0.4491536
## PsiBA:(Intercept) 0.3422651 0.1521157 0.0441182 0.6404119
## pA:(Intercept) 3.3385667 1.1198838 1.1435943 5.5335390
## pB:(Intercept) 2.3240806 0.3540616 1.6301198 3.0180414
## rA:(Intercept) 2.3477257 0.4065485 1.5508906 3.1445608
## rBA:(Intercept) 1.9753221 0.3632630 1.2633266 2.6873175
## rBa:(Intercept) 1.7333485 1.0822361 -0.3878342 3.8545312
##
##
## 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
summary(mod1.occind)
## 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.0734976 -0.4491536
## PsiBA:(Intercept) 0.3422651 0.1521157 0.0441182 0.6404119
## pA:(Intercept) 3.3385667 1.1198838 1.1435943 5.5335390
## pB:(Intercept) 2.3240806 0.3540616 1.6301198 3.0180414
## rA:(Intercept) 2.3477257 0.4065485 1.5508906 3.1445608
## rBA:(Intercept) 1.9753221 0.3632630 1.2633266 2.6873175
## rBa:(Intercept) 1.7333485 1.0822361 -0.3878342 3.8545312
##
##
## 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
names(mod1.occind)
## [1] "data" "model" "title"
## [4] "model.name" "links" "mixtures"
## [7] "call" "parameters" "time.intervals"
## [10] "number.of.groups" "group.labels" "nocc"
## [13] "begin.time" "covariates" "fixed"
## [16] "design.matrix" "pims" "design.data"
## [19] "strata.labels" "mlogit.list" "profile.int"
## [22] "simplify" "model.parameters" "results"
## [25] "output"
# Probability of occupancy of A
psiA = mod1.occind$results$real[1,]
psiA
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.3183585 0.0345629 0.2547385 0.389562
# If species occupy sites independently, then
# psiBA = psiBa
# i.e. regardlesss of occupancy by A, occupancy of B is
# the same
# Pr(Occupancy of B | A present)
psiBA = mod1.occind$results$real[2,]
psiBA
## estimate se lcl ucl fixed note
## PsiBA g1 a0 t1 0.5847406 0.0369366 0.5110278 0.6548466
# Pr(Occupancy of B | A absent)
psiBa = mod1.occind$results$real[2,]
psiBa
## estimate se lcl ucl fixed note
## PsiBA g1 a0 t1 0.5847406 0.0369366 0.5110278 0.6548466
# We can compute the marginal estimate of occupancy of B
psiB = psiA[1,1]*psiBA[1,1] +
(1-psiA[1,1])*psiBa[1,1]
psiB
## [1] 0.5847406
# RMark can generate this too, with confindence limits
mod1.occind$results$derived$`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.5847406 0.0369366 0.5110278 0.6548466
# Look at Species Interaction Factor (SIF), which is
# similar to odds ratio of occupancy of B|A to B|a.
# SIF has the form:
# SIF = psi(A and B)/(psi(A)*psi(B))
# = psi(B|A)*psi(A)/(psi(A)*(psi(A)*psi(B|A) + (1-psi(A))*psi(B|a)))
# If there is no interaction between species then this is 1
# If less than 1, then occur less often than expected if independent
# If greater than 1, then occur more often than expected if independent
mod1.occind$results$derived$`Species Interaction Factor`
## estimate se lcl ucl
## 1 1 0 1 1
# look at odds ratio of occupancy of B|A to B|a.
# If there is no interaction between species then this is 1
# Note, no confidence intervals (But RPrence analysis has them)
(psiBA[1,1]/(1-psiBA[1,1])) /
(psiBa[1,1]/(1-psiBa[1,1]))
## [1] 1
# For example, the marginal estimates are and prob if indendent are
c(psiA[1,1],psiB, psiA[1,1]*psiB)
## [1] 0.3183585 0.5847406 0.1861571
# Actual estimation of occupancy of both species
psiBA[1,1]*psiA[1,1]
## [1] 0.1861571
# Probability of detection of species A if alone
pA = mod1.occind$results$real[3,]
pA
## estimate se lcl ucl fixed note
## pA g1 a0 t1 0.9657284 0.0370648 0.7583389 0.9960636
# Probability of detection of species B if alone
pB = mod1.occind$results$real[4,]
pB
## estimate se lcl ucl fixed note
## pB g1 a0 t1 0.9108518 0.0287501 0.8361861 0.9533826
# Probability of detection when both species present
# rA - just A; rBA B given A also detected; rBa B given A not detected
# Again, if species do not interact, then
# rBA = rBa
rA = mod1.occind$results$real[5,]
rA
## estimate se lcl ucl fixed note
## rA g1 a0 t1 0.9127533 0.0323754 0.8250423 0.9586939
rBA = mod1.occind$results$real[6,]
rBA
## estimate se lcl ucl fixed note
## rBA g1 a0 t1 0.8781816 0.0388614 0.7795982 0.9362741
rBa = mod1.occind$results$real[7,]
rBa
## estimate se lcl ucl fixed note
## rBa g1 a0 t1 0.8498402 0.1381061 0.4042388 0.9792559
# similarly, the odds ratio of detection of B|A vs B|a
# rho < 1, B is harder to detect if A detected compared to detecting B if A not detected
# Note: RMark doesn't automatically produce confidence interaval for rho
# RPresence does though
(rBA[1,1]/(1-rBA[1,1])) /
(rBa[1,1]/(1-rBa[1,1]))
## [1] 1.273761
# Create AICc table comparing all 4 models
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model npar AICc
## 3 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB()rA(~1)rBA(~1)rBa() 6 667.3158
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8 670.6256
## 2 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 7 678.5194
## DeltaAICc weight Deviance
## 3 0.000000 0.83695807 654.8385
## 1 3.309793 0.15995255 653.7980
## 2 11.203607 0.00308938 663.8794
rm(model.set)
###########################################################################
# Look at the effect of temp
mod.psi.temp = RMark::mark(trout.data, ddl=trout.ddl,
model="2SpecConOccup",
model.parameters=list(
PsiA =list(formula= ~temp.std),
PsiBA =list(formula= ~temp.std),
PsiBa =list(formula= ~temp.std),
pA =list(formula= ~1, share=TRUE),
rA =list(formula= ~1),
rBA =list(formula= ~1, share=TRUE)
))
##
## 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.1535588 0.2448024 -1.6333715 -0.6737461
## PsiA:temp.std 2.5263632 0.4057295 1.7311334 3.3215930
## PsiBA:(Intercept) 3.8249049 1.1487518 1.5733513 6.0764584
## PsiBA:temp.std -2.0174865 0.6754512 -3.3413709 -0.6936021
## PsiBa:(Intercept) 0.4976977 0.2753908 -0.0420682 1.0374637
## PsiBa:temp.std 0.9257767 0.3689922 0.2025519 1.6490014
## pA:(Intercept) 2.4436670 0.3379555 1.7812743 3.1060598
## rA:(Intercept) 2.3495750 0.4072936 1.5512795 3.1478705
## rBA:(Intercept) 1.9541313 0.3471787 1.2736611 2.6346016
##
##
## Real Parameter PsiA
## 1
## 0.2398397
##
##
## 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
summary(mod.psi.temp)
## 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.1535588 0.2448024 -1.6333715 -0.6737461
## PsiA:temp.std 2.5263632 0.4057295 1.7311334 3.3215930
## PsiBA:(Intercept) 3.8249049 1.1487518 1.5733513 6.0764584
## PsiBA:temp.std -2.0174865 0.6754512 -3.3413709 -0.6936021
## PsiBa:(Intercept) 0.4976977 0.2753908 -0.0420682 1.0374637
## PsiBa:temp.std 0.9257767 0.3689922 0.2025519 1.6490014
## pA:(Intercept) 2.4436670 0.3379555 1.7812743 3.1060598
## rA:(Intercept) 2.3495750 0.4072936 1.5512795 3.1478705
## rBA:(Intercept) 1.9541313 0.3471787 1.2736611 2.6346016
##
##
## Real Parameter PsiA
## 1
## 0.2398397
##
##
## 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
names(mod.psi.temp)
## [1] "data" "model" "title"
## [4] "model.name" "links" "mixtures"
## [7] "call" "parameters" "time.intervals"
## [10] "number.of.groups" "group.labels" "nocc"
## [13] "begin.time" "covariates" "fixed"
## [16] "design.matrix" "pims" "design.data"
## [19] "strata.labels" "mlogit.list" "profile.int"
## [22] "simplify" "model.parameters" "results"
## [25] "output"
# Probability of occupancy of A, B|A, anb B|a now depends on occupancy (on the logit scale)
# We need to look at the beta values
mod.psi.temp$results$beta[1:6,]
## estimate se lcl ucl
## PsiA:(Intercept) -1.1535588 0.2448024 -1.6333715 -0.6737461
## PsiA:temp.std 2.5263632 0.4057295 1.7311334 3.3215930
## PsiBA:(Intercept) 3.8249049 1.1487518 1.5733513 6.0764584
## PsiBA:temp.std -2.0174865 0.6754512 -3.3413709 -0.6936021
## PsiBa:(Intercept) 0.4976977 0.2753908 -0.0420682 1.0374637
## PsiBa:temp.std 0.9257767 0.3689922 0.2025519 1.6490014
str(mod.psi.temp$results$beta$estimate[1:6])
## num [1:6] -1.154 2.526 3.825 -2.017 0.498 ...
range(input.history$temp)
## [1] 3.375832 11.561497
predictPSI <- data.frame(temp=seq(3,12,.1))
predictPSI$temp.std <- (predictPSI$temp - temp.mean)/temp.sd
predictPSI$psiA.logit <- mod.psi.temp$results$beta$estimate[1] +
mod.psi.temp$results$beta$estimate[2]*predictPSI$temp.std
predictPSI$psiBA.logit <- mod.psi.temp$results$beta$estimate[3] +
mod.psi.temp$results$beta$estimate[4]*predictPSI$temp.std
predictPSI$psiBa.logit <- mod.psi.temp$results$beta$estimate[5] +
mod.psi.temp$results$beta$estimate[6]*predictPSI$temp.std
expit <- function (x) {1/(1+exp(-x))}
predictPSI$psiA <- expit(predictPSI$psiA.logit)
predictPSI$psiBA <- expit(predictPSI$psiBA.logit)
predictPSI$psiBa <- expit(predictPSI$psiBa.logit)
# We can compute the marginal estimate of occupancy of B
predictPSI$psiB <- predictPSI$psiA * predictPSI$psiBA+
(1-predictPSI$psiA) * predictPSI$psiBa
predictPSI$psiAandB <- predictPSI$psiA * predictPSI$psiBA
predictPSI[1:3,]
## temp temp.std psiA.logit psiBA.logit psiBa.logit psiA psiBA
## 1 3.0 -2.012314 -6.237395 7.884721 -1.365256 0.001951127 0.9996237
## 2 3.1 -1.961570 -6.109197 7.782346 -1.318278 0.002217407 0.9995831
## 3 3.2 -1.910826 -5.980998 7.679970 -1.271300 0.002519937 0.9995382
## psiBa psiB psiAandB
## 1 0.2033874 0.2049410 0.001950393
## 2 0.2111049 0.2128533 0.002216483
## 3 0.2190348 0.2210016 0.002518773
plotdata <- reshape2::melt(predictPSI,
id.vars=c("temp","temp.std"),
measure.vars=c("psiA","psiB","psiAandB"),
variable.name="Species",
value.name='P.Occupancy')
species.by.temp <- ggplot2::ggplot(dat=plotdata, aes(x=temp, y=P.Occupancy, color=Species))+
ggtitle("Occupancy as a function of temp\nSpecies A = BKTR; Species B=BLTR")+
geom_line()
species.by.temp
ggsave(plot=species.by.temp, file='species-by-temp.png',
h=4, w=6, units="in", dpi=300)
# Probability of occupancy of A
psiA = mod.psi.temp$results$real[1,]
psiA
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.2398397 0.0446315 0.163369 0.3376585
# If species occupy sites independently, then
# psiBA = psiBa
# i.e. regardlesss of occupancy by A, occupancy of B is
# the same
# Pr(Occupancy of B | A present)
psiBA = mod.psi.temp$results$real[2,]
psiBA
## estimate se lcl ucl fixed note
## PsiBA g1 a0 t1 0.9786455 0.0240072 0.8282608 0.997709
# Pr(Occupancy of B | A absent)
psiBa = mod.psi.temp$results$real[3,]
psiBa
## estimate se lcl ucl fixed note
## PsiBa g1 a0 t1 0.6219181 0.0647543 0.4894845 0.7383603
# We can compute the marginal estimate of occupancy of B
psiB = psiA[1,1]*psiBA[1,1] +
(1-psiA[1,1])*psiBa[1,1]
psiB
## [1] 0.7074755
# RMark can generate this too, with confindence limits
mod.psi.temp$results$derived$`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.7074755 0.05219491 0.5960015 0.7985855
# Look at Species Interaction Factor (SIF), which is
# similar to odds ratio of occupancy of B|A to B|a.
# SIF has the form:
# SIF = psi(A and B)/(psi(A)*psi(B))
# = psi(B|A)*psi(A)/(psi(A)*(psi(A)*psi(B|A) + (1-psi(A))*psi(B|a)))
# If there is no interaction between species then this is 1
# If less than 1, then occur less often than expected if independent
# If greater than 1, then occur more often than expected if independent
mod.psi.temp$results$derived$`Species Interaction Factor`
## estimate se lcl ucl
## 1 1.383292 0.1048246 1.177836 1.588749
# look at odds ratio of occupancy of B|A to B|a.
# If there is no interaction between species then this is 1
# Note, no confidence intervals (But RPrence analysis has them)
(psiBA[1,1]/(1-psiBA[1,1])) /
(psiBa[1,1]/(1-psiBa[1,1]))
## [1] 27.86049
# For example, the marginal estimates are and prob if indendent are
c(psiA[1,1],psiB, psiA[1,1]*psiB)
## [1] 0.2398397 0.7074755 0.1696807
# Actual estimation of occupancy of both species
psiBA[1,1]*psiA[1,1]
## [1] 0.234718
# Probability of detection of species A if alone
pA = mod.psi.temp$results$real[4,]
pA
## estimate se lcl ucl fixed note
## pA g1 a0 t1 0.9200971 0.024846 0.8558541 0.957142
# Probability of detection of species B if alone
pB = mod.psi.temp$results$real[4,]
pB
## estimate se lcl ucl fixed note
## pA g1 a0 t1 0.9200971 0.024846 0.8558541 0.957142
# Probability of detection when both species present
# rA - just A; rBA B given A also detected; rBa B given A not detected
# Again, if species do not interact, then
# rBA = rBa
rA = mod.psi.temp$results$real[5,]
rA
## estimate se lcl ucl fixed note
## rA g1 a0 t1 0.9129004 0.0323852 0.8250985 0.9588247
rBA = mod.psi.temp$results$real[6,]
rBA
## estimate se lcl ucl fixed note
## rBA g1 a0 t1 0.8758964 0.037739 0.7813688 0.9330555
rBa = mod.psi.temp$results$real[6,]
rBa
## estimate se lcl ucl fixed note
## rBA g1 a0 t1 0.8758964 0.037739 0.7813688 0.9330555
# similarly, the odds ratio of detection of B|A vs B|a
# rho < 1, B is harder to detect if A detected compared to detecting B if A not detected
# Note: RMark doesn't automatically produce confidence interaval for rho
# RPresence does though
(rBA[1,1]/(1-rBA[1,1])) /
(rBa[1,1]/(1-rBa[1,1]))
## [1] 1
###########################################################################
# Create AICc table comparing all 4 models
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model
## 1 PsiA(~temp.std)PsiBA(~temp.std)PsiBa(~temp.std)pA(~1)pB()rA(~1)rBA(~1)rBa()
## 4 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB()rA(~1)rBA(~1)rBa()
## 2 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
## 1 9 540.4340 0.0000 1 521.3936
## 4 6 667.3158 126.8818 0 654.8385
## 2 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)