#Analysis of butterfly data set
#3 seasons, 2 sampling occasions per season
#Interested in occupancy of CLLESSU
#Using RPresence
#2018-11-26 submitted by Neil Faught
library(car)
library(readxl)
library(RPresence)
## Warning: package 'RPresence' was built under R version 3.5.0
library(ggplot2)
library(reshape)
# Get the RMark additional functions
source(file.path("..","..","..","AdditionalFunctions","RPresence.additional.functions.R"))
#Read in data
bfly<-read.csv("../Butterfly.csv")
#Data is in long format, need to convert to wide format for RMark
#Want one row for each transect
head(bfly)
## Field Border Transect Treatment Date Visit AMLA AMSN BLSW BUCK
## 1 181 181E 1810 diskcontrol 6/3/2008 1 0 0 0 0
## 2 181 181E 18100 diskcontrol 6/3/2008 1 0 0 0 0
## 3 181 181E 181000 diskcontrol 6/3/2008 1 0 1 0 0
## 4 181 181N 181N1 diskcontrol 6/3/2008 1 0 0 0 0
## 5 181 181N 181N2 diskcontrol 6/3/2008 1 1 0 0 0
## 6 181 181N 181N3 diskcontrol 6/3/2008 1 0 0 0 0
## CHSK CLLESSU CLSK DESK DUSK ETBL ETSW EUSK FISK GISW GPHA GRHA GUFR HAEM
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 3 0 1 0 0 0 0 0 0 0 0 0 0 0 1
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## HODW LIYE MONA ORSU PECR PISW RBHA RSPP SACH SICH SLOR SOCL SOSK SPSW
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 1 1 0 0 0 0 0 0 0 0 0 0
## 3 0 0 1 0 1 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## SSSK SWSK TESK VAFR VICE
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 0 0 0
bfly<-subset(bfly,select=c(Field,Border,Transect,Treatment,Date,Visit,CLLESSU))
head(bfly)
## Field Border Transect Treatment Date Visit CLLESSU
## 1 181 181E 1810 diskcontrol 6/3/2008 1 0
## 2 181 181E 18100 diskcontrol 6/3/2008 1 0
## 3 181 181E 181000 diskcontrol 6/3/2008 1 1
## 4 181 181N 181N1 diskcontrol 6/3/2008 1 0
## 5 181 181N 181N2 diskcontrol 6/3/2008 1 0
## 6 181 181N 181N3 diskcontrol 6/3/2008 1 0
bfly$Transect=as.factor(bfly$Transect)
#Convert to wide format
junk<-melt(bfly,id.var=c("Transect","Visit"),measure.var="CLLESSU")
j2=cast(junk,Transect ~ Visit)
head(j2)
## Transect 1 2 3 4 5 6
## 1 1810 0 0 1 0 0 0
## 2 18100 0 0 0 1 0 0
## 3 181000 1 0 0 2 1 1
## 4 181N1 0 0 0 0 0 0
## 5 181N2 0 0 0 0 1 0
## 6 181N3 0 0 0 0 0 1
#Note how there are counts in these columns, need to convert to simply 0/1 for
#present/absent
j2$`1`=ifelse(j2$`1`>0,1,0)
j2$`2`=ifelse(j2$`2`>0,1,0)
j2$`3`=ifelse(j2$`3`>0,1,0)
j2$`4`=ifelse(j2$`4`>0,1,0)
j2$`5`=ifelse(j2$`5`>0,1,0)
j2$`6`=ifelse(j2$`6`>0,1,0)
head(j2)
## Transect 1 2 3 4 5 6
## 1 1810 0 0 1 0 0 0
## 2 18100 0 0 0 1 0 0
## 3 181000 1 0 0 1 1 1
## 4 181N1 0 0 0 0 0 0
## 5 181N2 0 0 0 0 1 0
## 6 181N3 0 0 0 0 0 1
input.history = j2[,2:7]
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 129
ncol(input.history)
## [1] 6
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 0
#Is the treatment a site or survey level covariate (i.e. does each transect only
#have a single treatment applied to it over the course of the 6 sampling occasions?)
jcov = melt(bfly,id.var=c("Transect","Visit"),measure.var="Treatment")
jcov2=cast(jcov,Transect ~ Visit)
head(jcov2)
## Transect 1 2 3 4 5
## 1 1810 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 2 18100 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 3 181000 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 4 181N1 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 5 181N2 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6 181N3 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6
## 1 diskcontrol
## 2 diskcontrol
## 3 diskcontrol
## 4 diskcontrol
## 5 diskcontrol
## 6 diskcontrol
jcov2$test = ifelse(jcov2$`1`==jcov2$`2` & jcov2$`1`==jcov2$`3` &
jcov2$`1`==jcov2$`4` & jcov2$`1`==jcov2$`5` &
jcov2$`1`==jcov2$`6`,0,1)
head(jcov2)
## Transect 1 2 3 4 5
## 1 1810 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 2 18100 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 3 181000 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 4 181N1 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 5 181N2 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6 181N3 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6 test
## 1 diskcontrol 0
## 2 diskcontrol 0
## 3 diskcontrol 0
## 4 diskcontrol 0
## 5 diskcontrol 0
## 6 diskcontrol 0
#Looks like this is a site level covariate
sum(jcov2$test)
## [1] 0
#Add site covariates to input history
site.covar = data.frame(Treatment = jcov2$`1`)
head(site.covar)
## Treatment
## 1810 diskcontrol
## 18100 diskcontrol
## 181000 diskcontrol
## 181N1 diskcontrol
## 181N2 diskcontrol
## 181N3 diskcontrol
# Number of visits in each season
Nvisits.per.season <- rep(2,3)
# Create the *.pao file
bfly.pao <- RPresence::createPao(input.history,
nsurveyseason=Nvisits.per.season,
unitcov=site.covar,
title='Nighingale SSMS')
bfly.pao
## $nunits
## [1] 129
##
## $nsurveys
## [1] 6
##
## $nseasons
## [1] 3
##
## $nmethods
## [1] 1
##
## $det.data
## 1 2 3 4 5 6
## 1 0 0 1 0 0 0
## 2 0 0 0 1 0 0
## 3 1 0 0 1 1 1
## 4 0 0 0 0 0 0
## 5 0 0 0 0 1 0
## 6 0 0 0 0 0 1
## 7 0 1 0 1 1 1
## 8 0 1 0 0 1 1
## 9 0 1 0 1 1 1
## 10 0 0 0 1 0 0
## 11 0 0 0 1 0 0
## 12 1 0 0 0 0 0
## 13 0 0 0 1 0 1
## 14 0 0 1 0 1 0
## 15 1 0 0 0 0 0
## 16 0 0 0 0 0 0
## 17 0 0 0 1 1 1
## 18 0 0 0 0 1 1
## 19 1 0 0 0 1 1
## 20 1 0 0 0 1 1
## 21 1 0 0 0 0 0
## 22 1 0 0 0 0 0
## 23 1 0 1 1 0 0
## 24 1 0 0 0 1 0
## 25 1 1 0 1 1 0
## 26 1 0 0 0 1 0
## 27 1 0 0 0 1 0
## 28 1 0 0 1 1 0
## 29 1 1 0 1 1 0
## 30 1 0 0 1 1 0
## 31 1 0 0 0 1 1
## 32 0 0 0 1 0 1
## 33 1 0 0 0 0 0
## 34 1 1 0 1 1 1
## 35 0 0 1 0 1 1
## 36 0 0 0 0 1 1
## 37 1 0 0 1 0 1
## 38 1 0 0 0 1 1
## 39 1 0 0 0 1 1
## 40 0 0 0 0 0 1
## 41 0 0 0 0 0 1
## 42 0 0 1 0 0 1
## 43 0 0 0 1 0 0
## 44 1 0 0 0 0 0
## 45 0 1 0 0 0 0
## 46 0 0 0 1 0 0
## 47 1 0 0 0 1 0
## 48 1 0 0 0 1 0
## 49 1 0 0 0 0 0
## 50 0 0 0 0 0 0
## 51 0 0 0 0 0 1
## 52 1 0 0 1 1 1
## 53 1 0 0 0 0 1
## 54 0 0 0 0 0 1
## 55 0 0 1 0 1 1
## 56 0 0 0 1 0 0
## 57 1 0 1 1 1 1
## 58 1 0 1 1 1 1
## 59 1 0 1 1 1 1
## 60 1 0 1 1 1 1
## 61 1 1 0 0 1 1
## 62 0 0 0 0 1 1
## 63 1 0 1 1 1 1
## 64 0 0 0 0 0 0
## 65 1 0 0 0 1 0
## 66 0 0 0 0 0 0
## 67 1 1 0 1 0 1
## 68 1 1 0 1 1 1
## 69 1 0 1 1 1 1
## 70 0 0 0 0 0 0
## 71 1 0 0 0 0 0
## 72 0 0 0 0 0 0
## 73 1 0 0 1 1 1
## 74 1 0 0 0 1 1
## 75 0 0 0 1 0 1
## 76 1 0 0 0 1 1
## 77 0 0 0 0 0 1
## 78 1 0 0 0 0 1
## 79 0 0 0 0 1 1
## 80 0 1 1 1 1 1
## 81 1 0 0 0 1 1
## 82 1 0 0 0 0 1
## 83 1 0 0 0 0 0
## 84 0 0 0 0 0 1
## 85 0 0 1 1 1 1
## 86 0 0 0 1 1 1
## 87 0 0 0 1 1 1
## 88 0 0 0 0 0 1
## 89 0 1 0 0 1 0
## 90 0 1 1 0 0 1
## 91 1 0 0 1 1 1
## 92 0 0 0 0 0 1
## 93 0 0 0 0 0 1
## 94 0 0 0 0 0 0
## 95 0 0 0 0 0 0
## 96 1 0 0 0 0 0
## 97 0 0 0 0 0 0
## 98 0 0 0 0 0 1
## 99 0 0 0 0 0 0
## 100 0 1 0 1 0 0
## 101 1 1 0 0 0 0
## 102 1 0 0 0 0 0
## 103 0 0 1 0 0 1
## 104 1 0 0 1 0 1
## 105 1 0 1 0 1 1
## 106 1 1 1 1 1 1
## 107 1 1 1 1 1 1
## 108 1 1 1 0 1 1
## 109 1 0 1 0 1 1
## 110 1 0 0 1 1 0
## 111 1 0 0 1 0 1
## 112 1 1 1 0 0 1
## 113 1 1 1 0 0 1
## 114 1 1 1 0 0 1
## 115 1 0 1 1 1 1
## 116 1 1 1 1 1 1
## 117 1 0 1 1 1 1
## 118 1 0 0 0 0 1
## 119 1 0 1 0 0 1
## 120 1 0 0 1 0 0
## 121 0 0 0 0 0 1
## 122 0 1 0 0 0 0
## 123 0 0 0 0 0 0
## 124 0 0 0 0 0 0
## 125 0 0 0 0 0 0
## 126 0 0 0 0 0 0
## 127 0 1 1 1 1 1
## 128 0 1 0 1 1 1
## 129 0 0 1 1 0 1
##
## $nunitcov
## [1] 1
##
## $unitcov
## Treatment
## 1810 diskcontrol
## 18100 diskcontrol
## 181000 diskcontrol
## 181N1 diskcontrol
## 181N2 diskcontrol
## 181N3 diskcontrol
## 181SW1 disk
## 181SW2 disk
## 181SW3 disk
## 183NE1 burn
## 183NE2 burn
## 183NE3 burn
## 1850 fieldcontrol
## 18500 fieldcontrol
## 185000 fieldcontrol
## 185W1 fieldcontrol
## 185W2 fieldcontrol
## 185W3 fieldcontrol
## 1870 diskcontrol
## 18700 diskcontrol
## 187000 diskcontrol
## 187N1 disk
## 187N2 disk
## 187N3 disk
## 187S1 diskcontrol
## 187S2 diskcontrol
## 187S3 diskcontrol
## 187W1 diskcontrol
## 187W2 diskcontrol
## 187W3 diskcontrol
## 189N1 fieldcontrol
## 189N2 fieldcontrol
## 189N3 fieldcontrol
## 189SE1 fieldcontrol
## 189SE2 fieldcontrol
## 189SE3 fieldcontrol
## 189W1 fieldcontrol
## 189W2 fieldcontrol
## 189W3 fieldcontrol
## 191S1 burncontrol
## 191S2 burncontrol
## 191S3 burncontrol
## 193N1 fieldcontrol
## 193N2 fieldcontrol
## 193N3 fieldcontrol
## 193SE1 fieldcontrol
## 193SE2 fieldcontrol
## 193SE3 fieldcontrol
## 193W1 fieldcontrol
## 193W2 fieldcontrol
## 193W3 fieldcontrol
## 197ESE1 diskcontrol
## 197ESE2 diskcontrol
## 197ESE3 diskcontrol
## 197MDNW1 burncontrol
## 197MDNW2 burncontrol
## 197MDNW3 burncontrol
## 197MDSE1 diskcontrol
## 197MDSE2 diskcontrol
## 197MDSE3 diskcontrol
## 197NE1 disk
## 197NE2 disk
## 197NE3 disk
## 197NW1 burn
## 197NW2 burn
## 197NW3 burn
## 197SPUR1 diskcontrol
## 197SPUR2 diskcontrol
## 197SPUR3 diskcontrol
## 197W1 burncontrol
## 197W2 burncontrol
## 197W3 burncontrol
## 197WSE1 burncontrol
## 197WSE2 burncontrol
## 197WSE3 burncontrol
## 198NE1 burncontrol
## 198NE2 burncontrol
## 198NE3 burncontrol
## 198NW1 burn
## 198NW2 burn
## 198NW3 burn
## 198S1 burncontrol
## 198S2 burncontrol
## 198S3 burncontrol
## 199NW1 fieldcontrol
## 199NW2 fieldcontrol
## 199NW3 fieldcontrol
## 199SE1 fieldcontrol
## 199SE2 fieldcontrol
## 199SE3 fieldcontrol
## 199SW1 fieldcontrol
## 199SW2 fieldcontrol
## 199SW3 fieldcontrol
## 201NE1 diskcontrol
## 201NE2 diskcontrol
## 201NE3 diskcontrol
## 201NW1 diskcontrol
## 201NW2 diskcontrol
## 201NW3 diskcontrol
## 201S1 disk
## 201S2 disk
## 201S3 disk
## 2210 burncontrol
## 22100 burncontrol
## 221000 burncontrol
## 221N1 burn
## 221N2 burn
## 221N3 burn
## 221S1 burncontrol
## 221S2 burncontrol
## 221S3 burncontrol
## 223N1 diskcontrol
## 223N2 diskcontrol
## 223N3 diskcontrol
## 223SE1 disk
## 223SE2 disk
## 223SE3 disk
## 223W1 diskcontrol
## 223W2 diskcontrol
## 223W3 diskcontrol
## 2400 burncontrol
## 24000 burncontrol
## 240000 burncontrol
## 240S1 burn
## 240S2 burn
## 240S3 burn
## 240W1 burncontrol
## 240W2 burncontrol
## 240W3 burncontrol
##
## $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 1
## 41 1
## 42 1
## 43 1
## 44 1
## 45 1
## 46 1
## 47 1
## 48 1
## 49 1
## 50 1
## 51 1
## 52 1
## 53 1
## 54 1
## 55 1
## 56 1
## 57 1
## 58 1
## 59 1
## 60 1
## 61 1
## 62 1
## 63 1
## 64 1
## 65 1
## 66 1
## 67 1
## 68 1
## 69 1
## 70 1
## 71 1
## 72 1
## 73 1
## 74 1
## 75 1
## 76 1
## 77 1
## 78 1
## 79 1
## 80 1
## 81 1
## 82 1
## 83 1
## 84 1
## 85 1
## 86 1
## 87 1
## 88 1
## 89 1
## 90 1
## 91 1
## 92 1
## 93 1
## 94 1
## 95 1
## 96 1
## 97 1
## 98 1
## 99 1
## 100 1
## 101 1
## 102 1
## 103 1
## 104 1
## 105 1
## 106 1
## 107 1
## 108 1
## 109 1
## 110 1
## 111 1
## 112 1
## 113 1
## 114 1
## 115 1
## 116 1
## 117 1
## 118 1
## 119 1
## 120 1
## 121 1
## 122 1
## 123 1
## 124 1
## 125 1
## 126 1
## 127 1
## 128 1
## 129 1
## 130 2
## 131 2
## 132 2
## 133 2
## 134 2
## 135 2
## 136 2
## 137 2
## 138 2
## 139 2
## 140 2
## 141 2
## 142 2
## 143 2
## 144 2
## 145 2
## 146 2
## 147 2
## 148 2
## 149 2
## 150 2
## 151 2
## 152 2
## 153 2
## 154 2
## 155 2
## 156 2
## 157 2
## 158 2
## 159 2
## 160 2
## 161 2
## 162 2
## 163 2
## 164 2
## 165 2
## 166 2
## 167 2
## 168 2
## 169 2
## 170 2
## 171 2
## 172 2
## 173 2
## 174 2
## 175 2
## 176 2
## 177 2
## 178 2
## 179 2
## 180 2
## 181 2
## 182 2
## 183 2
## 184 2
## 185 2
## 186 2
## 187 2
## 188 2
## 189 2
## 190 2
## 191 2
## 192 2
## 193 2
## 194 2
## 195 2
## 196 2
## 197 2
## 198 2
## 199 2
## 200 2
## 201 2
## 202 2
## 203 2
## 204 2
## 205 2
## 206 2
## 207 2
## 208 2
## 209 2
## 210 2
## 211 2
## 212 2
## 213 2
## 214 2
## 215 2
## 216 2
## 217 2
## 218 2
## 219 2
## 220 2
## 221 2
## 222 2
## 223 2
## 224 2
## 225 2
## 226 2
## 227 2
## 228 2
## 229 2
## 230 2
## 231 2
## 232 2
## 233 2
## 234 2
## 235 2
## 236 2
## 237 2
## 238 2
## 239 2
## 240 2
## 241 2
## 242 2
## 243 2
## 244 2
## 245 2
## 246 2
## 247 2
## 248 2
## 249 2
## 250 2
## 251 2
## 252 2
## 253 2
## 254 2
## 255 2
## 256 2
## 257 2
## 258 2
## 259 3
## 260 3
## 261 3
## 262 3
## 263 3
## 264 3
## 265 3
## 266 3
## 267 3
## 268 3
## 269 3
## 270 3
## 271 3
## 272 3
## 273 3
## 274 3
## 275 3
## 276 3
## 277 3
## 278 3
## 279 3
## 280 3
## 281 3
## 282 3
## 283 3
## 284 3
## 285 3
## 286 3
## 287 3
## 288 3
## 289 3
## 290 3
## 291 3
## 292 3
## 293 3
## 294 3
## 295 3
## 296 3
## 297 3
## 298 3
## 299 3
## 300 3
## 301 3
## 302 3
## 303 3
## 304 3
## 305 3
## 306 3
## 307 3
## 308 3
## 309 3
## 310 3
## 311 3
## 312 3
## 313 3
## 314 3
## 315 3
## 316 3
## 317 3
## 318 3
## 319 3
## 320 3
## 321 3
## 322 3
## 323 3
## 324 3
## 325 3
## 326 3
## 327 3
## 328 3
## 329 3
## 330 3
## 331 3
## 332 3
## 333 3
## 334 3
## 335 3
## 336 3
## 337 3
## 338 3
## 339 3
## 340 3
## 341 3
## 342 3
## 343 3
## 344 3
## 345 3
## 346 3
## 347 3
## 348 3
## 349 3
## 350 3
## 351 3
## 352 3
## 353 3
## 354 3
## 355 3
## 356 3
## 357 3
## 358 3
## 359 3
## 360 3
## 361 3
## 362 3
## 363 3
## 364 3
## 365 3
## 366 3
## 367 3
## 368 3
## 369 3
## 370 3
## 371 3
## 372 3
## 373 3
## 374 3
## 375 3
## 376 3
## 377 3
## 378 3
## 379 3
## 380 3
## 381 3
## 382 3
## 383 3
## 384 3
## 385 3
## 386 3
## 387 3
## 388 4
## 389 4
## 390 4
## 391 4
## 392 4
## 393 4
## 394 4
## 395 4
## 396 4
## 397 4
## 398 4
## 399 4
## 400 4
## 401 4
## 402 4
## 403 4
## 404 4
## 405 4
## 406 4
## 407 4
## 408 4
## 409 4
## 410 4
## 411 4
## 412 4
## 413 4
## 414 4
## 415 4
## 416 4
## 417 4
## 418 4
## 419 4
## 420 4
## 421 4
## 422 4
## 423 4
## 424 4
## 425 4
## 426 4
## 427 4
## 428 4
## 429 4
## 430 4
## 431 4
## 432 4
## 433 4
## 434 4
## 435 4
## 436 4
## 437 4
## 438 4
## 439 4
## 440 4
## 441 4
## 442 4
## 443 4
## 444 4
## 445 4
## 446 4
## 447 4
## 448 4
## 449 4
## 450 4
## 451 4
## 452 4
## 453 4
## 454 4
## 455 4
## 456 4
## 457 4
## 458 4
## 459 4
## 460 4
## 461 4
## 462 4
## 463 4
## 464 4
## 465 4
## 466 4
## 467 4
## 468 4
## 469 4
## 470 4
## 471 4
## 472 4
## 473 4
## 474 4
## 475 4
## 476 4
## 477 4
## 478 4
## 479 4
## 480 4
## 481 4
## 482 4
## 483 4
## 484 4
## 485 4
## 486 4
## 487 4
## 488 4
## 489 4
## 490 4
## 491 4
## 492 4
## 493 4
## 494 4
## 495 4
## 496 4
## 497 4
## 498 4
## 499 4
## 500 4
## 501 4
## 502 4
## 503 4
## 504 4
## 505 4
## 506 4
## 507 4
## 508 4
## 509 4
## 510 4
## 511 4
## 512 4
## 513 4
## 514 4
## 515 4
## 516 4
## 517 5
## 518 5
## 519 5
## 520 5
## 521 5
## 522 5
## 523 5
## 524 5
## 525 5
## 526 5
## 527 5
## 528 5
## 529 5
## 530 5
## 531 5
## 532 5
## 533 5
## 534 5
## 535 5
## 536 5
## 537 5
## 538 5
## 539 5
## 540 5
## 541 5
## 542 5
## 543 5
## 544 5
## 545 5
## 546 5
## 547 5
## 548 5
## 549 5
## 550 5
## 551 5
## 552 5
## 553 5
## 554 5
## 555 5
## 556 5
## 557 5
## 558 5
## 559 5
## 560 5
## 561 5
## 562 5
## 563 5
## 564 5
## 565 5
## 566 5
## 567 5
## 568 5
## 569 5
## 570 5
## 571 5
## 572 5
## 573 5
## 574 5
## 575 5
## 576 5
## 577 5
## 578 5
## 579 5
## 580 5
## 581 5
## 582 5
## 583 5
## 584 5
## 585 5
## 586 5
## 587 5
## 588 5
## 589 5
## 590 5
## 591 5
## 592 5
## 593 5
## 594 5
## 595 5
## 596 5
## 597 5
## 598 5
## 599 5
## 600 5
## 601 5
## 602 5
## 603 5
## 604 5
## 605 5
## 606 5
## 607 5
## 608 5
## 609 5
## 610 5
## 611 5
## 612 5
## 613 5
## 614 5
## 615 5
## 616 5
## 617 5
## 618 5
## 619 5
## 620 5
## 621 5
## 622 5
## 623 5
## 624 5
## 625 5
## 626 5
## 627 5
## 628 5
## 629 5
## 630 5
## 631 5
## 632 5
## 633 5
## 634 5
## 635 5
## 636 5
## 637 5
## 638 5
## 639 5
## 640 5
## 641 5
## 642 5
## 643 5
## 644 5
## 645 5
## 646 6
## 647 6
## 648 6
## 649 6
## 650 6
## 651 6
## 652 6
## 653 6
## 654 6
## 655 6
## 656 6
## 657 6
## 658 6
## 659 6
## 660 6
## 661 6
## 662 6
## 663 6
## 664 6
## 665 6
## 666 6
## 667 6
## 668 6
## 669 6
## 670 6
## 671 6
## 672 6
## 673 6
## 674 6
## 675 6
## 676 6
## 677 6
## 678 6
## 679 6
## 680 6
## 681 6
## 682 6
## 683 6
## 684 6
## 685 6
## 686 6
## 687 6
## 688 6
## 689 6
## 690 6
## 691 6
## 692 6
## 693 6
## 694 6
## 695 6
## 696 6
## 697 6
## 698 6
## 699 6
## 700 6
## 701 6
## 702 6
## 703 6
## 704 6
## 705 6
## 706 6
## 707 6
## 708 6
## 709 6
## 710 6
## 711 6
## 712 6
## 713 6
## 714 6
## 715 6
## 716 6
## 717 6
## 718 6
## 719 6
## 720 6
## 721 6
## 722 6
## 723 6
## 724 6
## 725 6
## 726 6
## 727 6
## 728 6
## 729 6
## 730 6
## 731 6
## 732 6
## 733 6
## 734 6
## 735 6
## 736 6
## 737 6
## 738 6
## 739 6
## 740 6
## 741 6
## 742 6
## 743 6
## 744 6
## 745 6
## 746 6
## 747 6
## 748 6
## 749 6
## 750 6
## 751 6
## 752 6
## 753 6
## 754 6
## 755 6
## 756 6
## 757 6
## 758 6
## 759 6
## 760 6
## 761 6
## 762 6
## 763 6
## 764 6
## 765 6
## 766 6
## 767 6
## 768 6
## 769 6
## 770 6
## 771 6
## 772 6
## 773 6
## 774 6
##
## $nsurveyseason
## [1] 2 2 2
##
## $title
## [1] "Nighingale 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" "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"
## [78] "unit78" "unit79" "unit80" "unit81" "unit82" "unit83" "unit84"
## [85] "unit85" "unit86" "unit87" "unit88" "unit89" "unit90" "unit91"
## [92] "unit92" "unit93" "unit94" "unit95" "unit96" "unit97" "unit98"
## [99] "unit99" "unit100" "unit101" "unit102" "unit103" "unit104" "unit105"
## [106] "unit106" "unit107" "unit108" "unit109" "unit110" "unit111" "unit112"
## [113] "unit113" "unit114" "unit115" "unit116" "unit117" "unit118" "unit119"
## [120] "unit120" "unit121" "unit122" "unit123" "unit124" "unit125" "unit126"
## [127] "unit127" "unit128" "unit129"
##
## $surveynames
## [1] "1-1" "1-2" "2-1" "2-2" "3-1" "3-2"
##
## $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 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
## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## attr(,"class")
## [1] "pao"
summary(bfly.pao)
## paoname=pres.pao
## title=Nighingale SSMS
## Naive occ=0.8837209
## nunits nsurveys nseasons nsurveyseason nmethods
## "129" "6" "3" "2,2,2" "1"
## nunitcov nsurvcov
## "1" "1"
## unit covariates : Treatment
## survey covariates: SURVEY
#Fit Null Model
mod.fit <- RPresence::occMod(
model=list(psi~1, gamma~1, epsilon~1, p~1),
type="do.1", data=bfly.pao)
## PRESENCE Version 2.12.21.
summary(mod.fit)
## Model name=psi()gamma()epsilon()p()
## AIC=1024.4158
## -2*log-likelihood=1016.4158
## num. par=4
mod.fit <- RPresence.add.derived(mod.fit)
## Loading required package: plyr
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:reshape':
##
## rename, round_any
# Estimate of initial occupance
mod.fit$real$psi[1,]
## est se lower_0.95 upper_0.95
## unit1_1 0.771483 0.06066435 0.6323413 0.8688848
# 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.8012322 0.04891314 0.6882716 0.8803739
## unit1_3 0.8152842 0.05846394 0.6734490 0.9042712
# 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.7714830 0.06066435 0.6323413 0.8688848
## unit1_2 0.8012322 0.04891314 0.6882716 0.8803739
## unit1_3 0.8152842 0.05846394 0.6734490 0.9042712
# 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.09082806 0.05287451 0.0276922 0.259492
## epsilon2_unit1 0.09082806 0.05287451 0.0276922 0.259492
# 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.436823 0.1099428 0.2441637 0.65064
## gamma2_unit1 0.436823 0.1099428 0.2441637 0.65064
# 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.4970695 0.03295812 0.4328782 0.5613576
## p2_unit1 0.4970695 0.03295812 0.4328782 0.5613576
## p3_unit1 0.4970695 0.03295812 0.4328782 0.5613576
## p4_unit1 0.4970695 0.03295812 0.4328782 0.5613576
## p5_unit1 0.4970695 0.03295812 0.4328782 0.5613576
## p6_unit1 0.4970695 0.03295812 0.4328782 0.5613576
# 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.038561 NA NA NA
## unit1_2 1.017538 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.038561 NA NA NA
## unit1_2 1.017538 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.194000 NA NA NA
## unit1_2 1.094946 NA NA NA
# Plot each season's predicted occupancy
plotdata = mod.fit$derived$all_psi[ grepl('unit1_', row.names(mod.fit$derived$all_psi)),]
plotdata$Season = c(1:3)
head(plotdata)
## est se lower_0.95 upper_0.95 Season
## unit1_1 0.7714830 0.06066435 0.6323413 0.8688848 1
## unit1_2 0.8012322 0.04891314 0.6882716 0.8803739 2
## unit1_3 0.8152842 0.05846394 0.6734490 0.9042712 3
ggplot(data=plotdata, aes(x=Season,y=est))+
ggtitle("Estimated occupancy 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:3)
