# Single Species, Multi Season Occupancy analyais
# Salamanders over four seasons with covariates on elevation and stream type.
# We will use only the last 2 years (unsure why)
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
# RPresence package
library(readxl)
library(RPresence)
## Warning: package 'RPresence' was built under R version 3.5.0
library(ggplot2)
# Get the RPResence additional functions
source(file.path("..","..","..","AdditionalFunctions","Rpresence.additional.functions.R"))
# 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(file.path("..","SalamanderMultiSeason.xls"),
sheet="Detections",
skip=2, na="-",
col_names=FALSE)
head(input.data)
## # A tibble: 6 x 22
## X__1 X__2 X__3 X__4 X__5 X__6 X__7 X__8 X__9 X__10 X__11 X__12
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00 1.00 0 1.00 0 1.00 1.00 0 1.00 1.00 1.00 0
## 2 2.00 0 0 0 0 0 0 0 0 0 0 0
## 3 3.00 0 1.00 0 1.00 1.00 1.00 0 0 0 1.00 1.00
## 4 4.00 1.00 1.00 1.00 1.00 0 0 0 1.00 1.00 1.00 1.00
## 5 5.00 1.00 1.00 1.00 1.00 0 1.00 0 1.00 0 0 0
## 6 6.00 0 0 0 0 0 0 0 0 1.00 0 0
## # ... with 10 more variables: X__13 <dbl>, X__14 <dbl>, X__15 <dbl>, X__16
## # <dbl>, X__17 <dbl>, X__18 <dbl>, X__19 <dbl>, X__20 <lgl>, X__21
## # <dbl>, X__22 <dbl>
input.history <- input.data[, 10:19]
head(input.history)
## # A tibble: 6 x 10
## X__10 X__11 X__12 X__13 X__14 X__15 X__16 X__17 X__18 X__19
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00 1.00 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 1.00 1.00 0 0 0 0 1.00 0 0
## 4 1.00 1.00 1.00 1.00 0 1.00 1.00 1.00 1.00 0
## 5 0 0 0 0 1.00 0 0 0 0 0
## 6 1.00 0 0 0 0 0 0 1.00 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] 10
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 0
site.covar <- data.frame(Site=1:nrow(input.data),
Elevation=unlist(input.data[,21]),
Prox=car::recode(unlist(input.data[,22]),
"1='Yes'; 0='No';"))
row.names(site.covar) <- NULL
head(site.covar)
## Site Elevation Prox
## 1 1 841.248 Yes
## 2 2 853.440 No
## 3 3 780.288 Yes
## 4 4 707.136 Yes
## 5 5 670.560 No
## 6 6 670.560 Yes
# Number of visits in each season
Nvisits.per.season <- c(5,5)
# Create the *.pao file
salamander.pao <- RPresence::createPao(input.history,
nsurveyseason=Nvisits.per.season,
unitcov=site.covar,
title='salamander SSMS')
salamander.pao
## $nunits
## [1] 39
##
## $nsurveys
## [1] 10
##
## $nseasons
## [1] 2
##
## $nmethods
## [1] 1
##
## $det.data
## # A tibble: 39 x 10
## X__10 X__11 X__12 X__13 X__14 X__15 X__16 X__17 X__18 X__19
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.00 1.00 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 3 0 1.00 1.00 0 0 0 0 1.00 0 0
## 4 1.00 1.00 1.00 1.00 0 1.00 1.00 1.00 1.00 0
## 5 0 0 0 0 1.00 0 0 0 0 0
## 6 1.00 0 0 0 0 0 0 1.00 0 0
## 7 0 0 1.00 1.00 0 0 0 1.00 0 0
## 8 0 0 0 1.00 0 0 0 0 0 0
## 9 0 0 0 0 0 0 0 0 0 0
## 10 0 0 0 0 0 0 0 0 0 0
## # ... with 29 more rows
##
## $nunitcov
## [1] 3
##
## $unitcov
## Site Elevation Prox
## 1 1 841.248 Yes
## 2 2 853.440 No
## 3 3 780.288 Yes
## 4 4 707.136 Yes
## 5 5 670.560 No
## 6 6 670.560 Yes
## 7 7 682.752 Yes
## 8 8 694.944 Yes
## 9 9 768.096 No
## 10 10 768.096 No
## 11 11 743.712 Yes
## 12 12 774.192 No
## 13 13 804.672 No
## 14 14 822.960 No
## 15 15 835.152 No
## 16 16 768.096 No
## 17 17 768.096 No
## 18 18 804.672 No
## 19 19 762.000 Yes
## 20 20 847.344 No
## 21 21 865.632 No
## 22 22 902.208 No
## 23 23 950.976 No
## 24 24 1060.704 No
## 25 25 981.456 Yes
## 26 26 987.552 No
## 27 27 999.744 No
## 28 28 963.168 No
## 29 29 975.360 Yes
## 30 30 969.264 Yes
## 31 31 914.400 No
## 32 32 963.168 Yes
## 33 33 963.168 No
## 34 34 914.400 Yes
## 35 35 792.480 Yes
## 36 36 841.248 Yes
## 37 37 841.248 No
## 38 38 865.632 No
## 39 39 902.208 No
##
## $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
## 196 6
## 197 6
## 198 6
## 199 6
## 200 6
## 201 6
## 202 6
## 203 6
## 204 6
## 205 6
## 206 6
## 207 6
## 208 6
## 209 6
## 210 6
## 211 6
## 212 6
## 213 6
## 214 6
## 215 6
## 216 6
## 217 6
## 218 6
## 219 6
## 220 6
## 221 6
## 222 6
## 223 6
## 224 6
## 225 6
## 226 6
## 227 6
## 228 6
## 229 6
## 230 6
## 231 6
## 232 6
## 233 6
## 234 6
## 235 7
## 236 7
## 237 7
## 238 7
## 239 7
## 240 7
## 241 7
## 242 7
## 243 7
## 244 7
## 245 7
## 246 7
## 247 7
## 248 7
## 249 7
## 250 7
## 251 7
## 252 7
## 253 7
## 254 7
## 255 7
## 256 7
## 257 7
## 258 7
## 259 7
## 260 7
## 261 7
## 262 7
## 263 7
## 264 7
## 265 7
## 266 7
## 267 7
## 268 7
## 269 7
## 270 7
## 271 7
## 272 7
## 273 7
## 274 8
## 275 8
## 276 8
## 277 8
## 278 8
## 279 8
## 280 8
## 281 8
## 282 8
## 283 8
## 284 8
## 285 8
## 286 8
## 287 8
## 288 8
## 289 8
## 290 8
## 291 8
## 292 8
## 293 8
## 294 8
## 295 8
## 296 8
## 297 8
## 298 8
## 299 8
## 300 8
## 301 8
## 302 8
## 303 8
## 304 8
## 305 8
## 306 8
## 307 8
## 308 8
## 309 8
## 310 8
## 311 8
## 312 8
## 313 9
## 314 9
## 315 9
## 316 9
## 317 9
## 318 9
## 319 9
## 320 9
## 321 9
## 322 9
## 323 9
## 324 9
## 325 9
## 326 9
## 327 9
## 328 9
## 329 9
## 330 9
## 331 9
## 332 9
## 333 9
## 334 9
## 335 9
## 336 9
## 337 9
## 338 9
## 339 9
## 340 9
## 341 9
## 342 9
## 343 9
## 344 9
## 345 9
## 346 9
## 347 9
## 348 9
## 349 9
## 350 9
## 351 9
## 352 10
## 353 10
## 354 10
## 355 10
## 356 10
## 357 10
## 358 10
## 359 10
## 360 10
## 361 10
## 362 10
## 363 10
## 364 10
## 365 10
## 366 10
## 367 10
## 368 10
## 369 10
## 370 10
## 371 10
## 372 10
## 373 10
## 374 10
## 375 10
## 376 10
## 377 10
## 378 10
## 379 10
## 380 10
## 381 10
## 382 10
## 383 10
## 384 10
## 385 10
## 386 10
## 387 10
## 388 10
## 389 10
## 390 10
##
## $nsurveyseason
## [1] 5 5
##
## $title
## [1] "salamander SSMS"
##
## $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" "2-1" "2-2" "2-3" "2-4" "2-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 models.
# We start with the first formulation in terms of psi, gamma, epsilon, p (do.1_)
# Note that formula DO NOT HAVE AN = SIGN
mod.fit<- RPresence::occMod(
model=list(psi~1, gamma~1, epsilon~1, p~1),
type="do.1", data=salamander.pao)
## PRESENCE Version 2.12.21.
## Warning in file.remove(outfile): cannot remove file 'del_16768_81425.out',
## reason 'Permission denied'
mod.fit <- RPresence.add.derived(mod.fit) # combine the psi; add the lambda
## Loading required package: plyr
summary(mod.fit)
## Model name=psi()gamma()epsilon()p()
## AIC=324.773
## -2*log-likelihood=316.773
## num. par=4
names(mod.fit)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "aic" "npar" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
names(mod.fit$real)
## [1] "psi" "gamma" "epsilon" "p" "theta" "th0pi"
# Estimate of initial occupancy
mod.fit$real$psi[1,]
## est se lower_0.95 upper_0.95
## unit1_1 0.5313125 0.1026873 0.3356346 0.7178126
# Derived parameters - estimated occupancy for each unit in years 2....
names(mod.fit$derived)
## [1] "psi" "all_psi" "lambda" "lambdap"
mod.fit$derived$psi[ grepl('unit1_', row.names(mod.fit$derived$psi)),]
## est se lower_0.95 upper_0.95
## unit1_2 0.5625662 0.1038638 0.359935 0.7462693
# Additional derived parameters - all of the psi stacked together
mod.fit$derived$all_psi[ grepl('unit1_', row.names(mod.fit$derived$all_psi)),]
## est se lower_0.95 upper_0.95
## unit1_1 0.5313125 0.1026873 0.3356346 0.7178126
## unit1_2 0.5625662 0.1038638 0.3599350 0.7462693
# Estimate of local extinction probability for each unit
mod.fit$real$epsilon[ seq(1, by=nrow(input.history), length.out=length(Nvisits.per.season)-1),]
## est se lower_0.95 upper_0.95
## epsilon1_unit1 0.1396047 0.145713 0.01482902 0.6362385
# Estimate of local colonization probability for each unit
mod.fit$real$gamma[ seq(1, by=nrow(input.history), length.out=length(Nvisits.per.season)-1),]
## est se lower_0.95 upper_0.95
## gamma1_unit1 0.2249419 0.1514079 0.05024835 0.614207
# Estimate of probability of detection at each time point for each unit
mod.fit$real$p[ grepl('_unit1$', row.names(mod.fit$real$p)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p2_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p3_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p4_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p5_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p6_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p7_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p8_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p9_unit1 0.2906618 0.04186267 0.2158314 0.3789006
## p10_unit1 0.2906618 0.04186267 0.2158314 0.3789006
# Derived parameters - estimated occupancy for each unit in years 2....
names(mod.fit$derived)
## [1] "psi" "all_psi" "lambda" "lambdap"
mod.fit$derived$lambda[ grepl('unit1_', row.names(mod.fit$derived$lambda)),]
## est se lower_0.95 upper_0.95
## unit1_1 1.058824 NA NA NA
# Get the change in occupancy
# Not yet possible to estimate the se of these values. May have to use bootstrapping.
mod.fit$derived$lambda [grepl('unit1_', row.names(mod.fit$derived$lambda), fixed=TRUE),]
## est se lower_0.95 upper_0.95
## unit1_1 1.058824 NA NA NA
mod.fit$derived$lambdap[grepl('unit1_', row.names(mod.fit$derived$lambdap), fixed=TRUE),]
## est se lower_0.95 upper_0.95
## unit1_1 1.134474 NA NA NA