# Analysis of Hawaii Elepaio nesting data using the logistic exposure model
# 2019-05-01 CJS Initial code
# Paul C Banko, Kelly A Jaenecke, Robert W Peck, Kevin W Brinck,
# Increased nesting success of Hawaii Elepaio in response to the removal
# of invasive black rats,
# The Condor: Ornithological Applications, Volume 121, Issue 1, 1 February 2019, duz003,
# https://doi.org/10.1093/condor/duz003
# Data from
# Banko, P.C., Jaenecke, K.A., Peck, R.W., and
# Brinck, K.W., 2018, Hawaii Volcanoes National Park Elepaio
# nest monitoring and black rat mark recapture data 2015–2017:
# US Geological Survey data release, https://doi.org/10.5066/P93TOM58
# In Hawai‘i and other oceanic islands with few native land mammals,
# black rats (Rattus rattus) are among the most damaging invasive
# vertebrate species to native forest bird populations and habitats,
# due to their arboreal behavior and generalist foraging habitats and habitat use.
# We evaluated the nesting response of Hawai‘i ‘Elepaio (Chasiempis sandwichensis; Monarchidae),
# a generalist insectivore, to the removal of black rats using rodenticide
# in a before-after-control-impact study in high and low, mesic montane habitat
# recovering from long-term damage from introduced ungulates and weeds.
# We monitored nesting and rat activity during 2015–2016 before applying diphacinone
# bait in 2017 to remove rats from two 700 x 700-m treatment plots that were paired
# with two non-treatment plots. We continued monitoring through July 2017.
# This data release consists of three tabular datasets including nest monitoring data,
# rat capture data, and spatial data for all field plots.
library(AICcmodavg)
library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
library(plyr)
library(readxl)
source(file.path("..","..", "logistic-exposure-model.R"))
# The dataframe must contain the following fields with the following names
#
# FirstFound: day the nest was first found
# LastPresent: last day that a chick was present in the nest
# LastChecked: last day the nest was checked
# Fate: fate of the nest; 0=hatch an
# Freq: number of nests with this data
#
# In this example, the multiple visits to a nest have been collapsed
# to a single record for each nest.
# In more complex examples, you may have multple records per nest
# as shown in the mallard example.
#
haeldata <- read.csv(file.path("..","HAEL_NestMonitoring_HAVO_2015-2017.csv"), header=TRUE, as.is=TRUE, strip.white=TRUE)
head(haeldata)
## Nest.ID year site Fate FirstFoundCD LastPresentCD LastCheckedCD AgeFound
## 1 2 2017 KKI 0 40 69 69 5
## 2 3 2017 KKI 0 47 82 82 1
## 3 11 2017 LKE 0 59 93 93 1
## 4 15 2017 LKE 0 66 99 99 2
## 5 16 2017 LKE 0 69 100 100 4
## 6 17 2017 LKE 0 73 105 105 3
## tmt elev FirstFound LastPresent LastChecked AgeDay1 cause
## 1 impact low 1 30 30 5 fledge
## 2 impact low 8 43 43 -6 fledge
## 3 impact high 20 54 54 -18 fledge
## 4 impact high 27 60 60 -24 fledge
## 5 impact high 30 61 61 -25 fledge
## 6 impact high 34 66 66 -30 fledge
# We expand the data to generate the effective sample size
haeldata2 <- expand.nest.data(haeldata)
head(haeldata2)
## Nest.ID year site Fate FirstFoundCD LastPresentCD LastCheckedCD AgeFound
## 1 2 2017 KKI 0 40 69 69 5
## 2 2 2017 KKI 0 40 69 69 5
## 3 2 2017 KKI 0 40 69 69 5
## 4 2 2017 KKI 0 40 69 69 5
## 5 2 2017 KKI 0 40 69 69 5
## 6 2 2017 KKI 0 40 69 69 5
## tmt elev FirstFound LastPresent LastChecked AgeDay1 cause Day
## 1 impact low 1 30 30 5 fledge 1
## 2 impact low 1 30 30 5 fledge 2
## 3 impact low 1 30 30 5 fledge 3
## 4 impact low 1 30 30 5 fledge 4
## 5 impact low 1 30 30 5 fledge 5
## 6 impact low 1 30 30 5 fledge 6
## Exposure Survive NestAge
## 1 1 1 5
## 2 1 1 6
## 3 1 1 7
## 4 1 1 8
## 5 1 1 9
## 6 1 1 10
rm(haeldata)
# create factor variable for year
xtabs(~year, data=haeldata2, exclude=NULL, na.action=na.pass)
## year
## 2015 2016 2017
## 985 879 1653
haeldata2$yearF <- factor(haeldata2$year)
# Set up the set of model to fit
model.list.csv <- textConnection(
" S
~NestAge + tmt
~ tmt
~ tmt + elev
")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list$model.number <- 1:nrow(model.list)
model.list
## S model.number
## 1 ~NestAge + tmt 1
## 2 ~ tmt 2
## 3 ~ tmt + elev 3
model.fits <- plyr::dlply(model.list, c("model.number","S"), function(x,input.data, input.ddl){
cat("\n\n***** Starting ", unlist(x), "\n")
fit <- glm(formula=as.formula(paste("Survive", eval(x$S))),
family=binomial(link=logexp(input.data$Exposure)),
data=input.data)
fit
},input.data=haeldata2)
##
##
## ***** Starting ~NestAge + tmt 1
##
##
## ***** Starting ~ tmt 2
##
##
## ***** Starting ~ tmt + elev 3
# examine results for a particular model
model.number <-3
names(model.fits[[model.number]])
## [1] "coefficients" "residuals" "fitted.values"
## [4] "effects" "R" "rank"
## [7] "qr" "family" "linear.predictors"
## [10] "deviance" "aic" "null.deviance"
## [13] "iter" "weights" "prior.weights"
## [16] "df.residual" "df.null" "y"
## [19] "converged" "boundary" "model"
## [22] "call" "formula" "terms"
## [25] "data" "offset" "control"
## [28] "method" "contrasts" "xlevels"
model.fits[[model.number]]$formula
## Survive ~ tmt + elev
## <environment: 0x7fc4a1cda0e8>
summary(model.fits[[model.number]])
##
## Call:
## glm(formula = as.formula(paste("Survive", eval(x$S))), family = binomial(link = logexp(input.data$Exposure)),
## data = input.data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8274 0.2110 0.2588 0.2834 0.2834
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3796 0.1415 23.889 <2e-16 ***
## tmtimpact 0.5991 0.2379 2.519 0.0118 *
## elevlow -0.1851 0.1883 -0.983 0.3257
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 999.31 on 3516 degrees of freedom
## Residual deviance: 745.04 on 3514 degrees of freedom
## AIC: 751.04
##
## Number of Fisher Scoring iterations: 7
# Model comparision and averaging
# collect models and make AICc table
AICcmodavg::aictab(model.fits)
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## 1.~NestAge + tmt 3 746.14 0.00 0.81 0.81 -370.06
## 2.~ tmt 2 749.97 3.83 0.12 0.93 -372.98
## 3.~ tmt + elev 3 751.04 4.91 0.07 1.00 -372.52
# results differ slight from paper because of approximation used for
# NestAge in final intervals where the mid-point is used, but the individual
# ages are used in MARK during likelihood constructions. All of the results
# for the other models match perfectly.