# Multi state single season
# Nichols et al (2007) Ecology 88:1395-1400 looked at breeding and
# non-breeding California Spotted Owls in s = 54 sites over k = 5
# surveys.
# Hand-fed mice con rmed occupancy; breeding status only
# confrmed when owl took mouse to nest.
# . = site not visited.
# 0 = owl not detected.
# 1 = owl detected, but no detection of young.
# 2 = owl detected, along with evidence of young.
# Territories that were expected to be occupied were selectively monitored
# and so psi is expected to be close to 1 and not really of interest.
# See MacKenzie et al, Section 5.7.1 for more details
# 2018-10-01 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
# General model averaging framework
library(ggplot2)
library(plyr)
library(readxl)
library(reshape2)
library(RPresence)
# get the data
occ.data <- readxl::read_excel(file.path("..","CalSpottedOwl.xls"),
sheet="RawData",
skip=10, na='.')
head(occ.data)
## # A tibble: 6 x 6
## Site `1` `2` `3` `4` `5`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 NA 2 2 1
## 2 2 0 1 1 NA NA
## 3 3 1 NA 1 NA NA
## 4 4 1 NA 1 NA NA
## 5 5 NA 1 2 1 2
## 6 6 1 1 1 1 NA
input.history <- occ.data[,-1]
input.history[1:5,]
## # A tibble: 5 x 5
## `1` `2` `3` `4` `5`
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 NA 2 2 1
## 2 0 1 1 NA NA
## 3 1 NA 1 NA NA
## 4 1 NA 1 NA NA
## 5 NA 1 2 1 2
Nsites <- nrow(input.history)
Nvisits<- ncol(input.history)
# For biological reasons, we believe that errors in identifying states
# varies between the first 2 visits and the last 3 visit (early bs late
# breeding season). We create a site-time covariates.
# The format is the "long" format ordered by visit major order,i.e.
# stack the columns of the matrix of period covariates
#
# The covariate matrix would be arranged originally as
# 1 2 3 4 5
# site1 [E E L L L]
# site2 [E E L L L]
# : [: : : : :]
# siteN [E E L L L]
#
# and this is stacked by columns
period <- data.frame(Site=rep(1:Nsites, Nvisits),
Visit=rep(1:Nvisits, each=Nsites),
Per =rep(c("E","E","L","L","L"), each=Nsites))
period[1:10,]
## Site Visit Per
## 1 1 1 E
## 2 2 1 E
## 3 3 1 E
## 4 4 1 E
## 5 5 1 E
## 6 6 1 E
## 7 7 1 E
## 8 8 1 E
## 9 9 1 E
## 10 10 1 E
# create the pao file
owl.pao <- createPao(input.history,
survcov=period,
title="owl multi state single season")
summary(owl.pao)
## paoname=pres.pao
## title=owl multi state single season
## Naive occ=0.8703704
## naiveR =0.3518519
## nunits nsurveys nseasons nsurveyseason nmethods
## "54" "5" "1" "5" "1"
## nunitcov nsurvcov
## "1" "4"
## unit covariates : TEMP
## survey covariates: Site Visit Per SURVEY
# Use the multi-state model, but only for a single season
# define the list of models to fit
# Notice the commas between the column and the placement of the quotes
# psi occupancy probability,
# r probability of being in the second state, conditional on the unit being occupied.
# p detection probability. .
# delta probability of detecting the second state in a survey, conditional on the species being detected in the survey.
model.list.csv <- textConnection("
psi, r, p, delta
~1, ~1, ~1, ~1
~1, ~1, ~SURVEY, ~1
~1, ~1, ~STATE, ~1
~1, ~1, ~1, ~Per
~1, ~1, ~SURVEY,~Per
~1, ~1, ~STATE, ~Per")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## psi r p delta
## 1 ~1 ~1 ~1 ~1
## 2 ~1 ~1 ~SURVEY ~1
## 3 ~1 ~1 ~STATE ~1
## 4 ~1 ~1 ~1 ~Per
## 5 ~1 ~1 ~SURVEY ~Per
## 6 ~1 ~1 ~STATE ~Per
# 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("r" ,x$r )),
as.formula(paste("p" ,x$p )),
as.formula(paste("delta",x$delta))),
data=detect.pao,type="do.ms.2")
fit
},detect.pao=owl.pao)
##
##
## ***** Starting ~1 ~1 ~1 ~1
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~1 ~1 ~SURVEY ~1
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~1 ~1 ~STATE ~1
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~1 ~1 ~1 ~Per
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 2.62 signifcant digits.
##
##
##
## ***** Starting ~1 ~1 ~SURVEY ~Per
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 0.39 signifcant digits.
##
##
##
## ***** Starting ~1 ~1 ~STATE ~Per
## PRESENCE Version 2.12.21.
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 3.25 signifcant digits.
# 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] "warnings"
model.fits[[check.model]]$beta
## $psi
## est se
## A1_psi0.psi0.int 3.67294 1.563724
##
## $psi.VC
## [1] 2.445234
##
## $r
## est se
## A2_R0.R0.int 0.218015 0.45513
##
## $r.VC
## [1] 0.207143
##
## $delta
## est se
## C1_delta(1-1).delta.int -0.31056 0.345561
##
## $delta.VC
## [1] 0.119413
##
## $p
## est se
## B1_p1(1-1).p1.int 1.071073 0.19103
##
## $p.VC
## [1] 0.036492
##
## $VC
## A1 A2 B1 C1
## A1 2.445234 0.000000 -0.092154 0.000000
## A2 0.000000 0.207143 0.000000 -0.082485
## B1 -0.092154 0.000000 0.036492 0.000000
## C1 0.000000 -0.082485 0.000000 0.119413
names(model.fits[[check.model]]$real)
## [1] "psi" "Cpsi0" "Cpsi1" "Cpsi2" "r" "CR0" "CR1" "CR2"
## [9] "delta" "p1" "p2"
model.fits[[check.model]]$real$psi[1:5,]
## est se lower upper
## unit1 0.9752 0.0378 0.6475 0.9988
## unit2 0.9752 0.0378 0.6475 0.9988
## unit3 0.9752 0.0378 0.6475 0.9988
## unit4 0.9752 0.0378 0.6475 0.9988
## unit5 0.9752 0.0378 0.6475 0.9988
model.fits[[check.model]]$real$Cpsi0[1:5,]
## data frame with 0 columns and 5 rows
model.fits[[check.model]]$real$Cpsi1[1:5,]
## data frame with 0 columns and 5 rows
model.fits[[check.model]]$real$Cpsi2[1:5,]
## data frame with 0 columns and 5 rows
model.fits[[check.model]]$real$r[1:5,]
## est se lower upper
## unit1 0.5543 0.1124 0.3376 0.7521
## unit2 0.5543 0.1124 0.3376 0.7521
## unit3 0.5543 0.1124 0.3376 0.7521
## unit4 0.5543 0.1124 0.3376 0.7521
## unit5 0.5543 0.1124 0.3376 0.7521
model.fits[[check.model]]$real$CR0[1:5,]
## data frame with 0 columns and 5 rows
model.fits[[check.model]]$real$CR1[1:5,]
## data frame with 0 columns and 5 rows
model.fits[[check.model]]$real$CR2[1:5,]
## data frame with 0 columns and 5 rows
model.fits[[check.model]]$real$delta[1:5,]
## est se lower upper
## delta(1-1)1unit1 0.423 0.0843 0.2713 0.5907
## delta(1-1)2unit2 0.423 0.0843 0.2713 0.5907
## delta(1-1)3unit3 0.423 0.0843 0.2713 0.5907
## delta(1-1)4unit4 0.423 0.0843 0.2713 0.5907
## delta(1-1)5unit5 0.423 0.0843 0.2713 0.5907
model.fits[[check.model]]$real$p1[1:5,]
## est se lower upper
## p1(1-1)1unit1 0.7448 0.0363 0.6674 0.8093
## p1(1-1)2unit2 0.7448 0.0363 0.6674 0.8093
## p1(1-1)3unit3 0.7448 0.0363 0.6674 0.8093
## p1(1-1)4unit4 0.7448 0.0363 0.6674 0.8093
## p1(1-1)5unit5 0.7448 0.0363 0.6674 0.8093
model.fits[[check.model]]$real$p2[1:5,]
## est se lower upper
## p2(1-1)1unit1 0.7448 0.0363 0.6674 0.8093
## p2(1-1)2unit2 0.7448 0.0363 0.6674 0.8093
## p2(1-1)3unit3 0.7448 0.0363 0.6674 0.8093
## p2(1-1)4unit4 0.7448 0.0363 0.6674 0.8093
## p2(1-1)5unit5 0.7448 0.0363 0.6674 0.8093
names(model.fits[[check.model]]$derived)
## NULL
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model AIC neg2ll npar warn.conv warn.VC
## 5 psi()r()p(SURVEY)delta(Per) 279.3550 261.3550 9 0.39 1
## 4 psi()r()p()delta(Per) 288.3492 278.3492 5 2.62 1
## 6 psi()r()p(STATE)delta(Per) 290.2400 278.2400 6 3.25 1
## 2 psi()r()p(SURVEY)delta() 326.3843 310.3843 8 0.00 0
## 1 psi()r()p()delta() 335.3784 327.3784 4 0.00 0
## 3 psi()r()p(STATE)delta() 337.3627 327.3627 5 0.00 0
## DAIC modlike wgt
## 5 0.0000 1.0000 0.9848
## 4 8.9942 0.0111 0.0110
## 6 10.8850 0.0043 0.0043
## 2 47.0293 0.0000 0.0000
## 1 56.0234 0.0000 0.0000
## 3 58.0077 0.0000 0.0000
names(aic.table)
## [1] "table" "models" "ess"
RPresence::modAvg(aic.table, param="psi")[1,]
## est se lower_0.95 upper_0.95
## unit1 0.9788445 0.03249207 0.681173 0.999003
RPresence::modAvg(aic.table, param="r") [1,]
## est se lower_0.95 upper_0.95
## unit1 0.4494685 0.08021052 0.3019253 0.6064727
RPresence::modAvg(aic.table, param="p1")[seq(1,by=Nsites, length.out=Nvisits),]
## est se lower_0.95 upper_0.95
## p1(1-1)1unit1 0.5814721 0.08402199 0.4138997 0.7321398
## p1(1-2)1unit1 0.8225428 0.06823416 0.6496395 0.9205532
## p1(1-3)1unit1 0.9196407 0.04828519 0.7607676 0.9762945
## p1(1-4)1unit1 0.6976745 0.07687768 0.5304329 0.8250024
## p1(1-5)1unit1 0.6214536 0.11255304 0.3912629 0.8074386
RPresence::modAvg(aic.table, param="p2")[seq(1,by=Nsites, length.out=Nvisits),]
## est se lower_0.95 upper_0.95
## p2(1-1)1unit1 0.5815911 0.08426813 0.4135243 0.7326347
## p2(1-2)1unit1 0.8226617 0.06811742 0.6500722 0.9205333
## p2(1-3)1unit1 0.9197596 0.04787949 0.7627009 0.9761219
## p2(1-4)1unit1 0.6977934 0.07696727 0.5303263 0.8252270
## p2(1-5)1unit1 0.6215725 0.11269472 0.3910747 0.8077185
RPresence::modAvg(aic.table, param="delta")[seq(1,by=Nsites, length.out=Nvisits),]
## est se lower_0.95 upper_0.95
## delta(1-1)1unit1 2.593953e-11 NaN NaN NaN
## delta(1-2)1unit1 2.593953e-11 NaN NaN NaN
## delta(1-3)1unit1 8.683107e-01 0.07269637 0.6547659 0.9581997
## delta(1-4)1unit1 8.683107e-01 0.07269637 0.6547659 0.9581997
## delta(1-5)1unit1 8.683107e-01 0.07269637 0.6547659 0.9581997