# Analysis of killdeer data illustrating logistic exposure models
# Incorportating an age variable
# 2019-05-01 CJS Initial code
# This is the killdeer data that ships with RMark
# I added arbitrary nest ages.
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
# AgeDay1: age of nest at day 1 of study (can be negative)
#
# 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.
#
killdata <- readxl::read_excel(file.path("..","Killdeer.xlsx"),
sheet="killdeer-age")
killdata <- as.data.frame(killdata) # sometimes doesn't play nice with tibbles
head(killdata)
## id FirstFound LastPresent LastChecked Fate Freq AgeDay1
## 1 /*A*/ 1 9 9 0 1 0
## 2 /*B*/ 5 5 9 1 1 -2
## 3 /*C*/ 5 40 40 0 1 -3
## 4 /*D*/ 9 32 32 0 1 -4
## 5 /*E*/ 7 8 8 0 1 -4
## 6 /*F*/ 3 15 15 0 1 1
killdata2 <- expand.nest.data(killdata)
head(killdata2)
## id FirstFound LastPresent LastChecked Fate Freq AgeDay1 Day Exposure
## 1 /*A*/ 1 9 9 0 1 0 1 1
## 2 /*A*/ 1 9 9 0 1 0 2 1
## 3 /*A*/ 1 9 9 0 1 0 3 1
## 4 /*A*/ 1 9 9 0 1 0 4 1
## 5 /*A*/ 1 9 9 0 1 0 5 1
## 6 /*A*/ 1 9 9 0 1 0 6 1
## Survive NestAge
## 1 1 0
## 2 1 1
## 3 1 2
## 4 1 3
## 5 1 4
## 6 1 5
print(killdata2[,c(1,2,3,4,7,8,11)], row.names=FALSE)
## id FirstFound LastPresent LastChecked AgeDay1 Day NestAge
## /*A*/ 1 9 9 0 1.0 0.0
## /*A*/ 1 9 9 0 2.0 1.0
## /*A*/ 1 9 9 0 3.0 2.0
## /*A*/ 1 9 9 0 4.0 3.0
## /*A*/ 1 9 9 0 5.0 4.0
## /*A*/ 1 9 9 0 6.0 5.0
## /*A*/ 1 9 9 0 7.0 6.0
## /*A*/ 1 9 9 0 8.0 7.0
## /*B*/ 5 5 9 -2 7.0 4.0
## /*C*/ 5 40 40 -3 5.0 1.0
## /*C*/ 5 40 40 -3 6.0 2.0
## /*C*/ 5 40 40 -3 7.0 3.0
## /*C*/ 5 40 40 -3 8.0 4.0
## /*C*/ 5 40 40 -3 9.0 5.0
## /*C*/ 5 40 40 -3 10.0 6.0
## /*C*/ 5 40 40 -3 11.0 7.0
## /*C*/ 5 40 40 -3 12.0 8.0
## /*C*/ 5 40 40 -3 13.0 9.0
## /*C*/ 5 40 40 -3 14.0 10.0
## /*C*/ 5 40 40 -3 15.0 11.0
## /*C*/ 5 40 40 -3 16.0 12.0
## /*C*/ 5 40 40 -3 17.0 13.0
## /*C*/ 5 40 40 -3 18.0 14.0
## /*C*/ 5 40 40 -3 19.0 15.0
## /*C*/ 5 40 40 -3 20.0 16.0
## /*C*/ 5 40 40 -3 21.0 17.0
## /*C*/ 5 40 40 -3 22.0 18.0
## /*C*/ 5 40 40 -3 23.0 19.0
## /*C*/ 5 40 40 -3 24.0 20.0
## /*C*/ 5 40 40 -3 25.0 21.0
## /*C*/ 5 40 40 -3 26.0 22.0
## /*C*/ 5 40 40 -3 27.0 23.0
## /*C*/ 5 40 40 -3 28.0 24.0
## /*C*/ 5 40 40 -3 29.0 25.0
## /*C*/ 5 40 40 -3 30.0 26.0
## /*C*/ 5 40 40 -3 31.0 27.0
## /*C*/ 5 40 40 -3 32.0 28.0
## /*C*/ 5 40 40 -3 33.0 29.0
## /*C*/ 5 40 40 -3 34.0 30.0
## /*C*/ 5 40 40 -3 35.0 31.0
## /*C*/ 5 40 40 -3 36.0 32.0
## /*C*/ 5 40 40 -3 37.0 33.0
## /*C*/ 5 40 40 -3 38.0 34.0
## /*C*/ 5 40 40 -3 39.0 35.0
## /*D*/ 9 32 32 -4 9.0 4.0
## /*D*/ 9 32 32 -4 10.0 5.0
## /*D*/ 9 32 32 -4 11.0 6.0
## /*D*/ 9 32 32 -4 12.0 7.0
## /*D*/ 9 32 32 -4 13.0 8.0
## /*D*/ 9 32 32 -4 14.0 9.0
## /*D*/ 9 32 32 -4 15.0 10.0
## /*D*/ 9 32 32 -4 16.0 11.0
## /*D*/ 9 32 32 -4 17.0 12.0
## /*D*/ 9 32 32 -4 18.0 13.0
## /*D*/ 9 32 32 -4 19.0 14.0
## /*D*/ 9 32 32 -4 20.0 15.0
## /*D*/ 9 32 32 -4 21.0 16.0
## /*D*/ 9 32 32 -4 22.0 17.0
## /*D*/ 9 32 32 -4 23.0 18.0
## /*D*/ 9 32 32 -4 24.0 19.0
## /*D*/ 9 32 32 -4 25.0 20.0
## /*D*/ 9 32 32 -4 26.0 21.0
## /*D*/ 9 32 32 -4 27.0 22.0
## /*D*/ 9 32 32 -4 28.0 23.0
## /*D*/ 9 32 32 -4 29.0 24.0
## /*D*/ 9 32 32 -4 30.0 25.0
## /*D*/ 9 32 32 -4 31.0 26.0
## /*E*/ 7 8 8 -4 7.0 2.0
## /*F*/ 3 15 15 1 3.0 3.0
## /*F*/ 3 15 15 1 4.0 4.0
## /*F*/ 3 15 15 1 5.0 5.0
## /*F*/ 3 15 15 1 6.0 6.0
## /*F*/ 3 15 15 1 7.0 7.0
## /*F*/ 3 15 15 1 8.0 8.0
## /*F*/ 3 15 15 1 9.0 9.0
## /*F*/ 3 15 15 1 10.0 10.0
## /*F*/ 3 15 15 1 11.0 11.0
## /*F*/ 3 15 15 1 12.0 12.0
## /*F*/ 3 15 15 1 13.0 13.0
## /*F*/ 3 15 15 1 14.0 14.0
## /*G*/ 8 32 32 -7 8.0 0.0
## /*G*/ 8 32 32 -7 9.0 1.0
## /*G*/ 8 32 32 -7 10.0 2.0
## /*G*/ 8 32 32 -7 11.0 3.0
## /*G*/ 8 32 32 -7 12.0 4.0
## /*G*/ 8 32 32 -7 13.0 5.0
## /*G*/ 8 32 32 -7 14.0 6.0
## /*G*/ 8 32 32 -7 15.0 7.0
## /*G*/ 8 32 32 -7 16.0 8.0
## /*G*/ 8 32 32 -7 17.0 9.0
## /*G*/ 8 32 32 -7 18.0 10.0
## /*G*/ 8 32 32 -7 19.0 11.0
## /*G*/ 8 32 32 -7 20.0 12.0
## /*G*/ 8 32 32 -7 21.0 13.0
## /*G*/ 8 32 32 -7 22.0 14.0
## /*G*/ 8 32 32 -7 23.0 15.0
## /*G*/ 8 32 32 -7 24.0 16.0
## /*G*/ 8 32 32 -7 25.0 17.0
## /*G*/ 8 32 32 -7 26.0 18.0
## /*G*/ 8 32 32 -7 27.0 19.0
## /*G*/ 8 32 32 -7 28.0 20.0
## /*G*/ 8 32 32 -7 29.0 21.0
## /*G*/ 8 32 32 -7 30.0 22.0
## /*G*/ 8 32 32 -7 31.0 23.0
## /*H*/ 14 14 16 -10 15.0 4.0
## /*I*/ 8 14 14 -7 8.0 0.0
## /*I*/ 8 14 14 -7 9.0 1.0
## /*I*/ 8 14 14 -7 10.0 2.0
## /*I*/ 8 14 14 -7 11.0 3.0
## /*I*/ 8 14 14 -7 12.0 4.0
## /*I*/ 8 14 14 -7 13.0 5.0
## /*J*/ 13 14 14 -12 13.0 0.0
## /*K*/ 14 33 33 -13 14.0 0.0
## /*K*/ 14 33 33 -13 15.0 1.0
## /*K*/ 14 33 33 -13 16.0 2.0
## /*K*/ 14 33 33 -13 17.0 3.0
## /*K*/ 14 33 33 -13 18.0 4.0
## /*K*/ 14 33 33 -13 19.0 5.0
## /*K*/ 14 33 33 -13 20.0 6.0
## /*K*/ 14 33 33 -13 21.0 7.0
## /*K*/ 14 33 33 -13 22.0 8.0
## /*K*/ 14 33 33 -13 23.0 9.0
## /*K*/ 14 33 33 -13 24.0 10.0
## /*K*/ 14 33 33 -13 25.0 11.0
## /*K*/ 14 33 33 -13 26.0 12.0
## /*K*/ 14 33 33 -13 27.0 13.0
## /*K*/ 14 33 33 -13 28.0 14.0
## /*K*/ 14 33 33 -13 29.0 15.0
## /*K*/ 14 33 33 -13 30.0 16.0
## /*K*/ 14 33 33 -13 31.0 17.0
## /*K*/ 14 33 33 -13 32.0 18.0
## /*L*/ 15 37 37 -10 15.0 4.0
## /*L*/ 15 37 37 -10 16.0 5.0
## /*L*/ 15 37 37 -10 17.0 6.0
## /*L*/ 15 37 37 -10 18.0 7.0
## /*L*/ 15 37 37 -10 19.0 8.0
## /*L*/ 15 37 37 -10 20.0 9.0
## /*L*/ 15 37 37 -10 21.0 10.0
## /*L*/ 15 37 37 -10 22.0 11.0
## /*L*/ 15 37 37 -10 23.0 12.0
## /*L*/ 15 37 37 -10 24.0 13.0
## /*L*/ 15 37 37 -10 25.0 14.0
## /*L*/ 15 37 37 -10 26.0 15.0
## /*L*/ 15 37 37 -10 27.0 16.0
## /*L*/ 15 37 37 -10 28.0 17.0
## /*L*/ 15 37 37 -10 29.0 18.0
## /*L*/ 15 37 37 -10 30.0 19.0
## /*L*/ 15 37 37 -10 31.0 20.0
## /*L*/ 15 37 37 -10 32.0 21.0
## /*L*/ 15 37 37 -10 33.0 22.0
## /*L*/ 15 37 37 -10 34.0 23.0
## /*L*/ 15 37 37 -10 35.0 24.0
## /*L*/ 15 37 37 -10 36.0 25.0
## /*M*/ 16 37 40 -11 16.0 4.0
## /*M*/ 16 37 40 -11 17.0 5.0
## /*M*/ 16 37 40 -11 18.0 6.0
## /*M*/ 16 37 40 -11 19.0 7.0
## /*M*/ 16 37 40 -11 20.0 8.0
## /*M*/ 16 37 40 -11 21.0 9.0
## /*M*/ 16 37 40 -11 22.0 10.0
## /*M*/ 16 37 40 -11 23.0 11.0
## /*M*/ 16 37 40 -11 24.0 12.0
## /*M*/ 16 37 40 -11 25.0 13.0
## /*M*/ 16 37 40 -11 26.0 14.0
## /*M*/ 16 37 40 -11 27.0 15.0
## /*M*/ 16 37 40 -11 28.0 16.0
## /*M*/ 16 37 40 -11 29.0 17.0
## /*M*/ 16 37 40 -11 30.0 18.0
## /*M*/ 16 37 40 -11 31.0 19.0
## /*M*/ 16 37 40 -11 32.0 20.0
## /*M*/ 16 37 40 -11 33.0 21.0
## /*M*/ 16 37 40 -11 34.0 22.0
## /*M*/ 16 37 40 -11 35.0 23.0
## /*M*/ 16 37 40 -11 36.0 24.0
## /*M*/ 16 37 40 -11 38.5 26.5
## /*N*/ 16 28 32 -11 16.0 4.0
## /*N*/ 16 28 32 -11 17.0 5.0
## /*N*/ 16 28 32 -11 18.0 6.0
## /*N*/ 16 28 32 -11 19.0 7.0
## /*N*/ 16 28 32 -11 20.0 8.0
## /*N*/ 16 28 32 -11 21.0 9.0
## /*N*/ 16 28 32 -11 22.0 10.0
## /*N*/ 16 28 32 -11 23.0 11.0
## /*N*/ 16 28 32 -11 24.0 12.0
## /*N*/ 16 28 32 -11 25.0 13.0
## /*N*/ 16 28 32 -11 26.0 14.0
## /*N*/ 16 28 32 -11 27.0 15.0
## /*N*/ 16 28 32 -11 30.0 18.0
## /*O*/ 16 17 17 -12 16.0 3.0
## /*P*/ 21 28 33 -15 21.0 5.0
## /*P*/ 21 28 33 -15 22.0 6.0
## /*P*/ 21 28 33 -15 23.0 7.0
## /*P*/ 21 28 33 -15 24.0 8.0
## /*P*/ 21 28 33 -15 25.0 9.0
## /*P*/ 21 28 33 -15 26.0 10.0
## /*P*/ 21 28 33 -15 27.0 11.0
## /*P*/ 21 28 33 -15 30.5 14.5
## /*Q*/ 23 33 34 -17 23.0 5.0
## /*Q*/ 23 33 34 -17 24.0 6.0
## /*Q*/ 23 33 34 -17 25.0 7.0
## /*Q*/ 23 33 34 -17 26.0 8.0
## /*Q*/ 23 33 34 -17 27.0 9.0
## /*Q*/ 23 33 34 -17 28.0 10.0
## /*Q*/ 23 33 34 -17 29.0 11.0
## /*Q*/ 23 33 34 -17 30.0 12.0
## /*Q*/ 23 33 34 -17 31.0 13.0
## /*Q*/ 23 33 34 -17 32.0 14.0
## /*Q*/ 23 33 34 -17 33.5 15.5
## /*R*/ 27 29 29 -23 27.0 3.0
## /*R*/ 27 29 29 -23 28.0 4.0
# Fit a particular model
# This is a model with S a function of nest age
mod.age <- glm(Survive~NestAge,
family=binomial(link=logexp(killdata2$Exposure)),
data=killdata2)
summary(mod.age)
##
## Call:
## glm(formula = Survive ~ NestAge, family = binomial(link = logexp(killdata2$Exposure)),
## data = killdata2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6533 0.2151 0.2300 0.2499 0.3140
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.89120 0.79883 4.871 1.11e-06 ***
## NestAge -0.02589 0.04996 -0.518 0.604
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 54.491 on 209 degrees of freedom
## Residual deviance: 42.257 on 208 degrees of freedom
## AIC: 46.257
##
## Number of Fisher Scoring iterations: 7
# predict surival as a function of nest age
pred.data <- data.frame(NestAge=1:10)
logit.dsr.pred.age <- predict(mod.age, newdata=pred.data, se.fit=TRUE)
# put these together in a data frame
dsr.age <- cbind(pred.data, logit.dsr=logit.dsr.pred.age$fit, logit.dsr.se=logit.dsr.pred.age$se.fit)
head(dsr.age)
## NestAge logit.dsr logit.dsr.se
## 1 1 3.865309 0.7563272
## 2 2 3.839419 0.7147924
## 3 3 3.813528 0.6744015
## 4 4 3.787638 0.6353727
## 5 5 3.761747 0.5979726
## 6 6 3.735857 0.5625264
dsr.age$dsr <- expit(dsr.age$logit.dsr)
dsr.age$dsr.se <- dsr.age$logit.dsr.se* dsr.age$dsr * (1-dsr.age$dsr)
head(dsr.age)
## NestAge logit.dsr logit.dsr.se dsr dsr.se
## 1 1 3.865309 0.7563272 0.9794737 0.01520592
## 2 2 3.839419 0.7147924 0.9789467 0.01473193
## 3 3 3.813528 0.6744015 0.9784064 0.01424829
## 4 4 3.787638 0.6353727 0.9778526 0.01376021
## 5 5 3.761747 0.5979726 0.9772849 0.01327448
## 6 6 3.735857 0.5625264 0.9767030 0.01279988
# plot this in the usual way
killdeer.age.le <- ggplot(data=dsr.age, aes(x=NestAge, y=dsr, group=1))+
ggtitle("DSR as a function of nest age", subtitle="Logistic exposure model")+
geom_line()+
geom_ribbon(aes(ymin=expit(logit.dsr-1.96*logit.dsr.se),
ymax=expit(logit.dsr+1.96*logit.dsr.se)), alpha=.1)+
ylim(0.5, 1)
killdeer.age.le
ggsave(killdeer.age.le,
file=file.path("..","..","..","..","MyStuff","Images","killdeer-age-le.png"),h=4, w=6, units="in", dpi=300)
#_----------------------------------------------------------------------
# Compare to the null model
mod.null <- glm(Survive~1,
family=binomial(link=logexp(killdata2$Exposure)),
data=killdata2)
AICcmodavg::aictab(list(mod.age, mod.null))
## Warning in aictab.AICglm.lm(list(mod.age, mod.null)):
## Model names have been supplied automatically in the table
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod2 1 44.53 0.00 0.71 0.71 -21.26
## Mod1 2 46.31 1.79 0.29 1.00 -21.13