# Blue Ridge Salamander - with missing data
# Notice how the data are read in with the missing values.
# Single Species Single Season Occupancy
# Fitting a single model
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 *.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 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] 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 a model
# Note that formula DO NOT HAVE AN = SIGN
mod.fit <- RPresence::occMod(model=list(psi~1, p~1),
type="so",
data=salamander.pao)
## PRESENCE Version 2.12.18.
summary(mod.fit)
## Model name=psi()p()
## AIC=141.6129
## -2*log-likelihood=137.6129
## num. par=2
# Look the objects returned in more details
names(mod.fit)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
# look at estimated occupancy probability. RPresence gives for EACH site in case it depends on covariates
mod.fit$beta$psi # on the logit scale
## est se
## A1_psi -0.022647 0.46368
mod.fit$real$psi # on the regular 0-1 scale for each site
## 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
## psi_unit6 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit7 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit8 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit9 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit10 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit11 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit12 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit13 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit14 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit15 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit16 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit17 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit18 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit19 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit20 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit21 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit22 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit23 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit24 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit25 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit26 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit27 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit28 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit29 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit30 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit31 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit32 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit33 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit34 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit35 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit36 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit37 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit38 0.4943385 0.1159054 0.2826317 0.7080952
## psi_unit39 0.4943385 0.1159054 0.2826317 0.7080952
mod.fit$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
# look at the estimated probability of detection. It gives an estimate for every site at very visit
mod.fit$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
# extract the detection probabilities for each visit. The row names have the unit names
mod.fit.p <- mod.fit$real$p[grepl("unit1$", row.names(mod.fit$real$p)),]
mod.fit.p
## est se lower_0.95 upper_0.95
## p1_unit1 0.2799979 0.06686437 0.1687471 0.4269244
## p2_unit1 0.2799979 0.06686437 0.1687471 0.4269244
## p3_unit1 0.2799979 0.06686437 0.1687471 0.4269244
## p4_unit1 0.2799979 0.06686437 0.1687471 0.4269244
## p5_unit1 0.2799979 0.06686437 0.1687471 0.4269244
# Look at the posterior probability of detection
names(mod.fit$derived)
## [1] "psi_c"
mod.fit$derived$psi_c
## 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
## unit11 1.0000000 0.0000000 1.00000000 1.0000000
## unit12 1.0000000 0.0000000 1.00000000 1.0000000
## unit13 1.0000000 0.0000000 1.00000000 1.0000000
## unit14 1.0000000 0.0000000 1.00000000 1.0000000
## unit15 0.2080612 0.1196824 0.05950888 0.5217297
## unit16 0.1590715 0.1077521 0.03754541 0.4784210
## unit17 0.1590715 0.1077521 0.03754541 0.4784210
## unit18 0.1590715 0.1077521 0.03754541 0.4784210
## unit19 0.1590715 0.1077521 0.03754541 0.4784210
## unit20 0.1590715 0.1077521 0.03754541 0.4784210
## unit21 0.2080612 0.1196824 0.05950888 0.5217297
## unit22 0.1590715 0.1077521 0.03754541 0.4784210
## unit23 0.1590715 0.1077521 0.03754541 0.4784210
## unit24 0.1590715 0.1077521 0.03754541 0.4784210
## unit25 0.1590715 0.1077521 0.03754541 0.4784210
## unit26 0.2080612 0.1196824 0.05950888 0.5217297
## unit27 0.1590715 0.1077521 0.03754541 0.4784210
## unit28 0.1590715 0.1077521 0.03754541 0.4784210
## unit29 0.1590715 0.1077521 0.03754541 0.4784210
## unit30 0.1590715 0.1077521 0.03754541 0.4784210
## unit31 0.1590715 0.1077521 0.03754541 0.4784210
## unit32 0.2080612 0.1196824 0.05950888 0.5217297
## unit33 0.1590715 0.1077521 0.03754541 0.4784210
## unit34 0.1590715 0.1077521 0.03754541 0.4784210
## unit35 0.1590715 0.1077521 0.03754541 0.4784210
## unit36 1.0000000 0.0000000 1.00000000 1.0000000
## unit37 1.0000000 0.0000000 1.00000000 1.0000000
## unit38 1.0000000 0.0000000 1.00000000 1.0000000
## unit39 1.0000000 0.0000000 1.00000000 1.0000000
# alternatively
RPresence::print_one_site_estimates(mod.fit, site = 1)
## psi()p()
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.4943385 0.11590537 0.2826317 0.7080952
## p_p1_unit1 0.2799979 0.06686437 0.1687471 0.4269244
fitted(mod.fit, param="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
fitted(mod.fit, param="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