# Multi species single season with a list of models
# Using RPresence
# Co-occurence of spotted (SO) and barred owls (BO)
# 2020-06-24 CJS real$rho and real$nu now in derived$
# 2018-12-13 code contributed by Neil Faught
library(car)
## Loading required package: carData
library(ggplot2)
library(readxl)
library(reshape2)
library(RPresence)
# get the data
SO.data <- readxl::read_excel(file.path("..","SpottedAndBarredOwls.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "B11:K161")
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * ...
BO.data <- readxl::read_excel(file.path("..","SpottedAndBarredOwls.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "M11:V161")
## 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.data <- 2*BO.data + SO.data
input.data
## ...1 ...2 ...3 ...4 ...5 ...6 ...7 ...8 ...9 ...10
## 1 0 1 1 1 NA NA NA NA NA NA
## 2 1 1 1 1 NA NA NA NA NA NA
## 3 0 3 2 0 NA NA NA NA NA NA
## 4 1 1 1 1 NA NA NA NA NA NA
## 5 0 0 0 1 NA NA NA NA NA NA
## 6 0 1 1 NA NA NA NA NA NA NA
## 7 1 1 1 NA NA NA NA NA NA NA
## 8 0 0 0 0 NA NA NA NA NA NA
## 9 0 2 2 2 NA NA NA NA NA NA
## 10 1 1 1 NA NA NA NA NA NA NA
## 11 1 1 1 NA NA NA NA NA NA NA
## 12 1 1 1 1 1 0 1 NA NA NA
## 13 0 0 0 0 NA NA NA NA NA NA
## 14 1 1 1 1 1 NA NA NA NA NA
## 15 1 0 0 0 0 0 0 1 1 NA
## 16 1 1 1 0 1 1 NA NA NA NA
## 17 1 0 0 0 NA NA NA NA NA NA
## 18 0 0 0 0 NA NA NA NA NA NA
## 19 1 1 1 1 1 1 1 1 1 1
## 20 0 0 0 0 NA NA NA NA NA NA
## 21 1 1 1 NA NA NA NA NA NA NA
## 22 1 1 NA NA NA NA NA NA NA NA
## 23 1 0 0 1 0 0 1 0 1 NA
## 24 0 0 0 0 NA NA NA NA NA NA
## 25 0 0 0 0 NA NA NA NA NA NA
## 26 0 1 1 1 NA NA NA NA NA NA
## 27 0 0 0 0 NA NA NA NA NA NA
## 28 3 1 0 0 0 0 NA NA NA NA
## 29 0 0 0 0 NA NA NA NA NA NA
## 30 0 1 0 1 1 NA NA NA NA NA
## 31 0 0 0 0 NA NA NA NA NA NA
## 32 0 0 1 1 1 NA NA NA NA NA
## 33 0 0 1 1 0 1 1 0 0 NA
## 34 0 0 0 0 NA NA NA NA NA NA
## 35 1 0 0 0 1 1 1 NA NA NA
## 36 1 1 0 1 1 NA NA NA NA NA
## 37 0 1 0 0 1 NA NA NA NA NA
## 38 1 1 1 NA NA NA NA NA NA NA
## 39 1 1 1 0 1 1 NA NA NA NA
## 40 0 1 0 1 1 1 NA NA NA NA
## 41 0 0 0 0 NA NA NA NA NA NA
## 42 0 0 1 0 0 NA NA NA NA NA
## 43 0 1 1 1 1 1 1 1 1 1
## 44 1 1 0 1 1 NA NA NA NA NA
## 45 2 0 0 1 0 NA NA NA NA NA
## 46 0 0 2 0 2 NA NA NA NA NA
## 47 1 1 1 1 NA NA NA NA NA NA
## 48 0 1 1 1 NA NA NA NA NA NA
## 49 1 1 3 NA NA NA NA NA NA NA
## 50 1 0 0 0 2 0 0 0 NA NA
## 51 0 0 3 NA NA NA NA NA NA NA
## 52 1 1 0 2 1 NA NA NA NA NA
## 53 1 1 1 1 1 1 NA NA NA NA
## 54 0 0 3 2 0 2 1 0 1 0
## 55 0 1 1 1 1 1 NA NA NA NA
## 56 1 0 2 0 NA NA NA NA NA NA
## 57 0 0 0 0 NA NA NA NA NA NA
## 58 1 1 1 0 0 1 1 NA NA NA
## 59 0 0 0 0 NA NA NA NA NA NA
## 60 1 0 1 1 0 0 NA NA NA NA
## 61 0 1 0 2 NA NA NA NA NA NA
## 62 0 1 1 1 NA NA NA NA NA NA
## 63 0 1 1 0 0 NA NA NA NA NA
## 64 2 0 0 0 NA NA NA NA NA NA
## 65 2 1 1 2 2 NA NA NA NA NA
## 66 0 1 1 NA NA NA NA NA NA NA
## 67 2 2 2 2 0 NA NA NA NA NA
## 68 2 0 0 2 NA NA NA NA NA NA
## 69 0 0 0 0 NA NA NA NA NA NA
## 70 1 1 1 1 1 1 1 1 NA NA
## 71 0 0 0 0 NA NA NA NA NA NA
## 72 1 0 1 1 1 1 0 1 0 1
## 73 0 2 2 2 0 NA NA NA NA NA
## 74 1 1 0 0 2 1 0 0 3 2
## 75 0 1 0 0 1 NA NA NA NA NA
## 76 2 2 2 2 NA NA NA NA NA NA
## 77 0 0 0 0 NA NA NA NA NA NA
## 78 1 0 0 0 0 0 NA NA NA NA
## 79 0 0 0 0 NA NA NA NA NA NA
## 80 0 0 0 0 NA NA NA NA NA NA
## 81 2 0 0 0 NA NA NA NA NA NA
## 82 0 0 0 0 NA NA NA NA NA NA
## 83 0 2 0 0 NA NA NA NA NA NA
## 84 0 0 2 2 NA NA NA NA NA NA
## 85 1 1 0 1 1 1 1 NA NA NA
## 86 1 0 0 0 NA NA NA NA NA NA
## 87 0 0 0 1 NA NA NA NA NA NA
## 88 1 1 1 1 1 NA NA NA NA NA
## 89 0 2 NA NA NA NA NA NA NA NA
## 90 0 0 1 0 NA NA NA NA NA NA
## 91 0 2 0 0 NA NA NA NA NA NA
## 92 1 0 1 1 1 1 NA NA NA NA
## 93 1 1 1 1 NA NA NA NA NA NA
## 94 0 0 0 0 1 1 1 NA NA NA
## 95 1 0 0 0 1 1 1 1 NA NA
## 96 0 1 0 0 NA NA NA NA NA NA
## 97 0 0 0 1 1 NA NA NA NA NA
## 98 1 1 1 NA NA NA NA NA NA NA
## 99 2 0 0 0 NA NA NA NA NA NA
## 100 2 1 1 1 1 NA NA NA NA NA
## 101 0 3 0 0 NA NA NA NA NA NA
## 102 0 0 1 1 0 NA NA NA NA NA
## 103 0 1 1 NA NA NA NA NA NA NA
## 104 1 1 1 1 NA NA NA NA NA NA
## 105 1 0 1 2 0 0 NA NA NA NA
## 106 1 0 0 1 0 1 0 0 1 1
## 107 1 1 1 NA NA NA NA NA NA NA
## 108 1 0 0 0 1 0 1 NA NA NA
## 109 0 2 0 0 NA NA NA NA NA NA
## 110 0 1 0 1 1 NA NA NA NA NA
## 111 0 0 0 0 NA NA NA NA NA NA
## 112 0 0 1 1 1 1 1 1 NA NA
## 113 2 0 0 0 NA NA NA NA NA NA
## 114 0 0 1 1 0 0 0 NA NA NA
## 115 0 2 3 2 NA NA NA NA NA NA
## 116 0 0 0 2 NA NA NA NA NA NA
## 117 0 1 1 1 NA NA NA NA NA NA
## 118 0 1 1 3 1 0 0 NA NA NA
## 119 3 1 1 1 1 1 NA NA NA NA
## 120 1 1 0 0 1 NA NA NA NA NA
## 121 0 1 1 0 0 NA NA NA NA NA
## 122 0 0 0 0 NA NA NA NA NA NA
## 123 0 1 0 0 NA NA NA NA NA NA
## 124 0 0 0 0 NA NA NA NA NA NA
## 125 0 0 0 0 NA NA NA NA NA NA
## 126 1 1 1 1 0 0 1 NA NA NA
## 127 1 0 1 1 1 NA NA NA NA NA
## 128 0 0 0 1 1 NA NA NA NA NA
## 129 1 1 1 NA NA NA NA NA NA NA
## 130 0 0 0 2 NA NA NA NA NA NA
## 131 1 0 2 3 1 NA NA NA NA NA
## 132 0 0 0 0 NA NA NA NA NA NA
## 133 1 0 0 0 0 0 1 1 0 NA
## 134 1 1 1 1 1 1 NA NA NA NA
## 135 0 2 0 0 NA NA NA NA NA NA
## 136 0 0 0 0 0 NA NA NA NA NA
## 137 0 2 2 0 NA NA NA NA NA NA
## 138 0 0 0 0 NA NA NA NA NA NA
## 139 1 0 1 1 NA NA NA NA NA NA
## 140 0 0 0 1 NA NA NA NA NA NA
## 141 0 0 0 0 NA NA NA NA NA NA
## 142 0 2 0 0 NA NA NA NA NA NA
## 143 1 1 NA NA NA NA NA NA NA NA
## 144 1 0 1 1 1 0 3 NA NA NA
## 145 0 1 0 1 0 NA NA NA NA NA
## 146 0 0 0 0 0 NA NA NA NA NA
## 147 0 0 0 0 NA NA NA NA NA NA
## 148 0 0 0 0 NA NA NA NA NA NA
## 149 0 1 1 1 1 1 1 NA NA NA
## 150 0 NA NA NA NA NA NA NA NA NA
## 151 0 1 1 0 1 1 NA NA NA NA
input.history = input.data
# Read in survey level covariates
night = readxl::read_excel(file.path("..","SpottedAndBarredOwls.xls"),
sheet="Nite Covariate", na='-',
col_names=FALSE,
range = "B11:K161")
## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * ...
head(night)
## # A tibble: 6 x 10
## ...1 ...2 ...3 ...4 ...5 ...6 ...7 ...8 ...9 ...10
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 NA NA NA NA NA NA
## 2 0 0 1 0 NA NA NA NA NA NA
## 3 0 0 0 0 NA NA NA NA NA NA
## 4 0 0 1 1 NA NA NA NA NA NA
## 5 0 0 0 0 NA NA NA NA NA NA
## 6 0 0 0 NA NA NA NA NA NA NA
# Convert to a survey covariate. You need to stack the data by columns
survey.cov <- data.frame(Site=1:nrow(night),
visit=as.factor(rep(1:ncol(night), each=nrow(night))),
night =unlist(night), stringsAsFactors=FALSE)
head(survey.cov)
## Site visit night
## ...11 1 1 0
## ...12 2 1 0
## ...13 3 1 0
## ...14 4 1 0
## ...15 5 1 0
## ...16 6 1 0
# create the pao file
owl.pao <- createPao(input.history,
survcov = survey.cov,
title="Owl multi species - co-occurance")
summary(owl.pao)
## paoname=pres.pao
## title=Owl multi species - co-occurance
## Naive occ=0.794702
## naiveR =0.218543
## nunits nsurveys nseasons nsurveyseason nmethods
## "151" "10" "1" "10" "1"
## nunitcov nsurvcov
## "1" "4"
## unit covariates : TEMP
## survey covariates: Site visit night 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, ~SP+INT
~SP+INT_o+SP:INT_o+ INT_d, ~SP
~SP+INT_o+SP:INT_o, ~SP
~SP+INT_o + SP:INT_o + SP:INT_o + night, ~SP+INT")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p psi
## 1 ~SP+INT_o+SP:INT_o + INT_d ~SP+INT
## 2 ~SP+INT_o+SP:INT_o ~SP+INT
## 3 ~SP+INT_o+SP:INT_o+ INT_d ~SP
## 4 ~SP+INT_o+SP:INT_o ~SP
## 5 ~SP+INT_o + SP:INT_o + SP:INT_o + night ~SP+INT
# 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=owl.pao)
##
##
## ***** Starting ~SP+INT_o+SP:INT_o + INT_d ~SP+INT
##
##
## ***** Starting ~SP+INT_o+SP:INT_o ~SP+INT
##
##
## ***** Starting ~SP+INT_o+SP:INT_o+ INT_d ~SP
##
##
## ***** Starting ~SP+INT_o+SP:INT_o ~SP
##
##
## ***** Starting ~SP+INT_o + SP:INT_o + SP:INT_o + night ~SP+INT
# Look the output from a specific model
check.model <- 5
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.948291 0.227282
## A2_psiBA -1.693274 0.349993
## A3_psiBa 0.590229 0.603968
##
## $psi.VC
## A1 A2 A3
## A1 0.051657 -0.031756 -0.033266
## A2 -0.031756 0.122495 -0.099582
## A3 -0.033266 -0.099582 0.364778
##
## $p
## est se
## B1_pA[1] 0.399719 0.156349
## B2_pB[1] -0.932362 0.368743
## B3_rA[1] -1.030705 0.301301
## B4_pA[1].night_pA 2.502748 0.447166
## B5_rBA[1] 0.013165 0.461707
##
## $p.VC
## B1 B2 B3 B4 B5
## B1 0.024445 -0.019773 -0.031302 -0.018429 0.010380
## B2 -0.019773 0.135971 0.008861 0.013772 -0.115840
## B3 -0.031302 0.008861 0.090782 0.017022 -0.043933
## B4 -0.018429 0.013772 0.017022 0.199957 -0.002903
## B5 0.010380 -0.115840 -0.043933 -0.002903 0.213173
##
## $VC
## A1 A2 A3 B1 B2 B3 B4
## A1 0.051657 -0.031756 -0.033266 0.001515 0.013239 -0.020699 0.000283
## A2 -0.031756 0.122495 -0.099582 0.020882 -0.014485 -0.026489 -0.019685
## A3 -0.033266 -0.099582 0.364778 -0.029796 -0.070343 0.074887 0.025999
## B1 0.001515 0.020882 -0.029796 0.024445 -0.019773 -0.031302 -0.018429
## B2 0.013239 -0.014485 -0.070343 -0.019773 0.135971 0.008861 0.013772
## B3 -0.020699 -0.026489 0.074887 -0.031302 0.008861 0.090782 0.017022
## B4 0.000283 -0.019685 0.025999 -0.018429 0.013772 0.017022 0.199957
## B5 -0.001005 -0.018514 0.066979 0.010380 -0.115840 -0.043933 -0.002903
## B5
## A1 -0.001005
## A2 -0.018514
## A3 0.066979
## B1 0.010380
## B2 -0.115840
## B3 -0.043933
## B4 -0.002903
## B5 0.213173
model.fits[[check.model]]$dmat
## $psi
## a1 a2 a3
## psiA "1" "0" "0"
## psiBA "1" "1" "0"
## psiBa "1" "1" "1"
##
## $p
## b1 b2 b3 b4 b5
## pA[1] "1" "0" "0" "night_pA" "0"
## pA[2] "1" "0" "0" "night_pA" "0"
## pA[3] "1" "0" "0" "night_pA" "0"
## pA[4] "1" "0" "0" "night_pA" "0"
## pA[5] "1" "0" "0" "night_pA" "0"
## pA[6] "1" "0" "0" "night_pA" "0"
## pA[7] "1" "0" "0" "night_pA" "0"
## pA[8] "1" "0" "0" "night_pA" "0"
## pA[9] "1" "0" "0" "night_pA" "0"
## pA[10] "1" "0" "0" "night_pA" "0"
## pB[1] "1" "1" "0" "night_pB" "0"
## pB[2] "1" "1" "0" "night_pB" "0"
## pB[3] "1" "1" "0" "night_pB" "0"
## pB[4] "1" "1" "0" "night_pB" "0"
## pB[5] "1" "1" "0" "night_pB" "0"
## pB[6] "1" "1" "0" "night_pB" "0"
## pB[7] "1" "1" "0" "night_pB" "0"
## pB[8] "1" "1" "0" "night_pB" "0"
## pB[9] "1" "1" "0" "night_pB" "0"
## pB[10] "1" "1" "0" "night_pB" "0"
## rA[1] "1" "0" "1" "night_rA" "0"
## rA[2] "1" "0" "1" "night_rA" "0"
## rA[3] "1" "0" "1" "night_rA" "0"
## rA[4] "1" "0" "1" "night_rA" "0"
## rA[5] "1" "0" "1" "night_rA" "0"
## rA[6] "1" "0" "1" "night_rA" "0"
## rA[7] "1" "0" "1" "night_rA" "0"
## rA[8] "1" "0" "1" "night_rA" "0"
## rA[9] "1" "0" "1" "night_rA" "0"
## rA[10] "1" "0" "1" "night_rA" "0"
## rBA[1] "1" "1" "1" "night_rBA" "1"
## rBA[2] "1" "1" "1" "night_rBA" "1"
## rBA[3] "1" "1" "1" "night_rBA" "1"
## rBA[4] "1" "1" "1" "night_rBA" "1"
## rBA[5] "1" "1" "1" "night_rBA" "1"
## rBA[6] "1" "1" "1" "night_rBA" "1"
## rBA[7] "1" "1" "1" "night_rBA" "1"
## rBA[8] "1" "1" "1" "night_rBA" "1"
## rBA[9] "1" "1" "1" "night_rBA" "1"
## rBA[10] "1" "1" "1" "night_rBA" "1"
## rBa[1] "1" "1" "1" "night_rBa" "1"
## rBa[2] "1" "1" "1" "night_rBa" "1"
## rBa[3] "1" "1" "1" "night_rBa" "1"
## rBa[4] "1" "1" "1" "night_rBa" "1"
## rBa[5] "1" "1" "1" "night_rBa" "1"
## rBa[6] "1" "1" "1" "night_rBa" "1"
## rBa[7] "1" "1" "1" "night_rBa" "1"
## rBa[8] "1" "1" "1" "night_rBa" "1"
## rBa[9] "1" "1" "1" "night_rBa" "1"
## rBa[10] "1" "1" "1" "night_rBa" "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.7208 0.0457 0.6231 0.8012
## unit2 0.7208 0.0457 0.6231 0.8012
## unit3 0.7208 0.0457 0.6231 0.8012
## unit4 0.7208 0.0457 0.6231 0.8012
## unit5 0.7208 0.0457 0.6231 0.8012
model.fits[[check.model]]$real$psiBA[1:5,]
## est se lower upper
## unit1 0.3219 0.0726 0.1983 0.4768
## unit2 0.3219 0.0726 0.1983 0.4768
## unit3 0.3219 0.0726 0.1983 0.4768
## unit4 0.3219 0.0726 0.1983 0.4768
## unit5 0.3219 0.0726 0.1983 0.4768
model.fits[[check.model]]$real$psiBa[1:5,]
## est se lower upper
## unit1 0.4614 0.1138 0.2588 0.6776
## unit2 0.4614 0.1138 0.2588 0.6776
## unit3 0.4614 0.1138 0.2588 0.6776
## unit4 0.4614 0.1138 0.2588 0.6776
## unit5 0.4614 0.1138 0.2588 0.6776
model.fits[[check.model]]$derived$nu[1:5,]
## est se lower upper
## unit1 0.5542004 0.3347196 0.1696518 1.810403
## unit2 0.5542004 0.3347196 0.1696518 1.810403
## unit3 0.5542004 0.3347196 0.1696518 1.810403
## unit4 0.5542004 0.3347196 0.1696518 1.810403
## unit5 0.5542004 0.3347196 0.1696518 1.810403
model.fits[[check.model]]$real$pA[1:5,]
## est se lower upper
## unit1_1-1 0.5986 0.0376 0.5233 0.6696
## unit2_1-1 0.5986 0.0376 0.5233 0.6696
## unit3_1-1 0.5986 0.0376 0.5233 0.6696
## unit4_1-1 0.5986 0.0376 0.5233 0.6696
## unit5_1-1 0.5986 0.0376 0.5233 0.6696
model.fits[[check.model]]$real$pB[1:5,]
## est se lower upper
## unit1_1-1 0.3699 0.081 0.229 0.5371
## unit2_1-1 0.3699 0.081 0.229 0.5371
## unit3_1-1 0.3699 0.081 0.229 0.5371
## unit4_1-1 0.3699 0.081 0.229 0.5371
## unit5_1-1 0.3699 0.081 0.229 0.5371
model.fits[[check.model]]$real$rA[1:5,]
## est se lower upper
## unit1_1-1 0.3473 0.052 0.2534 0.4548
## unit2_1-1 0.3473 0.052 0.2534 0.4548
## unit3_1-1 0.3473 0.052 0.2534 0.4548
## unit4_1-1 0.3473 0.052 0.2534 0.4548
## unit5_1-1 0.3473 0.052 0.2534 0.4548
model.fits[[check.model]]$real$rBA[1:5,]
## est se lower upper
## unit1_1-1 0.1751 0.0411 0.1083 0.2706
## unit2_1-1 0.1751 0.0411 0.1083 0.2706
## unit3_1-1 0.1751 0.0411 0.1083 0.2706
## unit4_1-1 0.1751 0.0411 0.1083 0.2706
## unit5_1-1 0.1751 0.0411 0.1083 0.2706
model.fits[[check.model]]$real$rBa[1:5,]
## est se lower upper
## unit1_1-1 0.1751 0.0411 0.1083 0.2706
## unit2_1-1 0.1751 0.0411 0.1083 0.2706
## unit3_1-1 0.1751 0.0411 0.1083 0.2706
## unit4_1-1 0.1751 0.0411 0.1083 0.2706
## unit5_1-1 0.1751 0.0411 0.1083 0.2706
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 AIC
## 5 psi(SP P INT)p(SP P INT_o P SP T INT_o P SP T INT_o P night) 1228.112
## 4 psi(SP)p(SP P INT_o P SP T INT_o) 1270.952
## 2 psi(SP P INT)p(SP P INT_o P SP T INT_o) 1272.201
## 3 psi(SP)p(SP P INT_o P SP T INT_o P INT_d) 1272.847
## 1 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d) 1273.980
## neg2ll npar warn.conv warn.VC DAIC modlike wgt
## 5 1212.112 8 7 0 0.0000 1 1
## 4 1258.952 6 7 0 42.8395 0 0
## 2 1258.201 7 7 0 44.0892 0 0
## 3 1258.847 7 7 0 44.7345 0 0
## 1 1257.980 8 7 0 45.8681 0 0