# Single Species Single Season Occupancy
# Weta Example
# RPresence package
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
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.history <- readxl::read_excel(file.path("..","weta.xls"),
sheet="detection_histories",
na="-",
col_names=FALSE) # notice no column names in row 1 of data file.
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 72
ncol(input.history)
## [1] 5
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 98
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 0 NA
## 2 0 0 0 0 NA
## 3 0 0 0 1 NA
## 4 0 0 0 0 NA
## 5 0 0 0 0 NA
## 6 0 0 0 0 NA
# Get the site level covariates
site_covar <- readxl::read_excel(file.path("..","weta.xls"),
sheet="site_covar",
na="-",
col_names=TRUE) # notice col_names in row 1 of table.
head(site_covar)
## # A tibble: 6 x 3
## Browsed Unbrowsed BrowseCat
## <dbl> <dbl> <chr>
## 1 1 0 b
## 2 1 0 b
## 3 1 0 b
## 4 0 1 n
## 5 1 0 b
## 6 0 1 n
# Create an alternate site level covariate that is a categorical variable rather
# than indicator variables
site_covar$BrowCat <- site_covar$BrowseCat
xtabs(~BrowCat, data=site_covar,exclude=NULL, na.action=na.pass)
## BrowCat
## b n
## 35 37
head(site_covar)
## # A tibble: 6 x 4
## Browsed Unbrowsed BrowseCat BrowCat
## <dbl> <dbl> <chr> <chr>
## 1 1 0 b b
## 2 1 0 b b
## 3 1 0 b b
## 4 0 1 n n
## 5 1 0 b b
## 6 0 1 n n
# Get the individual covariates.
obs1 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs1",
na="-",
col_names=FALSE)
obs2 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs2",
na="-",
col_names=FALSE)
obs3 <- readxl::read_excel(file.path("..","weta.xls"),
sheet="Obs3",
na="-",
col_names=FALSE)
Obs <- obs1*1 + obs2*2 + obs3*3
head(Obs)
## X__1 X__2 X__3 X__4 X__5
## 1 1 3 2 3 NA
## 2 1 3 2 3 NA
## 3 1 3 2 3 NA
## 4 1 3 2 3 NA
## 5 1 3 2 3 NA
## 6 1 3 2 3 NA
# Observational covariate needs to be "stacked" so that sites1...siteS for survey occastion 1
# are then followed by covariate at survey occastion 2 for sites1...siteS, etc
survey.cov <- data.frame(site=rep(1:nrow(input.history) , ncol(input.history)),
visit=as.character(rep(1:ncol(input.history), each=nrow(input.history))), # notice we make a character
obs1 =as.vector(unlist(obs1)),
obs2 =as.vector(unlist(obs2)),
obs3 =as.vector(unlist(obs3)),
Obs =paste("O",as.vector(unlist(Obs)),sep=""), # notice we make a character string
stringsAsFactors=FALSE)
survey.cov$Obs[ grepl("NA",survey.cov$Obs)] <- NA
head(survey.cov[,c("visit","site","Obs")],n=100)
## visit site Obs
## 1 1 1 O1
## 2 1 2 O1
## 3 1 3 O1
## 4 1 4 O1
## 5 1 5 O1
## 6 1 6 O1
## 7 1 7 O1
## 8 1 8 O1
## 9 1 9 O1
## 10 1 10 O1
## 11 1 11 O1
## 12 1 12 O1
## 13 1 13 O1
## 14 1 14 O1
## 15 1 15 O1
## 16 1 16 <NA>
## 17 1 17 <NA>
## 18 1 18 <NA>
## 19 1 19 <NA>
## 20 1 20 O1
## 21 1 21 O1
## 22 1 22 O1
## 23 1 23 O1
## 24 1 24 O1
## 25 1 25 O1
## 26 1 26 O2
## 27 1 27 O2
## 28 1 28 O2
## 29 1 29 O2
## 30 1 30 O2
## 31 1 31 O2
## 32 1 32 O2
## 33 1 33 O2
## 34 1 34 O2
## 35 1 35 O2
## 36 1 36 O2
## 37 1 37 O2
## 38 1 38 <NA>
## 39 1 39 <NA>
## 40 1 40 <NA>
## 41 1 41 <NA>
## 42 1 42 <NA>
## 43 1 43 <NA>
## 44 1 44 <NA>
## 45 1 45 <NA>
## 46 1 46 <NA>
## 47 1 47 <NA>
## 48 1 48 O3
## 49 1 49 O3
## 50 1 50 O3
## 51 1 51 O3
## 52 1 52 O3
## 53 1 53 O3
## 54 1 54 O3
## 55 1 55 O3
## 56 1 56 O3
## 57 1 57 O3
## 58 1 58 O3
## 59 1 59 O3
## 60 1 60 O3
## 61 1 61 O3
## 62 1 62 O3
## 63 1 63 <NA>
## 64 1 64 <NA>
## 65 1 65 <NA>
## 66 1 66 <NA>
## 67 1 67 <NA>
## 68 1 68 O3
## 69 1 69 O3
## 70 1 70 O3
## 71 1 71 O3
## 72 1 72 O3
## 73 2 1 O3
## 74 2 2 O3
## 75 2 3 O3
## 76 2 4 O3
## 77 2 5 O3
## 78 2 6 O3
## 79 2 7 O3
## 80 2 8 O3
## 81 2 9 O3
## 82 2 10 O3
## 83 2 11 O3
## 84 2 12 O3
## 85 2 13 O3
## 86 2 14 O3
## 87 2 15 O3
## 88 2 16 O3
## 89 2 17 O3
## 90 2 18 O3
## 91 2 19 O3
## 92 2 20 O3
## 93 2 21 <NA>
## 94 2 22 <NA>
## 95 2 23 <NA>
## 96 2 24 <NA>
## 97 2 25 <NA>
## 98 2 26 O1
## 99 2 27 O1
## 100 2 28 O1
str(survey.cov)
## 'data.frame': 360 obs. of 6 variables:
## $ site : int 1 2 3 4 5 6 7 8 9 10 ...
## $ visit: chr "1" "1" "1" "1" ...
## $ obs1 : num 1 1 1 1 1 1 1 1 1 1 ...
## $ obs2 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ obs3 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Obs : chr "O1" "O1" "O1" "O1" ...
# check that missing values in history and observer covariates align
select <- is.na(as.vector(unlist(input.history)))
survey.cov[select,]
## site visit obs1 obs2 obs3 Obs
## 16 16 1 NA NA NA <NA>
## 17 17 1 NA NA NA <NA>
## 18 18 1 NA NA NA <NA>
## 19 19 1 NA NA NA <NA>
## 38 38 1 NA NA NA <NA>
## 39 39 1 NA NA NA <NA>
## 40 40 1 NA NA NA <NA>
## 41 41 1 NA NA NA <NA>
## 42 42 1 NA NA NA <NA>
## 43 43 1 NA NA NA <NA>
## 44 44 1 NA NA NA <NA>
## 45 45 1 NA NA NA <NA>
## 46 46 1 NA NA NA <NA>
## 47 47 1 NA NA NA <NA>
## 63 63 1 NA NA NA <NA>
## 64 64 1 NA NA NA <NA>
## 65 65 1 NA NA NA <NA>
## 66 66 1 NA NA NA <NA>
## 67 67 1 NA NA NA <NA>
## 93 21 2 NA NA NA <NA>
## 94 22 2 NA NA NA <NA>
## 95 23 2 NA NA NA <NA>
## 96 24 2 NA NA NA <NA>
## 97 25 2 NA NA NA <NA>
## 115 43 2 NA NA NA <NA>
## 116 44 2 NA NA NA <NA>
## 117 45 2 NA NA NA <NA>
## 118 46 2 NA NA NA <NA>
## 119 47 2 NA NA NA <NA>
## 140 68 2 NA NA NA <NA>
## 141 69 2 NA NA NA <NA>
## 142 70 2 NA NA NA <NA>
## 143 71 2 NA NA NA <NA>
## 144 72 2 NA NA NA <NA>
## 155 11 3 NA NA NA <NA>
## 156 12 3 NA NA NA <NA>
## 157 13 3 NA NA NA <NA>
## 158 14 3 NA NA NA <NA>
## 159 15 3 NA NA NA <NA>
## 182 38 3 NA NA NA <NA>
## 183 39 3 NA NA NA <NA>
## 202 58 3 NA NA NA <NA>
## 203 59 3 NA NA NA <NA>
## 204 60 3 NA NA NA <NA>
## 205 61 3 NA NA NA <NA>
## 206 62 3 NA NA NA <NA>
## 232 16 4 NA NA NA <NA>
## 233 17 4 NA NA NA <NA>
## 234 18 4 NA NA NA <NA>
## 235 19 4 NA NA NA <NA>
## 236 20 4 NA NA NA <NA>
## 237 21 4 NA NA NA <NA>
## 238 22 4 NA NA NA <NA>
## 239 23 4 NA NA NA <NA>
## 240 24 4 NA NA NA <NA>
## 241 25 4 NA NA NA <NA>
## 242 26 4 NA NA NA <NA>
## 243 27 4 NA NA NA <NA>
## 244 28 4 NA NA NA <NA>
## 245 29 4 NA NA NA <NA>
## 279 63 4 NA NA NA <NA>
## 280 64 4 NA NA NA <NA>
## 281 65 4 NA NA NA <NA>
## 282 66 4 NA NA NA <NA>
## 283 67 4 NA NA NA <NA>
## 284 68 4 NA NA NA <NA>
## 285 69 4 NA NA NA <NA>
## 286 70 4 NA NA NA <NA>
## 287 71 4 NA NA NA <NA>
## 288 72 4 NA NA NA <NA>
## 289 1 5 NA NA NA <NA>
## 290 2 5 NA NA NA <NA>
## 291 3 5 NA NA NA <NA>
## 292 4 5 NA NA NA <NA>
## 293 5 5 NA NA NA <NA>
## 294 6 5 NA NA NA <NA>
## 295 7 5 NA NA NA <NA>
## 296 8 5 NA NA NA <NA>
## 297 9 5 NA NA NA <NA>
## 298 10 5 NA NA NA <NA>
## 314 26 5 NA NA NA <NA>
## 315 27 5 NA NA NA <NA>
## 316 28 5 NA NA NA <NA>
## 317 29 5 NA NA NA <NA>
## 318 30 5 NA NA NA <NA>
## 319 31 5 NA NA NA <NA>
## 320 32 5 NA NA NA <NA>
## 321 33 5 NA NA NA <NA>
## 336 48 5 NA NA NA <NA>
## 337 49 5 NA NA NA <NA>
## 338 50 5 NA NA NA <NA>
## 339 51 5 NA NA NA <NA>
## 340 52 5 NA NA NA <NA>
## 341 53 5 NA NA NA <NA>
## 342 54 5 NA NA NA <NA>
## 343 55 5 NA NA NA <NA>
## 344 56 5 NA NA NA <NA>
## 345 57 5 NA NA NA <NA>
sum(is.na(survey.cov[!select,]))
## [1] 0
# Create the *.pao file
weta.pao <- RPresence::createPao(input.history,
unitcov=site_covar,
survcov=survey.cov,
title='weta SSSS')
weta.pao
## $nunits
## [1] 72
##
## $nsurveys
## [1] 5
##
## $nseasons
## [1] 1
##
## $nmethods
## [1] 1
##
## $det.data
## # A tibble: 72 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 NA
## 2 0 0 0 0 NA
## 3 0 0 0 1 NA
## 4 0 0 0 0 NA
## 5 0 0 0 0 NA
## 6 0 0 0 0 NA
## 7 0 0 0 0 NA
## 8 0 0 0 0 NA
## 9 0 0 0 0 NA
## 10 1 1 0 0 NA
## # ... with 62 more rows
##
## $nunitcov
## [1] 4
##
## $unitcov
## # A tibble: 72 x 4
## Browsed Unbrowsed BrowseCat BrowCat
## <dbl> <dbl> <chr> <chr>
## 1 1 0 b b
## 2 1 0 b b
## 3 1 0 b b
## 4 0 1 n n
## 5 1 0 b b
## 6 0 1 n n
## 7 0 1 n n
## 8 0 1 n n
## 9 0 1 n n
## 10 1 0 b b
## # ... with 62 more rows
##
## $nsurvcov
## [1] 7
##
## $survcov
## site visit obs1 obs2 obs3 Obs SURVEY
## 1 1 1 1 0 0 O1 1
## 2 2 1 1 0 0 O1 1
## 3 3 1 1 0 0 O1 1
## 4 4 1 1 0 0 O1 1
## 5 5 1 1 0 0 O1 1
## 6 6 1 1 0 0 O1 1
## 7 7 1 1 0 0 O1 1
## 8 8 1 1 0 0 O1 1
## 9 9 1 1 0 0 O1 1
## 10 10 1 1 0 0 O1 1
## 11 11 1 1 0 0 O1 1
## 12 12 1 1 0 0 O1 1
## 13 13 1 1 0 0 O1 1
## 14 14 1 1 0 0 O1 1
## 15 15 1 1 0 0 O1 1
## 16 16 1 NA NA NA <NA> 1
## 17 17 1 NA NA NA <NA> 1
## 18 18 1 NA NA NA <NA> 1
## 19 19 1 NA NA NA <NA> 1
## 20 20 1 1 0 0 O1 1
## 21 21 1 1 0 0 O1 1
## 22 22 1 1 0 0 O1 1
## 23 23 1 1 0 0 O1 1
## 24 24 1 1 0 0 O1 1
## 25 25 1 1 0 0 O1 1
## 26 26 1 0 1 0 O2 1
## 27 27 1 0 1 0 O2 1
## 28 28 1 0 1 0 O2 1
## 29 29 1 0 1 0 O2 1
## 30 30 1 0 1 0 O2 1
## 31 31 1 0 1 0 O2 1
## 32 32 1 0 1 0 O2 1
## 33 33 1 0 1 0 O2 1
## 34 34 1 0 1 0 O2 1
## 35 35 1 0 1 0 O2 1
## 36 36 1 0 1 0 O2 1
## 37 37 1 0 1 0 O2 1
## 38 38 1 NA NA NA <NA> 1
## 39 39 1 NA NA NA <NA> 1
## 40 40 1 NA NA NA <NA> 1
## 41 41 1 NA NA NA <NA> 1
## 42 42 1 NA NA NA <NA> 1
## 43 43 1 NA NA NA <NA> 1
## 44 44 1 NA NA NA <NA> 1
## 45 45 1 NA NA NA <NA> 1
## 46 46 1 NA NA NA <NA> 1
## 47 47 1 NA NA NA <NA> 1
## 48 48 1 0 0 1 O3 1
## 49 49 1 0 0 1 O3 1
## 50 50 1 0 0 1 O3 1
## 51 51 1 0 0 1 O3 1
## 52 52 1 0 0 1 O3 1
## 53 53 1 0 0 1 O3 1
## 54 54 1 0 0 1 O3 1
## 55 55 1 0 0 1 O3 1
## 56 56 1 0 0 1 O3 1
## 57 57 1 0 0 1 O3 1
## 58 58 1 0 0 1 O3 1
## 59 59 1 0 0 1 O3 1
## 60 60 1 0 0 1 O3 1
## 61 61 1 0 0 1 O3 1
## 62 62 1 0 0 1 O3 1
## 63 63 1 NA NA NA <NA> 1
## 64 64 1 NA NA NA <NA> 1
## 65 65 1 NA NA NA <NA> 1
## 66 66 1 NA NA NA <NA> 1
## 67 67 1 NA NA NA <NA> 1
## 68 68 1 0 0 1 O3 1
## 69 69 1 0 0 1 O3 1
## 70 70 1 0 0 1 O3 1
## 71 71 1 0 0 1 O3 1
## 72 72 1 0 0 1 O3 1
## 73 1 2 0 0 1 O3 2
## 74 2 2 0 0 1 O3 2
## 75 3 2 0 0 1 O3 2
## 76 4 2 0 0 1 O3 2
## 77 5 2 0 0 1 O3 2
## 78 6 2 0 0 1 O3 2
## 79 7 2 0 0 1 O3 2
## 80 8 2 0 0 1 O3 2
## 81 9 2 0 0 1 O3 2
## 82 10 2 0 0 1 O3 2
## 83 11 2 0 0 1 O3 2
## 84 12 2 0 0 1 O3 2
## 85 13 2 0 0 1 O3 2
## 86 14 2 0 0 1 O3 2
## 87 15 2 0 0 1 O3 2
## 88 16 2 0 0 1 O3 2
## 89 17 2 0 0 1 O3 2
## 90 18 2 0 0 1 O3 2
## 91 19 2 0 0 1 O3 2
## 92 20 2 0 0 1 O3 2
## 93 21 2 NA NA NA <NA> 2
## 94 22 2 NA NA NA <NA> 2
## 95 23 2 NA NA NA <NA> 2
## 96 24 2 NA NA NA <NA> 2
## 97 25 2 NA NA NA <NA> 2
## 98 26 2 1 0 0 O1 2
## 99 27 2 1 0 0 O1 2
## 100 28 2 1 0 0 O1 2
## 101 29 2 1 0 0 O1 2
## 102 30 2 1 0 0 O1 2
## 103 31 2 1 0 0 O1 2
## 104 32 2 1 0 0 O1 2
## 105 33 2 1 0 0 O1 2
## 106 34 2 1 0 0 O1 2
## 107 35 2 1 0 0 O1 2
## 108 36 2 1 0 0 O1 2
## 109 37 2 1 0 0 O1 2
## 110 38 2 1 0 0 O1 2
## 111 39 2 1 0 0 O1 2
## 112 40 2 1 0 0 O1 2
## 113 41 2 1 0 0 O1 2
## 114 42 2 1 0 0 O1 2
## 115 43 2 NA NA NA <NA> 2
## 116 44 2 NA NA NA <NA> 2
## 117 45 2 NA NA NA <NA> 2
## 118 46 2 NA NA NA <NA> 2
## 119 47 2 NA NA NA <NA> 2
## 120 48 2 0 1 0 O2 2
## 121 49 2 0 1 0 O2 2
## 122 50 2 0 1 0 O2 2
## 123 51 2 0 1 0 O2 2
## 124 52 2 0 1 0 O2 2
## 125 53 2 0 1 0 O2 2
## 126 54 2 0 1 0 O2 2
## 127 55 2 0 1 0 O2 2
## 128 56 2 0 1 0 O2 2
## 129 57 2 0 1 0 O2 2
## 130 58 2 0 1 0 O2 2
## 131 59 2 0 1 0 O2 2
## 132 60 2 0 1 0 O2 2
## 133 61 2 0 1 0 O2 2
## 134 62 2 0 1 0 O2 2
## 135 63 2 0 1 0 O2 2
## 136 64 2 0 1 0 O2 2
## 137 65 2 0 1 0 O2 2
## 138 66 2 0 1 0 O2 2
## 139 67 2 0 1 0 O2 2
## 140 68 2 NA NA NA <NA> 2
## 141 69 2 NA NA NA <NA> 2
## 142 70 2 NA NA NA <NA> 2
## 143 71 2 NA NA NA <NA> 2
## 144 72 2 NA NA NA <NA> 2
## 145 1 3 0 1 0 O2 3
## 146 2 3 0 1 0 O2 3
## 147 3 3 0 1 0 O2 3
## 148 4 3 0 1 0 O2 3
## 149 5 3 0 1 0 O2 3
## 150 6 3 0 1 0 O2 3
## 151 7 3 0 1 0 O2 3
## 152 8 3 0 1 0 O2 3
## 153 9 3 0 1 0 O2 3
## 154 10 3 0 1 0 O2 3
## 155 11 3 NA NA NA <NA> 3
## 156 12 3 NA NA NA <NA> 3
## 157 13 3 NA NA NA <NA> 3
## 158 14 3 NA NA NA <NA> 3
## 159 15 3 NA NA NA <NA> 3
## 160 16 3 0 1 0 O2 3
## 161 17 3 0 1 0 O2 3
## 162 18 3 0 1 0 O2 3
## 163 19 3 0 1 0 O2 3
## 164 20 3 0 1 0 O2 3
## 165 21 3 0 1 0 O2 3
## 166 22 3 0 1 0 O2 3
## 167 23 3 0 1 0 O2 3
## 168 24 3 0 1 0 O2 3
## 169 25 3 0 1 0 O2 3
## 170 26 3 0 0 1 O3 3
## 171 27 3 0 0 1 O3 3
## 172 28 3 0 0 1 O3 3
## 173 29 3 0 0 1 O3 3
## 174 30 3 0 0 1 O3 3
## 175 31 3 0 0 1 O3 3
## 176 32 3 0 0 1 O3 3
## 177 33 3 0 0 1 O3 3
## 178 34 3 0 0 1 O3 3
## 179 35 3 0 0 1 O3 3
## 180 36 3 0 0 1 O3 3
## 181 37 3 0 0 1 O3 3
## 182 38 3 NA NA NA <NA> 3
## 183 39 3 NA NA NA <NA> 3
## 184 40 3 0 0 1 O3 3
## 185 41 3 0 0 1 O3 3
## 186 42 3 0 0 1 O3 3
## 187 43 3 0 0 1 O3 3
## 188 44 3 0 0 1 O3 3
## 189 45 3 0 0 1 O3 3
## 190 46 3 0 0 1 O3 3
## 191 47 3 0 0 1 O3 3
## 192 48 3 1 0 0 O1 3
## 193 49 3 1 0 0 O1 3
## 194 50 3 1 0 0 O1 3
## 195 51 3 1 0 0 O1 3
## 196 52 3 1 0 0 O1 3
## 197 53 3 1 0 0 O1 3
## 198 54 3 1 0 0 O1 3
## 199 55 3 1 0 0 O1 3
## 200 56 3 1 0 0 O1 3
## 201 57 3 1 0 0 O1 3
## 202 58 3 NA NA NA <NA> 3
## 203 59 3 NA NA NA <NA> 3
## 204 60 3 NA NA NA <NA> 3
## 205 61 3 NA NA NA <NA> 3
## 206 62 3 NA NA NA <NA> 3
## 207 63 3 1 0 0 O1 3
## 208 64 3 1 0 0 O1 3
## 209 65 3 1 0 0 O1 3
## 210 66 3 1 0 0 O1 3
## 211 67 3 1 0 0 O1 3
## 212 68 3 1 0 0 O1 3
## 213 69 3 1 0 0 O1 3
## 214 70 3 1 0 0 O1 3
## 215 71 3 1 0 0 O1 3
## 216 72 3 1 0 0 O1 3
## 217 1 4 0 0 1 O3 4
## 218 2 4 0 0 1 O3 4
## 219 3 4 0 0 1 O3 4
## 220 4 4 0 0 1 O3 4
## 221 5 4 0 0 1 O3 4
## 222 6 4 0 0 1 O3 4
## 223 7 4 0 0 1 O3 4
## 224 8 4 0 0 1 O3 4
## 225 9 4 0 0 1 O3 4
## 226 10 4 0 0 1 O3 4
## 227 11 4 0 0 1 O3 4
## 228 12 4 0 0 1 O3 4
## 229 13 4 0 0 1 O3 4
## 230 14 4 0 0 1 O3 4
## 231 15 4 0 0 1 O3 4
## 232 16 4 NA NA NA <NA> 4
## 233 17 4 NA NA NA <NA> 4
## 234 18 4 NA NA NA <NA> 4
## 235 19 4 NA NA NA <NA> 4
## 236 20 4 NA NA NA <NA> 4
## 237 21 4 NA NA NA <NA> 4
## 238 22 4 NA NA NA <NA> 4
## 239 23 4 NA NA NA <NA> 4
## 240 24 4 NA NA NA <NA> 4
## 241 25 4 NA NA NA <NA> 4
## 242 26 4 NA NA NA <NA> 4
## 243 27 4 NA NA NA <NA> 4
## 244 28 4 NA NA NA <NA> 4
## 245 29 4 NA NA NA <NA> 4
## 246 30 4 1 0 0 O1 4
## 247 31 4 1 0 0 O1 4
## 248 32 4 1 0 0 O1 4
## 249 33 4 1 0 0 O1 4
## 250 34 4 1 0 0 O1 4
## 251 35 4 1 0 0 O1 4
## 252 36 4 1 0 0 O1 4
## 253 37 4 1 0 0 O1 4
## 254 38 4 1 0 0 O1 4
## 255 39 4 1 0 0 O1 4
## 256 40 4 1 0 0 O1 4
## 257 41 4 1 0 0 O1 4
## 258 42 4 1 0 0 O1 4
## 259 43 4 1 0 0 O1 4
## 260 44 4 1 0 0 O1 4
## 261 45 4 1 0 0 O1 4
## 262 46 4 1 0 0 O1 4
## 263 47 4 1 0 0 O1 4
## 264 48 4 0 1 0 O2 4
## 265 49 4 0 1 0 O2 4
## 266 50 4 0 1 0 O2 4
## 267 51 4 0 1 0 O2 4
## 268 52 4 0 1 0 O2 4
## 269 53 4 0 1 0 O2 4
## 270 54 4 0 1 0 O2 4
## 271 55 4 0 1 0 O2 4
## 272 56 4 0 1 0 O2 4
## 273 57 4 0 1 0 O2 4
## 274 58 4 0 1 0 O2 4
## 275 59 4 0 1 0 O2 4
## 276 60 4 0 1 0 O2 4
## 277 61 4 0 1 0 O2 4
## 278 62 4 0 1 0 O2 4
## 279 63 4 NA NA NA <NA> 4
## 280 64 4 NA NA NA <NA> 4
## 281 65 4 NA NA NA <NA> 4
## 282 66 4 NA NA NA <NA> 4
## 283 67 4 NA NA NA <NA> 4
## 284 68 4 NA NA NA <NA> 4
## 285 69 4 NA NA NA <NA> 4
## 286 70 4 NA NA NA <NA> 4
## 287 71 4 NA NA NA <NA> 4
## 288 72 4 NA NA NA <NA> 4
## 289 1 5 NA NA NA <NA> 5
## 290 2 5 NA NA NA <NA> 5
## 291 3 5 NA NA NA <NA> 5
## 292 4 5 NA NA NA <NA> 5
## 293 5 5 NA NA NA <NA> 5
## 294 6 5 NA NA NA <NA> 5
## 295 7 5 NA NA NA <NA> 5
## 296 8 5 NA NA NA <NA> 5
## 297 9 5 NA NA NA <NA> 5
## 298 10 5 NA NA NA <NA> 5
## 299 11 5 1 0 0 O1 5
## 300 12 5 1 0 0 O1 5
## 301 13 5 1 0 0 O1 5
## 302 14 5 1 0 0 O1 5
## 303 15 5 1 0 0 O1 5
## 304 16 5 1 0 0 O1 5
## 305 17 5 1 0 0 O1 5
## 306 18 5 1 0 0 O1 5
## 307 19 5 1 0 0 O1 5
## 308 20 5 1 0 0 O1 5
## 309 21 5 1 0 0 O1 5
## 310 22 5 1 0 0 O1 5
## 311 23 5 1 0 0 O1 5
## 312 24 5 1 0 0 O1 5
## 313 25 5 1 0 0 O1 5
## 314 26 5 NA NA NA <NA> 5
## 315 27 5 NA NA NA <NA> 5
## 316 28 5 NA NA NA <NA> 5
## 317 29 5 NA NA NA <NA> 5
## 318 30 5 NA NA NA <NA> 5
## 319 31 5 NA NA NA <NA> 5
## 320 32 5 NA NA NA <NA> 5
## 321 33 5 NA NA NA <NA> 5
## 322 34 5 0 1 0 O2 5
## 323 35 5 0 1 0 O2 5
## 324 36 5 0 1 0 O2 5
## 325 37 5 0 1 0 O2 5
## 326 38 5 0 1 0 O2 5
## 327 39 5 0 1 0 O2 5
## 328 40 5 0 1 0 O2 5
## 329 41 5 0 1 0 O2 5
## 330 42 5 0 1 0 O2 5
## 331 43 5 0 1 0 O2 5
## 332 44 5 0 1 0 O2 5
## 333 45 5 0 1 0 O2 5
## 334 46 5 0 1 0 O2 5
## 335 47 5 0 1 0 O2 5
## 336 48 5 NA NA NA <NA> 5
## 337 49 5 NA NA NA <NA> 5
## 338 50 5 NA NA NA <NA> 5
## 339 51 5 NA NA NA <NA> 5
## 340 52 5 NA NA NA <NA> 5
## 341 53 5 NA NA NA <NA> 5
## 342 54 5 NA NA NA <NA> 5
## 343 55 5 NA NA NA <NA> 5
## 344 56 5 NA NA NA <NA> 5
## 345 57 5 NA NA NA <NA> 5
## 346 58 5 0 0 1 O3 5
## 347 59 5 0 0 1 O3 5
## 348 60 5 0 0 1 O3 5
## 349 61 5 0 0 1 O3 5
## 350 62 5 0 0 1 O3 5
## 351 63 5 0 0 1 O3 5
## 352 64 5 0 0 1 O3 5
## 353 65 5 0 0 1 O3 5
## 354 66 5 0 0 1 O3 5
## 355 67 5 0 0 1 O3 5
## 356 68 5 0 0 1 O3 5
## 357 69 5 0 0 1 O3 5
## 358 70 5 0 0 1 O3 5
## 359 71 5 0 0 1 O3 5
## 360 72 5 0 0 1 O3 5
##
## $nsurveyseason
## [1] 5
##
## $title
## [1] "weta 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"
##
## $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 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
##
## attr(,"class")
## [1] "pao"
# Fit some models.
# Note that formula DO NOT HAVE AN = SIGN
mod.pdot <- RPresence::occMod(model=list(psi~1, p~1), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pdot)
## Model name=psi()p()
## AIC=265.7872
## -2*log-likelihood=261.7872
## num. par=2
# look at estimated occupancy probability. RPresence gives for EACH site in case it depends on covariates
mod.pdot.psi <-mod.pdot$real$psi[1,] # occupancy probability
mod.pdot.psi
## est se lower_0.95 upper_0.95
## psi_unit1 0.6165958 0.08854195 0.4356219 0.770157
# look at the estimated probability of detection. It gives an estimate for every site at very visit
mod.pdot.p <- mod.pdot$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
mod.pdot.p
## est se lower_0.95 upper_0.95
## p1_unit1 0.3493651 0.05372869 0.2525413 0.4604435
## p2_unit1 0.3493651 0.05372869 0.2525413 0.4604435
## p3_unit1 0.3493651 0.05372869 0.2525413 0.4604435
## p4_unit1 0.3493651 0.05372869 0.2525413 0.4604435
## p5_unit1 0.3493651 0.05372869 0.2525413 0.4604435
# alternatively
RPresence::print_one_site_estimates(mod.pdot, site = 1)
## psi()p()
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.6165958 0.08854195 0.4356219 0.7701570
## p_p1_unit1 0.3493651 0.05372869 0.2525413 0.4604435
# Impact of browse.
# Use the two indicator variables or the categorized column
# This give identical results
# Cell means models.
# We need to EXCLUDE the intercept if we are using the two indicator variables.
mod.pdot.psiB.1 <- RPresence::occMod(model=list(psi~-1+Browsed+Unbrowsed, p~1), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pdot.psiB.1)
## Model name=psi(-1 P Browsed P Unbrowsed)p()
## AIC=264.2643
## -2*log-likelihood=258.2643
## num. par=3
names(mod.pdot.psiB.1)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.pdot.psiB.1.psi <-mod.pdot.psiB.1$real$psi # occupancy probability
mod.pdot.psiB.1.psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7593708 0.1198231 0.4660490 0.9194190
## psi_unit2 0.7593708 0.1198231 0.4660490 0.9194190
## psi_unit3 0.7593708 0.1198231 0.4660490 0.9194190
## psi_unit4 0.4809982 0.1078853 0.2843322 0.6837339
## psi_unit5 0.7593708 0.1198231 0.4660490 0.9194190
# This is the logit occupancy for each browse category
mod.pdot.psiB.1$dmat$psi
## a1 a2
## psi "psi.Browsed" "psi.Unbrowsed"
mod.pdot.psiB.1$beta$psi
## est se
## A1_psi.psi.Browsed 1.149233 0.655750
## A2_psi.psi.Unbrowsed -0.076044 0.432165
mod.pdot.psiB.1.oddsratio.browse <- exp( sum(c(1,-1)*mod.pdot.psiB.1$beta$psi))
mod.pdot.psiB.1.oddsratio.browse
## [1] 4.258266
# Alternate way for cell means models.
mod.pdot.psiB.4 <- RPresence::occMod(model=list(psi~-1+BrowCat, p~1), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pdot.psiB.4)
## Model name=psi(-1 P BrowCat)p()
## AIC=264.2643
## -2*log-likelihood=258.2643
## num. par=3
names(mod.pdot.psiB.4)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.pdot.psiB.4.psi <-mod.pdot.psiB.4$real$psi # occupancy probability
mod.pdot.psiB.4.psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7593708 0.1198231 0.4660490 0.9194190
## psi_unit2 0.7593708 0.1198231 0.4660490 0.9194190
## psi_unit3 0.7593708 0.1198231 0.4660490 0.9194190
## psi_unit4 0.4809982 0.1078853 0.2843322 0.6837339
## psi_unit5 0.7593708 0.1198231 0.4660490 0.9194190
# This is the logit occupancy for each browse category
mod.pdot.psiB.4$dmat$psi
## a1 a2
## psi "psi.BrowCatb" "psi.BrowCatn"
mod.pdot.psiB.4$beta$psi
## est se
## A1_psi.psi.BrowCatb 1.149233 0.655750
## A2_psi.psi.BrowCatn -0.076044 0.432165
mod.pdot.psiB.4.oddsratio.browse <- exp( sum(c(1,-1)*mod.pdot.psiB.4$beta$psi))
mod.pdot.psiB.4.oddsratio.browse
## [1] 4.258266
# We can use just a single Browser indicator which is the cell effects model.
mod.pdot.psiB.2 <- RPresence::occMod(model=list(psi~Browsed, p~1), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pdot.psiB.2)
## Model name=psi(Browsed)p()
## AIC=264.2643
## -2*log-likelihood=258.2643
## num. par=3
names(mod.pdot.psiB.2)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.pdot.psiB.2.psi <-mod.pdot.psiB.2$real$psi # occupancy probability
mod.pdot.psiB.2.psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7593710 0.1198227 0.4660500 0.9194189
## psi_unit2 0.7593710 0.1198227 0.4660500 0.9194189
## psi_unit3 0.7593710 0.1198227 0.4660500 0.9194189
## psi_unit4 0.4809982 0.1078871 0.2843294 0.6837368
## psi_unit5 0.7593710 0.1198227 0.4660500 0.9194189
# This is the logit occupancy the baseline and the log(odds ratio) of the two
# browse classes
mod.pdot.psiB.2$dmat$psi
## a1 a2
## psi "1" "psi.Browsed"
mod.pdot.psiB.2$beta$psi
## est se
## A1_psi -0.076044 0.432172
## A2_psi.psi.Browsed 1.225278 0.720562
mod.pdot.psiB.2.oddsratio.browse <- exp( sum(c(0,1)*mod.pdot.psiB.2$beta$psi))
mod.pdot.psiB.2.oddsratio.browse
## [1] 6.999509
# We can use just a categorical browser variable (preferred).
mod.pdot.psiB.3 <- RPresence::occMod(model=list(psi~BrowCat, p~1), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pdot.psiB.3)
## Model name=psi(BrowCat)p()
## AIC=264.2643
## -2*log-likelihood=258.2643
## num. par=3
names(mod.pdot.psiB.3)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.pdot.psiB.3.psi <-mod.pdot.psiB.3$real$psi # occupancy probability
mod.pdot.psiB.3.psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7593708 0.1197860 0.4661479 0.9193896
## psi_unit2 0.7593708 0.1197860 0.4661479 0.9193896
## psi_unit3 0.7593708 0.1197860 0.4661479 0.9193896
## psi_unit4 0.4809979 0.1078882 0.2843274 0.6837386
## psi_unit5 0.7593708 0.1197860 0.4661479 0.9193896
# This is the logit occupancy the baseline and the log(odds ratio) of the two
# browse classes
mod.pdot.psiB.3$dmat$psi
## a1 a2
## psi "1" "psi.BrowCatn"
mod.pdot.psiB.3$beta$psi
## est se
## A1_psi 1.149233 0.655547
## A2_psi.psi.BrowCatn -1.225278 0.720332
mod.pdot.psiB.3.oddsratio.browse <- exp( sum(c(0,-1)*mod.pdot.psiB.3$beta$psi))
mod.pdot.psiB.3.oddsratio.browse
## [1] 1.656896
# Impact of browse on detectability as well?.
mod.pB.psiB <- RPresence::occMod(model=list(psi~BrowCat, p~BrowCat), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pB.psiB)
## Model name=psi(BrowCat)p(BrowCat)
## AIC=266.2606
## -2*log-likelihood=258.2606
## num. par=4
names(mod.pB.psiB)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.pB.psiB.psi <-mod.pB.psiB$real$psi # occupancy probability
mod.pB.psiB.psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7565126 0.1276939 0.4439733 0.9236045
## psi_unit2 0.7565126 0.1276939 0.4439733 0.9236045
## psi_unit3 0.7565126 0.1276939 0.4439733 0.9236045
## psi_unit4 0.4839101 0.1190770 0.2691584 0.7047733
## psi_unit5 0.7565126 0.1276939 0.4439733 0.9236045
mod.pB.psiB.p <-mod.pB.psiB$real$p # detection probability
mod.pB.psiB.p[1:5,]
## est se lower_0.95 upper_0.95
## p1_unit1 0.3518920 0.06891346 0.2309473 0.4953761
## p1_unit2 0.3518920 0.06891346 0.2309473 0.4953761
## p1_unit3 0.3518920 0.06891346 0.2309473 0.4953761
## p1_unit4 0.3451663 0.08600807 0.2000208 0.5263392
## p1_unit5 0.3518920 0.06891346 0.2309473 0.4953761
# This is the logit occupancy for each browse category
mod.pB.psiB$dmat$psi
## a1 a2
## psi "1" "psi.BrowCatn"
mod.pB.psiB$beta$psi
## est se
## A1_psi 1.133654 0.693230
## A2_psi.psi.BrowCatn -1.198036 0.841271
mod.pB.psiB.oddsratio.browse <- exp( sum(c(0,-1)*mod.pB.psiB$beta$psi))
mod.pB.psiB.oddsratio.browse
## [1] 1.4287
# Model where p(t) varies across survey occasions but psi(browse)
# Not sure why SURVEY key work no longer works for time varying factor
# but we have our own covariate (visit). Be sure to declare "visit" as a character or as a actor
#
mod.pt.psiB <- RPresence::occMod(model=list(psi~BrowCat, p~visit), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pt.psiB)
## Model name=psi(BrowCat)p(visit)
## AIC=259.4368
## -2*log-likelihood=245.4368
## num. par=7
mod.pt.psiB$real$psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7699149 0.1232454 0.4611357 0.9290003
## psi_unit2 0.7699149 0.1232454 0.4611357 0.9290003
## psi_unit3 0.7699149 0.1232454 0.4611357 0.9290003
## psi_unit4 0.4931709 0.1116897 0.2884113 0.7002475
## psi_unit5 0.7699149 0.1232454 0.4611357 0.9290003
mod.pt.psiB$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.3520562 0.09836946 0.18920071 0.5585270
## p2_unit1 0.3175294 0.08921332 0.17192709 0.5104317
## p3_unit1 0.1694830 0.06707655 0.07424149 0.3417958
## p4_unit1 0.3115815 0.08815339 0.16822898 0.5031897
## p5_unit1 0.5917252 0.10433600 0.38334739 0.7716353
print_one_site_estimates(mod.pt.psiB, site = 1)
## psi(BrowCat)p(visit)
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.7699149 0.12324545 0.4611357 0.9290003
## p_p1_unit1 0.3520562 0.09836946 0.1892007 0.5585270
fitted(mod.pt.psiB, param="psi")[1,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7699149 0.1232454 0.4611357 0.9290003
# Model where p varies by observer but constant over time psi(browse)
#
mod.pO.psiB <- RPresence::occMod(model=list(psi~BrowCat, p~Obs), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pt.psiB)
## Model name=psi(BrowCat)p(visit)
## AIC=259.4368
## -2*log-likelihood=245.4368
## num. par=7
mod.pO.psiB$real$psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7535936 0.1163803 0.4723969 0.9126370
## psi_unit2 0.7535936 0.1163803 0.4723969 0.9126370
## psi_unit3 0.7535936 0.1163803 0.4723969 0.9126370
## psi_unit4 0.4846713 0.1084314 0.2865457 0.6877352
## psi_unit5 0.7535936 0.1163803 0.4723969 0.9126370
mod.pO.psiB$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.2235210 0.06272749 0.1241580 0.3689096
## p2_unit1 0.4462478 0.08080860 0.2980129 0.6047011
## p3_unit1 0.3786094 0.07998710 0.2383368 0.5426237
## p4_unit1 0.4462478 0.08080860 0.2980129 0.6047011
## p5_unit1 NA NA NA NA
print_one_site_estimates(mod.pO.psiB, site = 1)
## psi(BrowCat)p(Obs)
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.7535936 0.11638027 0.4723969 0.9126370
## p_p1_unit1 0.2235210 0.06272749 0.1241580 0.3689096
fitted(mod.pO.psiB, param="psi")[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7535936 0.1163803 0.4723969 0.9126370
## psi_unit2 0.7535936 0.1163803 0.4723969 0.9126370
## psi_unit3 0.7535936 0.1163803 0.4723969 0.9126370
## psi_unit4 0.4846713 0.1084314 0.2865457 0.6877352
## psi_unit5 0.7535936 0.1163803 0.4723969 0.9126370
# Model where p varies by observer + visit but constant over time psi(browse)
# besure that obs and visit are declared as character
#
mod.pOpV.psiB <- RPresence::occMod(model=list(psi~BrowCat, p~Obs+visit), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
summary(mod.pOpV.psiB)
## Model name=psi(BrowCat)p(Obs P visit)
## AIC=257.5994
## -2*log-likelihood=239.5994
## num. par=9
mod.pOpV.psiB$real$psi[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7673332 0.1206126 0.4673488 0.9253540
## psi_unit2 0.7673332 0.1206126 0.4673488 0.9253540
## psi_unit3 0.7673332 0.1206126 0.4673488 0.9253540
## psi_unit4 0.5059982 0.1148981 0.2938164 0.7160411
## psi_unit5 0.7673332 0.1206126 0.4673488 0.9253540
mod.pOpV.psiB$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.2142123 0.08829610 0.08884418 0.4325123
## p2_unit1 0.4050194 0.11960241 0.20466934 0.6429461
## p3_unit1 0.1799850 0.08137467 0.06932630 0.3927389
## p4_unit1 0.4254702 0.12515648 0.21351570 0.6688854
## p5_unit1 NA NA NA NA
print_one_site_estimates(mod.pOpV.psiB, site = 1)
## psi(BrowCat)p(Obs P visit)
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.7673332 0.1206126 0.46734885 0.9253540
## p_p1_unit1 0.2142123 0.0882961 0.08884418 0.4325123
fitted(mod.pOpV.psiB, param="psi")[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7673332 0.1206126 0.4673488 0.9253540
## psi_unit2 0.7673332 0.1206126 0.4673488 0.9253540
## psi_unit3 0.7673332 0.1206126 0.4673488 0.9253540
## psi_unit4 0.5059982 0.1148981 0.2938164 0.7160411
## psi_unit5 0.7673332 0.1206126 0.4673488 0.9253540
# Other models
mod.pOpVpB.psiB <- RPresence::occMod(model=list(psi~BrowCat, p~Obs+visit+BrowCat), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
mod.pOpVpB.psi. <- RPresence::occMod(model=list(psi~1, p~Obs+visit+BrowCat), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
mod.pOpVpB.psi. <- RPresence::occMod(model=list(psi~1, p~Obs+visit), type="so", data=weta.pao)
## PRESENCE Version 2.12.21.
# Model averaging
models<-list(mod.pdot,
mod.pdot.psiB.1,
mod.pB.psiB,
mod.pt.psiB,
mod.pO.psiB,
mod.pOpV.psiB,
mod.pOpVpB.psiB,
mod.pOpVpB.psi.,
mod.pOpVpB.psi.
)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll warn.conv
## 1 psi(BrowCat)p(Obs P visit) 0.00 0.3138 9 239.60 0
## 2 psi()p(Obs P visit) 0.95 0.1954 8 242.55 0
## 3 psi()p(Obs P visit) 0.95 0.1954 8 242.55 0
## 4 psi(BrowCat)p(visit) 1.84 0.1252 7 245.44 0
## 5 psi(BrowCat)p(Obs P visit P BrowCat) 2.00 0.1155 10 239.60 0
## 6 psi(BrowCat)p(Obs) 4.44 0.0340 5 252.04 0
## 7 psi(-1 P Browsed P Unbrowsed)p() 6.66 0.0112 3 258.26 0
## 8 psi()p() 8.19 0.0052 2 261.79 0
## 9 psi(BrowCat)p(BrowCat) 8.66 0.0041 4 258.26 0
## warn.VC
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0
## 7 0
## 8 0
## 9 0
RPresence::modAvg(results, param="psi")[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7170024 0.1280240 0.4238522 0.8971786
## psi_unit2 0.7170024 0.1280240 0.4238522 0.8971786
## psi_unit3 0.7170024 0.1280240 0.4238522 0.8971786
## psi_unit4 0.5566424 0.1280841 0.3122397 0.7763918
## psi_unit5 0.7170024 0.1280240 0.4238522 0.8971786