# Analysis of mallard data using logistic exposure model
# How to fit a nest level continous covariate (Robel height
# 2019-05-01 CJS Initial code
# This is the mallard data that ships with RMark
library(AICcmodavg)
library(emmeans)
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
#
# Also contains the following fields
# Robel - Measurement of Robel pole of visibility of nest
# PpnGrass - proportion of grassland cover on the 10.4 km2 study site t
# AgeFound - Age of nest when found
# AgeDay1 - Age of nest on day 1 of study (can be negative)
# Habitat - N=Native; P=Planted; W=Wetland; R=roadside right of way
#
# 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.
#
malldata <- readxl::read_excel(file.path("..","mallard.xlsx"),
sheet="mallard")
head(malldata)
## # A tibble: 6 x 10
## FirstFound LastPresent LastChecked Fate Freq Robel PpnGrass AgeFound
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 73 89 89 0 1 6 0.800 13
## 2 63 90 90 0 1 3.75 0.668 5
## 3 70 70 76 1 1 3.75 0.800 13
## 4 63 81 85 1 1 3.12 0.668 6
## 5 61 61 66 1 1 4.5 0.668 4
## 6 57 57 61 1 1 4.25 0.668 1
## # … with 2 more variables: AgeDay1 <dbl>, Habitat <chr>
malldata2 <- expand.nest.data(malldata)
rm(malldata)
malldata2$Habitat <- factor(malldata2$Habitat) # RMark wants categorical variable to be factors
# Fit a particular model
# This is a model with S varying by robel height
mod.rob <- glm(Survive~Robel,
family=binomial(link=logexp(malldata2$Exposure)),
data=malldata2)
summary(mod.rob)
##
## Call:
## glm(formula = Survive ~ Robel, family = binomial(link = logexp(malldata2$Exposure)),
## data = malldata2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4814 0.3070 0.3106 0.3142 0.3232
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.90884 0.16827 17.287 <2e-16 ***
## Robel 0.02727 0.04499 0.606 0.544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2494.8 on 6125 degrees of freedom
## Residual deviance: 1566.8 on 6124 degrees of freedom
## AIC: 1570.8
##
## Number of Fisher Scoring iterations: 7
# You can make predictions just as before
pred.data <- data.frame(Robel=seq(min(malldata2$Robel), max(malldata2$Robel), length.out=50))
logit.dsr.pred.rob <- predict(mod.rob, newdata=pred.data, se.fit=TRUE)
# put these together in a data frame
dsr.rob <- cbind(pred.data, logit.dsr=logit.dsr.pred.rob$fit, logit.dsr.se=logit.dsr.pred.rob$se.fit)
dsr.rob$dsr <- expit(dsr.rob$logit.dsr)
dsr.rob$dsr.se <- dsr.rob$logit.dsr.se* dsr.rob$dsr * (1-dsr.rob$dsr)
dsr.rob
## Robel logit.dsr logit.dsr.se dsr dsr.se
## 1 0.6250000 2.925881 0.14202616 0.9491111 0.006859751
## 2 0.8010204 2.930681 0.13477235 0.9493425 0.006481385
## 3 0.9770408 2.935482 0.12759770 0.9495728 0.006109926
## 4 1.1530612 2.940282 0.12051636 0.9498022 0.005745981
## 5 1.3290816 2.945082 0.11354576 0.9500305 0.005390305
## 6 1.5051020 2.949882 0.10670764 0.9502579 0.005043838
## 7 1.6811224 2.954682 0.10002914 0.9504843 0.004707761
## 8 1.8571429 2.959482 0.09354448 0.9507097 0.004383564
## 9 2.0331633 2.964282 0.08729685 0.9509342 0.004073129
## 10 2.2091837 2.969082 0.08134089 0.9511577 0.003778835
## 11 2.3852041 2.973882 0.07574545 0.9513802 0.003503677
## 12 2.5612245 2.978682 0.07059628 0.9516017 0.003251374
## 13 2.7372449 2.983482 0.06599794 0.9518223 0.003026441
## 14 2.9132653 2.988282 0.06207296 0.9520420 0.002834132
## 15 3.0892857 2.993083 0.05895596 0.9522606 0.002680157
## 16 3.2653061 2.997883 0.05678019 0.9524784 0.002570060
## 17 3.4413265 3.002683 0.05565614 0.9526952 0.002508259
## 18 3.6173469 3.007483 0.05564756 0.9529110 0.002496995
## 19 3.7933673 3.012283 0.05675498 0.9531260 0.002535635
## 20 3.9693878 3.017083 0.05891549 0.9533399 0.002620732
## 21 4.1454082 3.021883 0.06201913 0.9535530 0.002746807
## 22 4.3214286 3.026683 0.06593284 0.9537651 0.002907454
## 23 4.4974490 3.031483 0.07052190 0.9539763 0.003096298
## 24 4.6734694 3.036283 0.07566352 0.9541866 0.003307593
## 25 4.8494898 3.041083 0.08125286 0.9543960 0.003536471
## 26 5.0255102 3.045884 0.08720388 0.9546045 0.003778960
## 27 5.2015306 3.050684 0.09344752 0.9548120 0.004031888
## 28 5.3775510 3.055484 0.09992892 0.9550187 0.004292746
## 29 5.5535714 3.060284 0.10660474 0.9552244 0.004559561
## 30 5.7295918 3.065084 0.11344066 0.9554293 0.004830775
## 31 5.9056122 3.069884 0.12040941 0.9556333 0.005105159
## 32 6.0816327 3.074684 0.12748922 0.9558363 0.005381733
## 33 6.2576531 3.079484 0.13466256 0.9560385 0.005659716
## 34 6.4336735 3.084284 0.14191525 0.9562398 0.005938477
## 35 6.6096939 3.089084 0.14923573 0.9564402 0.006217506
## 36 6.7857143 3.093884 0.15661450 0.9566398 0.006496387
## 37 6.9617347 3.098685 0.16404367 0.9568385 0.006774779
## 38 7.1377551 3.103485 0.17151672 0.9570363 0.007052401
## 39 7.3137755 3.108285 0.17902813 0.9572332 0.007329020
## 40 7.4897959 3.113085 0.18657328 0.9574293 0.007604441
## 41 7.6658163 3.117885 0.19414824 0.9576245 0.007878503
## 42 7.8418367 3.122685 0.20174964 0.9578188 0.008151071
## 43 8.0178571 3.127485 0.20937461 0.9580123 0.008422029
## 44 8.1938776 3.132285 0.21702066 0.9582050 0.008691281
## 45 8.3698980 3.137085 0.22468564 0.9583968 0.008958746
## 46 8.5459184 3.141885 0.23236768 0.9585878 0.009224356
## 47 8.7219388 3.146685 0.24006514 0.9587779 0.009488052
## 48 8.8979592 3.151485 0.24777658 0.9589672 0.009749784
## 49 9.0739796 3.156286 0.25550074 0.9591557 0.010009513
## 50 9.2500000 3.161086 0.26323649 0.9593433 0.010267203
# we plot the results
dsr.rob.plot <- ggplot(data=dsr.rob, aes(x=Robel, y=dsr))+
ggtitle("DSR by Robel height for mallards", subtitle="Logistic exposure model")+
geom_line()+
geom_ribbon(aes(ymin=expit(logit.dsr-1.95*logit.dsr.se),
ymax=expit(logit.dsr+1.95*logit.dsr.se)),alpha=0.1)+
ylab("DSR (95% ci)")
dsr.rob.plot
ggsave(dsr.rob.plot,
file=file.path("..","..","..","..","MyStuff","Images","mallard-dsr-robel-le.png"), h=4, w=6, units="in", dpi=300)
# fit a null models
mod.null <- glm(Survive~1,
family=binomial(link=logexp(malldata2$Exposure)),
data=malldata2)
summary(mod.null)
##
## Call:
## glm(formula = Survive ~ 1, family = binomial(link = logexp(malldata2$Exposure)),
## data = malldata2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4714 0.3109 0.3109 0.3109 0.3109
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.00565 0.05551 54.15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2494.8 on 6125 degrees of freedom
## Residual deviance: 1567.1 on 6125 degrees of freedom
## AIC: 1569.1
##
## Number of Fisher Scoring iterations: 7
AICcmodavg::aictab( list(mod.rob, mod.null))
## Warning in aictab.AICglm.lm(list(mod.rob, 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 1569.12 0.00 0.7 0.7 -783.56
## Mod1 2 1570.77 1.66 0.3 1.0 -783.39