# Multi species single season with a list of models
# 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 and real$rho 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)
# define the list of models to fit
# Notice the commas between the column and the placement of the quotes
model.list.csv <- textConnection("
p, psi
~SP+INT_o+SP:INT_o + INT_d, ~SP+INT
~SP+INT_o+SP:INT_o + INT_d, ~SP
~SP+INT_o+SP:INT_o , ~SP+INT
~SP+INT_o+SP:INT_o, ~SP + Std..Elevation+SP:Std..Elevation+INT+INT:Std..Elevation")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p
## 1 ~SP+INT_o+SP:INT_o + INT_d
## 2 ~SP+INT_o+SP:INT_o + INT_d
## 3 ~SP+INT_o+SP:INT_o
## 4 ~SP+INT_o+SP:INT_o
## psi
## 1 ~SP+INT
## 2 ~SP
## 3 ~SP+INT
## 4 ~SP + Std..Elevation+SP:Std..Elevation+INT+INT:Std..Elevation
# fit the models
model.fits <- plyr::alply(model.list, 1, function(x,detect.pao){
cat("\n\n***** Starting ", unlist(x), "\n")
fit <- RPresence::occMod(model=list(as.formula(paste("psi",x$psi)),
as.formula(paste("p" ,x$p ))),
data=detect.pao,type="so.2sp.1", randint=10)
fit
},detect.pao=salamander.pao)
##
##
## ***** Starting ~SP+INT_o+SP:INT_o + INT_d ~SP+INT
##
##
## ***** Starting ~SP+INT_o+SP:INT_o + INT_d ~SP
##
##
## ***** Starting ~SP+INT_o+SP:INT_o ~SP+INT
##
##
## ***** Starting ~SP+INT_o+SP:INT_o ~SP + Std..Elevation+SP:Std..Elevation+INT+INT:Std..Elevation
# Look the output from a specific model
check.model <- 4
names(model.fits[[check.model]])
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "warnings"
model.fits[[check.model]]$beta
## $psi
## est se
## A1_psiA 0.372485 0.244388
## A2_psiBA -0.984296 0.451688
## A3_psiA.Std..Elevation_psiA -0.686944 0.275267
## A4_psiBa 4.503822 4.326472
## A5_psiBA.SP2:Std..Elevation_psiBA 3.142486 0.810764
## A6_psiBa.Std..Elevation:INT2_psiBa 16.043128 18.128542
##
## $psi.VC
## A1 A2 A3 A4 A5 A6
## A1 0.059726 -0.056538 -0.010194 0.046398 0.012918 0.291997
## A2 -0.056538 0.204022 0.011789 -0.132702 -0.050124 0.190452
## A3 -0.010194 0.011789 0.075772 -0.001861 -0.073781 0.014501
## A4 0.046398 -0.132702 -0.001861 18.718359 -0.014198 73.306481
## A5 0.012918 -0.050124 -0.073781 -0.014198 0.657338 -0.966796
## A6 0.291997 0.190452 0.014501 73.306481 -0.966796 328.644030
##
## $p
## est se
## B1_pA[1] 0.155253 0.168438
## B2_pB[1] 2.285572 0.461505
## B3_rA[1] -0.297429 0.315181
## B4_rBA[1] -1.949537 0.536616
##
## $p.VC
## B1 B2 B3 B4
## B1 0.028371 -0.027789 -0.029443 0.029098
## B2 -0.027789 0.212987 -0.004657 -0.186931
## B3 -0.029443 -0.004657 0.099339 -0.065900
## B4 0.029098 -0.186931 -0.065900 0.287957
##
## $VC
## A1 A2 A3 A4 A5 A6 B1
## A1 0.059726 -0.056538 -0.010194 0.046398 0.012918 0.291997 -0.003685
## A2 -0.056538 0.204022 0.011789 -0.132702 -0.050124 0.190452 0.004126
## A3 -0.010194 0.011789 0.075772 -0.001861 -0.073781 0.014501 0.005789
## A4 0.046398 -0.132702 -0.001861 18.718359 -0.014198 73.306481 -0.008110
## A5 0.012918 -0.050124 -0.073781 -0.014198 0.657338 -0.966796 -0.008045
## A6 0.291997 0.190452 0.014501 73.306481 -0.966796 328.644030 -0.008416
## B1 -0.003685 0.004126 0.005789 -0.008110 -0.008045 -0.008416 0.028371
## B2 0.019087 -0.003256 0.004419 0.202552 0.002082 1.379862 -0.027789
## B3 -0.004297 -0.002741 -0.012126 -0.075764 0.013603 -0.589965 -0.029443
## B4 -0.011609 -0.004460 0.002432 -0.135966 -0.010059 -0.929339 0.029098
## B2 B3 B4
## A1 0.019087 -0.004297 -0.011609
## A2 -0.003256 -0.002741 -0.004460
## A3 0.004419 -0.012126 0.002432
## A4 0.202552 -0.075764 -0.135966
## A5 0.002082 0.013603 -0.010059
## A6 1.379862 -0.589965 -0.929339
## B1 -0.027789 -0.029443 0.029098
## B2 0.212987 -0.004657 -0.186931
## B3 -0.004657 0.099339 -0.065900
## B4 -0.186931 -0.065900 0.287957
model.fits[[check.model]]$dmat
## $psi
## a1 a2 a3 a4 a5
## psiA "1" "0" "Std..Elevation_psiA" "0" "0"
## psiBA "1" "1" "Std..Elevation_psiBA" "0" "SP2:Std..Elevation_psiBA"
## psiBa "1" "1" "Std..Elevation_psiBa" "1" "SP2:Std..Elevation_psiBa"
## a6
## psiA "0"
## psiBA "0"
## psiBa "Std..Elevation:INT2_psiBa"
##
## $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(model.fits[[check.model]]$real)
## [1] "psiA" "psiBA" "psiBa" "pA" "pB" "rA" "rBA" "rBa"
model.fits[[check.model]]$real$psiA[1:5,]
## est se lower upper
## unit1 0.8640 0.0791 0.6292 0.9596
## unit2 0.8560 0.0800 0.6250 0.9549
## unit3 0.8560 0.0800 0.6250 0.9549
## unit4 0.8560 0.0800 0.6250 0.9549
## unit5 0.8476 0.0806 0.6208 0.9497
model.fits[[check.model]]$real$psiBA[1:5,]
## est se lower upper
## unit1 0.0028 0.0048 1e-04 0.0766
## unit2 0.0035 0.0058 1e-04 0.0837
## unit3 0.0035 0.0058 1e-04 0.0837
## unit4 0.0035 0.0058 1e-04 0.0837
## unit5 0.0044 0.0070 2e-04 0.0914
model.fits[[check.model]]$real$psiBa[1:5,]
## est se lower upper
## unit1 0 0 0 1
## unit2 0 0 0 1
## unit3 0 0 0 1
## unit4 0 0 0 1
## unit5 0 0 0 1
model.fits[[check.model]]$derived$nu[1:5,]
## est se lower upper
## unit1 1.036587e+13 3.622318e+14 1.862559e-17 5.769013e+42
## unit2 2.200140e+12 7.303357e+13 1.219893e-16 3.968065e+40
## unit3 2.200140e+12 7.303357e+13 1.219893e-16 3.968065e+40
## unit4 2.200140e+12 7.303357e+13 1.219893e-16 3.968065e+40
## unit5 4.669762e+11 1.468429e+13 7.986638e-16 2.730395e+38
model.fits[[check.model]]$real$pA[1:5,]
## est se lower upper
## unit1_1-1 0.5387 0.0419 0.4564 0.619
## unit2_1-1 0.5387 0.0419 0.4564 0.619
## unit3_1-1 0.5387 0.0419 0.4564 0.619
## unit4_1-1 0.5387 0.0419 0.4564 0.619
## unit5_1-1 0.5387 0.0419 0.4564 0.619
model.fits[[check.model]]$real$pB[1:5,]
## est se lower upper
## unit1_1-1 0.9199 0.0318 0.8315 0.9639
## unit2_1-1 0.9199 0.0318 0.8315 0.9639
## unit3_1-1 0.9199 0.0318 0.8315 0.9639
## unit4_1-1 0.9199 0.0318 0.8315 0.9639
## unit5_1-1 0.9199 0.0318 0.8315 0.9639
model.fits[[check.model]]$real$rA[1:5,]
## est se lower upper
## unit1_1-1 0.4645 0.0653 0.3416 0.5919
## unit2_1-1 0.4645 0.0653 0.3416 0.5919
## unit3_1-1 0.4645 0.0653 0.3416 0.5919
## unit4_1-1 0.4645 0.0653 0.3416 0.5919
## unit5_1-1 0.4645 0.0653 0.3416 0.5919
model.fits[[check.model]]$real$rBA[1:5,]
## est se lower upper
## unit1_1-1 0.5483 0.0593 0.4315 0.66
## unit2_1-1 0.5483 0.0593 0.4315 0.66
## unit3_1-1 0.5483 0.0593 0.4315 0.66
## unit4_1-1 0.5483 0.0593 0.4315 0.66
## unit5_1-1 0.5483 0.0593 0.4315 0.66
model.fits[[check.model]]$real$rBa[1:5,]
## est se lower upper
## unit1_1-1 0.5483 0.0593 0.4315 0.66
## unit2_1-1 0.5483 0.0593 0.4315 0.66
## unit3_1-1 0.5483 0.0593 0.4315 0.66
## unit4_1-1 0.5483 0.0593 0.4315 0.66
## unit5_1-1 0.5483 0.0593 0.4315 0.66
model.fits[[check.model]]$derived$rho[1:5,]
## est se lower upper
## unit1_1-1 1 0 1 1
## unit2_1-1 1 0 1 1
## unit3_1-1 1 0 1 1
## unit4_1-1 1 0 1 1
## unit5_1-1 1 0 1 1
names(model.fits[[check.model]]$derived)
## [1] "nu" "rho"
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model
## 4 psi(SP P Std..Elevation P SP T Std..Elevation P INT P INT T Std..Elevation)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)
## 1 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d)
## 2 psi(SP)p(SP P INT_o P SP T INT_o P INT_d)
## AIC neg2ll npar warn.conv warn.VC DAIC modlike wgt
## 4 687.4714 667.4714 10 7 0 0.0000 1 1
## 3 750.6240 736.6240 7 7 0 63.1526 0 0
## 1 751.8665 735.8665 8 7 0 64.3951 0 0
## 2 760.3286 746.3286 7 7 0 72.8572 0 0
names(aic.table)
## [1] "table" "models" "ess"
ma.psiA<- RPresence::modAvg(aic.table, param="psiA")
ma.psiA$parameter <- "psiA"
ma.psiA <- cbind(ma.psiA, site.covar)
ma.psiBA<- RPresence::modAvg(aic.table, param="psiBA")
ma.psiBA$parameter <- "psiBA"
ma.psiBA <- cbind(ma.psiBA, site.covar)
ma.psiBa<- RPresence::modAvg(aic.table, param="psiBa")
ma.psiBa$parameter <- "psiBa"
ma.psiBa <- cbind(ma.psiBa, site.covar)
ma.psiB <- site.covar
ma.psiB$parameter <- 'psiB'
ma.psiB$est <- ma.psiBA$est * ma.psiA$est +
ma.psiBa$est * (1-ma.psiA$est)
# too hard to compute the se
psi.all <- plyr::rbind.fill(ma.psiA, ma.psiBA, ma.psiBa, ma.psiB)
head(psi.all)
## est se lower_0.95 upper_0.95 parameter Elevation..m. Std..Elevation
## 1 0.8640 0.0791 0.6293739 0.9596239 psiA 502.92 -2.148793
## 2 0.8560 0.0800 0.6249035 0.9549764 psiA 518.16 -2.052179
## 3 0.8560 0.0800 0.6249035 0.9549764 psiA 518.16 -2.052179
## 4 0.8560 0.0800 0.6249035 0.9549764 psiA 518.16 -2.052179
## 5 0.8476 0.0806 0.6208022 0.9497337 psiA 533.40 -1.955564
## 6 0.8388 0.0811 0.6162744 0.9440060 psiA 548.64 -1.858950
ggplot(data=psi.all, aes(x=Elevation..m., y=est))+
ggtitle("Estimated occupancy parameters vs. elevation")+
geom_point()+
geom_ribbon( aes(ymin=lower_0.95, ymax=upper_0.95),alpha=0.2)+
facet_wrap(~parameter, ncol=2)
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
# 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.
RPresence::modAvg(aic.table, param="pA")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.5387 0.0419 0.4562721 0.6190646
RPresence::modAvg(aic.table, param="pB")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.9199 0.0318 0.8313355 0.9639751
RPresence::modAvg(aic.table, param="rA")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.4645 0.0653 0.3414643 0.5920139
RPresence::modAvg(aic.table, param="rBA")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.5483 0.0593 0.4315627 0.6599535
RPresence::modAvg(aic.table, param="rBa")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.5483 0.0593 0.4315627 0.6599535
# 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)
## `geom_smooth()` using formula 'y ~ x'
# 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)
## `geom_smooth()` using formula 'y ~ x'
