# 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 several models using RPresence
library(car)
## Loading required package: carData
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("../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 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 1 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 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 1 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
## Cro Cross2 Cross3
## <dbl> <dbl> <dbl>
## 1 2.87 2.61 2.73
## 2 0.88 2.54 1.16
## 3 1.5 0.96 1.52
## 4 0.26 1.5 1.13
## 5 0.89 2.59 2.69
## 6 0.67 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
## Cro1 1 1 2.87
## Cro2 2 1 0.88
## Cro3 3 1 1.50
## Cro4 4 1 0.26
## Cro5 5 1 0.89
## Cro6 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 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 1 0
## 7 0 0 0
## 8 0 1 1
## 9 1 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
## Cro1 1 1 2.87 1
## Cro2 2 1 0.88 1
## Cro3 3 1 1.50 1
## Cro4 4 1 0.26 1
## Cro5 5 1 0.89 1
## Cro6 6 1 0.67 1
## Cro7 7 1 2.74 1
## Cro8 8 1 1.65 1
## Cro9 9 1 0.93 1
## Cro10 10 1 1.87 1
## Cro11 11 1 1.19 1
## Cro12 12 1 0.80 1
## Cro13 13 1 1.99 1
## Cro14 14 1 0.78 1
## Cro15 15 1 2.50 1
## Cro16 16 1 0.51 1
## Cro17 17 1 0.84 1
## Cro18 18 1 2.77 1
## Cro19 19 1 2.18 1
## Cro20 20 1 2.95 1
## Cro21 21 1 0.88 1
## Cro22 22 1 1.13 1
## Cro23 23 1 2.34 1
## Cro24 24 1 0.82 1
## Cro25 25 1 2.82 1
## Cro26 26 1 2.27 1
## Cro27 27 1 1.87 1
## Cro28 28 1 1.31 1
## Cro29 29 1 2.56 1
## Cro30 30 1 2.13 1
## Cro31 31 1 1.49 1
## Cro32 32 1 1.62 1
## Cro33 33 1 2.65 1
## Cro34 34 1 0.78 1
## Cro35 35 1 2.43 1
## Cro36 36 1 0.93 1
## Cro37 37 1 1.56 1
## Cro38 38 1 0.71 1
## Cro39 39 1 2.42 1
## Cro40 40 1 1.38 1
## Cro41 41 1 0.59 1
## Cro42 42 1 0.34 1
## Cro43 43 1 0.26 1
## Cro44 44 1 2.37 1
## Cro45 45 1 1.73 1
## Cro46 46 1 1.93 1
## Cro47 47 1 0.33 1
## Cro48 48 1 2.07 1
## Cro49 49 1 0.87 1
## Cro50 50 1 2.50 1
## Cro51 51 1 0.47 1
## Cro52 52 1 1.85 1
## Cro53 53 1 0.25 1
## Cro54 54 1 1.28 1
## Cro55 55 1 1.68 1
## Cro56 56 1 2.48 1
## Cro57 57 1 0.38 1
## Cro58 58 1 2.34 1
## Cro59 59 1 1.48 1
## Cro60 60 1 1.83 1
## Cro61 61 1 1.21 1
## Cro62 62 1 2.85 1
## Cro63 63 1 0.27 1
## Cro64 64 1 1.36 1
## Cro65 65 1 0.60 1
## Cro66 66 1 1.23 1
## Cro67 67 1 2.15 1
## Cro68 68 1 2.87 1
## Cro69 69 1 1.17 1
## Cro70 70 1 0.82 1
## Cro71 71 1 2.83 1
## Cro72 72 1 0.71 1
## Cro73 73 1 0.95 1
## Cro74 74 1 1.67 1
## Cro75 75 1 1.43 1
## Cro76 76 1 1.73 1
## Cro77 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"
# 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
~1, ~Elevation
~Cross, ~1
~Cross, ~Elevation")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p psi
## 1 ~1 ~1
## 2 ~1 ~Elevation
## 3 ~Cross ~1
## 4 ~Cross ~Elevation
# 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=trout.pao)
##
##
## ***** Starting ~1 ~1
## PRESENCE Version 2.12.18.
##
##
## ***** Starting ~1 ~Elevation
## PRESENCE Version 2.12.18.
##
##
## ***** Starting ~Cross ~1
## PRESENCE Version 2.12.18.
##
##
## ***** Starting ~Cross ~Elevation
## 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.150051 0.264924
##
## $psi.VC
## [,1]
## [1,] 0.070185
##
## $p
## est se
## B1_p1 0.134018 0.247821
##
## $p.VC
## [,1]
## [1,] 0.061415
##
## $VC
## A1_psi B1_p1
## A1_psi 0.070185 -0.020671
## B1_p1 -0.020671 0.061415
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.4625575 0.06585972 0.3386551 0.5912636
## psi_unit2 0.4625575 0.06585972 0.3386551 0.5912636
## psi_unit3 0.4625575 0.06585972 0.3386551 0.5912636
## psi_unit4 0.4625575 0.06585972 0.3386551 0.5912636
## psi_unit5 0.4625575 0.06585972 0.3386551 0.5912636
model.fits[[check.model]]$real$p[1:5,]
## est se lower_0.95 upper_0.95
## p1_unit1 0.5334544 0.06167776 0.4129699 0.6501588
## p1_unit2 0.5334544 0.06167776 0.4129699 0.6501588
## p1_unit3 0.5334544 0.06167776 0.4129699 0.6501588
## p1_unit4 0.5334544 0.06167776 0.4129699 0.6501588
## p1_unit5 0.5334544 0.06167776 0.4129699 0.6501588
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 0.08037596 0.0400539 0.02933136 0.2017856
## unit2 1.00000000 0.0000000 1.00000000 1.0000000
## unit3 0.08037596 0.0400539 0.02933136 0.2017856
## unit4 0.08037596 0.0400539 0.02933136 0.2017856
## unit5 0.08037596 0.0400539 0.02933136 0.2017856
## unit6 1.00000000 0.0000000 1.00000000 1.0000000
## unit7 0.08037596 0.0400539 0.02933136 0.2017856
## unit8 1.00000000 0.0000000 1.00000000 1.0000000
## unit9 1.00000000 0.0000000 1.00000000 1.0000000
## unit10 0.08037596 0.0400539 0.02933136 0.2017856
tail(model.fits[[check.model]]$derived$psi_c)
## est se lower_0.95 upper_0.95
## unit72 0.08037596 0.0400539 0.02933136 0.2017856
## unit73 1.00000000 0.0000000 1.00000000 1.0000000
## unit74 1.00000000 0.0000000 1.00000000 1.0000000
## unit75 0.08037596 0.0400539 0.02933136 0.2017856
## unit76 0.08037596 0.0400539 0.02933136 0.2017856
## unit77 1.00000000 0.0000000 1.00000000 1.0000000
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model AIC neg2ll npar warn.conv warn.VC DAIC
## 4 psi(Elevation)p(Cross) 201.7903 193.7903 4 0 0 0.0000
## 2 psi(Elevation)p() 210.1568 204.1568 3 0 0 8.3665
## 3 psi()p(Cross) 226.1193 220.1193 3 0 0 24.3290
## 1 psi()p() 232.7886 228.7886 2 0 0 30.9983
## modlike wgt
## 4 1.0000 0.985
## 2 0.0152 0.015
## 3 0.0000 0.000
## 1 0.0000 0.000
names(aic.table)
## [1] "table" "models" "ess"
# plot occupancy as a function of elevation of stream
psi.ma <- RPresence::modAvg(aic.table, param="psi")[1:5,]
head(psi.ma)
## est se lower_0.95 upper_0.95
## psi_unit1 0.8713520 0.07429068 0.64886094 0.9612795
## psi_unit2 0.8281583 0.08388423 0.60284957 0.9386535
## psi_unit3 0.1925658 0.07345463 0.08632468 0.3757827
## psi_unit4 0.6625619 0.09519064 0.46014489 0.8189460
## psi_unit5 0.5179151 0.09050760 0.34552070 0.6861476
psi.ma$Site <- as.numeric(substring(row.names(psi.ma), 4+regexpr("unit",row.names(psi.ma), fixed=TRUE)))
plotdata <- merge(psi.ma, elevation.data)
head(plotdata)
## Site est se lower_0.95 upper_0.95 Elevation
## 1 1 0.8713520 0.07429068 0.64886094 0.9612795 4.26614
## 2 2 0.8281583 0.08388423 0.60284957 0.9386535 4.03882
## 3 3 0.1925658 0.07345463 0.08632468 0.3757827 2.03158
## 4 4 0.6625619 0.09519064 0.46014489 0.8189460 3.43917
## 5 5 0.5179151 0.09050760 0.34552070 0.6861476 3.03650
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")
# plot detection as a function of cross section
p.ma <- RPresence::modAvg(aic.table, param="p")
p.ma[grepl("unit1$", rownames(p.ma)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.2642356 0.09009285 0.1264934 0.4710791
## p2_unit1 0.3083741 0.08494278 0.1696196 0.4932161
## p3_unit1 0.2874741 0.08756870 0.1486038 0.4825638
head(p.ma)
## est se lower_0.95 upper_0.95
## p1_unit1 0.2642356 0.09009285 0.1264934 0.4710791
## p1_unit2 0.6550509 0.07045799 0.5075408 0.7777266
## p1_unit3 0.5303765 0.06123368 0.4109146 0.6464560
## p1_unit4 0.7610506 0.08195637 0.5683548 0.8851118
## p1_unit5 0.6531576 0.07023225 0.5063307 0.7756634
## p1_unit6 0.6935748 0.07510547 0.5310237 0.8189895
p.ma$Site <- as.numeric(substring(row.names(p.ma), 4+regexpr("unit", row.names(p.ma), fixed=TRUE)))
p.ma$visit<- as.character(substr(row.names(p.ma), 2, -1+regexpr("_", row.names(p.ma), fixed=TRUE)))
survey.cov$visit <- as.character(survey.cov$visit)
plotdata <- merge(p.ma, survey.cov)
head(plotdata)
## Site visit est se lower_0.95 upper_0.95 Cross
## 1 1 1 0.2642356 0.09009285 0.1264934 0.4710791 2.87
## 2 1 2 0.3083741 0.08494278 0.1696196 0.4932161 2.61
## 3 1 3 0.2874741 0.08756870 0.1486038 0.4825638 2.73
## 4 10 1 0.4529641 0.06527355 0.3307033 0.5811760 1.87
## 5 10 2 0.5962904 0.06425397 0.4667710 0.7136496 1.18
## 6 10 3 0.3526559 0.07868118 0.2170575 0.5170249 2.37
ggplot(data=plotdata, aes(x=Cross, y=est))+
ggtitle("Detection probability as a function of cross sectional width")+
geom_point()+
geom_ribbon(aes(ymin=lower_0.95, ymax=upper_0.95), alpha=.2)+
ylim(0,1)+
ylab("Estimated occupancy")
