# 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"
summary(salamander.pao)
## paoname=pres.pao
## title=salamander SSMS
## Naive occ=0.5897436
## nunits nsurveys nseasons nsurveyseason nmethods
## "39" "10" "2" "5,5" "1"
## nunitcov nsurvcov
## "3" "1"
## unit covariates : Site Elevation Prox
## survey covariates: SURVEY
# Define the models.
# model.type do.1 is dynamic occupancy first parameterization
# do.4 is dynamic occupancy 4th parameterization (random occupancy)
# Random occupancy are fit using type="do.4" in the call.
# Parameters are psi, p with gamma=1-epsilon enforced internally
model.list.csv <- textConnection("
p, psi, gamma, epsilon, model.type
~1, ~1, ~1, ~1, do.1
~1, ~Elevation, ~1, ~1, do.1
~1, ~Prox , ~Prox, ~1, do.1")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p psi gamma epsilon model.type
## 1 ~1 ~1 ~1 ~1 do.1
## 2 ~1 ~Elevation ~1 ~1 do.1
## 3 ~1 ~Prox ~Prox ~1 do.1
# fit the model
model.fits <- plyr::alply(model.list, 1, function(x,detect.pao){
cat("\n\n***** Starting ", unlist(x), "\n")
if(x$model.type == 'do.1'){
fit <- RPresence::occMod(model=list(as.formula(paste("psi",x$psi)),
as.formula(paste("p" ,x$p )),
as.formula(paste("gamma",x$gamma)),
as.formula(paste("epsilon",x$epsilon))),
data=detect.pao,type="do.1")
}
if(x$model.type == 'do.4'){
fit <- RPresence::occMod(model=list(as.formula(paste("psi",x$psi)),
as.formula(paste("p" ,x$p ))),
data=detect.pao,type="do.4")
}
fit <- RPresence.add.derived(fit)
fit
},detect.pao=salamander.pao)
##
##
## ***** Starting ~1 ~1 ~1 ~1 do.1
## PRESENCE Version 2.12.21.
## Warning in file.remove(outfile): cannot remove file 'del_18260_62899.out',
## reason 'Permission denied'
## Loading required package: plyr
##
##
## ***** Starting ~1 ~Elevation ~1 ~1 do.1
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~1 ~Prox ~Prox ~1 do.1
## PRESENCE Version 2.12.21.
# Look at output from a specified model
model.number <- 3
names(model.fits[[model.number]])
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "aic" "npar" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
names(model.fits[[model.number]]$real)
## [1] "psi" "gamma" "epsilon" "p" "theta" "th0pi"
model.fits[[model.number]]$beta
## $psi
## psi.coeff
## 1 -1.386639
## 2 5.477777
##
## $psi.VC
## A1 A2
## A1 0.330327 -0.26174
## A2 -0.261740 20.23872
##
## $gamma
## gamma.coeff
## 1 -0.955395
## 2 -21.397920
##
## $gamma.VC
## B1 B2
## B1 0.370072 -3.007130e-01
## B2 -0.300713 6.086494e+10
##
## $epsilon
## epsilon.coeff
## 1 -1.597178
##
## $epsilon.VC
## [1] 0.984154
##
## $p
## p.coeff
## 1 -0.840804
##
## $p.VC
## [1] 0.034062
##
## $VC
## A1 A2 B1 B2 C1 D1
## A1 0.330327 -0.261740 -0.033056 8.276000e-02 -0.017577 -0.012800
## A2 -0.261740 20.238725 0.107520 8.980581e+00 0.909270 -0.172859
## B1 -0.033056 0.107520 0.370072 -3.007130e-01 0.016423 -0.007535
## B2 0.082760 8.980581 -0.300713 6.086494e+10 1.153112 -0.112583
## C1 -0.017577 0.909270 0.016423 1.153112e+00 0.984154 0.046591
## D1 -0.012800 -0.172859 -0.007535 -1.125830e-01 0.046591 0.034062
names(model.fits[[model.number]]$derived)
## [1] "psi" "all_psi" "lambda" "lambdap"
model.fits[[model.number]]$derived$psi[1:10,]
## est se lower_0.95 upper_0.95
## unit1_2 0.8179474 0.1378051 0.4227961 0.9649841
## unit2_2 0.3885351 0.1090791 0.2053214 0.6097861
## unit3_2 0.8179474 0.1378051 0.4227961 0.9649841
## unit4_2 0.8179474 0.1378051 0.4227961 0.9649841
## unit5_2 0.3885351 0.1090791 0.2053214 0.6097861
## unit6_2 0.8179474 0.1378051 0.4227961 0.9649841
## unit7_2 0.8179474 0.1378051 0.4227961 0.9649841
## unit8_2 0.8179474 0.1378051 0.4227961 0.9649841
## unit9_2 0.3885351 0.1090791 0.2053214 0.6097861
## unit10_2 0.3885351 0.1090791 0.2053214 0.6097861
model.fits[[model.number]]$real$gamma[1:5,]
## est se lower_0.95 upper_0.95
## gamma1_unit1 1.959199e-10 4.833506e-05 0.0000000 1.0000000
## gamma1_unit2 2.778011e-01 1.220489e-01 0.1045454 0.5589546
## gamma1_unit3 1.959199e-10 4.833506e-05 0.0000000 1.0000000
## gamma1_unit4 1.959199e-10 4.833506e-05 0.0000000 1.0000000
## gamma1_unit5 2.778011e-01 1.220489e-01 0.1045454 0.5589546
model.fits[[model.number]]$real$epsilon[1:5,]
## est se lower_0.95 upper_0.95
## epsilon1_unit1 0.1683764 0.1389119 0.02815282 0.5859373
## epsilon1_unit2 0.1683764 0.1389119 0.02815282 0.5859373
## epsilon1_unit3 0.1683764 0.1389119 0.02815282 0.5859373
## epsilon1_unit4 0.1683764 0.1389119 0.02815282 0.5859373
## epsilon1_unit5 0.1683764 0.1389119 0.02815282 0.5859373
# Estimate of initial occupance
model.fits[[model.number]]$real$psi[grepl('unit1_', row.names(model.fits[[model.number]]$real$psi)),]
## est se lower_0.95 upper_0.95
## unit1_1 0.9835548 0.07241819 0.009156724 0.9999974
# Derived parameters - estimated occupancy for each unit in years 2....
names(model.fits[[model.number]]$derived)
## [1] "psi" "all_psi" "lambda" "lambdap"
model.fits[[model.number]]$derived$psi[ grepl('unit1_', row.names(model.fits[[model.number]]$derived$psi)),]
## est se lower_0.95 upper_0.95
## unit1_2 0.8179474 0.1378051 0.4227961 0.9649841
# Derived parameters - all of the psi stacked together
model.fits[[model.number]]$derived$all_psi[ grepl('unit1_', row.names(model.fits[[model.number]]$derived$all_psi)),]
## est se lower_0.95 upper_0.95
## unit1_1 0.9835548 0.07241819 0.009156724 0.9999974
## unit1_2 0.8179474 0.13780514 0.422796149 0.9649841
# Estimate of local extinction probability for each unit
model.fits[[model.number]]$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.1683764 0.1389119 0.02815282 0.5859373
# Estimate of local colonization probability for each unit
model.fits[[model.number]]$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 1.959199e-10 4.833506e-05 0 1
# Estimate of probability of detection at each time point for each unit
model.fits[[model.number]]$real$p[ grepl('unit1_', row.names(model.fits[[model.number]]$real$p), fixed=TRUE),]
## [1] est se lower_0.95 upper_0.95
## <0 rows> (or 0-length row.names)
# Get the change in occupancy
# Not yet possible to estimate the se of these values. May have to use bootstrapping.
model.fits[[model.number]]$derived$lambda [grepl('unit1_', row.names(model.fits[[model.number]]$derived$lambda), fixed=TRUE),]
## est se lower_0.95 upper_0.95
## unit1_1 0.8316236 NA NA NA
model.fits[[model.number]]$derived$lambdap[grepl('unit1_', row.names(model.fits[[model.number]]$derived$lambdap), fixed=TRUE),]
## est se lower_0.95 upper_0.95
## unit1_1 0.07512245 NA NA NA
# collect models and make AIC table
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model AIC neg2ll npar warn.conv
## 3 psi(Prox)p()gamma(Prox)epsilon() 308.6992 296.6992 6 0
## 1 psi()p()gamma()epsilon() 324.7730 316.7730 4 0
## 2 psi(Elevation)p()gamma()epsilon() 325.6731 315.6731 5 0
## warn.VC DAIC modlike wgt
## 3 0 0.0000 1e+00 0.9995
## 1 0 16.0738 3e-04 0.0003
## 2 0 16.9739 2e-04 0.0002
# model averaging in the usual way
# initial occupancy
RPresence::modAvg(aic.table, param="psi")[1:5,]
## est se lower_0.95 upper_0.95
## unit1_1 0.9833175 0.07322104 0.009268237 0.9999973
## unit2_1 0.2001200 0.09230403 0.074762127 0.4365073
## unit3_1 0.9833287 0.07314617 0.009303753 0.9999973
## unit4_1 0.9833415 0.07306949 0.009334694 0.9999973
## unit5_1 0.2001523 0.09241956 0.074687732 0.4368712
# model averaging of derived parameters such as the occupancy at each time step
ma_all_psi <- RPresence.modAvg.derived(aic.table, param="all_psi")
psi.est <- ma_all_psi[grepl('unit1_', row.names(ma_all_psi), fixed=TRUE),]
psi.est$Year <- as.numeric(substring(row.names(psi.est),1+regexpr("_",row.names(psi.est))))
psi.est$parameter <- 'psi'
psi.est
## est se lower_0.95 upper_0.95 Year parameter
## unit1_1 0.9833307 0.07262013 0.8409978 1.125664 1 psi
## unit1_2 0.8178207 0.13787825 0.5475843 1.088057 2 psi
# likely more interested in colonization and extinction probabilities
epsilon.ma <- RPresence::modAvg(aic.table, param="epsilon")
epsilon.ma <- epsilon.ma[grepl('unit1$', row.names(epsilon.ma)),]
epsilon.ma$Year <- as.numeric(substr(row.names(epsilon.ma),7+regexpr("epsilon",row.names(epsilon.ma)),-1+regexpr("_",row.names(epsilon.ma))))
epsilon.ma$parameter <- 'epsilon'
epsilon.ma
## est se lower_0.95 upper_0.95 Year parameter
## epsilon1_unit1 0.1683598 0.1389177 0.02814316 0.5859653 1 epsilon
gamma.ma <-RPresence::modAvg(aic.table, param="gamma")
gamma.ma <- gamma.ma[grepl('unit1$', row.names(gamma.ma)),]
gamma.ma$Year <- as.numeric(substr(row.names(gamma.ma),5+regexpr("gamma",row.names(gamma.ma)),-1+regexpr("_",row.names(gamma.ma))))
gamma.ma$parameter <- 'gamma'
gamma.ma
## est se lower_0.95 upper_0.95 Year
## gamma1_unit1 0.0001168308 0.006179183 1.102229e-49 1 1
## parameter
## gamma1_unit1 gamma
all.est <- rbind(psi.est, epsilon.ma, gamma.ma)
all.est$upper_0.95 = ifelse(all.est$upper_0.95 >1,1,all.est$upper_0.95)
ggplot(data=all.est, aes(x=Year,y=est, color=parameter))+
ggtitle("Estimated occupancy, extinction, colonization, over time")+
geom_point(position=position_dodge(w=0.2))+
geom_line(position=position_dodge(w=0.2))+
ylim(0,1)+
geom_errorbar(aes(ymin=lower_0.95, ymax=upper_0.95), width=.1,position=position_dodge(w=0.2))+
scale_x_continuous(breaks=1:10)
