# 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).
# Try and reverse the two species to see if we get the same answers
# 2020-06-24 CJS real$rho and real$nu now in derived$
# 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 <- 2*PG.data + PJ.data
input.history
## ...1 ...2 ...3 ...4 ...5
## 1 2 0 0 0 0
## 2 0 2 0 0 0
## 3 0 0 0 0 0
## 4 2 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 2 2 2
## 7 0 0 0 0 0
## 8 0 0 0 0 0
## 9 0 0 0 0 0
## 10 2 2 0 0 0
## 11 0 0 0 0 0
## 12 0 0 2 2 0
## 13 0 0 2 0 2
## 14 0 0 2 0 2
## 15 2 2 0 2 0
## 16 2 2 2 0 2
## 17 0 2 0 2 0
## 18 2 2 2 0 2
## 19 0 0 0 0 0
## 20 0 2 0 0 0
## 21 NA NA NA 2 0
## 22 0 0 0 2 2
## 23 0 2 2 2 2
## 24 0 0 0 0 0
## 25 2 2 2 3 3
## 26 3 2 2 2 2
## 27 2 2 2 2 2
## 28 0 0 0 0 0
## 29 0 0 0 0 0
## 30 0 0 0 0 0
## 31 0 2 2 0 2
## 32 1 0 1 0 0
## 33 2 2 0 0 0
## 34 0 2 0 0 0
## 35 2 2 2 2 2
## 36 NA 2 2 2 3
## 37 NA 1 1 1 1
## 38 0 0 0 0 0
## 39 2 2 2 3 2
## 40 2 2 2 2 2
## 41 2 2 2 2 2
## 42 0 2 0 0 2
## 43 2 0 2 0 0
## 44 2 2 0 2 0
## 45 2 2 2 2 2
## 46 2 0 2 0 1
## 47 0 2 0 0 0
## 48 NA 1 1 1 1
## 49 1 1 3 1 0
## 50 2 0 0 0 0
## 51 NA 2 2 2 2
## 52 NA 2 2 0 2
## 53 1 1 1 1 1
## 54 2 3 3 3 3
## 55 0 0 2 2 0
## 56 1 0 1 1 1
## 57 NA 1 0 1 1
## 58 1 1 1 1 1
## 59 1 1 1 1 1
## 60 NA 0 0 1 0
## 61 3 1 1 1 1
## 62 NA 1 0 0 1
## 63 NA 1 1 0 1
## 64 0 0 2 1 3
## 65 0 0 1 3 1
## 66 1 1 1 3 1
## 67 1 1 1 0 1
## 68 1 1 1 1 1
## 69 0 2 2 2 2
## 70 0 0 1 1 3
## 71 1 1 1 1 1
## 72 1 1 1 1 1
## 73 1 1 1 1 1
## 74 1 0 1 1 1
## 75 NA 1 1 1 1
## 76 0 3 3 1 0
## 77 1 1 1 1 1
## 78 NA 1 1 1 1
## 79 NA 1 1 1 1
## 80 1 3 3 1 1
## 81 0 1 1 1 1
## 82 1 1 1 1 1
## 83 1 1 1 1 1
## 84 1 1 1 1 1
## 85 0 2 2 2 0
## 86 1 1 1 1 1
## 87 0 0 1 1 1
## 88 1 1 1 1 1
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.4659091
## 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
# If species 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, ~Std..Elevation+SP+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 ~Std..Elevation+SP+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 ~Std..Elevation+SP+SP:Std..Elevation+INT+INT:Std..Elevation
# Look the output from a specific model
check.model <- 1
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.073855 0.215447
## A2_psiBA -0.427138 0.420281
## A3_psiBa 1.613929 0.515902
##
## $psi.VC
## A1 A2 A3
## A1 0.046417 -0.044590 -0.002408
## A2 -0.044590 0.176637 -0.132734
## A3 -0.002408 -0.132734 0.266154
##
## $p
## est se
## B1_pA[1] 2.227318 0.443821
## B2_pB[1] -2.068689 0.471414
## B3_rA[1] -2.005373 0.550806
## B4_rBa[1] 0.427812 0.502024
## B5_rBA[1] 1.608733 0.597632
##
## $p.VC
## B1 B2 B3 B4 B5
## B1 0.196977 -0.195557 -0.217453 -0.053310 0.174595
## B2 -0.195557 0.222231 0.215919 0.050688 -0.201737
## B3 -0.217453 0.215919 0.303387 0.066623 -0.256515
## B4 -0.053310 0.050688 0.066623 0.252028 -0.143552
## B5 0.174595 -0.201737 -0.256515 -0.143552 0.357164
##
## $VC
## A1 A2 A3 B1 B2 B3 B4
## A1 0.046417 -0.044590 -0.002408 0.001284 -0.001318 -0.004024 -0.000575
## A2 -0.044590 0.176637 -0.132734 0.054204 -0.053485 -0.063284 -0.023709
## A3 -0.002408 -0.132734 0.266154 -0.056508 0.049375 0.070106 0.024960
## B1 0.001284 0.054204 -0.056508 0.196977 -0.195557 -0.217453 -0.053310
## B2 -0.001318 -0.053485 0.049375 -0.195557 0.222231 0.215919 0.050688
## B3 -0.004024 -0.063284 0.070106 -0.217453 0.215919 0.303387 0.066623
## B4 -0.000575 -0.023709 0.024960 -0.053310 0.050688 0.066623 0.252028
## B5 0.003674 0.036766 -0.036471 0.174595 -0.201737 -0.256515 -0.143552
## B5
## A1 0.003674
## A2 0.036766
## A3 -0.036471
## B1 0.174595
## B2 -0.201737
## B3 -0.256515
## B4 -0.143552
## B5 0.357164
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.4815 0.0538 0.3785 0.5862
## unit2 0.4815 0.0538 0.3785 0.5862
## unit3 0.4815 0.0538 0.3785 0.5862
## unit4 0.4815 0.0538 0.3785 0.5862
## unit5 0.4815 0.0538 0.3785 0.5862
model.fits[[check.model]]$real$psiBA[1:5,]
## est se lower upper
## unit1 0.3773 0.086 0.2283 0.5538
## unit2 0.3773 0.086 0.2283 0.5538
## unit3 0.3773 0.086 0.2283 0.5538
## unit4 0.3773 0.086 0.2283 0.5538
## unit5 0.3773 0.086 0.2283 0.5538
model.fits[[check.model]]$real$psiBa[1:5,]
## est se lower upper
## unit1 0.7527 0.0671 0.6004 0.8604
## unit2 0.7527 0.0671 0.6004 0.8604
## unit3 0.7527 0.0671 0.6004 0.8604
## unit4 0.7527 0.0671 0.6004 0.8604
## unit5 0.7527 0.0671 0.6004 0.8604
model.fits[[check.model]]$derived$nu[1:5,]
## est se lower upper
## unit1 0.1991038 0.1027179 0.07243286 0.5472975
## unit2 0.1991038 0.1027179 0.07243286 0.5472975
## unit3 0.1991038 0.1027179 0.07243286 0.5472975
## unit4 0.1991038 0.1027179 0.07243286 0.5472975
## unit5 0.1991038 0.1027179 0.07243286 0.5472975
model.fits[[check.model]]$real$pA[1:5,]
## est se lower upper
## unit1_1-1 0.9027 0.039 0.7953 0.9568
## unit2_1-1 0.9027 0.039 0.7953 0.9568
## unit3_1-1 0.9027 0.039 0.7953 0.9568
## unit4_1-1 0.9027 0.039 0.7953 0.9568
## unit5_1-1 0.9027 0.039 0.7953 0.9568
model.fits[[check.model]]$real$pB[1:5,]
## est se lower upper
## unit1_1-1 0.5396 0.0416 0.4576 0.6194
## unit2_1-1 0.5396 0.0416 0.4576 0.6194
## unit3_1-1 0.5396 0.0416 0.4576 0.6194
## unit4_1-1 0.5396 0.0416 0.4576 0.6194
## unit5_1-1 0.5396 0.0416 0.4576 0.6194
model.fits[[check.model]]$real$rA[1:5,]
## est se lower upper
## unit1_1-1 0.5553 0.0632 0.4306 0.6734
## unit2_1-1 0.5553 0.0632 0.4306 0.6734
## unit3_1-1 0.5553 0.0632 0.4306 0.6734
## unit4_1-1 0.5553 0.0632 0.4306 0.6734
## unit5_1-1 0.5553 0.0632 0.4306 0.6734
model.fits[[check.model]]$real$rBA[1:5,]
## est se lower upper
## unit1_1-1 0.4408 0.0848 0.2866 0.6073
## unit2_1-1 0.4408 0.0848 0.2866 0.6073
## unit3_1-1 0.4408 0.0848 0.2866 0.6073
## unit4_1-1 0.4408 0.0848 0.2866 0.6073
## unit5_1-1 0.4408 0.0848 0.2866 0.6073
model.fits[[check.model]]$real$rBa[1:5,]
## est se lower upper
## unit1_1-1 0.5473 0.1139 0.3294 0.7485
## unit2_1-1 0.5473 0.1139 0.3294 0.7485
## unit3_1-1 0.5473 0.1139 0.3294 0.7485
## unit4_1-1 0.5473 0.1139 0.3294 0.7485
## unit5_1-1 0.5473 0.1139 0.3294 0.7485
model.fits[[check.model]]$derived$rho[1:5,]
## est se lower upper
## unit1_1-1 0.651934 0.3272864 0.2437094 1.743954
## unit2_1-1 0.651934 0.3272864 0.2437094 1.743954
## unit3_1-1 0.651934 0.3272864 0.2437094 1.743954
## unit4_1-1 0.651934 0.3272864 0.2437094 1.743954
## unit5_1-1 0.651934 0.3272864 0.2437094 1.743954
names(model.fits[[check.model]]$derived)
## [1] "nu" "rho"
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model
## 4 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)
## 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 679.3718 659.3718 10 7 0 0.0000 1 1
## 3 750.6240 736.6240 7 7 0 71.2522 0 0
## 1 751.8665 735.8665 8 7 0 72.4947 0 0
## 2 760.3286 746.3286 7 7 0 80.9568 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.0005 0.0008 2.170571e-05 0.01139758 psiA 502.92 -2.148793
## 2 0.0007 0.0011 3.212542e-05 0.01504383 psiA 518.16 -2.052179
## 3 0.0007 0.0011 3.212542e-05 0.01504383 psiA 518.16 -2.052179
## 4 0.0007 0.0011 3.212542e-05 0.01504383 psiA 518.16 -2.052179
## 5 0.0009 0.0014 4.259229e-05 0.01869484 psiA 533.40 -1.955564
## 6 0.0013 0.0019 7.392529e-05 0.02240528 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
