# Multi species single season
# 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).
# 2020-06-24 CJS real$nu real$rho are now in derived$nu and derived$rho
# 2018-09-27 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library(ggplot2)
library(readxl)
library(reshape2)
library(RPresence)
# get the data
PG.data <- readxl::read_excel(file.path("..","Salamander2.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "B3:F90")
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
PJ.data <- readxl::read_excel(file.path("..","Salamander2.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "H3:L90")
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
# Convert history to compressed format
# 0=neither species detected,
# 1=only species A detected,
# 2=only species B detected,
# 3=both species detected
input.history <- PG.data + 2*PJ.data
input.history
## ...1 ...2 ...3 ...4 ...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
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$Std..Elevation <- (site.covar$Elevation..m. - elevation.mean)/elevation.std
Nsites <- nrow(input.history)
Nvisits<- ncol(input.history)
# create the pao file
salamander.pao <- createPao(input.history,
unitcov=site.covar,
title="Salamander multi species - co-occurance")
summary(salamander.pao)
## paoname=pres.pao
## title=Salamander multi species - co-occurance
## Naive occ=0.8636364
## naiveR =0.4204545
## nunits nsurveys nseasons nsurveyseason nmethods
## "88" "5" "1" "5" "1"
## nunitcov nsurvcov
## "2" "1"
## unit covariates : Elevation..m. Std..Elevation
## survey covariates: SURVEY
# 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
#
# The following special variables are available for modelling PSI
# SP species effect
# INT interaction of effect on presence of species B when species A was, or was not present
#
# Model for PSI.... impact on parameters
# psi~1 psiA=psiBA=psiBa (1 parameter)
# psi~SP psiA psiBA=psiBa (2 parameters)
# psi~SP+INT psiA psiBA psiBa (3 parameters)
#
# The following special variables are available for p
# SP species effect
# INT_o is a detection-level interaction where the OCCUPANCY of one species
# changes the detection probability of the other species
# INT_d is a detection-level interaction where DETECTION of one species changes the
# detection probability of the other species in the same survey.
#
# Model for p.... impact on parameters
# p~1 pA = pB = rA = rBA = rBa (1 parameter)
# p~SP pA=rA, pB=rBA=rBa (2 parameters)
# p~SP+INT_o pA=rA, pB, rBA=rBa (3 parameters)
# p~SP+INT_o+SP:INT_o pA, rA, pB, rBA=rBa (4 parameters)
# p~SP+INT_o+SP:INT_o+INT_d pA, rA, pB, rBA, rBa (5 parameters)
mod1 <- occMod(model=list(psi~SP+INT, p~SP+INT_o+SP:INT_o + INT_d),
data=salamander.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=751.8665
## -2*log-likelihood=735.8665
## num. par=8
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 7 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.5719 0.0572 0.458 0.6787
# 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.3177 0.0716 0.1959 0.4709
# Pr(Occupancy of B | AA absent)
mod1$real$psiBa[1,]
## est se lower upper
## unit1 0.7005 0.0778 0.5307 0.8287
# 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.4815767 0.07744536 0.3773616 0.5858611
# 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 0.1991038 0.1027362 0.07241979 0.5473963
# 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.5719000 0.4815767 0.2754137
# Actual estimation of occupancy of both species
mod1$real$psiBA [1,1]*mod1$real$psiA[1,1]
## [1] 0.1816926
# Probability of detection of species A if alone
mod1$real$pA[1,]
## est se lower upper
## unit1_1-1 0.5396 0.0416 0.4576 0.6194
# Probability of detection of species B if alone
mod1$real$pB[1,]
## est se lower upper
## unit1_1-1 0.9027 0.0391 0.795 0.9569
# 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$real$rA[1,]
## est se lower upper
## unit1_1-1 0.4882 0.0766 0.3433 0.635
mod1$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.5014 0.0847 0.341 0.6614
mod1$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.6067 0.092 0.4202 0.7665
# 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 0.651934 0.3272929 0.2437046 1.743988
# 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~SP+INT_o+INT_d+SP:INT_o),
data=salamander.pao,
type="so.2sp.1")
summary(mod1.occind)
## Model name=psi(SP)p(SP P INT_o P INT_d P SP T INT_o)
## AIC=760.3286
## -2*log-likelihood=746.3286
## num. par=7
## 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.5949 0.057 0.4801 0.7002
# 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.4867 0.0546 0.3819 0.5926
# Pr(Occupancy of B | AA absent)
mod1.occind$real$psiBa[1,]
## est se lower upper
## unit1 0.4867 0.0546 0.3819 0.5926
# 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.4867 0.0546 0.3819 0.5926
# 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.5949000 0.4867000 0.2895378
# Actual estimation of occupancy of both species
mod1.occind$real$psiBA [1,1]*mod1.occind$real$psiA[1,1]
## [1] 0.2895378
# Probability of detection of species A if alone
mod1.occind$real$pA[1,]
## est se lower upper
## unit1_1-1 0.5478 0.0409 0.4671 0.626
# Probability of detection of species B if alone
mod1.occind$real$pB[1,]
## est se lower upper
## unit1_1-1 0.9358 0.0318 0.8376 0.9763
# 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.4273 0.0682 0.3018 0.563
mod1.occind$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.4915 0.0861 0.3298 0.655
mod1.occind$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.5819 0.0817 0.4188 0.7289
# 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 0.6945563 0.3157696 0.2849144 1.69317
# is there any support for the indpendent occupancy model?
models<-list(mod1,
mod1.occind)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll
## 1 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d) 0.00 0.986 8 735.87
## 2 psi(SP)p(SP P INT_o P INT_d P SP T INT_o) 8.46 0.014 7 746.33
## warn.conv warn.VC
## 1 7 0
## 2 7 0
# mod1.detind is same as model 1, but we assume that detection probabilities
# of B do not depend if A has been detected or not i.e. r(B|A)=r(B|a)
mod1.detind <- occMod(model=list(psi~SP+INT, p~SP+INT_o+SP:INT_o),
data=salamander.pao,
type="so.2sp.1")
summary(mod1.detind)
## Model name=psi(SP P INT)p(SP P INT_o P SP T INT_o)
## AIC=750.624
## -2*log-likelihood=736.624
## num. par=7
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 7 signifcant digits.
mod1.detind$dmat
## $psi
## a1 a2 a3
## psiA "1" "0" "0"
## psiBA "1" "1" "0"
## psiBa "1" "1" "1"
##
## $p
## b1 b2 b3 b4
## pA[1] "1" "0" "0" "0"
## pA[2] "1" "0" "0" "0"
## pA[3] "1" "0" "0" "0"
## pA[4] "1" "0" "0" "0"
## pA[5] "1" "0" "0" "0"
## pB[1] "1" "1" "0" "0"
## pB[2] "1" "1" "0" "0"
## pB[3] "1" "1" "0" "0"
## pB[4] "1" "1" "0" "0"
## pB[5] "1" "1" "0" "0"
## rA[1] "1" "0" "1" "0"
## rA[2] "1" "0" "1" "0"
## rA[3] "1" "0" "1" "0"
## rA[4] "1" "0" "1" "0"
## rA[5] "1" "0" "1" "0"
## rBA[1] "1" "1" "1" "1"
## rBA[2] "1" "1" "1" "1"
## rBA[3] "1" "1" "1" "1"
## rBA[4] "1" "1" "1" "1"
## rBA[5] "1" "1" "1" "1"
## rBa[1] "1" "1" "1" "1"
## rBa[2] "1" "1" "1" "1"
## rBa[3] "1" "1" "1" "1"
## rBa[4] "1" "1" "1" "1"
## rBa[5] "1" "1" "1" "1"
names(mod1.detind$real)
## [1] "psiA" "psiBA" "psiBa" "pA" "pB" "rA" "rBA" "rBa"
# Probability of occupancy of A
mod1.detind$real$psiA [1,]
## est se lower upper
## unit1 0.5754 0.0568 0.4621 0.6813
# 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.detind$real$psiBA [1,]
## est se lower upper
## unit1 0.3224 0.0712 0.2007 0.4741
# Pr(Occupancy of B | AA absent)
mod1.detind$real$psiBa[1,]
## est se lower upper
## unit1 0.6979 0.0782 0.5277 0.8269
# We can compute the marginal estimate of occupancy of B
psiB <- mod1.detind$real$psiA * mod1.detind$real$psiBA+
(1-mod1.detind$real$psiA) * mod1.detind$real$psiBa
psiB[1,]
## est se lower upper
## unit1 0.4818373 0.0778024 0.3765933 0.5865374
# 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.detind$derived$nu[1,]
## est se lower upper
## unit1 0.205909 0.1053035 0.07557117 0.5610408
# For example, the marginal estimates are and prob if indendent are
c(mod1.detind$real$psiA[1,1],psiB[1,1], mod1.detind$real$psiA[1,1]*psiB[1,1])
## [1] 0.5754000 0.4818373 0.2772492
# Actual estimation of occupancy of both species
mod1.detind$real$psiBA [1,1]*mod1.detind$real$psiA[1,1]
## [1] 0.185509
# Probability of detection of species A if alone
mod1.detind$real$pA[1,]
## est se lower upper
## unit1_1-1 0.5406 0.0416 0.4587 0.6204
# Probability of detection of species B if alone
mod1.detind$real$pB[1,]
## est se lower upper
## unit1_1-1 0.9092 0.0372 0.8053 0.9604
# 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.detind$real$rA[1,]
## est se lower upper
## unit1_1-1 0.4776 0.073 0.3401 0.6186
mod1.detind$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.5504 0.062 0.4283 0.6666
mod1.detind$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.5504 0.062 0.4283 0.6666
# 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.detind$derived$rho[1,]
## est se lower upper
## unit1_1-1 1 0 1 1
# Now do model averaging
models<-list(mod1,
mod1.occind,
mod1.detind)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll
## 1 psi(SP P INT)p(SP P INT_o P SP T INT_o) 0.00 0.6472 7 736.62
## 2 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d) 1.24 0.3477 8 735.87
## 3 psi(SP)p(SP P INT_o P INT_d P SP T INT_o) 9.70 0.0051 7 746.33
## warn.conv warn.VC
## 1 7 0
## 2 7 0
## 3 7 0
# Look at the effect of elevtion
mod.psi.elevation <- occMod(model=list(
psi~Std..Elevation+SP+SP:Std..Elevation+INT+INT:Std..Elevation,
p ~SP+INT_o+SP:INT_o),
data=salamander.pao,
type="so.2sp.1")
summary(mod.psi.elevation)
## Model name=psi(Std..Elevation P SP P SP T Std..Elevation P INT P INT T Std..Elevation)p(SP P INT_o P SP T INT_o)
## AIC=687.4714
## -2*log-likelihood=667.4714
## num. par=10
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 7 signifcant digits.
names(mod.psi.elevation$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)
# Ww need to look at the beta values
mod.psi.elevation$beta$psi
## est se
## A1_psiA 0.372485 0.244388
## A2_psiA.Std..Elevation_psiA -0.686944 0.275267
## A3_psiBA -0.984296 0.451688
## A4_psiBa 4.503822 4.326472
## A5_psiBA.Std..Elevation:SP2_psiBA 3.142486 0.810764
## A6_psiBa.Std..Elevation:INT2_psiBa 16.043129 18.128544
str(mod.psi.elevation$beta$psi)
## 'data.frame': 6 obs. of 2 variables:
## $ est: num 0.372 -0.687 -0.984 4.504 3.142 ...
## $ se : num 0.244 0.275 0.452 4.326 0.811 ...
mod.psi.elevation$dmat$psi
## a1 a2 a3 a4 a5
## psiA "1" "Std..Elevation_psiA" "0" "0" "0"
## psiBA "1" "Std..Elevation_psiBA" "1" "0" "Std..Elevation:SP2_psiBA"
## psiBa "1" "Std..Elevation_psiBa" "1" "1" "Std..Elevation:SP2_psiBa"
## a6
## psiA "0"
## psiBA "0"
## psiBa "Std..Elevation:INT2_psiBa"
range(site.covar$Elevation..m.)
## [1] 502.920 1182.624
predictPSI <- data.frame(Elevation..m.=seq(400,1200,10))
predictPSI$Std..Elevation <- (predictPSI$Elevation..m. - elevation.mean)/elevation.std
predictPSI$psiA.logit <- mod.psi.elevation$beta$psi[1,]$est + mod.psi.elevation$beta$psi[2,]$est*predictPSI$Std..Elevation
predictPSI$psiBA.logit <- mod.psi.elevation$beta$psi[3,]$est + mod.psi.elevation$beta$psi[5,]$est*predictPSI$Std..Elevation
predictPSI$psiBa.logit <- mod.psi.elevation$beta$psi[4,]$est + mod.psi.elevation$beta$psi[6,]$est*predictPSI$Std..Elevation
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,]
## Elevation..m. Std..Elevation psiA.logit psiBA.logit psiBa.logit psiA
## 1 400 -2.801258 2.296792 -9.787209 -40.43711 0.9086110
## 2 410 -2.737862 2.253243 -9.587990 -39.42005 0.9049299
## 3 420 -2.674467 2.209694 -9.388771 -38.40300 0.9011167
## psiBA psiBa psiB psiAandB
## 1 5.616231e-05 2.744004e-18 5.102969e-05 5.102969e-05
## 2 6.854239e-05 7.587310e-18 6.202606e-05 6.202606e-05
## 3 8.365123e-05 2.097930e-17 7.537952e-05 7.537952e-05
plotdata <- reshape2::melt(predictPSI,
id.vars=c("Elevation..m.","Std..Elevation"),
measure.vars=c("psiA","psiB","psiAandB"),
variable.name="Species",
value.name='P.Occupancy')
species.by.elev <- ggplot2::ggplot(dat=plotdata, aes(x=Elevation..m., y=P.Occupancy, color=Species))+
ggtitle("Occupancy as a function of elevation")+
geom_line()
species.by.elev
ggsave(plot=species.by.elev, file='species-by-elev.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.elevation$real$psiBA [1,]
## est se lower upper
## unit1 0.0028 0.0048 1e-04 0.0766
# Pr(Occupancy of B | AA absent)
mod.psi.elevation$real$psiBa[1,]
## est se lower upper
## unit1 0 0 0 1
# We can compute the marginal estimate of occupancy of B
psiB <- mod.psi.elevation$real$psiA$est * mod.psi.elevation$real$psiBA$est+
(1-mod.psi.elevation$real$psiA$est) * mod.psi.elevation$real$psiBa$est
psiB
## [1] 0.00241920 0.00299600 0.00299600 0.00299600 0.00372944 0.00469728
## [7] 0.00738090 0.00923514 0.01151424 0.01435616 0.01780475 0.02735931
## [13] 0.02735931 0.02735931 0.03242532 0.03673139 0.03828870 0.04511520
## [19] 0.05305890 0.07281200 0.07281200 0.07566261 0.10071844 0.10071844
## [25] 0.10071844 0.10071844 0.11091562 0.11722969 0.12464875 0.12464875
## [31] 0.17561024 0.17561024 0.22295775 0.22295775 0.28904230 0.32743526
## [37] 0.32743526 0.36707971 0.47475504 0.50297926 0.58032601 0.58032601
## [43] 0.60563243 0.60563243 0.60563243 0.62700797 0.64628222 0.64628222
## [49] 0.68191548 0.68191548 0.77153365 0.77153365 0.77907214 0.77907214
## [55] 0.77907214 0.80794630 0.80794630 0.80794630 0.83455273 0.86984140
## [61] 0.88030888 0.89485984 0.89485984 0.89938798 0.89938798 0.89938798
## [67] 0.90798896 0.91600000 0.91600000 0.91600000 0.92341710 0.92341710
## [73] 0.93028553 0.93028553 0.93663964 0.94245284 0.94245284 0.96854480
## [79] 0.96854480 0.97933312 0.98332373 0.98332373 0.98657452 0.98921180
## [85] 0.28904230 0.99306307 0.99646372 0.99774684
# 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.elevation$derived$nu[1,]
## est se lower upper
## unit1 1.036589e+13 3.622326e+14 1.862546e-17 5.769077e+42
# For example, the marginal estimates are and prob if indendent are
cbind(mod.psi.elevation$real$psiA$est, psiB, mod.psi.elevation$real$psiA$est*psiB)
## psiB
## [1,] 0.8640 0.00241920 0.002090189
## [2,] 0.8560 0.00299600 0.002564576
## [3,] 0.8560 0.00299600 0.002564576
## [4,] 0.8560 0.00299600 0.002564576
## [5,] 0.8476 0.00372944 0.003161073
## [6,] 0.8388 0.00469728 0.003940078
## [7,] 0.8201 0.00738090 0.006053076
## [8,] 0.8101 0.00923514 0.007481387
## [9,] 0.7996 0.01151424 0.009206786
## [10,] 0.7888 0.01435616 0.011324139
## [11,] 0.7775 0.01780475 0.013843193
## [12,] 0.7537 0.02735931 0.020620712
## [13,] 0.7537 0.02735931 0.020620712
## [14,] 0.7537 0.02735931 0.020620712
## [15,] 0.7437 0.03242532 0.024114710
## [16,] 0.7361 0.03673139 0.027037976
## [17,] 0.7335 0.03828870 0.028084761
## [18,] 0.7230 0.04511520 0.032618290
## [19,] 0.7122 0.05305890 0.037788549
## [20,] 0.6900 0.07281200 0.050240280
## [21,] 0.6900 0.07281200 0.050240280
## [22,] 0.6871 0.07566261 0.051987779
## [23,] 0.6668 0.10071844 0.067159056
## [24,] 0.6668 0.10071844 0.067159056
## [25,] 0.6668 0.10071844 0.067159056
## [26,] 0.6668 0.10071844 0.067159056
## [27,] 0.6609 0.11091562 0.073304133
## [28,] 0.6579 0.11722969 0.077125413
## [29,] 0.6549 0.12464875 0.081632466
## [30,] 0.6549 0.12464875 0.081632466
## [31,] 0.6428 0.17561024 0.112882262
## [32,] 0.6428 0.17561024 0.112882262
## [33,] 0.6367 0.22295775 0.141957199
## [34,] 0.6367 0.22295775 0.141957199
## [35,] 0.6305 0.28904230 0.182241170
## [36,] 0.6274 0.32743526 0.205432882
## [37,] 0.6274 0.32743526 0.205432882
## [38,] 0.6243 0.36707971 0.229167863
## [39,] 0.6149 0.47475504 0.291926874
## [40,] 0.6118 0.50297926 0.307722711
## [41,] 0.5991 0.58032601 0.347673313
## [42,] 0.5991 0.58032601 0.347673313
## [43,] 0.5927 0.60563243 0.358958341
## [44,] 0.5927 0.60563243 0.358958341
## [45,] 0.5927 0.60563243 0.358958341
## [46,] 0.5863 0.62700797 0.367614773
## [47,] 0.5798 0.64628222 0.374714431
## [48,] 0.5798 0.64628222 0.374714431
## [49,] 0.5668 0.68191548 0.386509694
## [50,] 0.5668 0.68191548 0.386509694
## [51,] 0.5307 0.77153365 0.409452908
## [52,] 0.5307 0.77153365 0.409452908
## [53,] 0.5274 0.77907214 0.410882647
## [54,] 0.5274 0.77907214 0.410882647
## [55,] 0.5274 0.77907214 0.410882647
## [56,] 0.5142 0.80794630 0.415445987
## [57,] 0.5142 0.80794630 0.415445987
## [58,] 0.5142 0.80794630 0.415445987
## [59,] 0.5009 0.83455273 0.418027462
## [60,] 0.4810 0.86984140 0.418393713
## [61,] 0.4744 0.88030888 0.417618533
## [62,] 0.4644 0.89485984 0.415572910
## [63,] 0.4644 0.89485984 0.415572910
## [64,] 0.4611 0.89938798 0.414707798
## [65,] 0.4611 0.89938798 0.414707798
## [66,] 0.4611 0.89938798 0.414707798
## [67,] 0.4546 0.90798896 0.412771781
## [68,] 0.4480 0.91600000 0.410368000
## [69,] 0.4480 0.91600000 0.410368000
## [70,] 0.4480 0.91600000 0.410368000
## [71,] 0.4414 0.92341710 0.407596308
## [72,] 0.4414 0.92341710 0.407596308
## [73,] 0.4349 0.93028553 0.404581177
## [74,] 0.4349 0.93028553 0.404581177
## [75,] 0.4284 0.93663964 0.401256422
## [76,] 0.4219 0.94245284 0.397620853
## [77,] 0.4219 0.94245284 0.397620853
## [78,] 0.3836 0.96854480 0.371533785
## [79,] 0.3836 0.96854480 0.371533785
## [80,] 0.3588 0.97933312 0.351384723
## [81,] 0.3467 0.98332373 0.340918337
## [82,] 0.3467 0.98332373 0.340918337
## [83,] 0.3348 0.98657452 0.330305149
## [84,] 0.3230 0.98921180 0.319515411
## [85,] 0.6305 0.28904230 0.182241170
## [86,] 0.3003 0.99306307 0.298216840
## [87,] 0.2679 0.99646372 0.266952631
## [88,] 0.2476 0.99774684 0.247042118
# Actual estimation of occupancy of both species
mod.psi.elevation$real$psiBA$est *mod.psi.elevation$real$psiA$est
## [1] 0.00241920 0.00299600 0.00299600 0.00299600 0.00372944 0.00469728
## [7] 0.00738090 0.00923514 0.01151424 0.01435616 0.01780475 0.02735931
## [13] 0.02735931 0.02735931 0.03242532 0.03673139 0.03828870 0.04511520
## [19] 0.05305890 0.07265700 0.07265700 0.07544358 0.09788624 0.09788624
## [25] 0.09788624 0.09788624 0.10508310 0.10888245 0.11270829 0.11270829
## [31] 0.12913852 0.12913852 0.13790922 0.13790922 0.14696955 0.15164258
## [37] 0.15164258 0.15638715 0.17106518 0.17607604 0.19662462 0.19662462
## [43] 0.20708938 0.20708938 0.20708938 0.21769319 0.22826726 0.22826726
## [49] 0.24927864 0.24927864 0.30223365 0.30223365 0.30647214 0.30647214
## [55] 0.30647214 0.32214630 0.32214630 0.32214630 0.33545273 0.35084140
## [61] 0.35470888 0.35925984 0.35925984 0.36048798 0.36048798 0.36048798
## [67] 0.36258896 0.36400000 0.36400000 0.36400000 0.36481710 0.36481710
## [73] 0.36518553 0.36518553 0.36503964 0.36435284 0.36435284 0.35214480
## [79] 0.35214480 0.33813312 0.33002373 0.33002373 0.32137452 0.31221180
## [85] 0.14696955 0.29336307 0.26436372 0.24534684
# Probability of detection of species A if alone
mod.psi.elevation$real$pA[1,]
## est se lower upper
## unit1_1-1 0.5387 0.0419 0.4564 0.619
# Probability of detection of species B if alone
mod.psi.elevation$real$pB[1,]
## est se lower upper
## unit1_1-1 0.9199 0.0318 0.8315 0.9639
# 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.elevation$real$rA[1,]
## est se lower upper
## unit1_1-1 0.4645 0.0653 0.3416 0.5919
mod.psi.elevation$real$rBA[1,]
## est se lower upper
## unit1_1-1 0.5483 0.0593 0.4315 0.66
mod.psi.elevation$real$rBa[1,]
## est se lower upper
## unit1_1-1 0.5483 0.0593 0.4315 0.66
# 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.elevation$derived$rho[1,]
## est se lower upper
## unit1_1-1 1 0 1 1
# Model averaging
models<-list(mod1,
mod1.occind,
mod1.detind,
mod.psi.elevation)
results<-RPresence::createAicTable(models)
summary(results)
## Model
## 1 psi(Std..Elevation P SP P SP T Std..Elevation P INT P INT T Std..Elevation)p(SP P INT_o P SP T INT_o)
## 2 psi(SP P INT)p(SP P INT_o P SP T INT_o)
## 3 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d)
## 4 psi(SP)p(SP P INT_o P INT_d P SP T INT_o)
## DAIC wgt npar neg2ll warn.conv warn.VC
## 1 0.00 1 10 667.47 7 0
## 2 63.15 0 7 736.62 7 0
## 3 64.40 0 8 735.87 7 0
## 4 72.86 0 7 746.33 7 0
RPresence::modAvg(results, param="psiA")[1,]
## est se lower_0.95 upper_0.95
## unit1 0.864 0.0791 0.6293739 0.9596239