# Single Species Single Season Occupancy
# Salamander Example with model averaging with individual model fits
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="CompleteData",
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 0 0 0
## 3 0 1 0 0 0
## 4 1 1 1 1 0
## 5 0 0 1 0 0
## 6 0 0 1 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 0 0 0
## 3 0 1 0 0 0
## 4 1 1 1 1 0
## 5 0 0 1 0 0
## 6 0 0 1 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] 0
# Create the *.pao file
salamander.pao <- RPresence::createPao(input.history,
title='Salamander SSSS')
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 0 0 0
## 3 0 1 0 0 0
## 4 1 1 1 1 0
## 5 0 0 1 0 0
## 6 0 0 1 0 0
## 7 0 0 1 0 0
## 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] 1
##
## $survcov
## SURVEY
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
## 7 1
## 8 1
## 9 1
## 10 1
## 11 1
## 12 1
## 13 1
## 14 1
## 15 1
## 16 1
## 17 1
## 18 1
## 19 1
## 20 1
## 21 1
## 22 1
## 23 1
## 24 1
## 25 1
## 26 1
## 27 1
## 28 1
## 29 1
## 30 1
## 31 1
## 32 1
## 33 1
## 34 1
## 35 1
## 36 1
## 37 1
## 38 1
## 39 1
## 40 2
## 41 2
## 42 2
## 43 2
## 44 2
## 45 2
## 46 2
## 47 2
## 48 2
## 49 2
## 50 2
## 51 2
## 52 2
## 53 2
## 54 2
## 55 2
## 56 2
## 57 2
## 58 2
## 59 2
## 60 2
## 61 2
## 62 2
## 63 2
## 64 2
## 65 2
## 66 2
## 67 2
## 68 2
## 69 2
## 70 2
## 71 2
## 72 2
## 73 2
## 74 2
## 75 2
## 76 2
## 77 2
## 78 2
## 79 3
## 80 3
## 81 3
## 82 3
## 83 3
## 84 3
## 85 3
## 86 3
## 87 3
## 88 3
## 89 3
## 90 3
## 91 3
## 92 3
## 93 3
## 94 3
## 95 3
## 96 3
## 97 3
## 98 3
## 99 3
## 100 3
## 101 3
## 102 3
## 103 3
## 104 3
## 105 3
## 106 3
## 107 3
## 108 3
## 109 3
## 110 3
## 111 3
## 112 3
## 113 3
## 114 3
## 115 3
## 116 3
## 117 3
## 118 4
## 119 4
## 120 4
## 121 4
## 122 4
## 123 4
## 124 4
## 125 4
## 126 4
## 127 4
## 128 4
## 129 4
## 130 4
## 131 4
## 132 4
## 133 4
## 134 4
## 135 4
## 136 4
## 137 4
## 138 4
## 139 4
## 140 4
## 141 4
## 142 4
## 143 4
## 144 4
## 145 4
## 146 4
## 147 4
## 148 4
## 149 4
## 150 4
## 151 4
## 152 4
## 153 4
## 154 4
## 155 4
## 156 4
## 157 5
## 158 5
## 159 5
## 160 5
## 161 5
## 162 5
## 163 5
## 164 5
## 165 5
## 166 5
## 167 5
## 168 5
## 169 5
## 170 5
## 171 5
## 172 5
## 173 5
## 174 5
## 175 5
## 176 5
## 177 5
## 178 5
## 179 5
## 180 5
## 181 5
## 182 5
## 183 5
## 184 5
## 185 5
## 186 5
## 187 5
## 188 5
## 189 5
## 190 5
## 191 5
## 192 5
## 193 5
## 194 5
## 195 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"
# Fit some models.
# Note that formula DO NOT HAVE AN = SIGN
mod.pdot <- RPresence::occMod(model=list(psi~1, p~1),
type="so", data=salamander.pao)
## PRESENCE Version 2.12.18.
summary(mod.pdot)
## Model name=psi()p()
## AIC=165.7586
## -2*log-likelihood=161.7586
## num. par=2
names(mod.pdot)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.pdot$beta$psi
## est se
## A1_psi 0.383108 0.508589
# look at estimated occupancy probability. RPresence gives for EACH site in case it depends on covariates
mod.pdot$real$psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.5946225 0.1225937 0.3512137 0.7989789
## psi_unit2 0.5946225 0.1225937 0.3512137 0.7989789
## psi_unit3 0.5946225 0.1225937 0.3512137 0.7989789
## psi_unit4 0.5946225 0.1225937 0.3512137 0.7989789
## psi_unit5 0.5946225 0.1225937 0.3512137 0.7989789
mod.pdot.psi <-mod.pdot$real$psi[1,] # occupancy probability
mod.pdot.psi
## est se lower_0.95 upper_0.95
## psi_unit1 0.5946225 0.1225937 0.3512137 0.7989789
# look at the estimated probability of detection. It gives an estimate for every site at very visit
head(mod.pdot$real$p)
## est se lower_0.95 upper_0.95
## p1_unit1 0.258729 0.0576996 0.1621604 0.3862912
## p1_unit2 0.258729 0.0576996 0.1621604 0.3862912
## p1_unit3 0.258729 0.0576996 0.1621604 0.3862912
## p1_unit4 0.258729 0.0576996 0.1621604 0.3862912
## p1_unit5 0.258729 0.0576996 0.1621604 0.3862912
## p1_unit6 0.258729 0.0576996 0.1621604 0.3862912
mod.pdot.p <- mod.pdot$real$p[grepl("unit1$", row.names(mod.pdot$real$p)),]
mod.pdot.p
## est se lower_0.95 upper_0.95
## p1_unit1 0.258729 0.0576996 0.1621604 0.3862912
## p2_unit1 0.258729 0.0576996 0.1621604 0.3862912
## p3_unit1 0.258729 0.0576996 0.1621604 0.3862912
## p4_unit1 0.258729 0.0576996 0.1621604 0.3862912
## p5_unit1 0.258729 0.0576996 0.1621604 0.3862912
# Look at the posterior probability of detection
names(mod.pdot$derived)
## [1] "psi_c"
mod.pdot$derived$psi_c
## est se lower_0.95 upper_0.95
## unit1 1.0000000 0.0000000 1.00000000 1.0000
## unit2 1.0000000 0.0000000 1.00000000 1.0000
## unit3 1.0000000 0.0000000 1.00000000 1.0000
## unit4 1.0000000 0.0000000 1.00000000 1.0000
## unit5 1.0000000 0.0000000 1.00000000 1.0000
## unit6 1.0000000 0.0000000 1.00000000 1.0000
## unit7 1.0000000 0.0000000 1.00000000 1.0000
## unit8 1.0000000 0.0000000 1.00000000 1.0000
## unit9 1.0000000 0.0000000 1.00000000 1.0000
## unit10 1.0000000 0.0000000 1.00000000 1.0000
## unit11 1.0000000 0.0000000 1.00000000 1.0000
## unit12 1.0000000 0.0000000 1.00000000 1.0000
## unit13 1.0000000 0.0000000 1.00000000 1.0000
## unit14 1.0000000 0.0000000 1.00000000 1.0000
## unit15 0.2471562 0.1471721 0.06512707 0.6074
## unit16 0.2471562 0.1471721 0.06512707 0.6074
## unit17 0.2471562 0.1471721 0.06512707 0.6074
## unit18 0.2471562 0.1471721 0.06512707 0.6074
## unit19 0.2471562 0.1471721 0.06512707 0.6074
## unit20 0.2471562 0.1471721 0.06512707 0.6074
## unit21 0.2471562 0.1471721 0.06512707 0.6074
## unit22 0.2471562 0.1471721 0.06512707 0.6074
## unit23 0.2471562 0.1471721 0.06512707 0.6074
## unit24 0.2471562 0.1471721 0.06512707 0.6074
## unit25 0.2471562 0.1471721 0.06512707 0.6074
## unit26 0.2471562 0.1471721 0.06512707 0.6074
## unit27 0.2471562 0.1471721 0.06512707 0.6074
## unit28 0.2471562 0.1471721 0.06512707 0.6074
## unit29 0.2471562 0.1471721 0.06512707 0.6074
## unit30 0.2471562 0.1471721 0.06512707 0.6074
## unit31 0.2471562 0.1471721 0.06512707 0.6074
## unit32 0.2471562 0.1471721 0.06512707 0.6074
## unit33 0.2471562 0.1471721 0.06512707 0.6074
## unit34 0.2471562 0.1471721 0.06512707 0.6074
## unit35 0.2471562 0.1471721 0.06512707 0.6074
## unit36 1.0000000 0.0000000 1.00000000 1.0000
## unit37 1.0000000 0.0000000 1.00000000 1.0000
## unit38 1.0000000 0.0000000 1.00000000 1.0000
## unit39 1.0000000 0.0000000 1.00000000 1.0000
# alternatively
RPresence::print_one_site_estimates(mod.pdot, site = 1)
## psi()p()
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.5946225 0.1225937 0.3512137 0.7989789
## p_p1_unit1 0.2587290 0.0576996 0.1621604 0.3862912
# Model where p(t) varies across survey occasions
#
mod.pt <- RPresence::occMod(model=list(psi~1, p~SURVEY),
type="so", data=salamander.pao)
## PRESENCE Version 2.12.18.
summary(mod.pt)
## Model name=psi()p(SURVEY)
## AIC=167.7144
## -2*log-likelihood=155.7144
## num. par=6
mod.pt$real$psi[1,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.5798786 0.117562 0.3489783 0.7804144
mod.pt.p <- mod.pt$real$p[grepl("unit1$", row.names(mod.pt$real$p)),]
mod.pt.p
## est se lower_0.95 upper_0.95
## p1_unit1 0.1768716 0.08449753 0.06445107 0.4012766
## p2_unit1 0.1326537 0.07405373 0.04151959 0.3506443
## p3_unit1 0.3979613 0.11900344 0.19981015 0.6363478
## p4_unit1 0.3537434 0.11369960 0.17116105 0.5919834
## p5_unit1 0.2653074 0.10101771 0.11564619 0.4992989
print_one_site_estimates(mod.pt, site = 1)
## psi()p(SURVEY)
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.5798786 0.11756195 0.34897835 0.7804144
## p_p1_unit1 0.1768716 0.08449753 0.06445107 0.4012766
fitted(mod.pt, param="psi")[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.5798786 0.117562 0.3489783 0.7804144
## psi_unit2 0.5798786 0.117562 0.3489783 0.7804144
## psi_unit3 0.5798786 0.117562 0.3489783 0.7804144
## psi_unit4 0.5798786 0.117562 0.3489783 0.7804144
## psi_unit5 0.5798786 0.117562 0.3489783 0.7804144
# Fit a model with the detection probability equal in first 2 occasions and last 3 occasions.
# We need to define a survey covariate that has two levels
# 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
mod.pcustom <- RPresence::occMod(model=list(psi~1, p~d),
cov.list=list(p.cov=survey.cov),
type="so",
data=salamander.pao)
## PRESENCE Version 2.12.18.
mod.pcustom$real$psi[1,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.5831459 0.1186453 0.3495736 0.7845389
mod.pcustom$real$p[grepl("unit1$", row.names(mod.pt$real$p)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.1538956 0.05838078 0.07023254 0.3045734
## p2_unit1 0.1538956 0.05838078 0.07023254 0.3045734
## p3_unit1 0.3371048 0.07678979 0.20591460 0.4993209
## p4_unit1 0.3371048 0.07678979 0.20591460 0.4993209
## p5_unit1 0.3371048 0.07678979 0.20591460 0.4993209
# Model averaging
models<-list(mod.pdot,mod.pt,mod.pcustom)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll warn.conv warn.VC
## 1 psi()p(d) 0.00 0.760 3 156.82 0 0
## 2 psi()p() 2.94 0.175 2 161.76 0 0
## 3 psi()p(SURVEY) 4.90 0.066 6 155.71 0 0
RPresence::modAvg(results, param="psi")[1,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.5849351 0.1193595 0.3496677 0.786949