# multi season multi-state model
# Robust design, Multi State, Reproduction parameterization
# Using RPresence
# Visits to hawk territories to establish occupancy and breeding
# 6 years x up to 9 visits/year
# 2018-11-27 Code contributed by Neil Faught
library(car)
## Loading required package: carData
library(ggplot2)
library(readxl)
library(reshape2)
library(RPresence)
input.data <- read.csv(file.path("..","hawk.csv"),
as.is=TRUE, strip.white=TRUE, header=TRUE)
head(input.data)
## TerritoryN tsh V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16
## 1 1 9 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 2 2 3 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 3 10 NA NA NA NA NA NA NA NA NA 0 2 0 NA NA NA NA
## 4 4 2 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 5 7 NA NA NA NA NA NA NA NA NA 0 NA NA NA NA NA NA
## 6 6 7 0 0 NA NA NA NA NA NA NA 0 0 NA NA NA NA NA
## V17 V18 V19 V20 V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34
## 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 2 NA NA 0 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 NA NA 0 NA NA NA NA NA NA NA NA 0 1 NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## V35 V36 V37 V38 V39 V40 V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52
## 1 NA NA NA NA NA NA NA NA NA NA NA 0 0 NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## V53 V54
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
input.history = input.data[,-c(1,2)]
head(input.history)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20
## 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 0 1
## 3 NA NA NA NA NA NA NA NA NA 0 2 0 NA NA NA NA NA NA 0 NA
## 4 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA 0 NA NA NA NA NA NA NA NA NA NA
## 6 0 0 NA NA NA NA NA NA NA 0 0 NA NA NA NA NA NA NA NA NA
## V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37 V38
## 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA 0 1 NA NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## V39 V40 V41 V42 V43 V44 V45 V46 V47 V48 V49 V50 V51 V52 V53 V54
## 1 NA NA NA NA NA NA NA 0 0 NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# Create site level covariate data frame
unit.cov = input.data[,c(1,2)]
names(unit.cov) = c("Territory","tsh")
head(unit.cov)
## Territory tsh
## 1 1 9
## 2 2 3
## 3 3 10
## 4 4 2
## 5 5 7
## 6 6 7
# make a plot of the imputed occupancy
plotdata <- reshape2::melt(input.data,
id.vars=c("TerritoryN","tsh"),
value.name="ObsOcc",
variable.name="Visit")
plotdata$Visit <- as.numeric(substring(plotdata$Visit,2))
head(plotdata)
## TerritoryN tsh Visit ObsOcc
## 1 1 9 1 NA
## 2 2 3 1 NA
## 3 3 10 1 NA
## 4 4 2 1 1
## 5 5 7 1 NA
## 6 6 7 1 0
ggplot(data=plotdata, aes(x=Visit, y=TerritoryN, color=as.factor(ObsOcc), shape=as.factor(ObsOcc)))+
ggtitle("Observed occupancy state", subtitle="Not corrected for false negatives")+
geom_point()+
scale_shape_manual(name="Observed\nOccupancy", values=c(1, 16, 19))+
scale_color_manual(name="Observed\nOccupancy", values=c("black","green","red"))+
theme_bw()
## Warning: Removed 7275 rows containing missing values (geom_point).
# Number of visits in each season
Nvisits.per.season <- rep(6,9)
# create the pao file
hawk.pao <- createPao(input.history,
nsurveyseason=Nvisits.per.season,
unitcov=unit.cov,
title="hawk multi state multi season")
summary(hawk.pao)
## paoname=pres.pao
## title=hawk multi state multi season
## Naive occ=0.5878378
## naiveR =0.4189189
## nunits nsurveys nseasons
## "148" "54" "9"
## nsurveyseason nmethods nunitcov
## "6,6,6,6,6,6,6,6,6" "1" "2"
## nsurvcov
## "1"
## unit covariates : Territory tsh
## survey covariates: SURVEY
# 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
~DYN + PREV_STATE, ~DYN+ PREV_STATE, ~1, ~1
~tsh*DYN, ~tsh, ~1, ~1
~tsh*DYN, ~1, ~1, ~1")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## psi r p delta
## 1 ~DYN + PREV_STATE ~DYN+ PREV_STATE ~1 ~1
## 2 ~tsh*DYN ~tsh ~1 ~1
## 3 ~tsh*DYN ~1 ~1 ~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("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=hawk.pao)
##
##
## ***** Starting ~DYN + PREV_STATE ~DYN+ PREV_STATE ~1 ~1
## PRESENCE Version 2.12.21.
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 2.46 signifcant digits.
##
##
##
## ***** Starting ~tsh*DYN ~tsh ~1 ~1
## PRESENCE Version 2.12.21.
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
##
##
## ***** Starting ~tsh*DYN ~1 ~1 ~1
## PRESENCE Version 2.12.21.
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
## Warning in while (nchar(v[i + j]) > 10) {: the condition has length > 1 and
## only the first element will be used
##
## Warning: Numerical convergence may not have been reached. Parameter esimates converged to approximately 3.52 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 0.674366 0.384866
## A2_Cpsi0(1) 0.495718 13.120190
## A3_Cpsi0(1) -2.590096 13.107463
## A4_Cpsi1(1) 0.596189 13.182372
## A5_Cpsi2(1) 2.488158 13.214786
##
## $psi.VC
## A1 A2 A3 A4 A5
## A1 0.148122 -0.135525 -0.028498 -0.030227 0.018934
## A2 -0.135525 172.139391 -171.728020 -172.316132 -172.428238
## A3 -0.028498 -171.728020 171.805587 171.671123 172.069860
## A4 -0.030227 -172.316132 171.671123 173.774919 172.518305
## A5 0.018934 -172.428238 172.069860 172.518305 174.630571
##
## $r
## est se
## A6_R0 0.018733 0.457595
## A7_CR0(1) -4.703130 15.401274
## A8_CR0(1) -16.657607 33520.287418
## A9_CR1(1) 4.097648 15.334333
## A10_CR2(1) 7.844406 15.553951
##
## $r.VC
## A6 A7 A8 A9 A10
## A6 0.209393 -0.267493 9.438000e-02 0.00719 0.113369
## A7 -0.267493 237.199241 -2.486312e+02 -235.88360 -239.091765
## A8 0.094380 -248.631242 1.123610e+09 247.01247 250.752325
## A9 0.007190 -235.883605 2.470125e+02 235.14177 237.784263
## A10 0.113369 -239.091765 2.507523e+02 237.78426 241.925390
##
## $delta
## est se
## E1_delta(1-1) 1.70713 0.358556
##
## $delta.VC
## [1] 0.128563
##
## $p
## est se
## E2_p1(1-1) -0.453572 0.13515
##
## $p.VC
## [1] 0.018265
##
## $VC
## A1 A2 A3 A4 A5 A6
## A1 0.148122 -0.135525 -0.028498 -0.030227 0.018934 -0.015136
## A2 -0.135525 172.139391 -171.728020 -172.316132 -172.428238 -0.016641
## A3 -0.028498 -171.728020 171.805587 171.671123 172.069860 0.004444
## A4 -0.030227 -172.316132 171.671123 173.774919 172.518305 0.073908
## A5 0.018934 -172.428238 172.069860 172.518305 174.630571 0.030782
## A6 -0.015136 -0.016641 0.004444 0.073908 0.030782 0.209393
## A7 -0.025520 0.321495 -0.018013 -0.963495 -0.318392 -0.267493
## A8 -0.060364 11.720149 -12.060779 -10.703628 -12.251199 0.094380
## A9 0.010869 -0.143852 0.001670 0.515620 0.017086 0.007190
## A10 0.020507 -0.285535 0.002529 0.976538 0.117647 0.113369
## E1 -0.006755 0.030374 -0.001167 -0.048943 -0.039138 -0.040064
## E2 -0.014535 0.035424 -0.009239 -0.064624 -0.083274 -0.002052
## A7 A8 A9 A10 E1 E2
## A1 -0.025520 -6.036400e-02 0.010869 0.020507 -0.006755 -0.014535
## A2 0.321495 1.172015e+01 -0.143852 -0.285535 0.030374 0.035424
## A3 -0.018013 -1.206078e+01 0.001670 0.002529 -0.001167 -0.009239
## A4 -0.963495 -1.070363e+01 0.515620 0.976538 -0.048943 -0.064624
## A5 -0.318392 -1.225120e+01 0.017086 0.117647 -0.039138 -0.083274
## A6 -0.267493 9.438000e-02 0.007190 0.113369 -0.040064 -0.002052
## A7 237.199241 -2.486312e+02 -235.883605 -239.091765 0.267922 0.048676
## A8 -248.631242 1.123610e+09 247.012468 250.752325 -0.350027 -0.017453
## A9 -235.883605 2.470125e+02 235.141774 237.784263 -0.215691 -0.026061
## A10 -239.091765 2.507523e+02 237.784263 241.925390 -0.346475 -0.047978
## E1 0.267922 -3.500270e-01 -0.215691 -0.346475 0.128563 0.004449
## E2 0.048676 -1.745300e-02 -0.026061 -0.047978 0.004449 0.018265
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.6625 0.0861 0.48 0.8067
## unit2 0.6625 0.0861 0.48 0.8067
## unit3 0.6625 0.0861 0.48 0.8067
## unit4 0.6625 0.0861 0.48 0.8067
## unit5 0.6625 0.0861 0.48 0.8067
model.fits[[check.model]]$real$Cpsi0[1:5,]
## est se lower upper
## unit1_1 0.1947 0.0871 0.0752 0.4181
## unit2_1 0.1947 0.0871 0.0752 0.4181
## unit3_1 0.1947 0.0871 0.0752 0.4181
## unit4_1 0.1947 0.0871 0.0752 0.4181
## unit5_1 0.1947 0.0871 0.0752 0.4181
model.fits[[check.model]]$real$Cpsi1[1:5,]
## est se lower upper
## unit1_1 0.854 0.1307 0.4285 0.9786
## unit2_1 0.854 0.1307 0.4285 0.9786
## unit3_1 0.854 0.1307 0.4285 0.9786
## unit4_1 0.854 0.1307 0.4285 0.9786
## unit5_1 0.854 0.1307 0.4285 0.9786
model.fits[[check.model]]$real$Cpsi2[1:5,]
## est se lower upper
## unit1_1 0.9749 0.0331 0.7326 0.9982
## unit2_1 0.9749 0.0331 0.7326 0.9982
## unit3_1 0.9749 0.0331 0.7326 0.9982
## unit4_1 0.9749 0.0331 0.7326 0.9982
## unit5_1 0.9749 0.0331 0.7326 0.9982
model.fits[[check.model]]$real$r[1:5,]
## est se lower upper
## unit1 0.5047 0.1144 0.2936 0.7141
## unit2 0.5047 0.1144 0.2936 0.7141
## unit3 0.5047 0.1144 0.2936 0.7141
## unit4 0.5047 0.1144 0.2936 0.7141
## unit5 0.5047 0.1144 0.2936 0.7141
model.fits[[check.model]]$real$CR0[1:5,]
## est se lower upper
## unit1_1 0 0 0 1
## unit2_1 0 0 0 1
## unit3_1 0 0 0 1
## unit4_1 0 0 0 1
## unit5_1 0 0 0 1
model.fits[[check.model]]$real$CR1[1:5,]
## est se lower upper
## unit1_1 0.3574 0.1177 0.1692 0.6029
## unit2_1 0.3574 0.1177 0.1692 0.6029
## unit3_1 0.3574 0.1177 0.1692 0.6029
## unit4_1 0.3574 0.1177 0.1692 0.6029
## unit5_1 0.3574 0.1177 0.1692 0.6029
model.fits[[check.model]]$real$CR2[1:5,]
## est se lower upper
## unit1_1 0.9593 0.0358 0.796 0.993
## unit2_1 0.9593 0.0358 0.796 0.993
## unit3_1 0.9593 0.0358 0.796 0.993
## unit4_1 0.9593 0.0358 0.796 0.993
## unit5_1 0.9593 0.0358 0.796 0.993
model.fits[[check.model]]$real$delta[1:5,]
## est se lower upper
## delta(1-1)1unit1 0.8465 0.0466 0.7319 0.9176
## delta(1-2)1unit1 0.8465 0.0466 0.7319 0.9176
## delta(1-3)1unit1 0.8465 0.0466 0.7319 0.9176
## delta(1-4)1unit1 0.8465 0.0466 0.7319 0.9176
## delta(1-5)1unit1 0.8465 0.0466 0.7319 0.9176
model.fits[[check.model]]$real$p1[1:5,]
## est se lower upper
## p1(1-1)1unit1 0.3885 0.0321 0.3277 0.453
## p1(1-2)1unit1 0.3885 0.0321 0.3277 0.453
## p1(1-3)1unit1 0.3885 0.0321 0.3277 0.453
## p1(1-4)1unit1 0.3885 0.0321 0.3277 0.453
## p1(1-5)1unit1 0.3885 0.0321 0.3277 0.453
model.fits[[check.model]]$real$p2[1:5,]
## est se lower upper
## p2(1-1)1unit1 0.3885 0.0321 0.3277 0.453
## p2(1-2)1unit1 0.3885 0.0321 0.3277 0.453
## p2(1-3)1unit1 0.3885 0.0321 0.3277 0.453
## p2(1-4)1unit1 0.3885 0.0321 0.3277 0.453
## p2(1-5)1unit1 0.3885 0.0321 0.3277 0.453
names(model.fits[[check.model]]$derived)
## NULL
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model AIC neg2ll
## 1 psi(DYN P PREV_STATE)r(DYN P PREV_STATE)p()delta() 1095.602 1071.602
## 2 psi(tsh X DYN)r(tsh)p()delta() 1118.896 1102.896
## 3 psi(tsh X DYN)r()p()delta() 1123.727 1109.727
## npar warn.conv warn.VC DAIC modlike wgt
## 1 12 2.46 0 0.0000 1 1
## 2 8 0.00 0 23.2934 0 0
## 3 7 3.52 1 28.1245 0 0
names(aic.table)
## [1] "table" "models" "ess"