# Multi species single season
# Co-occurence of Bull and Brook trout
# Data provided by Parks Canada
# 2020-06-24 CJS real$rho and real$nu now in derived$
library(ggplot2)
library(plyr)
library(readxl)
library(reshape2)
library(RPresence)
# 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 W3 BKTR 1 1 10.2 8.93
## 2 1 W3 BKTR 1 2 10.2 8.93
## 3 1 W3 BLTR 1 1 10.2 8.93
## 4 1 W3 BLTR 1 2 10.2 8.93
## 5 2 W3 BKTR 1 1 9.41 10.4
## 6 2 W3 BKTR 1 2 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.
# If you change the coding, the program doesn't coverge well (perhaps a
# change in initial values may be helpful)
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.history <- history[,c("r1","r2")]
input.history[1:5,]
## r1 r2
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 2 2
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
Nsites <- nrow(input.history)
Nvisits<- ncol(input.history)
# create the pao file
trout.pao <- createPao(input.history,
unitcov=site.covar,
title="trout multi species - co-occurance")
summary(trout.pao)
## paoname=pres.pao
## title=trout multi species - co-occurance
## Naive occ=0.6612022
## naiveR =0.3770492
## nunits nsurveys nseasons nsurveyseason nmethods
## "183" "2" "1" "2" "1"
## nunitcov nsurvcov
## "5" "1"
## unit covariates : watershed temp discharge temp.std discharge.std
## survey covariates: SURVEY
# Use the psiAB parameterization
mod1 <- occMod(model=list(psi~SP+INT, p~SP+INT_o+SP:INT_o + INT_d),
data=trout.pao,
type="so.2sp.1")# param="PsiBA")
summary(mod1)
## Model name=psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d)
## AIC=669.798
## -2*log-likelihood=653.798
## num. par=8
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 5.88 signifcant digits.
names(mod1$real)
## [1] "psiA" "psiBA" "psiBa" "pA" "pB" "rA" "rBA" "rBa"
# Probability of occupancy of A
mod1$real$psiA [1,]
## est se lower upper
## unit1 0.319 0.0346 0.2552 0.3903
# 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)
mod1$real$psiBA [1,]
## est se lower upper
## unit1 0.7562 0.059 0.6236 0.853
# Pr(Occupancy of B | A absent)
mod1$real$psiBa[1,]
## est se lower upper
## unit1 0.5064 0.0453 0.4184 0.594
# We can compute the marginal estimate of occupancy of B
psiB <- mod1$real$psiA * mod1$real$psiBA+
(1-mod1$real$psiA) * mod1$real$psiBa
psiB[1,]
## est se lower upper
## unit1 0.5860862 0.04577402 0.470767 0.6950877
# look at odds ratio of occupancy of B|A to 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$derived$nu[1,]
## est se lower upper
## unit1 3.022557 1.112362 1.469296 6.217842
# For example, the marginal estimates are and prob if indendent are
c(mod1$real$psiA[1,1],psiB[1,1], mod1$real$psiA[1,1]*psiB[1,1])
## [1] 0.3190000 0.5860862 0.1869615
# Actual estimation of occupancy of both species
mod1$real$psiBA [1,1]*mod1$real$psiA[1,1]
## [1] 0.2412278
# Probability of detection of species A if alone
mod1$real$pA[1,]
## est se lower upper
## unit1_1-1 0.9693 0.0367 0.7374 0.9972
# Probability of detection of species B if alone
mod1$real$pB[1,]
## est se lower upper
## unit1_1-1 0.9142 0.0273 0.8434 0.9547
# 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
mod1$real$rA[1,]
## est se lower upper
## unit1_1-1 0.9098 0.034 0.8175 0.9579
mod1$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.8689 0.0431 0.7595 0.933
mod1$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.8388 0.1462 0.3845 0.9774
# 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
mod1$derived$rho[1,]
## est se lower upper
## unit1_1-1 1.273952 1.438312 0.1393528 11.64637
# The detection rates are all high and almost all the same.
# Fit a much simpler model
mod1.psimp <- occMod(model=list(psi~SP+INT, p~1),
data=trout.pao,
type="so.2sp.1")# param="PsiBA")
summary(mod1.psimp)
## Model name=psi(SP P INT)p()
## AIC=664.501
## -2*log-likelihood=656.501
## num. par=4
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 7 signifcant digits.
names(mod1.psimp$real)
## [1] "psiA" "psiBA" "psiBa" "pA" "pB" "rA" "rBA" "rBa"
# Probability of occupancy of A
mod1.psimp$real$psiA [1,]
## est se lower upper
## unit1 0.3198 0.0347 0.2558 0.3913
# 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)
mod1.psimp$real$psiBA [1,]
## est se lower upper
## unit1 0.748 0.058 0.6188 0.8444
# Pr(Occupancy of B | AA absent)
mod1.psimp$real$psiBa[1,]
## est se lower upper
## unit1 0.5075 0.0453 0.4194 0.5952
# We can compute the marginal estimate of occupancy of B
psiB <- mod1.psimp$real$psiA * mod1.psimp$real$psiBA+
(1-mod1.psimp$real$psiA) * mod1.psimp$real$psiBa
psiB[1,]
## est se lower upper
## unit1 0.5844119 0.04574069 0.4704065 0.692712
# look at odds ratio of occupancy of B|A to 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$derived$nu[1,]
## est se lower upper
## unit1 2.880391 1.031079 1.428057 5.80975
# For example, the marginal estimates are and prob if indendent are
c(mod1.psimp$real$psiA[1,1],psiB[1,1], mod1.psimp$real$psiA[1,1]*psiB[1,1])
## [1] 0.3198000 0.5844119 0.1868949
# Actual estimation of occupancy of both species
mod1.psimp$real$psiBA [1,1]*mod1.psimp$real$psiA[1,1]
## [1] 0.2392104
# Probability of detection of species A if alone
mod1.psimp$real$pA[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
# Probability of detection of species B if alone
mod1.psimp$real$pB[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
# Probability of detection when both species present
# rA - just A; rBA B give A; rBa B given A not detected
# Again, if species do not interact, the
# rBA = rBa
mod1.psimp$real$rA[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
mod1.psimp$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
mod1.psimp$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
# 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
mod1.psimp$derived$rho[1,]
## est se lower upper
## unit1_1-1 1 0 1 1
models<-list(mod1,
mod1.psimp)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll
## 1 psi(SP P INT)p() 0.0 0.934 4 656.5
## 2 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d) 5.3 0.066 8 653.8
## warn.conv warn.VC
## 1 7.00 0
## 2 5.88 0
# mod1.occind is same as model 1, but we assume that species occupancy is
# independent. This implies that a psiBA=psiBa and is found by not
# fitting the interaction term in the occupancy mode
mod1.occind <- occMod(model=list(psi~SP, p~1),
data=trout.pao,
type="so.2sp.1")
summary(mod1.occind)
## Model name=psi(SP)p()
## AIC=671.9923
## -2*log-likelihood=665.9923
## num. par=3
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 7 signifcant digits.
names(mod1.occind$real)
## [1] "psiA" "psiBA" "psiBa" "pA" "pB" "rA" "rBA" "rBa"
# Probability of occupancy of A
mod1.occind$real$psiA [1,]
## est se lower upper
## unit1 0.3198 0.0347 0.2558 0.3913
# 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)
mod1.occind$real$psiBA [1,]
## est se lower upper
## unit1 0.5844 0.0369 0.5108 0.6544
# Pr(Occupancy of B | AA absent)
mod1.occind$real$psiBa[1,]
## est se lower upper
## unit1 0.5844 0.0369 0.5108 0.6544
# We can compute the marginal estimate of occupancy of B
psiB <- mod1.occind$real$psiA * mod1.occind$real$psiBA+
(1-mod1.occind$real$psiA) * mod1.occind$real$psiBa
psiB[1,]
## est se lower upper
## unit1 0.5844 0.0369 0.5108 0.6544
# look at odds ratio of occupancy of B|A to 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$derived$nu[1,]
## est se lower upper
## unit1 1 0 1 1
# For example, the marginal estimates are and prob if indendent are
c(mod1.occind$real$psiA[1,1],psiB[1,1], mod1.occind$real$psiA[1,1]*psiB[1,1])
## [1] 0.3198000 0.5844000 0.1868911
# Actual estimation of occupancy of both species
mod1.occind$real$psiBA [1,1]*mod1.occind$real$psiA[1,1]
## [1] 0.1868911
# Probability of detection of species A if alone
mod1.occind$real$pA[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
# Probability of detection of species B if alone
mod1.occind$real$pB[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
# Probability of detection when both species present
# rA - just A; rBA B give A; rBa B given A not detected
# Again, if species do not interact, the
# rBA = rBa
mod1.occind$real$rA[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
mod1.occind$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
mod1.occind$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.906 0.0177 0.8651 0.9354
# 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
mod1.occind$derived$rho[1,]
## est se lower upper
## unit1_1-1 1 0 1 1
models<-list(mod1,
mod1.psimp,
mod1.occind)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll
## 1 psi(SP P INT)p() 0.00 0.914 4 656.50
## 2 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d) 5.30 0.065 8 653.80
## 3 psi(SP)p() 7.49 0.022 3 665.99
## warn.conv warn.VC
## 1 7.00 0
## 2 5.88 0
## 3 7.00 0
# Look at the effect of temp
mod.psi.temp <- occMod(model=list(
psi~temp.std+SP+SP:temp.std+INT+INT:temp.std,
p ~1),
data=trout.pao,
type="so.2sp.1")
summary(mod.psi.temp)
## Model name=psi(temp.std P SP P SP T temp.std P INT P INT T temp.std)p()
## AIC=536.4685
## -2*log-likelihood=522.4685
## num. par=7
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 7 signifcant digits.
names(mod.psi.temp$real)
## [1] "psiA" "psiBA" "psiBa" "pA" "pB" "rA" "rBA" "rBa"
# 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$beta$psi
## est se
## A1_psiA -1.146655 0.245757
## A2_psiA.temp.std_psiA 2.558980 0.419684
## A3_psiBA 4.943428 1.135612
## A4_psiBa -3.270303 1.153666
## A5_psiBA.temp.std:SP2_psiBA -4.576692 0.769018
## A6_psiBa.temp.std:INT2_psiBa 2.979576 0.753615
str(mod.psi.temp$beta$psi)
## 'data.frame': 6 obs. of 2 variables:
## $ est: num -1.15 2.56 4.94 -3.27 -4.58 ...
## $ se : num 0.246 0.42 1.136 1.154 0.769 ...
mod.psi.temp$dmat$psi
## a1 a2 a3 a4 a5 a6
## psiA "1" "temp.std_psiA" "0" "0" "0" "0"
## psiBA "1" "temp.std_psiBA" "1" "0" "temp.std:SP2_psiBA" "0"
## psiBa "1" "temp.std_psiBa" "1" "1" "temp.std:SP2_psiBa" "temp.std:INT2_psiBa"
range(site.covar$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$beta$psi[1,1] + mod.psi.temp$beta$psi[2,1]*predictPSI$temp.std
predictPSI$psiBA.logit <- mod.psi.temp$beta$psi[3,1] + mod.psi.temp$beta$psi[5,1]*predictPSI$temp.std
predictPSI$psiBa.logit <- mod.psi.temp$beta$psi[4,1] + mod.psi.temp$beta$psi[6,1]*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.296127 14.15317 -9.266146 0.001840039 0.9999993
## 2 3.1 -1.961570 -6.166273 13.92093 -9.114950 0.002094648 0.9999991
## 3 3.2 -1.910826 -6.036420 13.68869 -8.963753 0.002384402 0.9999989
## psiBa psiB psiAandB
## 1 9.456337e-05 0.001934427 0.001840038
## 2 1.099968e-04 0.002204412 0.002094646
## 3 1.279487e-04 0.002512043 0.002384399
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)
# 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)
mod.psi.temp$real$psiBA [1,]
## est se lower upper
## unit1 0.9977 0.0042 0.9241 0.9999
# Pr(Occupancy of B | A absent)
mod.psi.temp$real$psiBa[1,]
## est se lower upper
## unit1 0.364 0.0674 0.2445 0.5031
# We can compute the marginal estimate of occupancy of B
psiB <- mod.psi.temp$real$psiA * mod.psi.temp$real$psiBA+
(1-mod.psi.temp$real$psiA) * mod.psi.temp$real$psiBa
psiB[1,]
## est se lower upper
## unit1 0.3750264 0.06678696 0.2484417 0.5286355
# look at odds ratio of occupancy of B|A to 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$derived$nu[1,]
## est se lower upper
## unit1 757.0836 1393.534 20.52775 27921.98
# For example, the marginal estimates are and prob if indendent are
c(mod.psi.temp$real$psiA[1,1],psiB[1,1], mod.psi.temp$real$psiA[1,1]*psiB[1,1])
## [1] 0.017400000 0.375026380 0.006525459
# Actual estimation of occupancy of both species
mod.psi.temp$real$psiBA [1,1]*mod.psi.temp$real$psiA[1,1]
## [1] 0.01735998
# Probability of detection of species A if alone
mod.psi.temp$real$pA[1,]
## est se lower upper
## unit1_1-1 0.9053 0.0179 0.8638 0.9351
# Probability of detection of species B if alone
mod.psi.temp$real$pB[1,]
## est se lower upper
## unit1_1-1 0.9053 0.0179 0.8638 0.9351
# Probability of detection when both species present
# rA - just A; rBA B give A; rBa B given A not detected
# Again, if species do not interact, the
# rBA = rBa
mod.psi.temp$real$rA[1,]
## est se lower upper
## unit1_1-1 0.9053 0.0179 0.8638 0.9351
mod.psi.temp$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.9053 0.0179 0.8638 0.9351
mod.psi.temp$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.9053 0.0179 0.8638 0.9351
# 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
mod.psi.temp$derived$rho[1,]
## est se lower upper
## unit1_1-1 1 0 1 1
# Model averaging
models<-list(mod1,
mod1.occind,
mod1.psimp,
mod.psi.temp)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar
## 1 psi(temp.std P SP P SP T temp.std P INT P INT T temp.std)p() 0.00 1 7
## 2 psi(SP P INT)p() 128.03 0 4
## 3 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d) 133.33 0 8
## 4 psi(SP)p() 135.52 0 3
## neg2ll warn.conv warn.VC
## 1 522.47 7.00 0
## 2 656.50 7.00 0
## 3 653.80 5.88 0
## 4 665.99 7.00 0
RPresence::modAvg(results, param="psiA")[1,]
## est se lower_0.95 upper_0.95
## unit1 0.0174 0.0097 0.005790647 0.05108831