nso-anal.R
cschwarz
Fri Aug 17 16:52:11 2018
# single species multi season
# Northern Spotted Owl (Strix occidentalis caurina) in California.
# s=55 sites visited up to K=4 times per season between 1997 and 2001 Y=5).
# Detection probabilities relatively constant within years, but likely different among years.
# Using the unmarked package
# 2018-08-17 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.ca@gmail.com)# unmarked package
library(readxl)
library(unmarked)## Loading required package: reshape## Loading required package: lattice## Loading required package: parallel## Loading required package: Rcpplibrary(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.
# Each year must have the same number of visit. Pad with missing if this is not true.
input.history <- read.csv(file.path("..","NSO.csv"), header=FALSE, skip=2, na.strings="-")
input.history$V1 <- NULL # drop the site number
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)## [1] 55ncol(input.history)## [1] 40range(input.history, na.rm=TRUE)## [1] 0 1sum(is.na(input.history))## [1] 752# Five years with 8 visits. input.history must have same number of visits/year
Nsites = nrow(input.history)
Nyears = 5
Nvisits= ncol(input.history)/Nyears
# Create the site covariates
site.covar <- data.frame(SiteNum=1:Nsites)
head(site.covar)## SiteNum
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6# Create the Yearly site (site x time) covariates.
# yearlySiteCovs is a data frame with Nsites x Nyears rows which are in site-major, year-minor order.
# I.e. list covariates for Site 1 for years 1..., then Site 2 for years .... etc
yearsite.covar <- data.frame(SiteNum=rep(1:Nsites, each=Nyears),
Year =as.factor(rep(1:Nyears,Nsites)))
head(yearsite.covar)## SiteNum Year
## 1 1 1
## 2 1 2
## 3 1 3
## 4 1 4
## 5 1 5
## 6 2 1# Create the Observation covarates (site x time x visit)
# Observation covaraites is a data frame with Nsite x Ntime x Nvisit rows
# which are ordered by site-year-observation, so that a column of obsCovs corresponds to as.vector(t(y))
test.history <- matrix(1:40, nrow=4, ncol=10, byrow=TRUE)
test.history## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 2 3 4 5 6 7 8 9 10
## [2,] 11 12 13 14 15 16 17 18 19 20
## [3,] 21 22 23 24 25 26 27 28 29 30
## [4,] 31 32 33 34 35 36 37 38 39 40as.vector(t(test.history))## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40obs.covar <- data.frame( SiteNum = rep(1:Nsites, each=Nyears * Nvisits),
Year = as.factor(rep( rep(1:Nyears, each=Nvisits), Nsites)),
Visit = as.factor(rep(1:Nvisits, Nyears*Nsites)))
obs.covar[1:30,]## SiteNum Year Visit
## 1 1 1 1
## 2 1 1 2
## 3 1 1 3
## 4 1 1 4
## 5 1 1 5
## 6 1 1 6
## 7 1 1 7
## 8 1 1 8
## 9 1 2 1
## 10 1 2 2
## 11 1 2 3
## 12 1 2 4
## 13 1 2 5
## 14 1 2 6
## 15 1 2 7
## 16 1 2 8
## 17 1 3 1
## 18 1 3 2
## 19 1 3 3
## 20 1 3 4
## 21 1 3 5
## 22 1 3 6
## 23 1 3 7
## 24 1 3 8
## 25 1 4 1
## 26 1 4 2
## 27 1 4 3
## 28 1 4 4
## 29 1 4 5
## 30 1 4 6# Create the UMF file
nso.UMF <- unmarked::unmarkedMultFrame(input.history,
siteCovs=site.covar,
yearlySiteCovs=yearsite.covar,
obsCovs = obs.covar,
numPrimary=Nyears)
nso.UMF## Data frame representation of unmarkedFrame object.
## y.1 y.2 y.3 y.4 y.5 y.6 y.7 y.8 y.9 y.10 y.11 y.12 y.13 y.14 y.15 y.16
## 1 0 1 1 1 NA NA NA NA 0 1 0 NA NA NA NA NA
## 2 0 0 NA NA NA NA NA NA 0 0 0 NA NA NA NA NA
## 3 1 0 0 0 1 1 1 NA 0 1 0 0 NA NA NA NA
## 4 1 1 1 NA NA NA NA NA 1 0 NA NA NA NA NA NA
## 5 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 0 NA
## 6 1 NA NA NA NA NA NA NA 0 0 0 1 1 1 NA NA
## 7 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 0 NA
## 8 0 0 0 0 0 0 NA NA 0 0 1 0 0 0 0 0
## 9 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 NA NA
## 10 1 1 1 NA NA NA NA NA 0 0 0 0 0 0 0 NA
## 11 1 1 1 NA NA NA NA NA 1 0 0 1 0 1 1 1
## 12 0 0 0 0 1 NA NA NA 1 1 0 0 0 0 1 NA
## 13 0 1 1 NA NA NA NA NA 0 0 1 1 1 NA NA NA
## 14 1 1 NA NA NA NA NA NA 0 1 0 0 0 NA NA NA
## 15 1 0 0 0 0 1 0 1 0 1 1 1 1 1 1 1
## 16 1 0 0 0 NA NA NA NA 0 0 1 1 1 NA NA NA
## 17 0 1 1 NA NA NA NA NA 0 1 1 1 1 1 NA NA
## 18 0 0 0 0 0 NA NA NA 0 0 0 0 NA NA NA NA
## 19 0 0 0 0 1 NA NA NA 0 0 0 0 0 0 NA NA
## 20 0 0 0 0 1 1 NA NA 0 1 1 1 1 NA NA NA
## 21 0 0 0 1 0 1 NA NA 1 1 1 NA NA NA NA NA
## 22 1 1 NA NA NA NA NA NA 0 1 0 1 0 1 NA NA
## 23 1 1 1 1 1 0 NA NA 0 1 0 1 NA NA NA NA
## 24 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 0 0
## 25 0 0 0 0 0 NA NA NA 0 0 0 NA NA NA NA NA
## 26 0 1 NA NA NA NA NA NA 0 1 1 NA NA NA NA NA
## 27 0 1 0 0 0 0 0 NA 1 1 0 0 0 0 NA NA
## 28 0 0 0 0 0 NA NA NA 0 0 0 0 0 0 0 NA
## 29 1 1 1 1 1 1 1 1 1 1 NA NA NA NA NA NA
## 30 0 0 1 1 1 1 1 NA 1 0 NA NA NA NA NA NA
## 31 1 1 0 NA NA NA NA NA 0 0 0 0 0 NA NA NA
## 32 0 0 0 NA NA NA NA NA 0 0 0 0 0 NA NA NA
## 33 0 1 0 0 0 0 NA NA 0 0 1 NA NA NA NA NA
## 34 0 1 1 1 1 1 NA NA 1 1 1 NA NA NA NA NA
## 35 0 0 0 1 1 1 NA NA 1 1 1 1 0 NA NA NA
## 36 0 0 0 0 0 NA NA NA 0 0 0 0 0 0 0 NA
## 37 0 1 1 NA NA NA NA NA 0 0 1 NA NA NA NA NA
## 38 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 NA NA
## 39 0 1 1 1 1 NA NA NA 0 0 0 0 0 0 1 NA
## 40 0 1 0 1 1 1 0 1 0 0 0 NA NA NA NA NA
## 41 0 0 1 1 1 NA NA NA 0 0 0 1 0 NA NA NA
## 42 0 0 0 0 NA NA NA NA 0 0 0 0 0 NA NA NA
## 43 0 0 1 1 0 1 1 NA 0 0 0 1 1 1 NA NA
## 44 0 0 0 0 0 0 0 0 0 1 1 1 0 NA NA NA
## 45 1 1 1 1 1 1 NA NA 0 0 1 NA NA NA NA NA
## 46 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 0 NA
## 47 0 0 1 0 1 1 NA NA 1 0 1 0 0 0 0 NA
## 48 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 NA NA
## 49 0 NA NA NA NA NA NA NA 0 0 NA NA NA NA NA NA
## 50 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 51 0 0 0 0 0 NA NA NA 0 0 0 0 0 0 NA NA
## 52 0 NA NA NA NA NA NA NA 1 1 1 1 1 1 1 NA
## 53 0 0 0 0 0 0 NA NA 0 0 0 0 0 NA NA NA
## 54 0 0 0 0 0 0 NA NA 0 0 0 0 0 0 NA NA
## 55 1 0 1 NA NA NA NA NA 1 1 1 NA NA NA NA NA
## y.17 y.18 y.19 y.20 y.21 y.22 y.23 y.24 y.25 y.26 y.27 y.28 y.29 y.30
## 1 1 1 NA NA NA NA NA NA 0 1 1 NA NA NA
## 2 0 0 0 NA NA NA NA NA 0 0 0 1 1 0
## 3 0 0 0 1 NA NA NA NA 0 0 0 0 0 0
## 4 0 1 1 NA NA NA NA NA 0 0 0 0 0 1
## 5 0 0 0 0 0 0 NA NA 0 0 0 0 0 0
## 6 0 1 1 1 NA NA NA NA 0 1 1 1 1 NA
## 7 0 0 0 0 0 0 NA NA 0 0 0 0 0 NA
## 8 0 0 0 0 1 NA NA NA 1 0 1 1 1 NA
## 9 0 0 0 0 0 0 NA NA 0 0 0 0 0 0
## 10 0 0 0 0 0 NA NA NA 1 1 0 1 NA NA
## 11 1 1 NA NA NA NA NA NA 1 0 0 NA NA NA
## 12 0 0 1 1 NA NA NA NA 0 0 1 0 1 1
## 13 0 0 NA NA NA NA NA NA 1 0 1 NA NA NA
## 14 0 0 0 1 1 1 NA NA 1 0 1 0 NA NA
## 15 1 1 1 NA NA NA NA NA 0 0 0 1 1 1
## 16 0 0 0 0 0 0 0 NA 0 0 0 0 0 NA
## 17 0 0 0 0 0 0 NA NA 0 0 0 0 0 0
## 18 0 0 0 0 0 0 NA NA 0 0 0 0 0 0
## 19 0 0 0 0 0 NA NA NA 0 0 0 0 0 NA
## 20 1 0 NA NA NA NA NA NA 1 1 NA NA NA NA
## 21 1 1 1 1 1 1 0 0 0 0 0 0 0 NA
## 22 1 0 0 1 0 1 1 NA 0 1 1 1 NA NA
## 23 1 1 1 NA NA NA NA NA 1 1 NA NA NA NA
## 24 0 0 0 0 0 NA NA NA 0 0 0 0 0 NA
## 25 0 0 0 0 0 NA NA NA 0 0 0 0 0 0
## 26 0 1 0 0 NA NA NA NA 0 0 0 0 0 0
## 27 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## 28 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## 29 0 0 0 0 0 NA NA NA 0 0 0 0 0 0
## 30 1 0 1 1 1 1 NA NA 1 1 0 NA NA NA
## 31 0 0 0 0 0 NA NA NA 0 0 0 0 0 0
## 32 0 0 0 0 0 NA NA NA 0 0 0 0 0 0
## 33 0 1 NA NA NA NA NA NA 0 0 0 0 0 0
## 34 0 0 0 0 1 1 NA NA 1 1 1 1 NA NA
## 35 0 1 0 0 1 NA NA NA 1 0 0 1 0 NA
## 36 0 0 0 0 1 NA NA NA 0 0 0 0 1 NA
## 37 0 0 0 0 0 0 1 1 0 0 0 1 1 NA
## 38 0 0 0 0 0 0 0 NA 0 0 0 0 0 NA
## 39 1 0 1 1 NA NA NA NA 1 0 0 0 0 0
## 40 0 0 0 0 0 1 0 NA 1 1 0 NA NA NA
## 41 0 0 0 0 0 0 0 NA 0 1 1 0 NA NA
## 42 0 0 0 NA NA NA NA NA 0 0 1 0 0 0
## 43 0 1 0 0 0 NA NA NA 0 0 0 0 NA NA
## 44 0 0 1 0 NA NA NA NA 0 0 0 0 0 0
## 45 0 0 0 0 0 NA NA NA 0 0 0 0 0 NA
## 46 0 0 0 0 0 0 0 NA 0 0 0 0 0 0
## 47 0 0 0 0 1 0 0 NA 0 1 0 0 0 0
## 48 0 0 0 0 0 0 NA NA 0 0 1 1 NA NA
## 49 0 0 0 0 1 1 0 1 0 0 0 0 0 NA
## 50 0 0 0 0 0 0 NA NA 0 0 0 0 0 0
## 51 0 0 0 0 0 0 NA NA 0 0 0 0 0 0
## 52 1 0 0 NA NA NA NA NA 0 1 0 0 0 0
## 53 0 0 0 0 0 0 NA NA 0 0 0 0 1 1
## 54 0 0 0 0 0 0 NA NA 0 0 1 1 0 NA
## 55 0 1 0 NA NA NA NA NA 1 0 NA NA NA NA
## y.31 y.32 y.33 y.34 y.35 y.36 y.37 y.38 y.39 y.40 SiteNum SiteNum.1
## 1 NA NA 0 1 0 1 NA NA NA NA 1 1
## 2 NA NA 0 0 1 1 NA NA NA NA 2 2
## 3 NA NA 0 0 0 0 0 NA NA NA 3 3
## 4 0 0 0 0 0 0 1 1 NA NA 4 4
## 5 NA NA 0 1 1 1 NA NA NA NA 5 5
## 6 NA NA 0 1 1 1 1 1 NA NA 6 6
## 7 NA NA 0 0 0 0 0 0 NA NA 7 7
## 8 NA NA 0 0 0 0 0 0 NA NA 8 8
## 9 NA NA 0 0 0 0 NA NA NA NA 9 9
## 10 NA NA 1 1 1 1 1 1 NA NA 10 10
## 11 NA NA 0 0 0 0 1 1 NA NA 11 11
## 12 1 NA 1 1 0 1 0 NA NA NA 12 12
## 13 NA NA 1 1 1 1 NA NA NA NA 13 13
## 14 NA NA 0 1 0 0 0 0 NA NA 14 14
## 15 NA NA 0 0 0 1 1 1 NA NA 15 15
## 16 NA NA 0 1 1 NA NA NA NA NA 16 16
## 17 0 0 0 0 0 0 0 0 0 0 17 17
## 18 0 NA 0 0 0 0 0 NA NA NA 18 18
## 19 NA NA 0 0 0 0 0 0 NA NA 19 19
## 20 NA NA 1 1 1 1 NA NA NA NA 20 20
## 21 NA NA 0 0 0 0 0 0 NA NA 21 21
## 22 NA NA 0 0 1 1 1 NA NA NA 22 22
## 23 NA NA 1 1 1 NA NA NA NA NA 23 23
## 24 NA NA 0 0 0 0 0 NA NA NA 24 24
## 25 NA NA 0 0 0 0 0 0 NA NA 25 25
## 26 NA NA 0 1 1 1 1 1 1 1 26 26
## 27 0 NA 0 0 0 0 0 0 NA NA 27 27
## 28 0 0 0 0 0 0 1 1 1 NA 28 28
## 29 0 NA 0 0 0 0 0 0 NA NA 29 29
## 30 NA NA 0 0 0 0 1 1 0 1 30 30
## 31 0 NA 0 0 0 0 0 0 NA NA 31 31
## 32 NA NA 0 0 0 0 0 NA NA NA 32 32
## 33 NA NA 0 0 0 0 0 0 0 0 33 33
## 34 NA NA 0 0 1 1 NA NA NA NA 34 34
## 35 NA NA 0 0 0 0 1 1 NA NA 35 35
## 36 NA NA 0 0 0 0 0 NA NA NA 36 36
## 37 NA NA 0 0 0 0 1 1 NA NA 37 37
## 38 NA NA 0 0 0 0 0 0 0 NA 38 38
## 39 NA NA 0 0 0 0 0 1 NA NA 39 39
## 40 NA NA 0 0 1 0 0 NA NA NA 40 40
## 41 NA NA 1 1 0 0 NA NA NA NA 41 41
## 42 0 0 0 0 0 0 0 0 0 NA 42 42
## 43 NA NA 0 0 0 0 0 0 NA NA 43 43
## 44 0 0 1 0 1 1 1 1 1 NA 44 44
## 45 NA NA 0 0 0 0 0 0 NA NA 45 45
## 46 1 1 0 0 0 1 1 1 1 1 46 46
## 47 0 0 1 1 1 1 1 1 1 NA 47 47
## 48 NA NA 0 0 1 1 NA NA NA NA 48 48
## 49 NA NA 0 0 0 0 0 NA NA NA 49 49
## 50 NA NA 0 0 0 0 0 0 NA NA 50 50
## 51 NA NA 0 0 0 0 0 NA NA NA 51 51
## 52 0 0 0 0 1 0 0 1 0 1 52 52
## 53 0 0 0 0 0 0 0 0 0 1 53 53
## 54 NA NA 0 0 0 1 1 NA NA NA 54 54
## 55 NA NA 0 0 0 0 0 NA NA NA 55 55
## SiteNum.2 SiteNum.3 SiteNum.4 SiteNum.5 SiteNum.6 SiteNum.7 SiteNum.8
## 1 1 1 1 1 1 1 1
## 2 2 2 2 2 2 2 2
## 3 3 3 3 3 3 3 3
## 4 4 4 4 4 4 4 4
## 5 5 5 5 5 5 5 5
## 6 6 6 6 6 6 6 6
## 7 7 7 7 7 7 7 7
## 8 8 8 8 8 8 8 8
## 9 9 9 9 9 9 9 9
## 10 10 10 10 10 10 10 10
## 11 11 11 11 11 11 11 11
## 12 12 12 12 12 12 12 12
## 13 13 13 13 13 13 13 13
## 14 14 14 14 14 14 14 14
## 15 15 15 15 15 15 15 15
## 16 16 16 16 16 16 16 16
## 17 17 17 17 17 17 17 17
## 18 18 18 18 18 18 18 18
## 19 19 19 19 19 19 19 19
## 20 20 20 20 20 20 20 20
## 21 21 21 21 21 21 21 21
## 22 22 22 22 22 22 22 22
## 23 23 23 23 23 23 23 23
## 24 24 24 24 24 24 24 24
## 25 25 25 25 25 25 25 25
## 26 26 26 26 26 26 26 26
## 27 27 27 27 27 27 27 27
## 28 28 28 28 28 28 28 28
## 29 29 29 29 29 29 29 29
## 30 30 30 30 30 30 30 30
## 31 31 31 31 31 31 31 31
## 32 32 32 32 32 32 32 32
## 33 33 33 33 33 33 33 33
## 34 34 34 34 34 34 34 34
## 35 35 35 35 35 35 35 35
## 36 36 36 36 36 36 36 36
## 37 37 37 37 37 37 37 37
## 38 38 38 38 38 38 38 38
## 39 39 39 39 39 39 39 39
## 40 40 40 40 40 40 40 40
## 41 41 41 41 41 41 41 41
## 42 42 42 42 42 42 42 42
## 43 43 43 43 43 43 43 43
## 44 44 44 44 44 44 44 44
## 45 45 45 45 45 45 45 45
## 46 46 46 46 46 46 46 46
## 47 47 47 47 47 47 47 47
## 48 48 48 48 48 48 48 48
## 49 49 49 49 49 49 49 49
## 50 50 50 50 50 50 50 50
## 51 51 51 51 51 51 51 51
## 52 52 52 52 52 52 52 52
## 53 53 53 53 53 53 53 53
## 54 54 54 54 54 54 54 54
## 55 55 55 55 55 55 55 55
## SiteNum.9 SiteNum.10 SiteNum.11 SiteNum.12 SiteNum.13 SiteNum.14
## 1 1 1 1 1 1 1
## 2 2 2 2 2 2 2
## 3 3 3 3 3 3 3
## 4 4 4 4 4 4 4
## 5 5 5 5 5 5 5
## 6 6 6 6 6 6 6
## 7 7 7 7 7 7 7
## 8 8 8 8 8 8 8
## 9 9 9 9 9 9 9
## 10 10 10 10 10 10 10
## 11 11 11 11 11 11 11
## 12 12 12 12 12 12 12
## 13 13 13 13 13 13 13
## 14 14 14 14 14 14 14
## 15 15 15 15 15 15 15
## 16 16 16 16 16 16 16
## 17 17 17 17 17 17 17
## 18 18 18 18 18 18 18
## 19 19 19 19 19 19 19
## 20 20 20 20 20 20 20
## 21 21 21 21 21 21 21
## 22 22 22 22 22 22 22
## 23 23 23 23 23 23 23
## 24 24 24 24 24 24 24
## 25 25 25 25 25 25 25
## 26 26 26 26 26 26 26
## 27 27 27 27 27 27 27
## 28 28 28 28 28 28 28
## 29 29 29 29 29 29 29
## 30 30 30 30 30 30 30
## 31 31 31 31 31 31 31
## 32 32 32 32 32 32 32
## 33 33 33 33 33 33 33
## 34 34 34 34 34 34 34
## 35 35 35 35 35 35 35
## 36 36 36 36 36 36 36
## 37 37 37 37 37 37 37
## 38 38 38 38 38 38 38
## 39 39 39 39 39 39 39
## 40 40 40 40 40 40 40
## 41 41 41 41 41 41 41
## 42 42 42 42 42 42 42
## 43 43 43 43 43 43 43
## 44 44 44 44 44 44 44
## 45 45 45 45 45 45 45
## 46 46 46 46 46 46 46
## 47 47 47 47 47 47 47
## 48 48 48 48 48 48 48
## 49 49 49 49 49 49 49
## 50 50 50 50 50 50 50
## 51 51 51 51 51 51 51
## 52 52 52 52 52 52 52
## 53 53 53 53 53 53 53
## 54 54 54 54 54 54 54
## 55 55 55 55 55 55 55
## SiteNum.15 SiteNum.16 SiteNum.17 SiteNum.18 SiteNum.19 SiteNum.20
## 1 1 1 1 1 1 1
## 2 2 2 2 2 2 2
## 3 3 3 3 3 3 3
## 4 4 4 4 4 4 4
## 5 5 5 5 5 5 5
## 6 6 6 6 6 6 6
## 7 7 7 7 7 7 7
## 8 8 8 8 8 8 8
## 9 9 9 9 9 9 9
## 10 10 10 10 10 10 10
## 11 11 11 11 11 11 11
## 12 12 12 12 12 12 12
## 13 13 13 13 13 13 13
## 14 14 14 14 14 14 14
## 15 15 15 15 15 15 15
## 16 16 16 16 16 16 16
## 17 17 17 17 17 17 17
## 18 18 18 18 18 18 18
## 19 19 19 19 19 19 19
## 20 20 20 20 20 20 20
## 21 21 21 21 21 21 21
## 22 22 22 22 22 22 22
## 23 23 23 23 23 23 23
## 24 24 24 24 24 24 24
## 25 25 25 25 25 25 25
## 26 26 26 26 26 26 26
## 27 27 27 27 27 27 27
## 28 28 28 28 28 28 28
## 29 29 29 29 29 29 29
## 30 30 30 30 30 30 30
## 31 31 31 31 31 31 31
## 32 32 32 32 32 32 32
## 33 33 33 33 33 33 33
## 34 34 34 34 34 34 34
## 35 35 35 35 35 35 35
## 36 36 36 36 36 36 36
## 37 37 37 37 37 37 37
## 38 38 38 38 38 38 38
## 39 39 39 39 39 39 39
## 40 40 40 40 40 40 40
## 41 41 41 41 41 41 41
## 42 42 42 42 42 42 42
## 43 43 43 43 43 43 43
## 44 44 44 44 44 44 44
## 45 45 45 45 45 45 45
## 46 46 46 46 46 46 46
## 47 47 47 47 47 47 47
## 48 48 48 48 48 48 48
## 49 49 49 49 49 49 49
## 50 50 50 50 50 50 50
## 51 51 51 51 51 51 51
## 52 52 52 52 52 52 52
## 53 53 53 53 53 53 53
## 54 54 54 54 54 54 54
## 55 55 55 55 55 55 55
## SiteNum.21 SiteNum.22 SiteNum.23 SiteNum.24 SiteNum.25 SiteNum.26
## 1 1 1 1 1 1 1
## 2 2 2 2 2 2 2
## 3 3 3 3 3 3 3
## 4 4 4 4 4 4 4
## 5 5 5 5 5 5 5
## 6 6 6 6 6 6 6
## 7 7 7 7 7 7 7
## 8 8 8 8 8 8 8
## 9 9 9 9 9 9 9
## 10 10 10 10 10 10 10
## 11 11 11 11 11 11 11
## 12 12 12 12 12 12 12
## 13 13 13 13 13 13 13
## 14 14 14 14 14 14 14
## 15 15 15 15 15 15 15
## 16 16 16 16 16 16 16
## 17 17 17 17 17 17 17
## 18 18 18 18 18 18 18
## 19 19 19 19 19 19 19
## 20 20 20 20 20 20 20
## 21 21 21 21 21 21 21
## 22 22 22 22 22 22 22
## 23 23 23 23 23 23 23
## 24 24 24 24 24 24 24
## 25 25 25 25 25 25 25
## 26 26 26 26 26 26 26
## 27 27 27 27 27 27 27
## 28 28 28 28 28 28 28
## 29 29 29 29 29 29 29
## 30 30 30 30 30 30 30
## 31 31 31 31 31 31 31
## 32 32 32 32 32 32 32
## 33 33 33 33 33 33 33
## 34 34 34 34 34 34 34
## 35 35 35 35 35 35 35
## 36 36 36 36 36 36 36
## 37 37 37 37 37 37 37
## 38 38 38 38 38 38 38
## 39 39 39 39 39 39 39
## 40 40 40 40 40 40 40
## 41 41 41 41 41 41 41
## 42 42 42 42 42 42 42
## 43 43 43 43 43 43 43
## 44 44 44 44 44 44 44
## 45 45 45 45 45 45 45
## 46 46 46 46 46 46 46
## 47 47 47 47 47 47 47
## 48 48 48 48 48 48 48
## 49 49 49 49 49 49 49
## 50 50 50 50 50 50 50
## 51 51 51 51 51 51 51
## 52 52 52 52 52 52 52
## 53 53 53 53 53 53 53
## 54 54 54 54 54 54 54
## 55 55 55 55 55 55 55
## SiteNum.27 SiteNum.28 SiteNum.29 SiteNum.30 SiteNum.31 SiteNum.32
## 1 1 1 1 1 1 1
## 2 2 2 2 2 2 2
## 3 3 3 3 3 3 3
## 4 4 4 4 4 4 4
## 5 5 5 5 5 5 5
## 6 6 6 6 6 6 6
## 7 7 7 7 7 7 7
## 8 8 8 8 8 8 8
## 9 9 9 9 9 9 9
## 10 10 10 10 10 10 10
## 11 11 11 11 11 11 11
## 12 12 12 12 12 12 12
## 13 13 13 13 13 13 13
## 14 14 14 14 14 14 14
## 15 15 15 15 15 15 15
## 16 16 16 16 16 16 16
## 17 17 17 17 17 17 17
## 18 18 18 18 18 18 18
## 19 19 19 19 19 19 19
## 20 20 20 20 20 20 20
## 21 21 21 21 21 21 21
## 22 22 22 22 22 22 22
## 23 23 23 23 23 23 23
## 24 24 24 24 24 24 24
## 25 25 25 25 25 25 25
## 26 26 26 26 26 26 26
## 27 27 27 27 27 27 27
## 28 28 28 28 28 28 28
## 29 29 29 29 29 29 29
## 30 30 30 30 30 30 30
## 31 31 31 31 31 31 31
## 32 32 32 32 32 32 32
## 33 33 33 33 33 33 33
## 34 34 34 34 34 34 34
## 35 35 35 35 35 35 35
## 36 36 36 36 36 36 36
## 37 37 37 37 37 37 37
## 38 38 38 38 38 38 38
## 39 39 39 39 39 39 39
## 40 40 40 40 40 40 40
## 41 41 41 41 41 41 41
## 42 42 42 42 42 42 42
## 43 43 43 43 43 43 43
## 44 44 44 44 44 44 44
## 45 45 45 45 45 45 45
## 46 46 46 46 46 46 46
## 47 47 47 47 47 47 47
## 48 48 48 48 48 48 48
## 49 49 49 49 49 49 49
## 50 50 50 50 50 50 50
## 51 51 51 51 51 51 51
## 52 52 52 52 52 52 52
## 53 53 53 53 53 53 53
## 54 54 54 54 54 54 54
## 55 55 55 55 55 55 55
## SiteNum.33 SiteNum.34 SiteNum.35 SiteNum.36 SiteNum.37 SiteNum.38
## 1 1 1 1 1 1 1
## 2 2 2 2 2 2 2
## 3 3 3 3 3 3 3
## 4 4 4 4 4 4 4
## 5 5 5 5 5 5 5
## 6 6 6 6 6 6 6
## 7 7 7 7 7 7 7
## 8 8 8 8 8 8 8
## 9 9 9 9 9 9 9
## 10 10 10 10 10 10 10
## 11 11 11 11 11 11 11
## 12 12 12 12 12 12 12
## 13 13 13 13 13 13 13
## 14 14 14 14 14 14 14
## 15 15 15 15 15 15 15
## 16 16 16 16 16 16 16
## 17 17 17 17 17 17 17
## 18 18 18 18 18 18 18
## 19 19 19 19 19 19 19
## 20 20 20 20 20 20 20
## 21 21 21 21 21 21 21
## 22 22 22 22 22 22 22
## 23 23 23 23 23 23 23
## 24 24 24 24 24 24 24
## 25 25 25 25 25 25 25
## 26 26 26 26 26 26 26
## 27 27 27 27 27 27 27
## 28 28 28 28 28 28 28
## 29 29 29 29 29 29 29
## 30 30 30 30 30 30 30
## 31 31 31 31 31 31 31
## 32 32 32 32 32 32 32
## 33 33 33 33 33 33 33
## 34 34 34 34 34 34 34
## 35 35 35 35 35 35 35
## 36 36 36 36 36 36 36
## 37 37 37 37 37 37 37
## 38 38 38 38 38 38 38
## 39 39 39 39 39 39 39
## 40 40 40 40 40 40 40
## 41 41 41 41 41 41 41
## 42 42 42 42 42 42 42
## 43 43 43 43 43 43 43
## 44 44 44 44 44 44 44
## 45 45 45 45 45 45 45
## 46 46 46 46 46 46 46
## 47 47 47 47 47 47 47
## 48 48 48 48 48 48 48
## 49 49 49 49 49 49 49
## 50 50 50 50 50 50 50
## 51 51 51 51 51 51 51
## 52 52 52 52 52 52 52
## 53 53 53 53 53 53 53
## 54 54 54 54 54 54 54
## 55 55 55 55 55 55 55
## SiteNum.39 SiteNum.40 Year.1 Year.2 Year.3 Year.4 Year.5 Year.6 Year.7
## 1 1 1 1 1 1 1 1 1 1
## 2 2 2 1 1 1 1 1 1 1
## 3 3 3 1 1 1 1 1 1 1
## 4 4 4 1 1 1 1 1 1 1
## 5 5 5 1 1 1 1 1 1 1
## 6 6 6 1 1 1 1 1 1 1
## 7 7 7 1 1 1 1 1 1 1
## 8 8 8 1 1 1 1 1 1 1
## 9 9 9 1 1 1 1 1 1 1
## 10 10 10 1 1 1 1 1 1 1
## 11 11 11 1 1 1 1 1 1 1
## 12 12 12 1 1 1 1 1 1 1
## 13 13 13 1 1 1 1 1 1 1
## 14 14 14 1 1 1 1 1 1 1
## 15 15 15 1 1 1 1 1 1 1
## 16 16 16 1 1 1 1 1 1 1
## 17 17 17 1 1 1 1 1 1 1
## 18 18 18 1 1 1 1 1 1 1
## 19 19 19 1 1 1 1 1 1 1
## 20 20 20 1 1 1 1 1 1 1
## 21 21 21 1 1 1 1 1 1 1
## 22 22 22 1 1 1 1 1 1 1
## 23 23 23 1 1 1 1 1 1 1
## 24 24 24 1 1 1 1 1 1 1
## 25 25 25 1 1 1 1 1 1 1
## 26 26 26 1 1 1 1 1 1 1
## 27 27 27 1 1 1 1 1 1 1
## 28 28 28 1 1 1 1 1 1 1
## 29 29 29 1 1 1 1 1 1 1
## 30 30 30 1 1 1 1 1 1 1
## 31 31 31 1 1 1 1 1 1 1
## 32 32 32 1 1 1 1 1 1 1
## 33 33 33 1 1 1 1 1 1 1
## 34 34 34 1 1 1 1 1 1 1
## 35 35 35 1 1 1 1 1 1 1
## 36 36 36 1 1 1 1 1 1 1
## 37 37 37 1 1 1 1 1 1 1
## 38 38 38 1 1 1 1 1 1 1
## 39 39 39 1 1 1 1 1 1 1
## 40 40 40 1 1 1 1 1 1 1
## 41 41 41 1 1 1 1 1 1 1
## 42 42 42 1 1 1 1 1 1 1
## 43 43 43 1 1 1 1 1 1 1
## 44 44 44 1 1 1 1 1 1 1
## 45 45 45 1 1 1 1 1 1 1
## 46 46 46 1 1 1 1 1 1 1
## 47 47 47 1 1 1 1 1 1 1
## 48 48 48 1 1 1 1 1 1 1
## 49 49 49 1 1 1 1 1 1 1
## 50 50 50 1 1 1 1 1 1 1
## 51 51 51 1 1 1 1 1 1 1
## 52 52 52 1 1 1 1 1 1 1
## 53 53 53 1 1 1 1 1 1 1
## 54 54 54 1 1 1 1 1 1 1
## 55 55 55 1 1 1 1 1 1 1
## Year.8 Year.9 Year.10 Year.11 Year.12 Year.13 Year.14 Year.15 Year.16
## 1 1 2 2 2 2 2 2 2 2
## 2 1 2 2 2 2 2 2 2 2
## 3 1 2 2 2 2 2 2 2 2
## 4 1 2 2 2 2 2 2 2 2
## 5 1 2 2 2 2 2 2 2 2
## 6 1 2 2 2 2 2 2 2 2
## 7 1 2 2 2 2 2 2 2 2
## 8 1 2 2 2 2 2 2 2 2
## 9 1 2 2 2 2 2 2 2 2
## 10 1 2 2 2 2 2 2 2 2
## 11 1 2 2 2 2 2 2 2 2
## 12 1 2 2 2 2 2 2 2 2
## 13 1 2 2 2 2 2 2 2 2
## 14 1 2 2 2 2 2 2 2 2
## 15 1 2 2 2 2 2 2 2 2
## 16 1 2 2 2 2 2 2 2 2
## 17 1 2 2 2 2 2 2 2 2
## 18 1 2 2 2 2 2 2 2 2
## 19 1 2 2 2 2 2 2 2 2
## 20 1 2 2 2 2 2 2 2 2
## 21 1 2 2 2 2 2 2 2 2
## 22 1 2 2 2 2 2 2 2 2
## 23 1 2 2 2 2 2 2 2 2
## 24 1 2 2 2 2 2 2 2 2
## 25 1 2 2 2 2 2 2 2 2
## 26 1 2 2 2 2 2 2 2 2
## 27 1 2 2 2 2 2 2 2 2
## 28 1 2 2 2 2 2 2 2 2
## 29 1 2 2 2 2 2 2 2 2
## 30 1 2 2 2 2 2 2 2 2
## 31 1 2 2 2 2 2 2 2 2
## 32 1 2 2 2 2 2 2 2 2
## 33 1 2 2 2 2 2 2 2 2
## 34 1 2 2 2 2 2 2 2 2
## 35 1 2 2 2 2 2 2 2 2
## 36 1 2 2 2 2 2 2 2 2
## 37 1 2 2 2 2 2 2 2 2
## 38 1 2 2 2 2 2 2 2 2
## 39 1 2 2 2 2 2 2 2 2
## 40 1 2 2 2 2 2 2 2 2
## 41 1 2 2 2 2 2 2 2 2
## 42 1 2 2 2 2 2 2 2 2
## 43 1 2 2 2 2 2 2 2 2
## 44 1 2 2 2 2 2 2 2 2
## 45 1 2 2 2 2 2 2 2 2
## 46 1 2 2 2 2 2 2 2 2
## 47 1 2 2 2 2 2 2 2 2
## 48 1 2 2 2 2 2 2 2 2
## 49 1 2 2 2 2 2 2 2 2
## 50 1 2 2 2 2 2 2 2 2
## 51 1 2 2 2 2 2 2 2 2
## 52 1 2 2 2 2 2 2 2 2
## 53 1 2 2 2 2 2 2 2 2
## 54 1 2 2 2 2 2 2 2 2
## 55 1 2 2 2 2 2 2 2 2
## Year.17 Year.18 Year.19 Year.20 Year.21 Year.22 Year.23 Year.24 Year.25
## 1 3 3 3 3 3 3 3 3 4
## 2 3 3 3 3 3 3 3 3 4
## 3 3 3 3 3 3 3 3 3 4
## 4 3 3 3 3 3 3 3 3 4
## 5 3 3 3 3 3 3 3 3 4
## 6 3 3 3 3 3 3 3 3 4
## 7 3 3 3 3 3 3 3 3 4
## 8 3 3 3 3 3 3 3 3 4
## 9 3 3 3 3 3 3 3 3 4
## 10 3 3 3 3 3 3 3 3 4
## 11 3 3 3 3 3 3 3 3 4
## 12 3 3 3 3 3 3 3 3 4
## 13 3 3 3 3 3 3 3 3 4
## 14 3 3 3 3 3 3 3 3 4
## 15 3 3 3 3 3 3 3 3 4
## 16 3 3 3 3 3 3 3 3 4
## 17 3 3 3 3 3 3 3 3 4
## 18 3 3 3 3 3 3 3 3 4
## 19 3 3 3 3 3 3 3 3 4
## 20 3 3 3 3 3 3 3 3 4
## 21 3 3 3 3 3 3 3 3 4
## 22 3 3 3 3 3 3 3 3 4
## 23 3 3 3 3 3 3 3 3 4
## 24 3 3 3 3 3 3 3 3 4
## 25 3 3 3 3 3 3 3 3 4
## 26 3 3 3 3 3 3 3 3 4
## 27 3 3 3 3 3 3 3 3 4
## 28 3 3 3 3 3 3 3 3 4
## 29 3 3 3 3 3 3 3 3 4
## 30 3 3 3 3 3 3 3 3 4
## 31 3 3 3 3 3 3 3 3 4
## 32 3 3 3 3 3 3 3 3 4
## 33 3 3 3 3 3 3 3 3 4
## 34 3 3 3 3 3 3 3 3 4
## 35 3 3 3 3 3 3 3 3 4
## 36 3 3 3 3 3 3 3 3 4
## 37 3 3 3 3 3 3 3 3 4
## 38 3 3 3 3 3 3 3 3 4
## 39 3 3 3 3 3 3 3 3 4
## 40 3 3 3 3 3 3 3 3 4
## 41 3 3 3 3 3 3 3 3 4
## 42 3 3 3 3 3 3 3 3 4
## 43 3 3 3 3 3 3 3 3 4
## 44 3 3 3 3 3 3 3 3 4
## 45 3 3 3 3 3 3 3 3 4
## 46 3 3 3 3 3 3 3 3 4
## 47 3 3 3 3 3 3 3 3 4
## 48 3 3 3 3 3 3 3 3 4
## 49 3 3 3 3 3 3 3 3 4
## 50 3 3 3 3 3 3 3 3 4
## 51 3 3 3 3 3 3 3 3 4
## 52 3 3 3 3 3 3 3 3 4
## 53 3 3 3 3 3 3 3 3 4
## 54 3 3 3 3 3 3 3 3 4
## 55 3 3 3 3 3 3 3 3 4
## Year.26 Year.27 Year.28 Year.29 Year.30 Year.31 Year.32 Year.33 Year.34
## 1 4 4 4 4 4 4 4 5 5
## 2 4 4 4 4 4 4 4 5 5
## 3 4 4 4 4 4 4 4 5 5
## 4 4 4 4 4 4 4 4 5 5
## 5 4 4 4 4 4 4 4 5 5
## 6 4 4 4 4 4 4 4 5 5
## 7 4 4 4 4 4 4 4 5 5
## 8 4 4 4 4 4 4 4 5 5
## 9 4 4 4 4 4 4 4 5 5
## 10 4 4 4 4 4 4 4 5 5
## 11 4 4 4 4 4 4 4 5 5
## 12 4 4 4 4 4 4 4 5 5
## 13 4 4 4 4 4 4 4 5 5
## 14 4 4 4 4 4 4 4 5 5
## 15 4 4 4 4 4 4 4 5 5
## 16 4 4 4 4 4 4 4 5 5
## 17 4 4 4 4 4 4 4 5 5
## 18 4 4 4 4 4 4 4 5 5
## 19 4 4 4 4 4 4 4 5 5
## 20 4 4 4 4 4 4 4 5 5
## 21 4 4 4 4 4 4 4 5 5
## 22 4 4 4 4 4 4 4 5 5
## 23 4 4 4 4 4 4 4 5 5
## 24 4 4 4 4 4 4 4 5 5
## 25 4 4 4 4 4 4 4 5 5
## 26 4 4 4 4 4 4 4 5 5
## 27 4 4 4 4 4 4 4 5 5
## 28 4 4 4 4 4 4 4 5 5
## 29 4 4 4 4 4 4 4 5 5
## 30 4 4 4 4 4 4 4 5 5
## 31 4 4 4 4 4 4 4 5 5
## 32 4 4 4 4 4 4 4 5 5
## 33 4 4 4 4 4 4 4 5 5
## 34 4 4 4 4 4 4 4 5 5
## 35 4 4 4 4 4 4 4 5 5
## 36 4 4 4 4 4 4 4 5 5
## 37 4 4 4 4 4 4 4 5 5
## 38 4 4 4 4 4 4 4 5 5
## 39 4 4 4 4 4 4 4 5 5
## 40 4 4 4 4 4 4 4 5 5
## 41 4 4 4 4 4 4 4 5 5
## 42 4 4 4 4 4 4 4 5 5
## 43 4 4 4 4 4 4 4 5 5
## 44 4 4 4 4 4 4 4 5 5
## 45 4 4 4 4 4 4 4 5 5
## 46 4 4 4 4 4 4 4 5 5
## 47 4 4 4 4 4 4 4 5 5
## 48 4 4 4 4 4 4 4 5 5
## 49 4 4 4 4 4 4 4 5 5
## 50 4 4 4 4 4 4 4 5 5
## 51 4 4 4 4 4 4 4 5 5
## 52 4 4 4 4 4 4 4 5 5
## 53 4 4 4 4 4 4 4 5 5
## 54 4 4 4 4 4 4 4 5 5
## 55 4 4 4 4 4 4 4 5 5
## Year.35 Year.36 Year.37 Year.38 Year.39 Year.40 Visit.1 Visit.2 Visit.3
## 1 5 5 5 5 5 5 1 2 3
## 2 5 5 5 5 5 5 1 2 3
## 3 5 5 5 5 5 5 1 2 3
## 4 5 5 5 5 5 5 1 2 3
## 5 5 5 5 5 5 5 1 2 3
## 6 5 5 5 5 5 5 1 2 3
## 7 5 5 5 5 5 5 1 2 3
## 8 5 5 5 5 5 5 1 2 3
## 9 5 5 5 5 5 5 1 2 3
## 10 5 5 5 5 5 5 1 2 3
## 11 5 5 5 5 5 5 1 2 3
## 12 5 5 5 5 5 5 1 2 3
## 13 5 5 5 5 5 5 1 2 3
## 14 5 5 5 5 5 5 1 2 3
## 15 5 5 5 5 5 5 1 2 3
## 16 5 5 5 5 5 5 1 2 3
## 17 5 5 5 5 5 5 1 2 3
## 18 5 5 5 5 5 5 1 2 3
## 19 5 5 5 5 5 5 1 2 3
## 20 5 5 5 5 5 5 1 2 3
## 21 5 5 5 5 5 5 1 2 3
## 22 5 5 5 5 5 5 1 2 3
## 23 5 5 5 5 5 5 1 2 3
## 24 5 5 5 5 5 5 1 2 3
## 25 5 5 5 5 5 5 1 2 3
## 26 5 5 5 5 5 5 1 2 3
## 27 5 5 5 5 5 5 1 2 3
## 28 5 5 5 5 5 5 1 2 3
## 29 5 5 5 5 5 5 1 2 3
## 30 5 5 5 5 5 5 1 2 3
## 31 5 5 5 5 5 5 1 2 3
## 32 5 5 5 5 5 5 1 2 3
## 33 5 5 5 5 5 5 1 2 3
## 34 5 5 5 5 5 5 1 2 3
## 35 5 5 5 5 5 5 1 2 3
## 36 5 5 5 5 5 5 1 2 3
## 37 5 5 5 5 5 5 1 2 3
## 38 5 5 5 5 5 5 1 2 3
## 39 5 5 5 5 5 5 1 2 3
## 40 5 5 5 5 5 5 1 2 3
## 41 5 5 5 5 5 5 1 2 3
## 42 5 5 5 5 5 5 1 2 3
## 43 5 5 5 5 5 5 1 2 3
## 44 5 5 5 5 5 5 1 2 3
## 45 5 5 5 5 5 5 1 2 3
## 46 5 5 5 5 5 5 1 2 3
## 47 5 5 5 5 5 5 1 2 3
## 48 5 5 5 5 5 5 1 2 3
## 49 5 5 5 5 5 5 1 2 3
## 50 5 5 5 5 5 5 1 2 3
## 51 5 5 5 5 5 5 1 2 3
## 52 5 5 5 5 5 5 1 2 3
## 53 5 5 5 5 5 5 1 2 3
## 54 5 5 5 5 5 5 1 2 3
## 55 5 5 5 5 5 5 1 2 3
## Visit.4 Visit.5 Visit.6 Visit.7 Visit.8 Visit.9 Visit.10 Visit.11
## 1 4 5 6 7 8 1 2 3
## 2 4 5 6 7 8 1 2 3
## 3 4 5 6 7 8 1 2 3
## 4 4 5 6 7 8 1 2 3
## 5 4 5 6 7 8 1 2 3
## 6 4 5 6 7 8 1 2 3
## 7 4 5 6 7 8 1 2 3
## 8 4 5 6 7 8 1 2 3
## 9 4 5 6 7 8 1 2 3
## 10 4 5 6 7 8 1 2 3
## 11 4 5 6 7 8 1 2 3
## 12 4 5 6 7 8 1 2 3
## 13 4 5 6 7 8 1 2 3
## 14 4 5 6 7 8 1 2 3
## 15 4 5 6 7 8 1 2 3
## 16 4 5 6 7 8 1 2 3
## 17 4 5 6 7 8 1 2 3
## 18 4 5 6 7 8 1 2 3
## 19 4 5 6 7 8 1 2 3
## 20 4 5 6 7 8 1 2 3
## 21 4 5 6 7 8 1 2 3
## 22 4 5 6 7 8 1 2 3
## 23 4 5 6 7 8 1 2 3
## 24 4 5 6 7 8 1 2 3
## 25 4 5 6 7 8 1 2 3
## 26 4 5 6 7 8 1 2 3
## 27 4 5 6 7 8 1 2 3
## 28 4 5 6 7 8 1 2 3
## 29 4 5 6 7 8 1 2 3
## 30 4 5 6 7 8 1 2 3
## 31 4 5 6 7 8 1 2 3
## 32 4 5 6 7 8 1 2 3
## 33 4 5 6 7 8 1 2 3
## 34 4 5 6 7 8 1 2 3
## 35 4 5 6 7 8 1 2 3
## 36 4 5 6 7 8 1 2 3
## 37 4 5 6 7 8 1 2 3
## 38 4 5 6 7 8 1 2 3
## 39 4 5 6 7 8 1 2 3
## 40 4 5 6 7 8 1 2 3
## 41 4 5 6 7 8 1 2 3
## 42 4 5 6 7 8 1 2 3
## 43 4 5 6 7 8 1 2 3
## 44 4 5 6 7 8 1 2 3
## 45 4 5 6 7 8 1 2 3
## 46 4 5 6 7 8 1 2 3
## 47 4 5 6 7 8 1 2 3
## 48 4 5 6 7 8 1 2 3
## 49 4 5 6 7 8 1 2 3
## 50 4 5 6 7 8 1 2 3
## 51 4 5 6 7 8 1 2 3
## 52 4 5 6 7 8 1 2 3
## 53 4 5 6 7 8 1 2 3
## 54 4 5 6 7 8 1 2 3
## 55 4 5 6 7 8 1 2 3
## Visit.12 Visit.13 Visit.14 Visit.15 Visit.16 Visit.17 Visit.18 Visit.19
## 1 4 5 6 7 8 1 2 3
## 2 4 5 6 7 8 1 2 3
## 3 4 5 6 7 8 1 2 3
## 4 4 5 6 7 8 1 2 3
## 5 4 5 6 7 8 1 2 3
## 6 4 5 6 7 8 1 2 3
## 7 4 5 6 7 8 1 2 3
## 8 4 5 6 7 8 1 2 3
## 9 4 5 6 7 8 1 2 3
## 10 4 5 6 7 8 1 2 3
## 11 4 5 6 7 8 1 2 3
## 12 4 5 6 7 8 1 2 3
## 13 4 5 6 7 8 1 2 3
## 14 4 5 6 7 8 1 2 3
## 15 4 5 6 7 8 1 2 3
## 16 4 5 6 7 8 1 2 3
## 17 4 5 6 7 8 1 2 3
## 18 4 5 6 7 8 1 2 3
## 19 4 5 6 7 8 1 2 3
## 20 4 5 6 7 8 1 2 3
## 21 4 5 6 7 8 1 2 3
## 22 4 5 6 7 8 1 2 3
## 23 4 5 6 7 8 1 2 3
## 24 4 5 6 7 8 1 2 3
## 25 4 5 6 7 8 1 2 3
## 26 4 5 6 7 8 1 2 3
## 27 4 5 6 7 8 1 2 3
## 28 4 5 6 7 8 1 2 3
## 29 4 5 6 7 8 1 2 3
## 30 4 5 6 7 8 1 2 3
## 31 4 5 6 7 8 1 2 3
## 32 4 5 6 7 8 1 2 3
## 33 4 5 6 7 8 1 2 3
## 34 4 5 6 7 8 1 2 3
## 35 4 5 6 7 8 1 2 3
## 36 4 5 6 7 8 1 2 3
## 37 4 5 6 7 8 1 2 3
## 38 4 5 6 7 8 1 2 3
## 39 4 5 6 7 8 1 2 3
## 40 4 5 6 7 8 1 2 3
## 41 4 5 6 7 8 1 2 3
## 42 4 5 6 7 8 1 2 3
## 43 4 5 6 7 8 1 2 3
## 44 4 5 6 7 8 1 2 3
## 45 4 5 6 7 8 1 2 3
## 46 4 5 6 7 8 1 2 3
## 47 4 5 6 7 8 1 2 3
## 48 4 5 6 7 8 1 2 3
## 49 4 5 6 7 8 1 2 3
## 50 4 5 6 7 8 1 2 3
## 51 4 5 6 7 8 1 2 3
## 52 4 5 6 7 8 1 2 3
## 53 4 5 6 7 8 1 2 3
## 54 4 5 6 7 8 1 2 3
## 55 4 5 6 7 8 1 2 3
## Visit.20 Visit.21 Visit.22 Visit.23 Visit.24 Visit.25 Visit.26 Visit.27
## 1 4 5 6 7 8 1 2 3
## 2 4 5 6 7 8 1 2 3
## 3 4 5 6 7 8 1 2 3
## 4 4 5 6 7 8 1 2 3
## 5 4 5 6 7 8 1 2 3
## 6 4 5 6 7 8 1 2 3
## 7 4 5 6 7 8 1 2 3
## 8 4 5 6 7 8 1 2 3
## 9 4 5 6 7 8 1 2 3
## 10 4 5 6 7 8 1 2 3
## 11 4 5 6 7 8 1 2 3
## 12 4 5 6 7 8 1 2 3
## 13 4 5 6 7 8 1 2 3
## 14 4 5 6 7 8 1 2 3
## 15 4 5 6 7 8 1 2 3
## 16 4 5 6 7 8 1 2 3
## 17 4 5 6 7 8 1 2 3
## 18 4 5 6 7 8 1 2 3
## 19 4 5 6 7 8 1 2 3
## 20 4 5 6 7 8 1 2 3
## 21 4 5 6 7 8 1 2 3
## 22 4 5 6 7 8 1 2 3
## 23 4 5 6 7 8 1 2 3
## 24 4 5 6 7 8 1 2 3
## 25 4 5 6 7 8 1 2 3
## 26 4 5 6 7 8 1 2 3
## 27 4 5 6 7 8 1 2 3
## 28 4 5 6 7 8 1 2 3
## 29 4 5 6 7 8 1 2 3
## 30 4 5 6 7 8 1 2 3
## 31 4 5 6 7 8 1 2 3
## 32 4 5 6 7 8 1 2 3
## 33 4 5 6 7 8 1 2 3
## 34 4 5 6 7 8 1 2 3
## 35 4 5 6 7 8 1 2 3
## 36 4 5 6 7 8 1 2 3
## 37 4 5 6 7 8 1 2 3
## 38 4 5 6 7 8 1 2 3
## 39 4 5 6 7 8 1 2 3
## 40 4 5 6 7 8 1 2 3
## 41 4 5 6 7 8 1 2 3
## 42 4 5 6 7 8 1 2 3
## 43 4 5 6 7 8 1 2 3
## 44 4 5 6 7 8 1 2 3
## 45 4 5 6 7 8 1 2 3
## 46 4 5 6 7 8 1 2 3
## 47 4 5 6 7 8 1 2 3
## 48 4 5 6 7 8 1 2 3
## 49 4 5 6 7 8 1 2 3
## 50 4 5 6 7 8 1 2 3
## 51 4 5 6 7 8 1 2 3
## 52 4 5 6 7 8 1 2 3
## 53 4 5 6 7 8 1 2 3
## 54 4 5 6 7 8 1 2 3
## 55 4 5 6 7 8 1 2 3
## Visit.28 Visit.29 Visit.30 Visit.31 Visit.32 Visit.33 Visit.34 Visit.35
## 1 4 5 6 7 8 1 2 3
## 2 4 5 6 7 8 1 2 3
## 3 4 5 6 7 8 1 2 3
## 4 4 5 6 7 8 1 2 3
## 5 4 5 6 7 8 1 2 3
## 6 4 5 6 7 8 1 2 3
## 7 4 5 6 7 8 1 2 3
## 8 4 5 6 7 8 1 2 3
## 9 4 5 6 7 8 1 2 3
## 10 4 5 6 7 8 1 2 3
## 11 4 5 6 7 8 1 2 3
## 12 4 5 6 7 8 1 2 3
## 13 4 5 6 7 8 1 2 3
## 14 4 5 6 7 8 1 2 3
## 15 4 5 6 7 8 1 2 3
## 16 4 5 6 7 8 1 2 3
## 17 4 5 6 7 8 1 2 3
## 18 4 5 6 7 8 1 2 3
## 19 4 5 6 7 8 1 2 3
## 20 4 5 6 7 8 1 2 3
## 21 4 5 6 7 8 1 2 3
## 22 4 5 6 7 8 1 2 3
## 23 4 5 6 7 8 1 2 3
## 24 4 5 6 7 8 1 2 3
## 25 4 5 6 7 8 1 2 3
## 26 4 5 6 7 8 1 2 3
## 27 4 5 6 7 8 1 2 3
## 28 4 5 6 7 8 1 2 3
## 29 4 5 6 7 8 1 2 3
## 30 4 5 6 7 8 1 2 3
## 31 4 5 6 7 8 1 2 3
## 32 4 5 6 7 8 1 2 3
## 33 4 5 6 7 8 1 2 3
## 34 4 5 6 7 8 1 2 3
## 35 4 5 6 7 8 1 2 3
## 36 4 5 6 7 8 1 2 3
## 37 4 5 6 7 8 1 2 3
## 38 4 5 6 7 8 1 2 3
## 39 4 5 6 7 8 1 2 3
## 40 4 5 6 7 8 1 2 3
## 41 4 5 6 7 8 1 2 3
## 42 4 5 6 7 8 1 2 3
## 43 4 5 6 7 8 1 2 3
## 44 4 5 6 7 8 1 2 3
## 45 4 5 6 7 8 1 2 3
## 46 4 5 6 7 8 1 2 3
## 47 4 5 6 7 8 1 2 3
## 48 4 5 6 7 8 1 2 3
## 49 4 5 6 7 8 1 2 3
## 50 4 5 6 7 8 1 2 3
## 51 4 5 6 7 8 1 2 3
## 52 4 5 6 7 8 1 2 3
## 53 4 5 6 7 8 1 2 3
## 54 4 5 6 7 8 1 2 3
## 55 4 5 6 7 8 1 2 3
## Visit.36 Visit.37 Visit.38 Visit.39 Visit.40 SiteNum.1.1 SiteNum.2.1
## 1 4 5 6 7 8 1 1
## 2 4 5 6 7 8 2 2
## 3 4 5 6 7 8 3 3
## 4 4 5 6 7 8 4 4
## 5 4 5 6 7 8 5 5
## 6 4 5 6 7 8 6 6
## 7 4 5 6 7 8 7 7
## 8 4 5 6 7 8 8 8
## 9 4 5 6 7 8 9 9
## 10 4 5 6 7 8 10 10
## 11 4 5 6 7 8 11 11
## 12 4 5 6 7 8 12 12
## 13 4 5 6 7 8 13 13
## 14 4 5 6 7 8 14 14
## 15 4 5 6 7 8 15 15
## 16 4 5 6 7 8 16 16
## 17 4 5 6 7 8 17 17
## 18 4 5 6 7 8 18 18
## 19 4 5 6 7 8 19 19
## 20 4 5 6 7 8 20 20
## 21 4 5 6 7 8 21 21
## 22 4 5 6 7 8 22 22
## 23 4 5 6 7 8 23 23
## 24 4 5 6 7 8 24 24
## 25 4 5 6 7 8 25 25
## 26 4 5 6 7 8 26 26
## 27 4 5 6 7 8 27 27
## 28 4 5 6 7 8 28 28
## 29 4 5 6 7 8 29 29
## 30 4 5 6 7 8 30 30
## 31 4 5 6 7 8 31 31
## 32 4 5 6 7 8 32 32
## 33 4 5 6 7 8 33 33
## 34 4 5 6 7 8 34 34
## 35 4 5 6 7 8 35 35
## 36 4 5 6 7 8 36 36
## 37 4 5 6 7 8 37 37
## 38 4 5 6 7 8 38 38
## 39 4 5 6 7 8 39 39
## 40 4 5 6 7 8 40 40
## 41 4 5 6 7 8 41 41
## 42 4 5 6 7 8 42 42
## 43 4 5 6 7 8 43 43
## 44 4 5 6 7 8 44 44
## 45 4 5 6 7 8 45 45
## 46 4 5 6 7 8 46 46
## 47 4 5 6 7 8 47 47
## 48 4 5 6 7 8 48 48
## 49 4 5 6 7 8 49 49
## 50 4 5 6 7 8 50 50
## 51 4 5 6 7 8 51 51
## 52 4 5 6 7 8 52 52
## 53 4 5 6 7 8 53 53
## 54 4 5 6 7 8 54 54
## 55 4 5 6 7 8 55 55
## SiteNum.3.1 SiteNum.4.1 SiteNum.5.1 Year.1.1 Year.2.1 Year.3.1 Year.4.1
## 1 1 1 1 1 2 3 4
## 2 2 2 2 1 2 3 4
## 3 3 3 3 1 2 3 4
## 4 4 4 4 1 2 3 4
## 5 5 5 5 1 2 3 4
## 6 6 6 6 1 2 3 4
## 7 7 7 7 1 2 3 4
## 8 8 8 8 1 2 3 4
## 9 9 9 9 1 2 3 4
## 10 10 10 10 1 2 3 4
## 11 11 11 11 1 2 3 4
## 12 12 12 12 1 2 3 4
## 13 13 13 13 1 2 3 4
## 14 14 14 14 1 2 3 4
## 15 15 15 15 1 2 3 4
## 16 16 16 16 1 2 3 4
## 17 17 17 17 1 2 3 4
## 18 18 18 18 1 2 3 4
## 19 19 19 19 1 2 3 4
## 20 20 20 20 1 2 3 4
## 21 21 21 21 1 2 3 4
## 22 22 22 22 1 2 3 4
## 23 23 23 23 1 2 3 4
## 24 24 24 24 1 2 3 4
## 25 25 25 25 1 2 3 4
## 26 26 26 26 1 2 3 4
## 27 27 27 27 1 2 3 4
## 28 28 28 28 1 2 3 4
## 29 29 29 29 1 2 3 4
## 30 30 30 30 1 2 3 4
## 31 31 31 31 1 2 3 4
## 32 32 32 32 1 2 3 4
## 33 33 33 33 1 2 3 4
## 34 34 34 34 1 2 3 4
## 35 35 35 35 1 2 3 4
## 36 36 36 36 1 2 3 4
## 37 37 37 37 1 2 3 4
## 38 38 38 38 1 2 3 4
## 39 39 39 39 1 2 3 4
## 40 40 40 40 1 2 3 4
## 41 41 41 41 1 2 3 4
## 42 42 42 42 1 2 3 4
## 43 43 43 43 1 2 3 4
## 44 44 44 44 1 2 3 4
## 45 45 45 45 1 2 3 4
## 46 46 46 46 1 2 3 4
## 47 47 47 47 1 2 3 4
## 48 48 48 48 1 2 3 4
## 49 49 49 49 1 2 3 4
## 50 50 50 50 1 2 3 4
## 51 51 51 51 1 2 3 4
## 52 52 52 52 1 2 3 4
## 53 53 53 53 1 2 3 4
## 54 54 54 54 1 2 3 4
## 55 55 55 55 1 2 3 4
## Year.5.1
## 1 5
## 2 5
## 3 5
## 4 5
## 5 5
## 6 5
## 7 5
## 8 5
## 9 5
## 10 5
## 11 5
## 12 5
## 13 5
## 14 5
## 15 5
## 16 5
## 17 5
## 18 5
## 19 5
## 20 5
## 21 5
## 22 5
## 23 5
## 24 5
## 25 5
## 26 5
## 27 5
## 28 5
## 29 5
## 30 5
## 31 5
## 32 5
## 33 5
## 34 5
## 35 5
## 36 5
## 37 5
## 38 5
## 39 5
## 40 5
## 41 5
## 42 5
## 43 5
## 44 5
## 45 5
## 46 5
## 47 5
## 48 5
## 49 5
## 50 5
## 51 5
## 52 5
## 53 5
## 54 5
## 55 5summary(nso.UMF)## unmarkedFrame Object
##
## 55 sites
## Maximum number of observations per site: 40
## Mean number of observations per site: 26.33
## Number of primary survey periods: 5
## Number of secondary survey periods: 8
## Sites with at least one detection: 46
##
## Tabulation of y observations:
## 0 1 <NA>
## 1052 396 752
##
## Site-level covariates:
## SiteNum
## Min. : 1.0
## 1st Qu.:14.5
## Median :28.0
## Mean :28.0
## 3rd Qu.:41.5
## Max. :55.0
##
## Observation-level covariates:
## SiteNum Year Visit
## Min. : 1 1:440 1 :275
## 1st Qu.:14 2:440 2 :275
## Median :28 3:440 3 :275
## Mean :28 4:440 4 :275
## 3rd Qu.:42 5:440 5 :275
## Max. :55 6 :275
## (Other):550
##
## Yearly-site-level covariates:
## SiteNum Year
## Min. : 1 1:55
## 1st Qu.:14 2:55
## Median :28 3:55
## Mean :28 4:55
## 3rd Qu.:42 5:55
## Max. :55plot(nso.UMF)
# Fit some models.
# We start with the most basic model for psi, gamma, epsilon, p
# Note that formula DO HAVE AN = SIGN now
mod.psiDot.gDot.eDot.pDot <- unmarked::colext(
psiformula= ~1,
gammaformula = ~ 1,
epsilonformula = ~ 1,
pformula = ~ 1,
data=nso.UMF,
se=TRUE)
summary(mod.psiDot.gDot.eDot.pDot)##
## Call:
## unmarked::colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
## pformula = ~1, data = nso.UMF, se = TRUE)
##
## Initial (logit-scale):
## Estimate SE z P(>|z|)
## 0.537 0.289 1.86 0.063
##
## Colonization (logit-scale):
## Estimate SE z P(>|z|)
## -1.49 0.284 -5.23 1.65e-07
##
## Extinction (logit-scale):
## Estimate SE z P(>|z|)
## -1.73 0.26 -6.66 2.73e-11
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -0.021 0.0743 -0.283 0.777
##
## AIC: 1363.32
## Number of sites: 55
## optim convergence code: 0
## optim iterations: 23
## Bootstrap iterations: 0slotNames(mod.psiDot.gDot.eDot.pDot)## [1] "phi" "psiformula" "gamformula"
## [4] "epsformula" "detformula" "projected"
## [7] "projected.mean" "smoothed" "smoothed.mean"
## [10] "projected.mean.bsse" "smoothed.mean.bsse" "fitType"
## [13] "call" "formula" "data"
## [16] "sitesRemoved" "estimates" "AIC"
## [19] "opt" "negLogLike" "nllFun"
## [22] "bootstrapSamples" "covMatBS"# We get some estimates using predict in the usual way
names(mod.psiDot.gDot.eDot.pDot)## [1] "psi" "col" "ext" "det"newdata <- data.frame(intercept=1)
predict(mod.psiDot.gDot.eDot.pDot, type='psi', newdata=newdata )## Predicted SE lower upper
## 1 0.6311598 0.06726424 0.4927218 0.7509165predict(mod.psiDot.gDot.eDot.pDot, type='col', newdata=newdata )## Predicted SE lower upper
## 1 0.1841754 0.0427179 0.1145044 0.2827042predict(mod.psiDot.gDot.eDot.pDot, type='ext', newdata=newdata )## Predicted SE lower upper
## 1 0.1507459 0.03322975 0.09643335 0.2279317predict(mod.psiDot.gDot.eDot.pDot, type='det', newdata=newdata )## Predicted SE lower upper
## 1 0.4947425 0.01858223 0.4584141 0.5311265# estimates of psi moving ahead.
# There are two types.
# - population level occupancy, i.e. for the sampled population from which the sites were selected
# - actual site level occupancy, i.e. for the sites actually sampled.
# What is the occupancy for the population of ALL potential sites
# first index = not occupied/occupied, second index=year, last index=site
mod.psiDot.gDot.eDot.pDot@projected[,,1]## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.3688402 0.3960537 0.4141528 0.4261901 0.4341959
## [2,] 0.6311598 0.6039463 0.5858472 0.5738099 0.5658041mod.psiDot.gDot.eDot.pDot@projected[2,,1] ## [1] 0.6311598 0.6039463 0.5858472 0.5738099 0.5658041projected(mod.psiDot.gDot.eDot.pDot) # this gives the mean population occupancy over all sites## 1 2 3 4 5
## unoccupied 0.3688402 0.3960537 0.4141528 0.4261901 0.4341959
## occupied 0.6311598 0.6039463 0.5858472 0.5738099 0.5658041# What is the actual occupancy for these particular sites
mod.psiDot.gDot.eDot.pDot@smoothed[2,,1:2]## [,1] [,2]
## [1,] 1 0.1069976
## [2,] 1 0.0544300
## [3,] 1 0.1538519
## [4,] 1 1.0000000
## [5,] 1 1.0000000smoothed(mod.psiDot.gDot.eDot.pDot)[2,] # this gives the mean actual occupancy over all sites## 1 2 3 4 5
## 0.6311617 0.6084676 0.5553642 0.5756990 0.5738898# Obtain standard errors for these forecasts using parametric bootstrapping
psi.projected <- function(model, site) {
logit <- function(x) log(x/(1-x))
expit <- function(x) 1/(1+exp(-x))
psi.projected <- model@projected[2,,site]
# compute population growth
lambda <- exp(diff(log(psi.projected),1))
lambda.overall <- prod(lambda) # overall growth rate over entire set of seasons
lambda.prime <- exp(diff(logit(psi.projected),1))
lambda.prime.overall <- prod(lambda.prime) # overall growth rate on logit scale over entire set of season
c(psi.projected=psi.projected,
lambda= lambda,
lambda.overall=lambda.overall,
lambda.prime =lambda.prime,
lambda.prime.overall=lambda.prime.overall )
}
psi.projected(mod.psiDot.gDot.eDot.pDot, site=1)## psi.projected1 psi.projected2 psi.projected3
## 0.6311598 0.6039463 0.5858472
## psi.projected4 psi.projected5 lambda1
## 0.5738099 0.5658041 0.9568833
## lambda2 lambda3 lambda4
## 0.9700319 0.9794531 0.9860480
## lambda.overall lambda.prime1 lambda.prime2
## 0.8964514 0.8911343 0.9276400
## lambda.prime3 lambda.prime4 lambda.prime.overall
## 0.9517894 0.9678671 0.7615164# use bootstrapping to find standard errors.
# You likely want to do 1000 bootstrap samples
mod.psiDot.gDot.eDot.pDot.boot <- parboot(mod.psiDot.gDot.eDot.pDot, statistic=psi.projected, nsim=100, site=1)
slotNames(mod.psiDot.gDot.eDot.pDot.boot)## [1] "call" "t0" "t.star"mod.psiDot.gDot.eDot.pDot.boot@t0## psi.projected1 psi.projected2 psi.projected3
## 0.6311598 0.6039463 0.5858472
## psi.projected4 psi.projected5 lambda1
## 0.5738099 0.5658041 0.9568833
## lambda2 lambda3 lambda4
## 0.9700319 0.9794531 0.9860480
## lambda.overall lambda.prime1 lambda.prime2
## 0.8964514 0.8911343 0.9276400
## lambda.prime3 lambda.prime4 lambda.prime.overall
## 0.9517894 0.9678671 0.7615164cbind(est=mod.psiDot.gDot.eDot.pDot.boot@t0,
se=apply(mod.psiDot.gDot.eDot.pDot.boot@t.star, 2, sd),
t(apply(mod.psiDot.gDot.eDot.pDot.boot@t.star, 2, quantile, probs=c(0.025, 0.975))))## est se 2.5% 97.5%
## psi.projected1 0.6311598 0.05979749 0.5258870 0.7335006
## psi.projected2 0.6039463 0.04513659 0.5156419 0.7047745
## psi.projected3 0.5858472 0.04937316 0.4826215 0.6842715
## psi.projected4 0.5738099 0.05784739 0.4594802 0.6784031
## psi.projected5 0.5658041 0.06534429 0.4416468 0.6811876
## lambda1 0.9568833 0.05621871 0.8780643 1.0724001
## lambda2 0.9700319 0.03996887 0.8954816 1.0467187
## lambda3 0.9794531 0.02877949 0.9180826 1.0305584
## lambda4 0.9860480 0.02113901 0.9440630 1.0217682
## lambda.overall 0.8964514 0.13532331 0.6644115 1.1834126
## lambda.prime1 0.8911343 0.14084369 0.6725713 1.1765924
## lambda.prime2 0.9276400 0.09310669 0.7727288 1.1214741
## lambda.prime3 0.9517894 0.06360376 0.8408314 1.0824218
## lambda.prime4 0.9678671 0.04458527 0.8907813 1.0577095
## lambda.prime.overall 0.7615164 0.30532331 0.3830622 1.5183664# Fit some models.
# Model for ps(.)i, gamma(.), epsilon(.), p(year)
# Note that formula DO HAVE AN = SIGN now
mod.psiDot.gDot.eDot.pYear <- unmarked::colext(
psiformula= ~1,
gammaformula = ~ 1,
epsilonformula = ~ 1,
pformula = ~ Year,
data=nso.UMF,
se=TRUE)
summary(mod.psiDot.gDot.eDot.pYear)##
## Call:
## unmarked::colext(psiformula = ~1, gammaformula = ~1, epsilonformula = ~1,
## pformula = ~Year, data = nso.UMF, se = TRUE)
##
## Initial (logit-scale):
## Estimate SE z P(>|z|)
## 0.51 0.285 1.79 0.0741
##
## Colonization (logit-scale):
## Estimate SE z P(>|z|)
## -1.52 0.293 -5.2 1.95e-07
##
## Extinction (logit-scale):
## Estimate SE z P(>|z|)
## -1.8 0.269 -6.68 2.33e-11
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.366 0.163 2.247 0.024668
## Year2 -0.276 0.229 -1.205 0.228334
## Year3 -0.741 0.245 -3.023 0.002499
## Year4 -0.835 0.239 -3.501 0.000464
## Year5 -0.220 0.227 -0.969 0.332495
##
## AIC: 1353.528
## Number of sites: 55
## optim convergence code: 0
## optim iterations: 39
## Bootstrap iterations: 0slotNames(mod.psiDot.gDot.eDot.pYear)## [1] "phi" "psiformula" "gamformula"
## [4] "epsformula" "detformula" "projected"
## [7] "projected.mean" "smoothed" "smoothed.mean"
## [10] "projected.mean.bsse" "smoothed.mean.bsse" "fitType"
## [13] "call" "formula" "data"
## [16] "sitesRemoved" "estimates" "AIC"
## [19] "opt" "negLogLike" "nllFun"
## [22] "bootstrapSamples" "covMatBS"# We get some estimates using predict in the usual way
names(mod.psiDot.gDot.eDot.pYear)## [1] "psi" "col" "ext" "det"newdata <- data.frame(Year=as.factor(1:5))
predict(mod.psiDot.gDot.eDot.pYear, type='psi', newdata=newdata )## Predicted SE lower upper
## 1 0.6247288 0.06689193 0.4876139 0.7443861
## 2 0.6247288 0.06689193 0.4876139 0.7443861
## 3 0.6247288 0.06689193 0.4876139 0.7443861
## 4 0.6247288 0.06689193 0.4876139 0.7443861
## 5 0.6247288 0.06689193 0.4876139 0.7443861predict(mod.psiDot.gDot.eDot.pYear, type='col', newdata=newdata )## Predicted SE lower upper
## 1 0.1788409 0.04301406 0.1092631 0.2788544
## 2 0.1788409 0.04301406 0.1092631 0.2788544
## 3 0.1788409 0.04301406 0.1092631 0.2788544
## 4 0.1788409 0.04301406 0.1092631 0.2788544
## 5 0.1788409 0.04301406 0.1092631 0.2788544predict(mod.psiDot.gDot.eDot.pYear, type='ext', newdata=newdata )## Predicted SE lower upper
## 1 0.1423004 0.03280271 0.08922963 0.2193345
## 2 0.1423004 0.03280271 0.08922963 0.2193345
## 3 0.1423004 0.03280271 0.08922963 0.2193345
## 4 0.1423004 0.03280271 0.08922963 0.2193345
## 5 0.1423004 0.03280271 0.08922963 0.2193345predict(mod.psiDot.gDot.eDot.pYear, type='det', newdata=newdata )## Predicted SE lower upper
## 1 0.5905707 0.03942752 0.5116810 0.6650582
## 2 0.5224669 0.04047178 0.4432440 0.6005754
## 3 0.4074948 0.04418033 0.3245423 0.4960781
## 4 0.3848536 0.04125135 0.3077810 0.4681730
## 5 0.5365865 0.03916093 0.4595733 0.6118931# estimates of psi moving ahead.
# There are two types.
# - population level occupancy, i.e. for the sampled population from which the sites were selected
# - actual site level occupancy, i.e. for the sites actually sampled.
# What is the occupancy for the population of ALL potential sites
# first index = not occupied/occupied, second index=year, last index=site
mod.psiDot.gDot.eDot.pYear@projected[,,1]## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.3752712 0.3970566 0.4118457 0.4218855 0.4287011
## [2,] 0.6247288 0.6029434 0.5881543 0.5781145 0.5712989mod.psiDot.gDot.eDot.pYear@projected[2,,1] ## [1] 0.6247288 0.6029434 0.5881543 0.5781145 0.5712989projected(mod.psiDot.gDot.eDot.pYear) # this gives the mean population occupancy over all sites## 1 2 3 4 5
## unoccupied 0.3752712 0.3970566 0.4118457 0.4218855 0.4287011
## occupied 0.6247288 0.6029434 0.5881543 0.5781145 0.5712989# What is the actual occupancy for these particular sites
mod.psiDot.gDot.eDot.pYear@smoothed[2,,1:2]## [,1] [,2]
## [1,] 1 0.07252568
## [2,] 1 0.05016299
## [3,] 1 0.21254094
## [4,] 1 1.00000000
## [5,] 1 1.00000000smoothed(mod.psiDot.gDot.eDot.pYear)[2,] # this gives the mean actual occupancy over all sites## 1 2 3 4 5
## 0.6247320 0.6067439 0.5689525 0.5943133 0.5710464# Obtain standard errors for these forecasts using parametric bootstrapping
psi.projected <- function(model, site) {
logit <- function(x) log(x/(1-x))
expit <- function(x) 1/(1+exp(-x))
psi.projected <- model@projected[2,,site]
# compute population growth
lambda <- exp(diff(log(psi.projected),1))
lambda.overall <- prod(lambda) # overall growth rate over entire set of seasons
lambda.prime <- exp(diff(logit(psi.projected),1))
lambda.prime.overall <- prod(lambda.prime) # overall growth rate on logit scale over entire set of season
c(psi.projected=psi.projected,
lambda= lambda,
lambda.overall=lambda.overall,
lambda.prime =lambda.prime,
lambda.prime.overall=lambda.prime.overall )
}
psi.projected(mod.psiDot.gDot.eDot.pYear, site=1)## psi.projected1 psi.projected2 psi.projected3
## 0.6247288 0.6029434 0.5881543
## psi.projected4 psi.projected5 lambda1
## 0.5781145 0.5712989 0.9651283
## lambda2 lambda3 lambda4
## 0.9754717 0.9829301 0.9882107
## lambda.overall lambda.prime1 lambda.prime2
## 0.9144751 0.9121745 0.9404430
## lambda.prime3 lambda.prime4 lambda.prime.overall
## 0.9595389 0.9724999 0.8005023# use bootstrapping to find standard errors.
# You likely want to do 1000 bootstrap samples
mod.psiDot.gDot.eDot.pYear.boot <- parboot(mod.psiDot.gDot.eDot.pYear, statistic=psi.projected, nsim=100, site=1)
slotNames(mod.psiDot.gDot.eDot.pYear.boot)## [1] "call" "t0" "t.star"mod.psiDot.gDot.eDot.pYear.boot@t0## psi.projected1 psi.projected2 psi.projected3
## 0.6247288 0.6029434 0.5881543
## psi.projected4 psi.projected5 lambda1
## 0.5781145 0.5712989 0.9651283
## lambda2 lambda3 lambda4
## 0.9754717 0.9829301 0.9882107
## lambda.overall lambda.prime1 lambda.prime2
## 0.9144751 0.9121745 0.9404430
## lambda.prime3 lambda.prime4 lambda.prime.overall
## 0.9595389 0.9724999 0.8005023cbind(est=mod.psiDot.gDot.eDot.pYear.boot@t0,
se=apply(mod.psiDot.gDot.eDot.pYear.boot@t.star, 2, sd),
t(apply(mod.psiDot.gDot.eDot.pYear.boot@t.star, 2, quantile, probs=c(0.025, 0.975))))## est se 2.5% 97.5%
## psi.projected1 0.6247288 0.06618520 0.5236033 0.7549659
## psi.projected2 0.6029434 0.04888238 0.5222311 0.6968253
## psi.projected3 0.5881543 0.04856831 0.5078085 0.6846648
## psi.projected4 0.5781145 0.05433856 0.4868542 0.6796451
## psi.projected5 0.5712989 0.06051012 0.4641454 0.6808083
## lambda1 0.9651283 0.05335727 0.8840290 1.0824879
## lambda2 0.9754717 0.03798808 0.9092223 1.0503913
## lambda3 0.9829301 0.02785093 0.9307243 1.0314385
## lambda4 0.9882107 0.02088996 0.9447421 1.0191890
## lambda.overall 0.9144751 0.12980332 0.7096046 1.1971346
## lambda.prime1 0.9121745 0.13967722 0.6556645 1.2049649
## lambda.prime2 0.9404430 0.09248807 0.7796846 1.1328267
## lambda.prime3 0.9595389 0.06348076 0.8471215 1.0861330
## lambda.prime4 0.9724999 0.04468537 0.8926752 1.0545864
## lambda.prime.overall 0.8005023 0.30504392 0.3961410 1.5733302# Fit some models.
# Model for ps(.)i, gamma(Year), epsilon(Year), p(year)
# Note that formula DO HAVE AN = SIGN now
mod.psiDot.gYear.eYear.pYear <- unmarked::colext(
psiformula= ~1,
gammaformula = ~ Year,
epsilonformula = ~ Year,
pformula = ~ Year,
data=nso.UMF,
se=TRUE)
summary(mod.psiDot.gYear.eYear.pYear)##
## Call:
## unmarked::colext(psiformula = ~1, gammaformula = ~Year, epsilonformula = ~Year,
## pformula = ~Year, data = nso.UMF, se = TRUE)
##
## Initial (logit-scale):
## Estimate SE z P(>|z|)
## 0.53 0.286 1.86 0.0635
##
## Colonization (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -2.1157 0.769 -2.7495 0.00597
## Year2 -0.4834 1.282 -0.3770 0.70621
## Year3 1.6522 0.890 1.8563 0.06341
## Year4 0.0874 1.142 0.0766 0.93898
##
## Extinction (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -2.332 0.623 -3.741 0.000183
## Year2 0.466 0.856 0.544 0.586631
## Year3 1.174 0.785 1.495 0.135008
## Year4 0.328 0.860 0.381 0.703172
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.361 0.163 2.211 0.027056
## Year2 -0.284 0.229 -1.239 0.215283
## Year3 -0.706 0.246 -2.870 0.004100
## Year4 -0.833 0.245 -3.394 0.000688
## Year5 -0.217 0.227 -0.953 0.340556
##
## AIC: 1355.647
## Number of sites: 55
## optim convergence code: 0
## optim iterations: 42
## Bootstrap iterations: 0slotNames(mod.psiDot.gYear.eYear.pYear)## [1] "phi" "psiformula" "gamformula"
## [4] "epsformula" "detformula" "projected"
## [7] "projected.mean" "smoothed" "smoothed.mean"
## [10] "projected.mean.bsse" "smoothed.mean.bsse" "fitType"
## [13] "call" "formula" "data"
## [16] "sitesRemoved" "estimates" "AIC"
## [19] "opt" "negLogLike" "nllFun"
## [22] "bootstrapSamples" "covMatBS"# We get some estimates using predict in the usual way
names(mod.psiDot.gYear.eYear.pYear)## [1] "psi" "col" "ext" "det"newdata <- data.frame(Year=as.factor(1:5))
predict(mod.psiDot.gYear.eYear.pYear, type='psi', newdata=newdata )## Predicted SE lower upper
## 1 0.6295107 0.06661706 0.492573 0.7483716
## 2 0.6295107 0.06661706 0.492573 0.7483716
## 3 0.6295107 0.06661706 0.492573 0.7483716
## 4 0.6295107 0.06661706 0.492573 0.7483716
## 5 0.6295107 0.06661706 0.492573 0.7483716predict(mod.psiDot.gYear.eYear.pYear, type='det', newdata=newdata )## Predicted SE lower upper
## 1 0.5893367 0.03954536 0.5102411 0.6640684
## 2 0.5193155 0.04027475 0.4405526 0.5971299
## 3 0.4145860 0.04465565 0.3305593 0.5038946
## 4 0.3842423 0.04330240 0.3035834 0.4718126
## 5 0.5360474 0.03935367 0.4586698 0.6117275newdata <- data.frame(Year=as.factor(2:5))
predict(mod.psiDot.gYear.eYear.pYear, type='col', newdata=newdata )## Predicted SE lower upper
## 1 0.10758318 0.07387533 0.02598804 0.3526181
## 2 0.06919589 0.06423449 0.01041749 0.3442486
## 3 0.38616593 0.10594337 0.20760167 0.6016939
## 4 0.11626660 0.08667531 0.02456337 0.4073551predict(mod.psiDot.gYear.eYear.pYear, type='ext', newdata=newdata )## Predicted SE lower upper
## 1 0.0885431 0.05029493 0.02784069 0.2478544
## 2 0.1340054 0.06722738 0.04736179 0.3250679
## 3 0.2391095 0.08686877 0.10974106 0.4447893
## 4 0.1187982 0.06204191 0.04048728 0.3010541# estimates of psi moving ahead.
# There are two types.
# - population level occupancy, i.e. for the sampled population from which the sites were selected
# - actual site level occupancy, i.e. for the sites actually sampled.
# What is the occupancy for the population of ALL potential sites
# first index = not occupied/occupied, second index=year, last index=site
mod.psiDot.gYear.eYear.pYear@projected[,,1]## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.3704893 0.3863698 0.4418643 0.4046869 0.4283574
## [2,] 0.6295107 0.6136302 0.5581357 0.5953131 0.5716426mod.psiDot.gYear.eYear.pYear@projected[2,,1] ## [1] 0.6295107 0.6136302 0.5581357 0.5953131 0.5716426projected(mod.psiDot.gYear.eYear.pYear) # this gives the mean population occupancy over all sites## 1 2 3 4 5
## unoccupied 0.3704893 0.3863698 0.4418643 0.4046869 0.4283574
## occupied 0.6295107 0.6136302 0.5581357 0.5953131 0.5716426# What is the actual occupancy for these particular sites
mod.psiDot.gYear.eYear.pYear@smoothed[2,,1:2]## [,1] [,2]
## [1,] 1 0.04241646
## [2,] 1 0.02170370
## [3,] 1 0.04352670
## [4,] 1 1.00000000
## [5,] 1 1.00000000smoothed(mod.psiDot.gYear.eYear.pYear)[2,] # this gives the mean actual occupancy over all sites## 1 2 3 4 5
## 0.6295096 0.6136176 0.5581249 0.5953073 0.5716324# Obtain standard errors for these forecasts using parametric bootstrapping
psi.projected <- function(model, site) {
logit <- function(x) log(x/(1-x))
expit <- function(x) 1/(1+exp(-x))
psi.projected <- model@projected[2,,site]
# compute population growth
lambda <- exp(diff(log(psi.projected),1))
lambda.overall <- prod(lambda) # overall growth rate over entire set of seasons
lambda.prime <- exp(diff(logit(psi.projected),1))
lambda.prime.overall <- prod(lambda.prime) # overall growth rate on logit scale over entire set of season
c(psi.projected=psi.projected,
lambda= lambda,
lambda.overall=lambda.overall,
lambda.prime =lambda.prime,
lambda.prime.overall=lambda.prime.overall )
}
psi.projected(mod.psiDot.gYear.eYear.pYear, site=1)## psi.projected1 psi.projected2 psi.projected3
## 0.6295107 0.6136302 0.5581357
## psi.projected4 psi.projected5 lambda1
## 0.5953131 0.5716426 0.9747734
## lambda2 lambda3 lambda4
## 0.9095635 1.0666100 0.9602385
## lambda.overall lambda.prime1 lambda.prime2
## 0.9080745 0.9347087 0.7953298
## lambda.prime3 lambda.prime4 lambda.prime.overall
## 1.1645963 0.9071768 0.7853999# use bootstrapping to find standard errors.
# You likely want to do 1000 bootstrap samples
mod.psiDot.gYear.eYear.pYear.boot <- parboot(mod.psiDot.gYear.eYear.pYear, statistic=psi.projected, nsim=100, site=1)
slotNames(mod.psiDot.gYear.eYear.pYear.boot)## [1] "call" "t0" "t.star"mod.psiDot.gYear.eYear.pYear.boot@t0## psi.projected1 psi.projected2 psi.projected3
## 0.6295107 0.6136302 0.5581357
## psi.projected4 psi.projected5 lambda1
## 0.5953131 0.5716426 0.9747734
## lambda2 lambda3 lambda4
## 0.9095635 1.0666100 0.9602385
## lambda.overall lambda.prime1 lambda.prime2
## 0.9080745 0.9347087 0.7953298
## lambda.prime3 lambda.prime4 lambda.prime.overall
## 1.1645963 0.9071768 0.7853999cbind(est=mod.psiDot.gYear.eYear.pYear.boot@t0,
se=apply(mod.psiDot.gYear.eYear.pYear.boot@t.star, 2, sd),
t(apply(mod.psiDot.gYear.eYear.pYear.boot@t.star, 2, quantile, probs=c(0.025, 0.975))))## est se 2.5% 97.5%
## psi.projected1 0.6295107 0.07113092 0.4871646 0.7666727
## psi.projected2 0.6136302 0.07251127 0.4930900 0.7659561
## psi.projected3 0.5581357 0.08018531 0.4225859 0.7252539
## psi.projected4 0.5953131 0.07054091 0.4736349 0.7277024
## psi.projected5 0.5716426 0.06876007 0.4347291 0.7090547
## lambda1 0.9747734 0.07719613 0.8495519 1.1149371
## lambda2 0.9095635 0.09353434 0.7393862 1.0775197
## lambda3 1.0666100 0.16935440 0.8172481 1.4310365
## lambda4 0.9602385 0.09558229 0.7697383 1.1408923
## lambda.overall 0.9080745 0.13743002 0.6834549 1.1452100
## lambda.prime1 0.9347087 0.20719415 0.6568265 1.4201007
## lambda.prime2 0.7953298 0.19843372 0.4621719 1.1892182
## lambda.prime3 1.1645963 0.48875323 0.5577297 2.2408422
## lambda.prime4 0.9071768 0.21569521 0.5313751 1.3156314
## lambda.prime.overall 0.7853999 0.30720095 0.3648621 1.4070747# Model averaging
models.u <-unmarked::fitList(
mod.psiDot.gDot.eDot.pDot,
mod.psiDot.gDot.eDot.pYear,
mod.psiDot.gYear.eYear.pYear)## Warning in unmarked::fitList(mod.psiDot.gDot.eDot.pDot,
## mod.psiDot.gDot.eDot.pYear, : Your list was unnamed, so model names were
## added as object namesaic.table.u <- unmarked::modSel(models.u)
aic.table.u## nPars AIC delta AICwt cumltvWt
## mod.psiDot.gDot.eDot.pYear 8 1353.53 0.00 0.7385 0.74
## mod.psiDot.gYear.eYear.pYear 14 1355.65 2.12 0.2560 0.99
## mod.psiDot.gDot.eDot.pDot 4 1363.32 9.79 0.0055 1.00# Get model averaged estimates of occupancy
predict(models.u, type="psi")[1,]## Predicted SE lower upper
## 1 0.6259883 0.0668573 0.4889115 0.7454423# Get model averaged estimates of detection.
# Notice these are one big list of nsites x nvisits long
# with the detection probabilities for each site
cbind(obs.covar, predict(models.u, type="det"))[1:10,]## SiteNum Year Visit Predicted SE lower upper
## 1 1 1 1 0.5897257 0.03978154 0.5110183 0.6640653
## 2 1 1 2 0.5897257 0.03978154 0.5110183 0.6640653
## 3 1 1 3 0.5897257 0.03978154 0.5110183 0.6640653
## 4 1 1 4 0.5897257 0.03978154 0.5110183 0.6640653
## 5 1 1 5 0.5897257 0.03978154 0.5110183 0.6640653
## 6 1 1 6 0.5897257 0.03978154 0.5110183 0.6640653
## 7 1 1 7 0.5897257 0.03978154 0.5110183 0.6640653
## 8 1 1 8 0.5897257 0.03978154 0.5110183 0.6640653
## 9 1 2 1 0.5215072 0.04040145 0.4426389 0.5993100
## 10 1 2 2 0.5215072 0.04040145 0.4426389 0.5993100# Not possible to model average projected psi's directly.
# You need to do this on your own (groan)