# Single Species Single Season Occupancy
# Weta Example with model averaging with a list of models to fit
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library(plyr)
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.
# 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"
# 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, ~BrowCat
~BrowCat, ~BrowCat
~visit, ~BrowCat
~Obs, ~BrowCat
~Obs+visit, ~BrowCat
~Obs+visit+BrowCat, ~BrowCat
~Obs+visit+BrowCat, ~1
~Obs+visit, ~1")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p psi
## 1 ~1 ~1
## 2 ~1 ~BrowCat
## 3 ~BrowCat ~BrowCat
## 4 ~visit ~BrowCat
## 5 ~Obs ~BrowCat
## 6 ~Obs+visit ~BrowCat
## 7 ~Obs+visit+BrowCat ~BrowCat
## 8 ~Obs+visit+BrowCat ~1
## 9 ~Obs+visit ~1
# 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=weta.pao)
##
##
## ***** Starting ~1 ~1
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~1 ~BrowCat
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~BrowCat ~BrowCat
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~visit ~BrowCat
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~Obs ~BrowCat
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~Obs+visit ~BrowCat
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~Obs+visit+BrowCat ~BrowCat
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~Obs+visit+BrowCat ~1
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~Obs+visit ~1
## PRESENCE Version 2.12.21.
# 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.475124 0.374534
##
## $psi.VC
## [,1]
## [1,] 0.140276
##
## $p
## est se
## B1_p1 -0.621831 0.236368
##
## $p.VC
## [,1]
## [1,] 0.05587
##
## $VC
## A1_psi B1_p1
## A1_psi 0.140276 -0.047963
## B1_p1 -0.047963 0.055870
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.6165958 0.08854195 0.4356219 0.770157
## psi_unit2 0.6165958 0.08854195 0.4356219 0.770157
## psi_unit3 0.6165958 0.08854195 0.4356219 0.770157
## psi_unit4 0.6165958 0.08854195 0.4356219 0.770157
## psi_unit5 0.6165958 0.08854195 0.4356219 0.770157
model.fits[[check.model]]$real$p[1:5,]
## est se lower_0.95 upper_0.95
## p1_unit1 0.3493651 0.05372869 0.2525413 0.4604435
## p1_unit2 0.3493651 0.05372869 0.2525413 0.4604435
## p1_unit3 0.3493651 0.05372869 0.2525413 0.4604435
## p1_unit4 0.3493651 0.05372869 0.2525413 0.4604435
## p1_unit5 0.3493651 0.05372869 0.2525413 0.4604435
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.2237227 0.107541 0.07887311 0.4923886
## unit2 0.2237227 0.107541 0.07887311 0.4923886
## unit3 1.0000000 0.000000 1.00000000 1.0000000
## unit4 0.2237227 0.107541 0.07887311 0.4923886
## unit5 0.2237227 0.107541 0.07887311 0.4923886
## unit6 0.2237227 0.107541 0.07887311 0.4923886
## unit7 0.2237227 0.107541 0.07887311 0.4923886
## unit8 0.2237227 0.107541 0.07887311 0.4923886
## unit9 0.2237227 0.107541 0.07887311 0.4923886
## unit10 1.0000000 0.000000 1.00000000 1.0000000
tail(model.fits[[check.model]]$derived$psi_c)
## est se lower_0.95 upper_0.95
## unit67 1.0000000 0.0000000 1.0000000 1.0000000
## unit68 1.0000000 0.0000000 1.0000000 1.0000000
## unit69 1.0000000 0.0000000 1.0000000 1.0000000
## unit70 0.3069758 0.1169478 0.1310483 0.5654055
## unit71 0.3069758 0.1169478 0.1310483 0.5654055
## unit72 1.0000000 0.0000000 1.0000000 1.0000000
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model AIC neg2ll npar warn.conv
## 6 psi(BrowCat)p(Obs P visit) 257.5994 239.5994 9 0
## 9 psi()p(Obs P visit) 258.5466 242.5466 8 0
## 8 psi()p(Obs P visit P BrowCat) 259.4147 241.4147 9 0
## 4 psi(BrowCat)p(visit) 259.4368 245.4368 7 0
## 7 psi(BrowCat)p(Obs P visit P BrowCat) 259.5993 239.5993 10 0
## 5 psi(BrowCat)p(Obs) 262.0421 252.0421 5 0
## 2 psi(BrowCat)p() 264.2643 258.2643 3 0
## 1 psi()p() 265.7872 261.7872 2 0
## 3 psi(BrowCat)p(BrowCat) 266.2606 258.2606 4 0
## warn.VC DAIC modlike wgt
## 6 0 0.0000 1.0000 0.3370
## 9 0 0.9472 0.6228 0.2099
## 8 0 1.8153 0.4035 0.1360
## 4 0 1.8374 0.3990 0.1345
## 7 0 1.9999 0.3679 0.1240
## 5 0 4.4427 0.1085 0.0366
## 2 0 6.6649 0.0357 0.0120
## 1 0 8.1878 0.0167 0.0056
## 3 0 8.6612 0.0132 0.0044
names(aic.table)
## [1] "table" "models" "ess"
RPresence::modAvg(aic.table, param="psi")[1:5,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.7274042 0.1270832 0.4317642 0.9035791
## psi_unit2 0.7274042 0.1270832 0.4317642 0.9035791
## psi_unit3 0.7274042 0.1270832 0.4317642 0.9035791
## psi_unit4 0.5551929 0.1330731 0.3027011 0.7820781
## psi_unit5 0.7274042 0.1270832 0.4317642 0.9035791
ma.p <- RPresence::modAvg(aic.table, param="p")
ma.p[grepl("unit1$", rownames(ma.p)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.2369357 0.10256635 0.09267160 0.4855878
## p2_unit1 0.3962845 0.11921908 0.19818838 0.6354603
## p3_unit1 0.1913661 0.09298252 0.06793497 0.4345124
## p4_unit1 0.4134718 0.12696411 0.20169365 0.6629525
## p5_unit1 NA NA NA NA