# Single Species Single Season Occupancy
# Salamander Example with model averaging with a list of models to fit
library(plyr)
library(readxl)
library(RPresence)
library(ggplot2)
# get the data read in
# Data for detections should be a data frame with rows corresponding to sites
# and columns to visits.
# The usual 1=detected; 0=not detected; NA=not visited is used.
input.data <- readxl::read_excel("../salamander.xls",
sheet="MissingData",
na="-",
col_names=FALSE) # notice no column names in row 1 of data file.
head(input.data)
## # A tibble: 6 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 NA 0 0
## 3 0 NA NA 0 0
## 4 1 NA NA 1 0
## 5 0 NA NA 0 0
## 6 0 0 NA 0 0
# Extract the history records
input.history <- input.data[, 1:5] # the history extracted
head(input.history)
## # A tibble: 6 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 NA 0 0
## 3 0 NA NA 0 0
## 4 1 NA NA 1 0
## 5 0 NA NA 0 0
## 6 0 0 NA 0 0
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 39
ncol(input.history)
## [1] 5
range(input.history, na.rm=TRUE) # check that all values are either 0 or 1
## [1] 0 1
sum(is.na(input.history)) # are there any missing values?
## [1] 14
# Create the survey covariates
# This covariate needs to be "stacked" so that sites1...site39 for survey occastion 1
# are then followed by covariate at survey occastion 2 for sites1...site39, etc
survey.cov <- data.frame(site=rep(1:nrow(input.history), ncol(input.history)),
visit=rep(1:ncol(input.history), each=nrow(input.history)),
d=rep( c("d1","d1","d2","d2","d2"), each=nrow(input.history)))
head(survey.cov)
## site visit d
## 1 1 1 d1
## 2 2 1 d1
## 3 3 1 d1
## 4 4 1 d1
## 5 5 1 d1
## 6 6 1 d1
survey.cov[c(1:4, 37:41),]
## site visit d
## 1 1 1 d1
## 2 2 1 d1
## 3 3 1 d1
## 4 4 1 d1
## 37 37 1 d1
## 38 38 1 d1
## 39 39 1 d1
## 40 1 2 d1
## 41 2 2 d1
# Create the *.pao file
salamander.pao <- RPresence::createPao(input.history,
title='Salamander SSSS',
survcov=survey.cov)
salamander.pao
## $nunits
## [1] 39
##
## $nsurveys
## [1] 5
##
## $nseasons
## [1] 1
##
## $nmethods
## [1] 1
##
## $det.data
## # A tibble: 39 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 NA 0 0
## 3 0 NA NA 0 0
## 4 1 NA NA 1 0
## 5 0 NA NA 0 0
## 6 0 0 NA 0 0
## 7 0 0 1 0 NA
## 8 0 0 1 0 0
## 9 0 0 1 0 0
## 10 1 0 0 0 0
## # ... with 29 more rows
##
## $nunitcov
## [1] 1
##
## $unitcov
## TEMP
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10
## 11 11
## 12 12
## 13 13
## 14 14
## 15 15
## 16 16
## 17 17
## 18 18
## 19 19
## 20 20
## 21 21
## 22 22
## 23 23
## 24 24
## 25 25
## 26 26
## 27 27
## 28 28
## 29 29
## 30 30
## 31 31
## 32 32
## 33 33
## 34 34
## 35 35
## 36 36
## 37 37
## 38 38
## 39 39
##
## $nsurvcov
## [1] 4
##
## $survcov
## site visit d SURVEY
## 1 1 1 d1 1
## 2 2 1 d1 1
## 3 3 1 d1 1
## 4 4 1 d1 1
## 5 5 1 d1 1
## 6 6 1 d1 1
## 7 7 1 d1 1
## 8 8 1 d1 1
## 9 9 1 d1 1
## 10 10 1 d1 1
## 11 11 1 d1 1
## 12 12 1 d1 1
## 13 13 1 d1 1
## 14 14 1 d1 1
## 15 15 1 d1 1
## 16 16 1 d1 1
## 17 17 1 d1 1
## 18 18 1 d1 1
## 19 19 1 d1 1
## 20 20 1 d1 1
## 21 21 1 d1 1
## 22 22 1 d1 1
## 23 23 1 d1 1
## 24 24 1 d1 1
## 25 25 1 d1 1
## 26 26 1 d1 1
## 27 27 1 d1 1
## 28 28 1 d1 1
## 29 29 1 d1 1
## 30 30 1 d1 1
## 31 31 1 d1 1
## 32 32 1 d1 1
## 33 33 1 d1 1
## 34 34 1 d1 1
## 35 35 1 d1 1
## 36 36 1 d1 1
## 37 37 1 d1 1
## 38 38 1 d1 1
## 39 39 1 d1 1
## 40 1 2 d1 2
## 41 2 2 d1 2
## 42 3 2 d1 2
## 43 4 2 d1 2
## 44 5 2 d1 2
## 45 6 2 d1 2
## 46 7 2 d1 2
## 47 8 2 d1 2
## 48 9 2 d1 2
## 49 10 2 d1 2
## 50 11 2 d1 2
## 51 12 2 d1 2
## 52 13 2 d1 2
## 53 14 2 d1 2
## 54 15 2 d1 2
## 55 16 2 d1 2
## 56 17 2 d1 2
## 57 18 2 d1 2
## 58 19 2 d1 2
## 59 20 2 d1 2
## 60 21 2 d1 2
## 61 22 2 d1 2
## 62 23 2 d1 2
## 63 24 2 d1 2
## 64 25 2 d1 2
## 65 26 2 d1 2
## 66 27 2 d1 2
## 67 28 2 d1 2
## 68 29 2 d1 2
## 69 30 2 d1 2
## 70 31 2 d1 2
## 71 32 2 d1 2
## 72 33 2 d1 2
## 73 34 2 d1 2
## 74 35 2 d1 2
## 75 36 2 d1 2
## 76 37 2 d1 2
## 77 38 2 d1 2
## 78 39 2 d1 2
## 79 1 3 d2 3
## 80 2 3 d2 3
## 81 3 3 d2 3
## 82 4 3 d2 3
## 83 5 3 d2 3
## 84 6 3 d2 3
## 85 7 3 d2 3
## 86 8 3 d2 3
## 87 9 3 d2 3
## 88 10 3 d2 3
## 89 11 3 d2 3
## 90 12 3 d2 3
## 91 13 3 d2 3
## 92 14 3 d2 3
## 93 15 3 d2 3
## 94 16 3 d2 3
## 95 17 3 d2 3
## 96 18 3 d2 3
## 97 19 3 d2 3
## 98 20 3 d2 3
## 99 21 3 d2 3
## 100 22 3 d2 3
## 101 23 3 d2 3
## 102 24 3 d2 3
## 103 25 3 d2 3
## 104 26 3 d2 3
## 105 27 3 d2 3
## 106 28 3 d2 3
## 107 29 3 d2 3
## 108 30 3 d2 3
## 109 31 3 d2 3
## 110 32 3 d2 3
## 111 33 3 d2 3
## 112 34 3 d2 3
## 113 35 3 d2 3
## 114 36 3 d2 3
## 115 37 3 d2 3
## 116 38 3 d2 3
## 117 39 3 d2 3
## 118 1 4 d2 4
## 119 2 4 d2 4
## 120 3 4 d2 4
## 121 4 4 d2 4
## 122 5 4 d2 4
## 123 6 4 d2 4
## 124 7 4 d2 4
## 125 8 4 d2 4
## 126 9 4 d2 4
## 127 10 4 d2 4
## 128 11 4 d2 4
## 129 12 4 d2 4
## 130 13 4 d2 4
## 131 14 4 d2 4
## 132 15 4 d2 4
## 133 16 4 d2 4
## 134 17 4 d2 4
## 135 18 4 d2 4
## 136 19 4 d2 4
## 137 20 4 d2 4
## 138 21 4 d2 4
## 139 22 4 d2 4
## 140 23 4 d2 4
## 141 24 4 d2 4
## 142 25 4 d2 4
## 143 26 4 d2 4
## 144 27 4 d2 4
## 145 28 4 d2 4
## 146 29 4 d2 4
## 147 30 4 d2 4
## 148 31 4 d2 4
## 149 32 4 d2 4
## 150 33 4 d2 4
## 151 34 4 d2 4
## 152 35 4 d2 4
## 153 36 4 d2 4
## 154 37 4 d2 4
## 155 38 4 d2 4
## 156 39 4 d2 4
## 157 1 5 d2 5
## 158 2 5 d2 5
## 159 3 5 d2 5
## 160 4 5 d2 5
## 161 5 5 d2 5
## 162 6 5 d2 5
## 163 7 5 d2 5
## 164 8 5 d2 5
## 165 9 5 d2 5
## 166 10 5 d2 5
## 167 11 5 d2 5
## 168 12 5 d2 5
## 169 13 5 d2 5
## 170 14 5 d2 5
## 171 15 5 d2 5
## 172 16 5 d2 5
## 173 17 5 d2 5
## 174 18 5 d2 5
## 175 19 5 d2 5
## 176 20 5 d2 5
## 177 21 5 d2 5
## 178 22 5 d2 5
## 179 23 5 d2 5
## 180 24 5 d2 5
## 181 25 5 d2 5
## 182 26 5 d2 5
## 183 27 5 d2 5
## 184 28 5 d2 5
## 185 29 5 d2 5
## 186 30 5 d2 5
## 187 31 5 d2 5
## 188 32 5 d2 5
## 189 33 5 d2 5
## 190 34 5 d2 5
## 191 35 5 d2 5
## 192 36 5 d2 5
## 193 37 5 d2 5
## 194 38 5 d2 5
## 195 39 5 d2 5
##
## $nsurveyseason
## [1] 5
##
## $title
## [1] "Salamander SSSS"
##
## $unitnames
## [1] "unit1" "unit2" "unit3" "unit4" "unit5" "unit6" "unit7"
## [8] "unit8" "unit9" "unit10" "unit11" "unit12" "unit13" "unit14"
## [15] "unit15" "unit16" "unit17" "unit18" "unit19" "unit20" "unit21"
## [22] "unit22" "unit23" "unit24" "unit25" "unit26" "unit27" "unit28"
## [29] "unit29" "unit30" "unit31" "unit32" "unit33" "unit34" "unit35"
## [36] "unit36" "unit37" "unit38" "unit39"
##
## $surveynames
## [1] "1-1" "1-2" "1-3" "1-4" "1-5"
##
## $paoname
## [1] "pres.pao"
##
## $frq
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1
##
## attr(,"class")
## [1] "pao"
summary(salamander.pao)
## paoname=pres.pao
## title=Salamander SSSS
## Naive occ=0.3846154
## nunits nsurveys nseasons nsurveyseason nmethods
## "39" "5" "1" "5" "1"
## nunitcov nsurvcov
## "1" "4"
## unit covariates : TEMP
## survey covariates: site visit d SURVEY
# 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
~1, ~1
~SURVEY, ~1
~d, ~1")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p psi
## 1 ~1 ~1
## 2 ~SURVEY ~1
## 3 ~d ~1
# 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")
fit
},detect.pao=salamander.pao)
##
##
## ***** Starting ~1 ~1
## PRESENCE Version 2.12.18.
##
##
## ***** Starting ~SURVEY ~1
## PRESENCE Version 2.12.18.
##
##
## ***** Starting ~d ~1
## PRESENCE Version 2.12.18.
# 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" "gof" "warnings" "version"
model.fits[[check.model]]$beta
## $psi
## est se
## A1_psi -0.022647 0.46368
##
## $psi.VC
## [,1]
## [1,] 0.215
##
## $p
## est se
## B1_p1 -0.944472 0.33167
##
## $p.VC
## [,1]
## [1,] 0.110005
##
## $VC
## A1_psi B1_p1
## A1_psi 0.215000 -0.077947
## B1_p1 -0.077947 0.110005
names(model.fits[[check.model]]$real)
## [1] "psi" "p"
model.fits[[check.model]]$real$psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit2 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit3 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit4 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit5 0.4943385 0.1159054 0.2826317 0.7080952
model.fits[[check.model]]$real$p[1:5,]
## est se lower_0.95 upper_0.95
## p1_unit1 0.2799979 0.06686437 0.1687471 0.4269244
## p1_unit2 0.2799979 0.06686437 0.1687471 0.4269244
## p1_unit3 0.2799979 0.06686437 0.1687471 0.4269244
## p1_unit4 0.2799979 0.06686437 0.1687471 0.4269244
## p1_unit5 0.2799979 0.06686437 0.1687471 0.4269244
names(model.fits[[check.model]]$derived)
## [1] "psi_c"
model.fits[[check.model]]$derived$psi_c[1:10,]
## est se lower_0.95 upper_0.95
## unit1 1.0000000 0.0000000 1.00000000 1.0000000
## unit2 1.0000000 0.0000000 1.00000000 1.0000000
## unit3 0.2673420 0.1274765 0.09247886 0.5664633
## unit4 1.0000000 0.0000000 1.00000000 1.0000000
## unit5 0.2673420 0.1274765 0.09247886 0.5664633
## unit6 0.2080612 0.1196824 0.05950888 0.5217297
## unit7 1.0000000 0.0000000 1.00000000 1.0000000
## unit8 1.0000000 0.0000000 1.00000000 1.0000000
## unit9 1.0000000 0.0000000 1.00000000 1.0000000
## unit10 1.0000000 0.0000000 1.00000000 1.0000000
tail(model.fits[[check.model]]$derived$psi_c)
## est se lower_0.95 upper_0.95
## unit34 0.1590715 0.1077521 0.03754541 0.478421
## unit35 0.1590715 0.1077521 0.03754541 0.478421
## unit36 1.0000000 0.0000000 1.00000000 1.000000
## unit37 1.0000000 0.0000000 1.00000000 1.000000
## unit38 1.0000000 0.0000000 1.00000000 1.000000
## unit39 1.0000000 0.0000000 1.00000000 1.000000
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model AIC neg2ll npar warn.conv warn.VC DAIC modlike
## 3 psi()p(d) 136.5309 130.5309 3 0 0 0.0000 1.0000
## 2 psi()p(SURVEY) 140.2942 128.2942 6 0 0 3.7633 0.1523
## 1 psi()p() 141.6129 137.6129 2 0 0 5.0820 0.0788
## wgt
## 3 0.8123
## 2 0.1237
## 1 0.0640
names(aic.table)
## [1] "table" "models" "ess"
RPresence::modAvg(aic.table, param="psi")[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.4789378 0.1100774 0.2791215 0.6857289
## psi_unit2 0.4789378 0.1100774 0.2791215 0.6857289
## psi_unit3 0.4789378 0.1100774 0.2791215 0.6857289
## psi_unit4 0.4789378 0.1100774 0.2791215 0.6857289
## psi_unit5 0.4789378 0.1100774 0.2791215 0.6857289
ma.p <- RPresence::modAvg(aic.table, param="p")
ma.p[grepl("unit1$", rownames(ma.p)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.1578821 0.07938321 0.05497670 0.3766382
## p2_unit1 0.1384016 0.07661587 0.04360894 0.3613859
## p3_unit1 0.3898950 0.10295813 0.21482764 0.5988222
## p4_unit1 0.3934898 0.10214168 0.21899934 0.6001702
## p5_unit1 0.3827237 0.10065927 0.21196846 0.5883367