# Multi state single season
# Coosa Bass collected via electrofishing from four 50-80m sections in streams
# at 54 sites in the Upper Coosa River basin.
# Objective was to evaluate the effect of stream flow variability on Coosa
# Bass reproduction. States are:
# Species not detected (state = 0)
# Adult Detected (state = 1)
# YOY present (state = 2)
# 2 covariates: stream link magnitude (standardized), and CV of stream flow
# during summer.
# Analysis in RPresence
# 2018-11-27 code submitted by Neil Faught
# General model averaging framework
library(ggplot2)
library(plyr)
library(readxl)
library(reshape2)
library(RPresence)
## Warning: package 'RPresence' was built under R version 3.5.0
# get the data
input.data <- readxl::read_excel(file.path("..","CoosaBass.xls"),
sheet="Sheet1",
skip=13, na='.')
head(input.data)
## # A tibble: 6 x 8
## Site `1` `2` `3` `4` X__1 Link Flow
## <dbl> <dbl> <dbl> <dbl> <dbl> <lgl> <dbl> <dbl>
## 1 1.00 1.00 1.00 1.00 1.00 NA 0.0500 2.73
## 2 2.00 1.00 1.00 1.00 2.00 NA 1.97 0.800
## 3 3.00 0 0 0 0 NA -0.980 1.40
## 4 4.00 2.00 2.00 0 0 NA 1.12 1.94
## 5 5.00 0 1.00 1.00 2.00 NA 0.410 1.85
## 6 6.00 0 1.00 1.00 1.00 NA 1.85 0.510
input.history <- input.data[,2:5]
input.history[1:5,]
## # A tibble: 5 x 4
## `1` `2` `3` `4`
## <dbl> <dbl> <dbl> <dbl>
## 1 1.00 1.00 1.00 1.00
## 2 1.00 1.00 1.00 2.00
## 3 0 0 0 0
## 4 2.00 2.00 0 0
## 5 0 1.00 1.00 2.00
nrow(input.history)
## [1] 54
Nsites <- nrow(input.history)
Nvisits<- ncol(input.history)
# Extract site level covariates from input data
unit.cov = input.data[,7:8]
names(unit.cov) = c("L","CV")
head(unit.cov)
## # A tibble: 6 x 2
## L CV
## <dbl> <dbl>
## 1 0.0500 2.73
## 2 1.97 0.800
## 3 -0.980 1.40
## 4 1.12 1.94
## 5 0.410 1.85
## 6 1.85 0.510
# create the pao file
bass.pao <- createPao(input.history,
unitcov = unit.cov,
title="bass multi state single season")
summary(bass.pao)
## paoname=pres.pao
## title=bass multi state single season
## Naive occ=0.7777778
## naiveR =0.4444444
## nunits nsurveys nseasons nsurveyseason nmethods
## "54" "4" "1" "4" "1"
## nunitcov nsurvcov
## "2" "1"
## unit covariates : L CV
## survey covariates: 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
~L, ~L, ~L, ~L
~CV, ~CV, ~CV, ~CV")
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 ~L ~L ~L ~L
## 3 ~CV ~CV ~CV ~CV
# 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=bass.pao)
##
##
## ***** Starting ~1 ~1 ~1 ~1
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~L ~L ~L ~L
## PRESENCE Version 2.12.21.
##
##
## ***** Starting ~CV ~CV ~CV ~CV
## PRESENCE Version 2.12.21.
## Warning in matrix(as.numeric(unlist(strsplit(gsub(" +", ",", gsub(".+:
## +", : NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
# Look the output from a specific model
check.model <- 3
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 1.321843 0.712218
## A2_psi0.psi0.CV -0.045662 0.424532
##
## $psi.VC
## A1 A2
## A1 0.507255 -0.268446
## A2 -0.268446 0.180227
##
## $r
## est se
## A3_R0.R0.int 4837.563 NA
## A4_R0.R0.CV -2292.491 NA
##
## $r.VC
## A3 A4
## A3 -293.5891 160.33726
## A4 160.3373 -66.73733
##
## $delta
## est se
## C1_delta(1-1).delta.int -0.611118 0.226437
## C2_delta(1-1).delta.CV 0.064839 0.220730
##
## $delta.VC
## C1 C2
## C1 0.051274 -0.022291
## C2 -0.022291 0.048722
##
## $p
## est se
## B1_p1(1-1).p1.int 1.558546 0.244315
## B2_p1(1-1).p1.CV -0.250366 0.214747
##
## $p.VC
## B1 B2
## B1 0.059690 -0.031473
## B2 -0.031473 0.046116
##
## $VC
## A1 A2 A3 A4 B1 B2
## A1 0.507255 -0.268446 0.000000 0.000000 -0.000632 0.000509
## A2 -0.268446 0.180227 0.000000 0.000000 0.000109 -0.000066
## A3 0.000000 0.000000 -293.589056 160.337264 0.000000 0.000000
## A4 0.000000 0.000000 160.337264 -66.737328 0.000000 0.000000
## B1 -0.000632 0.000109 0.000000 0.000000 0.059690 -0.031473
## B2 0.000509 -0.000066 0.000000 0.000000 -0.031473 0.046116
## C1 0.000000 0.000000 -0.006184 -0.002705 0.000000 0.000000
## C2 0.000000 0.000000 -0.033462 -0.014638 0.000000 0.000000
## C1 C2
## A1 0.000000 0.000000
## A2 0.000000 0.000000
## A3 -0.006184 -0.033462
## A4 -0.002705 -0.014638
## B1 0.000000 0.000000
## B2 0.000000 0.000000
## C1 0.051274 -0.022291
## C2 -0.022291 0.048722
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.7680 0.1105 0.4954 0.9178
## unit2 0.7834 0.0746 0.6045 0.8953
## unit3 0.7787 0.0569 0.6482 0.8704
## unit4 0.7744 0.0663 0.6200 0.8784
## unit5 0.7751 0.0631 0.6291 0.8750
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 0 0 0
## unit2 1 NA NA NA
## unit3 1 0 1 1
## unit4 1 0 1 1
## unit5 1 0 1 1
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.3525 0.0506 0.2607 0.4568
## delta(1-1)2unit2 0.3815 0.0921 0.2229 0.5701
## delta(1-1)3unit3 0.3375 0.0842 0.1958 0.5158
## delta(1-1)4unit4 0.3685 0.0582 0.2634 0.4878
## delta(1-1)5unit5 0.3579 0.0466 0.2724 0.4534
model.fits[[check.model]]$real$p1[1:5,]
## est se lower upper
## p1(1-1)1unit1 0.8243 0.0345 0.7464 0.8821
## p1(1-1)2unit2 0.7437 0.0645 0.5991 0.8493
## p1(1-1)3unit3 0.8586 0.0494 0.7323 0.9310
## p1(1-1)4unit4 0.7821 0.0370 0.7012 0.8460
## p1(1-1)5unit5 0.8109 0.0313 0.7419 0.8648
model.fits[[check.model]]$real$p2[1:5,]
## est se lower upper
## p2(1-1)1unit1 0.8243 0.0345 0.7464 0.8821
## p2(1-1)2unit2 0.7437 0.0645 0.5991 0.8493
## p2(1-1)3unit3 0.8586 0.0494 0.7323 0.9310
## p2(1-1)4unit4 0.7821 0.0370 0.7012 0.8460
## p2(1-1)5unit5 0.8109 0.0313 0.7419 0.8648
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
## 2 psi(L)r(L)p(L)delta(L) 360.7405 344.7405 8 0 0
## 3 psi(CV)r(CV)p(CV)delta(CV) 376.4100 360.4100 8 0 1
## 1 psi()r()p()delta() 389.1956 381.1956 4 0 0
## DAIC modlike wgt
## 2 0.0000 1e+00 0.9996
## 3 15.6695 4e-04 0.0004
## 1 28.4551 0e+00 0.0000
names(aic.table)
## [1] "table" "models" "ess"
# No need to do model averaging as nearly 100% of AIC weight is in one model