# Brook Trout
# This is from MARK demo files on occupancy.
# Collected via electrofishing three 50 m sections of streams at 77 sites
# in the Upper Chattachochee River basin.
# 77 streams 3 occasions, 4 covariates: elevation, cross sectional area each occasion.
# Single Species Single Season Occupancy
# Fitting a single model using RPresence
# 2018-12-02 Code contributed by Neil Faught
library(car)
library(readxl)
library(RPresence)
## Warning: package 'RPresence' was built under R version 3.5.0
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("../BrookTrout.xls",
sheet="DetectionHistory",
na="-",
col_names=FALSE) # notice no column names in row 1 of data file.
head(input.data)
## # A tibble: 6 x 3
## X__1 X__2 X__3
## <dbl> <dbl> <dbl>
## 1 0 0 0
## 2 1.00 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 1.00 0
# Extract the history records
input.history <- input.data # the history extracted
head(input.history)
## # A tibble: 6 x 3
## X__1 X__2 X__3
## <dbl> <dbl> <dbl>
## 1 0 0 0
## 2 1.00 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 1.00 0
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 77
ncol(input.history)
## [1] 3
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
# Get the elevation information
elevation.data <- readxl::read_excel("../BrookTrout.xls",
sheet="Elevation",
na="-",
col_names=TRUE)
elevation.data$Elevation <- elevation.data$Elevation/1000 # standardized it a bit
elevation.data$Site <- 1:nrow(elevation.data)
head(elevation.data)
## # A tibble: 6 x 2
## Elevation Site
## <dbl> <int>
## 1 4.27 1
## 2 4.04 2
## 3 2.03 3
## 4 3.44 4
## 5 3.04 5
## 6 2.05 6
# Get the cross sectional width
cross.data <- readxl::read_excel("../BrookTrout.xls",
sheet="CrossSectionWidth",
na="-",
col_names=TRUE)
head(cross.data)
## # A tibble: 6 x 3
## Cross1 Cross2 Cross3
## <dbl> <dbl> <dbl>
## 1 2.87 2.61 2.73
## 2 0.880 2.54 1.16
## 3 1.50 0.960 1.52
## 4 0.260 1.50 1.13
## 5 0.890 2.59 2.69
## 6 0.670 2.32 1.58
# Convert to a survey covariate. You need to stack the data by columns
survey.cov <- data.frame(Site=1:nrow(cross.data),
visit=as.factor(rep(1:ncol(cross.data), each=nrow(cross.data))),
Cross =unlist(cross.data), stringsAsFactors=FALSE)
head(survey.cov)
## Site visit Cross
## Cross11 1 1 2.87
## Cross12 2 1 0.88
## Cross13 3 1 1.50
## Cross14 4 1 0.26
## Cross15 5 1 0.89
## Cross16 6 1 0.67
# Create the *.pao file
trout.pao <- RPresence::createPao(input.history,
unitcov=elevation.data,
survcov=survey.cov,
title='Brook Trout SSSS')
trout.pao
## $nunits
## [1] 77
##
## $nsurveys
## [1] 3
##
## $nseasons
## [1] 1
##
## $nmethods
## [1] 1
##
## $det.data
## # A tibble: 77 x 3
## X__1 X__2 X__3
## <dbl> <dbl> <dbl>
## 1 0 0 0
## 2 1.00 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 1.00 0
## 7 0 0 0
## 8 0 1.00 1.00
## 9 1.00 0 0
## 10 0 0 0
## # ... with 67 more rows
##
## $nunitcov
## [1] 2
##
## $unitcov
## # A tibble: 77 x 2
## Elevation Site
## <dbl> <int>
## 1 4.27 1
## 2 4.04 2
## 3 2.03 3
## 4 3.44 4
## 5 3.04 5
## 6 2.05 6
## 7 1.82 7
## 8 3.19 8
## 9 4.41 9
## 10 3.40 10
## # ... with 67 more rows
##
## $nsurvcov
## [1] 4
##
## $survcov
## Site visit Cross SURVEY
## Cross11 1 1 2.87 1
## Cross12 2 1 0.88 1
## Cross13 3 1 1.50 1
## Cross14 4 1 0.26 1
## Cross15 5 1 0.89 1
## Cross16 6 1 0.67 1
## Cross17 7 1 2.74 1
## Cross18 8 1 1.65 1
## Cross19 9 1 0.93 1
## Cross110 10 1 1.87 1
## Cross111 11 1 1.19 1
## Cross112 12 1 0.80 1
## Cross113 13 1 1.99 1
## Cross114 14 1 0.78 1
## Cross115 15 1 2.50 1
## Cross116 16 1 0.51 1
## Cross117 17 1 0.84 1
## Cross118 18 1 2.77 1
## Cross119 19 1 2.18 1
## Cross120 20 1 2.95 1
## Cross121 21 1 0.88 1
## Cross122 22 1 1.13 1
## Cross123 23 1 2.34 1
## Cross124 24 1 0.82 1
## Cross125 25 1 2.82 1
## Cross126 26 1 2.27 1
## Cross127 27 1 1.87 1
## Cross128 28 1 1.31 1
## Cross129 29 1 2.56 1
## Cross130 30 1 2.13 1
## Cross131 31 1 1.49 1
## Cross132 32 1 1.62 1
## Cross133 33 1 2.65 1
## Cross134 34 1 0.78 1
## Cross135 35 1 2.43 1
## Cross136 36 1 0.93 1
## Cross137 37 1 1.56 1
## Cross138 38 1 0.71 1
## Cross139 39 1 2.42 1
## Cross140 40 1 1.38 1
## Cross141 41 1 0.59 1
## Cross142 42 1 0.34 1
## Cross143 43 1 0.26 1
## Cross144 44 1 2.37 1
## Cross145 45 1 1.73 1
## Cross146 46 1 1.93 1
## Cross147 47 1 0.33 1
## Cross148 48 1 2.07 1
## Cross149 49 1 0.87 1
## Cross150 50 1 2.50 1
## Cross151 51 1 0.47 1
## Cross152 52 1 1.85 1
## Cross153 53 1 0.25 1
## Cross154 54 1 1.28 1
## Cross155 55 1 1.68 1
## Cross156 56 1 2.48 1
## Cross157 57 1 0.38 1
## Cross158 58 1 2.34 1
## Cross159 59 1 1.48 1
## Cross160 60 1 1.83 1
## Cross161 61 1 1.21 1
## Cross162 62 1 2.85 1
## Cross163 63 1 0.27 1
## Cross164 64 1 1.36 1
## Cross165 65 1 0.60 1
## Cross166 66 1 1.23 1
## Cross167 67 1 2.15 1
## Cross168 68 1 2.87 1
## Cross169 69 1 1.17 1
## Cross170 70 1 0.82 1
## Cross171 71 1 2.83 1
## Cross172 72 1 0.71 1
## Cross173 73 1 0.95 1
## Cross174 74 1 1.67 1
## Cross175 75 1 1.43 1
## Cross176 76 1 1.73 1
## Cross177 77 1 0.52 1
## Cross21 1 2 2.61 2
## Cross22 2 2 2.54 2
## Cross23 3 2 0.96 2
## Cross24 4 2 1.50 2
## Cross25 5 2 2.59 2
## Cross26 6 2 2.32 2
## Cross27 7 2 1.13 2
## Cross28 8 2 1.78 2
## Cross29 9 2 0.86 2
## Cross210 10 2 1.18 2
## Cross211 11 2 2.83 2
## Cross212 12 2 2.89 2
## Cross213 13 2 1.16 2
## Cross214 14 2 2.49 2
## Cross215 15 2 1.51 2
## Cross216 16 2 0.92 2
## Cross217 17 2 2.31 2
## Cross218 18 2 0.66 2
## Cross219 19 2 2.80 2
## Cross220 20 2 1.58 2
## Cross221 21 2 2.42 2
## Cross222 22 2 0.48 2
## Cross223 23 2 2.44 2
## Cross224 24 2 1.88 2
## Cross225 25 2 1.21 2
## Cross226 26 2 0.95 2
## Cross227 27 2 1.97 2
## Cross228 28 2 2.06 2
## Cross229 29 2 1.74 2
## Cross230 30 2 1.71 2
## Cross231 31 2 1.33 2
## Cross232 32 2 1.25 2
## Cross233 33 2 0.31 2
## Cross234 34 2 1.02 2
## Cross235 35 2 2.18 2
## Cross236 36 2 1.12 2
## Cross237 37 2 0.95 2
## Cross238 38 2 0.62 2
## Cross239 39 2 2.78 2
## Cross240 40 2 2.44 2
## Cross241 41 2 0.75 2
## Cross242 42 2 1.89 2
## Cross243 43 2 1.18 2
## Cross244 44 2 1.64 2
## Cross245 45 2 2.57 2
## Cross246 46 2 2.37 2
## Cross247 47 2 1.44 2
## Cross248 48 2 1.24 2
## Cross249 49 2 0.26 2
## Cross250 50 2 0.79 2
## Cross251 51 2 1.72 2
## Cross252 52 2 2.10 2
## Cross253 53 2 2.26 2
## Cross254 54 2 1.39 2
## Cross255 55 2 2.18 2
## Cross256 56 2 1.25 2
## Cross257 57 2 2.50 2
## Cross258 58 2 2.67 2
## Cross259 59 2 2.51 2
## Cross260 60 2 2.83 2
## Cross261 61 2 2.30 2
## Cross262 62 2 0.35 2
## Cross263 63 2 2.11 2
## Cross264 64 2 0.94 2
## Cross265 65 2 1.05 2
## Cross266 66 2 1.04 2
## Cross267 67 2 0.74 2
## Cross268 68 2 1.96 2
## Cross269 69 2 1.05 2
## Cross270 70 2 1.41 2
## Cross271 71 2 0.93 2
## Cross272 72 2 1.13 2
## Cross273 73 2 1.54 2
## Cross274 74 2 2.24 2
## Cross275 75 2 1.85 2
## Cross276 76 2 1.10 2
## Cross277 77 2 2.19 2
## Cross31 1 3 2.73 3
## Cross32 2 3 1.16 3
## Cross33 3 3 1.52 3
## Cross34 4 3 1.13 3
## Cross35 5 3 2.69 3
## Cross36 6 3 1.58 3
## Cross37 7 3 1.73 3
## Cross38 8 3 1.42 3
## Cross39 9 3 1.47 3
## Cross310 10 3 2.37 3
## Cross311 11 3 2.89 3
## Cross312 12 3 0.77 3
## Cross313 13 3 1.51 3
## Cross314 14 3 1.35 3
## Cross315 15 3 2.21 3
## Cross316 16 3 0.78 3
## Cross317 17 3 0.39 3
## Cross318 18 3 0.43 3
## Cross319 19 3 2.57 3
## Cross320 20 3 2.43 3
## Cross321 21 3 2.70 3
## Cross322 22 3 0.72 3
## Cross323 23 3 0.86 3
## Cross324 24 3 0.48 3
## Cross325 25 3 0.88 3
## Cross326 26 3 0.85 3
## Cross327 27 3 2.34 3
## Cross328 28 3 2.52 3
## Cross329 29 3 3.00 3
## Cross330 30 3 2.06 3
## Cross331 31 3 2.51 3
## Cross332 32 3 1.38 3
## Cross333 33 3 1.25 3
## Cross334 34 3 1.21 3
## Cross335 35 3 1.44 3
## Cross336 36 3 0.26 3
## Cross337 37 3 2.51 3
## Cross338 38 3 1.35 3
## Cross339 39 3 1.74 3
## Cross340 40 3 2.02 3
## Cross341 41 3 0.46 3
## Cross342 42 3 0.35 3
## Cross343 43 3 2.76 3
## Cross344 44 3 0.78 3
## Cross345 45 3 2.41 3
## Cross346 46 3 1.22 3
## Cross347 47 3 0.28 3
## Cross348 48 3 0.32 3
## Cross349 49 3 2.97 3
## Cross350 50 3 0.51 3
## Cross351 51 3 2.67 3
## Cross352 52 3 0.92 3
## Cross353 53 3 1.42 3
## Cross354 54 3 0.85 3
## Cross355 55 3 2.03 3
## Cross356 56 3 0.53 3
## Cross357 57 3 0.58 3
## Cross358 58 3 0.72 3
## Cross359 59 3 1.77 3
## Cross360 60 3 1.40 3
## Cross361 61 3 2.76 3
## Cross362 62 3 2.44 3
## Cross363 63 3 2.73 3
## Cross364 64 3 1.36 3
## Cross365 65 3 1.25 3
## Cross366 66 3 2.92 3
## Cross367 67 3 2.71 3
## Cross368 68 3 0.38 3
## Cross369 69 3 1.39 3
## Cross370 70 3 0.86 3
## Cross371 71 3 0.78 3
## Cross372 72 3 2.90 3
## Cross373 73 3 1.05 3
## Cross374 74 3 2.90 3
## Cross375 75 3 1.62 3
## Cross376 76 3 2.25 3
## Cross377 77 3 0.64 3
##
## $nsurveyseason
## [1] 3
##
## $title
## [1] "Brook Trout 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" "unit40" "unit41" "unit42"
## [43] "unit43" "unit44" "unit45" "unit46" "unit47" "unit48" "unit49"
## [50] "unit50" "unit51" "unit52" "unit53" "unit54" "unit55" "unit56"
## [57] "unit57" "unit58" "unit59" "unit60" "unit61" "unit62" "unit63"
## [64] "unit64" "unit65" "unit66" "unit67" "unit68" "unit69" "unit70"
## [71] "unit71" "unit72" "unit73" "unit74" "unit75" "unit76" "unit77"
##
## $surveynames
## [1] "1-1" "1-2" "1-3"
##
## $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 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
## [71] 1 1 1 1 1 1 1
##
## attr(,"class")
## [1] "pao"
# Model where occupancy varies with elevation of stream
mod.elev <- RPresence::occMod(model=list(psi~Elevation, p~1),
type="so", data=trout.pao)
## PRESENCE Version 2.12.21.
summary(mod.elev)
## Model name=psi(Elevation)p()
## AIC=210.1568
## -2*log-likelihood=204.1568
## num. par=3
head(mod.elev$real$psi)
## est se lower_0.95 upper_0.95
## psi_unit1 0.8385830 0.08074740 0.61740368 0.9435824
## psi_unit2 0.7918434 0.08754743 0.57323101 0.9150643
## psi_unit3 0.1953217 0.07164442 0.09036233 0.3722981
## psi_unit4 0.6257368 0.09111721 0.43816822 0.7818603
## psi_unit5 0.4904815 0.08492637 0.33089709 0.6520298
## psi_unit6 0.1993520 0.07204241 0.09321430 0.3762035
mod.elev$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.5363987 0.06035689 0.4182753 0.650574
## p2_unit1 0.5363987 0.06035689 0.4182753 0.650574
## p3_unit1 0.5363987 0.06035689 0.4182753 0.650574
# plot occupancy as a function of elevation of stream
plotdata = data.frame(mod.elev$real$psi, Elevation = elevation.data$Elevation)
ggplot(data=plotdata, aes(x=Elevation, y=est))+
ggtitle("Occupancy as a function of elevation")+
geom_point()+
geom_ribbon(aes(ymin=lower_0.95, ymax=upper_0.95), alpha=0.2)+
ylim(0,1)+
ylab("Estimated occupancy")
