# Analysis of Redstart data illustrating logistic exposure models
# Reproduce Table 1b.
# 2019-05-01 CJS Initial code
# Sherry TW, Wilson S, Hunter S, Holmes RT (2015)
# Impacts of nest predators and weather on reproductive success and
# population limitation in a long-distance migratory songbird.
# Journal of Avian Biology 46(6): 559-569. https://doi.org/10.1111/jav.00536
# Data from
# Sherry TW, Wilson S, Hunter S, Holmes RT (2015)
# Data from: Impacts of nest predators and weather on reproductive success and
# population limitation in a long-distance migratory songbird.
# Dryad Digital Repository. https://doi.org/10.5061/dryad.73870
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.
#
reddata <- readxl::read_excel(file.path("..","Sherry.xlsx"),
sheet="NestData")
head(reddata)
## # A tibble: 6 x 11
## NestId FirstFound LastPresent LastChecked Fate Freq Height BaffleStatus
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 /*198… 11 11 18 1 1 16.9 N
## 2 /*198… 13 22 23 1 1 10.5 N
## 3 /*198… 7 26 26 0 1 4.8 N
## 4 /*198… 12 14 20 1 1 3.7 N
## 5 /*198… 6 20 25 1 1 10.7 N
## 6 /*198… 11 19 20 1 1 17.8 N
## # … with 3 more variables: DBH <dbl>, AgeDay1 <dbl>, Year <chr>
reddata <- as.data.frame(reddata)
# get the yearly date and merge with the nest data
annual.data <- readxl::read_excel(file.path("..","Sherry.xlsx"),
sheet="AnnualCovariates", skip=9)
head(annual.data)
## # A tibble: 6 x 7
## Year Density Predators MayTemp JuneTemp MayRain JuneRain
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1983 0.867 2.4 8.89 17.1 194. 35.3
## 2 1984 0.467 0.5 9.06 16.2 312. 180.
## 3 1985 0.589 3.4 10.6 13.2 107. 130
## 4 1986 0.439 6 12.2 13.4 128. 89.6
## 5 1987 0.372 1.7 10.9 15.4 91.6 200.
## 6 1988 0.361 0.3 12 14.4 97.8 70.2
dim(reddata)
## [1] 537 11
reddata <- merge(reddata, annual.data, all.x=TRUE)
dim(reddata)
## [1] 537 17
# any missing data?
sum(!complete.cases(reddata))
## [1] 0
# as noted in the paper, we remove nests with baffeling
dim(reddata)
## [1] 537 17
reddata <- reddata[ !(reddata$BaffleStatus=="Y"),]
dim(reddata)
## [1] 466 17
# expand the data
reddata2 <- expand.nest.data(reddata)
rm(reddata)
# create factor variable for year
reddata2$YearF <- factor(reddata2$Year)
# Set up the set of model to fit
model.list.csv <- textConnection(
" S
~DBH
~DBH+Predators+ MayTemp +JuneRain +NestAge
~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge
~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge
~DBH+Predators+ MayTemp +JuneRain +NestAge + Density
~DBH+Predators+ MayTemp +JuneRain
~DBH+Predators+ MayTemp +MayRain +NestAge
~DBH+Predators+ MayTemp +NestAge
~DBH+Predators+ MayTemp +MayRain +NestAge + Density
")
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
## 1 ~DBH
## 2 ~DBH+Predators+ MayTemp +JuneRain +NestAge
## 3 ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge
## 4 ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge
## 5 ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density
## 6 ~DBH+Predators+ MayTemp +JuneRain
## 7 ~DBH+Predators+ MayTemp +MayRain +NestAge
## 8 ~DBH+Predators+ MayTemp +NestAge
## 9 ~DBH+Predators+ MayTemp +MayRain +NestAge + Density
## model.number
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
model.fits <- plyr::dlply(model.list, c("S","model.number"), 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=reddata2)
##
##
## ***** Starting ~DBH 1
##
##
## ***** Starting ~DBH+Predators+ MayTemp +NestAge 8
##
##
## ***** Starting ~DBH+Predators+ MayTemp +JuneRain 6
##
##
## ***** Starting ~DBH+Predators+ MayTemp +JuneRain +NestAge 2
##
##
## ***** Starting ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density 5
##
##
## ***** Starting ~DBH+Predators+ MayTemp +MayRain +NestAge 7
##
##
## ***** Starting ~DBH+Predators+ MayTemp +MayRain +NestAge + Density 9
##
##
## ***** Starting ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge 3
##
##
## ***** Starting ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge 4
# examine individual model results
model.number <-1
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 ~ DBH
## <environment: 0x7fbf3bc6fdc8>
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
## -3.0828 0.2319 0.2908 0.3050 0.3192
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.930362 0.098788 29.663 < 2e-16 ***
## DBH 0.021201 0.005456 3.886 0.000102 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1760.6 on 5373 degrees of freedom
## Residual deviance: 1442.7 on 5372 degrees of freedom
## AIC: 1446.7
##
## Number of Fisher Scoring iterations: 6
# Model comparision and averaging
# collect models and make AICc table
AICcmodavg::aictab(model.fits)
##
## Model selection based on AICc:
##
## K
## ~DBH+Predators+ MayTemp +JuneRain +NestAge.2 6
## ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge.3 7
## ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge.4 7
## ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density.5 7
## ~DBH+Predators+ MayTemp +MayRain +NestAge.7 6
## ~DBH+Predators+ MayTemp +NestAge.8 5
## ~DBH+Predators+ MayTemp +MayRain +NestAge + Density.9 7
## ~DBH+Predators+ MayTemp +JuneRain.6 5
## ~DBH.1 2
## AICc
## ~DBH+Predators+ MayTemp +JuneRain +NestAge.2 1416.92
## ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge.3 1417.93
## ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge.4 1418.28
## ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density.5 1418.93
## ~DBH+Predators+ MayTemp +MayRain +NestAge.7 1422.22
## ~DBH+Predators+ MayTemp +NestAge.8 1422.44
## ~DBH+Predators+ MayTemp +MayRain +NestAge + Density.9 1422.48
## ~DBH+Predators+ MayTemp +JuneRain.6 1424.83
## ~DBH.1 1446.68
## Delta_AICc
## ~DBH+Predators+ MayTemp +JuneRain +NestAge.2 0.00
## ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge.3 1.01
## ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge.4 1.36
## ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density.5 2.01
## ~DBH+Predators+ MayTemp +MayRain +NestAge.7 5.30
## ~DBH+Predators+ MayTemp +NestAge.8 5.51
## ~DBH+Predators+ MayTemp +MayRain +NestAge + Density.9 5.55
## ~DBH+Predators+ MayTemp +JuneRain.6 7.90
## ~DBH.1 29.75
## AICcWt
## ~DBH+Predators+ MayTemp +JuneRain +NestAge.2 0.37
## ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge.3 0.22
## ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge.4 0.19
## ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density.5 0.14
## ~DBH+Predators+ MayTemp +MayRain +NestAge.7 0.03
## ~DBH+Predators+ MayTemp +NestAge.8 0.02
## ~DBH+Predators+ MayTemp +MayRain +NestAge + Density.9 0.02
## ~DBH+Predators+ MayTemp +JuneRain.6 0.01
## ~DBH.1 0.00
## Cum.Wt
## ~DBH+Predators+ MayTemp +JuneRain +NestAge.2 0.37
## ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge.3 0.60
## ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge.4 0.78
## ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density.5 0.92
## ~DBH+Predators+ MayTemp +MayRain +NestAge.7 0.95
## ~DBH+Predators+ MayTemp +NestAge.8 0.97
## ~DBH+Predators+ MayTemp +MayRain +NestAge + Density.9 0.99
## ~DBH+Predators+ MayTemp +JuneRain.6 1.00
## ~DBH.1 1.00
## LL
## ~DBH+Predators+ MayTemp +JuneRain +NestAge.2 -702.45
## ~DBH+Predators+ MayTemp +MayRain +JuneRain +NestAge.3 -701.96
## ~DBH+Predators+ MayTemp+JuneTemp +JuneRain +NestAge.4 -702.13
## ~DBH+Predators+ MayTemp +JuneRain +NestAge + Density.5 -702.45
## ~DBH+Predators+ MayTemp +MayRain +NestAge.7 -705.10
## ~DBH+Predators+ MayTemp +NestAge.8 -706.21
## ~DBH+Predators+ MayTemp +MayRain +NestAge + Density.9 -704.23
## ~DBH+Predators+ MayTemp +JuneRain.6 -707.41
## ~DBH.1 -721.34
# Estimates of beta from the top model
summary(model.fits[[2]])
##
## 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
## -3.0459 0.2078 0.2652 0.3139 0.4474
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.37284 0.66311 2.070 0.038422 *
## DBH 0.02233 0.00547 4.083 4.45e-05 ***
## Predators -0.13793 0.03155 -4.371 1.24e-05 ***
## MayTemp 0.22591 0.06532 3.458 0.000543 ***
## NestAge -0.03560 0.01152 -3.090 0.001999 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1760.6 on 5373 degrees of freedom
## Residual deviance: 1412.4 on 5369 degrees of freedom
## AIC: 1422.4
##
## Number of Fisher Scoring iterations: 7