Programmer: John Fieberg
R Code in support of Fieberg et al. Do Capture and Survey Methods Influence Whether Marked Animals are Representative of Unmarked Animals?
Goal: explore whether detection rates for individual animals changed as a function of time since they were captured.
rm(list = ls())
# Load libraries
library(lme4)
library(ggplot2)
library(nlme)
library(doBy)
library(mosaic)
First read in detection data…
sight.dat<-read.csv("../data/sightdat.csv")
Evaluate whether there is a trend in the age-distrubion relative to year of capture
boxplot(sight.dat$age~sight.dat$year, xlab="", ylab="Age distribution")
Now, get years since capture
sight.dat$yr.s.cap<-sight.dat$year-sight.dat$capyear
One observation was missing the date of capture, delete it:
sight.dat<-sight.dat[is.na(sight.dat$capyear)!=T,]
Get number of individuals by capture year (use median or mean, below - they give the same answer as they should)
capyrs<-summaryBy(capyear~id, data=sight.dat, FUN=c(median, mean))
table(capyrs$capyear.median)
##
## 2002 2003 2004 2005
## 9 29 14 15
Look at correlation between age and time since capture, and also the number of observations missing age
with(sight.dat[is.na(sight.dat$age)!=T,], cor.test(x=age, y=yr.s.cap))
##
## Pearson's product-moment correlation
##
## data: age and yr.s.cap
## t = 3.272, df = 140, p-value = 0.001345
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1065 0.4131
## sample estimates:
## cor
## 0.2665
table(is.na(sight.dat$age))
##
## FALSE TRUE
## 142 17
Net guns were used in 2002, darting in all other years
sight.dat$net<-0
sight.dat$net[sight.dat$capyear==2002]<-1
prop.table(table(sight.dat$net, sight.dat$observed), margin=1) # net gunners were more likely to be detected
##
## 0 1
## 0 0.5255 0.4745
## 1 0.3636 0.6364
chisq.test(table(sight.dat$net, sight.dat$observed)) # not significant, but small sample size
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(sight.dat$net, sight.dat$observed)
## X-squared = 1.393, df = 1, p-value = 0.2379
Look at simple table and plot of detection probabilities as a function of time since capture
temp.tab<-prop.table(table(sight.dat$observed, sight.dat$yr.s.cap), margin=2)
temp.tab
##
## 1 2 3 4 5
## 0 0.4630 0.4423 0.5417 0.6667 0.6000
## 1 0.5370 0.5577 0.4583 0.3333 0.4000
Note: a chi-squared test of Ho: detection probabilities were constant across years ends up being non-significant (but, the test is not geared towards alterantives where there is a monotonic increase/decrease). Also, the test does not adjust for any annual differences in detection that may be due to the distribution of animals in varying levels of cover/visual obstruction.
chisq.test(table(sight.dat$observed, sight.dat$yr.s.cap), simulate.p.value=T)
##
## Pearson's Chi-squared test with simulated p-value (based on 2000
## replicates)
##
## data: table(sight.dat$observed, sight.dat$yr.s.cap)
## X-squared = 4.016, df = NA, p-value = 0.4223
Detection models:
sight.dat$id<-as.factor(sight.dat$id)
fit.glmer1<-glmer(observed~voc+yr.s.cap + net + age + (1|id), family=binomial(), data=sight.dat, nAGQ=10 )
summary(fit.glmer1)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: observed ~ voc + yr.s.cap + net + age + (1 | id)
## Data: sight.dat
##
## AIC BIC logLik deviance df.resid
## 159.6 177.3 -73.8 147.6 136
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.404 -0.623 -0.291 0.602 3.842
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 142, groups: id, 61
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.17419 0.95464 4.37 1.2e-05 ***
## voc -0.04741 0.00862 -5.50 3.8e-08 ***
## yr.s.cap -0.45455 0.21247 -2.14 0.032 *
## net 0.30534 1.00237 0.30 0.761
## age -0.09947 0.07727 -1.29 0.198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) voc yr.s.c net
## voc -0.714
## yr.s.cap -0.386 0.159
## net 0.111 -0.067 -0.092
## age -0.663 0.221 -0.206 -0.128
fit.glmer2<-glmer(observed~voc+yr.s.cap + (1|id), family=binomial(), data=sight.dat, nAGQ=10 )
summary(fit.glmer2)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: observed ~ voc + yr.s.cap + (1 | id)
## Data: sight.dat
##
## AIC BIC logLik deviance df.resid
## 185.0 197.3 -88.5 177.0 155
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.090 -0.675 -0.346 0.720 3.209
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 159, groups: id, 67
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.02903 0.63528 4.77 1.9e-06 ***
## voc -0.04129 0.00749 -5.51 3.6e-08 ***
## yr.s.cap -0.35962 0.16256 -2.21 0.027 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) voc
## voc -0.783
## yr.s.cap -0.715 0.227
fit.glmer3<-glmer(observed~voc+ net +(1|id), family=binomial(), data=sight.dat, nAGQ=10 )
summary(fit.glmer3)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: observed ~ voc + net + (1 | id)
## Data: sight.dat
##
## AIC BIC logLik deviance df.resid
## 188.4 200.7 -90.2 180.4 155
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.584 -0.697 -0.389 0.693 2.332
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 159, groups: id, 67
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.99247 0.43476 4.58 4.6e-06 ***
## voc -0.03879 0.00712 -5.45 5.2e-08 ***
## net 0.68179 0.52503 1.30 0.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) voc
## voc -0.894
## net -0.138 -0.032
fit.glmer4<-glmer(observed~yr.s.cap+(1|id), family=binomial(), data=sight.dat, nAGQ=10 )
summary(fit.glmer4)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: observed ~ yr.s.cap + (1 | id)
## Data: sight.dat
##
## AIC BIC logLik deviance df.resid
## 223.5 232.7 -108.7 217.5 156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.146 -1.018 -0.713 0.983 1.403
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.3e-18 1.14e-09
## Number of obs: 159, groups: id, 67
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.510 0.347 1.47 0.142
## yr.s.cap -0.237 0.140 -1.69 0.091 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## yr.s.cap -0.887
fit.glmer5<-glmer(observed~ net +(1|id), family=binomial(), data=sight.dat, nAGQ=10 )
summary(fit.glmer5)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: observed ~ net + (1 | id)
## Data: sight.dat
##
## AIC BIC logLik deviance df.resid
## 224.4 233.6 -109.2 218.4 156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.32 -0.95 -0.95 1.05 1.05
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0 0
## Number of obs: 159, groups: id, 67
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.102 0.171 -0.60 0.55
## net 0.662 0.475 1.39 0.16
##
## Correlation of Fixed Effects:
## (Intr)
## net -0.360
# Can't compare models with and without age using AIC because of missing age values. Look at amount of missing data, below:
table(is.na(sight.dat$age))
##
## FALSE TRUE
## 142 17
For plotting purposes fit a model that assumes independence (since re variance for animal is ~0)
fit1<-glm(observed~voc+yr.s.cap, family=binomial(), data=sight.dat)
summary(fit1)
##
## Call:
## glm(formula = observed ~ voc + yr.s.cap, family = binomial(),
## data = sight.dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.171 -0.866 -0.475 0.913 2.202
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.02903 0.63529 4.77 1.9e-06 ***
## voc -0.04129 0.00749 -5.51 3.6e-08 ***
## yr.s.cap -0.35962 0.16256 -2.21 0.027 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 220.41 on 158 degrees of freedom
## Residual deviance: 177.05 on 156 degrees of freedom
## AIC: 183
##
## Number of Fisher Scoring iterations: 4
newdat<-data.frame(observe=1, voc=mean(sight.dat$voc), yr.s.cap=1:5)
p1<-predict(fit1, newdat, se=T)
Odds of detection per unit increase in years since capture (to report effect size)
exp(coef(fit1)[3])
## yr.s.cap
## 0.6979
exp(coef(fit1)[3]+ 1.96*sqrt(vcov(fit1)[3,3]))
## yr.s.cap
## 0.9598
exp(coef(fit1)[3]- 1.96*sqrt(vcov(fit1)[3,3]))
## yr.s.cap
## 0.5075
#Uncomment to produce single, 2 panel figure
#par(mfrow=c(1,2), bty="L", cex.lab=1.8, cex.axis=1.8, mar=c(6,5,4.1,2.1) )
plot(as.numeric(colnames(temp.tab)), temp.tab[2,], xlab="Years since capture", ylab="Detection Probability", ylim=c(0,.8), main="", xlim=c(0.5, 5.5), pch=16)
text(1:5, temp.tab[2,], paste("n = ", table(sight.dat$yr.s.cap), sep=""), pos=3, cex=1.4)
lines(1:5,plogis(p1$fit))
lines(1:5,plogis(p1$fit+1.96*p1$se.fit), lty=2)
lines(1:5,plogis(p1$fit-1.96*p1$se.fit), lty=2)
mtext(side=3, outer=F, cex=2, "A)")
boxplot(voc~yr.s.cap, data=sight.dat, ylab="Visual Obstruction", xlab="Years since capture")
# Alternative plot
densityplot(~voc, groups=as.factor(yr.s.cap), data=sight.dat, xlab="Visual Obstruction", auto.key=TRUE)
mtext(side=3, outer=F, cex=2, "B)")
## Error: plot.new has not been called yet
# Lets look at voc relative to time since capture
fit.voc1<-lmer(voc~yr.s.cap+(1|id), data=sight.dat)
summary(fit.voc1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: voc ~ yr.s.cap + (1 | id)
## Data: sight.dat
##
## REML criterion at convergence: 1507
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8157 -0.6759 0.0866 0.7728 1.6077
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 167 12.9
## Residual 673 25.9
## Number of obs: 159, groups: id, 67
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 54.484 5.070 10.75
## yr.s.cap -0.346 2.019 -0.17
##
## Correlation of Fixed Effects:
## (Intr)
## yr.s.cap -0.850
fit.voc2<-lmer(voc~age+yr.s.cap+(1|id), data=sight.dat)
summary(fit.voc2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: voc ~ age + yr.s.cap + (1 | id)
## Data: sight.dat
##
## REML criterion at convergence: 1344
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8098 -0.7108 0.0743 0.7737 1.5755
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 151 12.3
## Residual 697 26.4
## Number of obs: 142, groups: id, 61
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 57.178 8.718 6.56
## age -0.663 1.047 -0.63
## yr.s.cap 1.239 2.526 0.49
##
## Correlation of Fixed Effects:
## (Intr) age
## age -0.780
## yr.s.cap -0.262 -0.317
# Use lme to get approximate pvalue
fit.voc.lme<-lme(voc~yr.s.cap, random=list(id=~1), data=sight.dat)
summary(fit.voc.lme)
## Linear mixed-effects model fit by REML
## Data: sight.dat
## AIC BIC logLik
## 1515 1528 -753.7
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 12.93 25.93
##
## Fixed effects: voc ~ yr.s.cap
## Value Std.Error DF t-value p-value
## (Intercept) 54.48 5.070 91 10.747 0.0000
## yr.s.cap -0.35 2.019 91 -0.171 0.8645
## Correlation:
## (Intr)
## yr.s.cap -0.85
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.81565 -0.67586 0.08661 0.77280 1.60767
##
## Number of Observations: 159
## Number of Groups: 67
# Figure (for presentation)
plot(as.numeric(colnames(temp.tab)), temp.tab[2,], xlab="Years since capture", ylab="Detection Probability", ylim=c(0,.8), main="", xlim=c(0.5, 5.5), pch=16)
text(1:5, temp.tab[2,], paste("n = ", table(sight.dat$yr.s.cap), sep=""), pos=3, cex=1.4)
lines(1:5,plogis(p1$fit))
lines(1:5,plogis(p1$fit+1.96*p1$se.fit), lty=2)
lines(1:5,plogis(p1$fit-1.96*p1$se.fit), lty=2)
# Print out Session info
sessionInfo()
## R version 3.1.1 (2014-07-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] splines stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] knitr_1.6 mosaic_0.9.1-3 lattice_0.20-29 dplyr_0.4.1
## [5] car_2.0-21 doBy_4.5-10 MASS_7.3-33 survival_2.37-7
## [9] nlme_3.1-117 ggplot2_1.0.0 lme4_1.1-7 Rcpp_0.11.2
## [13] Matrix_1.1-4
##
## loaded via a namespace (and not attached):
## [1] assertthat_0.1 colorspace_1.2-4 DBI_0.3.0 digest_0.6.4
## [5] evaluate_0.5.5 formatR_1.0 ggdendro_0.1-14 grid_3.1.1
## [9] gridExtra_0.9.1 gtable_0.1.2 magrittr_1.0.1 markdown_0.7.4
## [13] mime_0.1.2 minqa_1.2.3 mosaicData_0.9.1 munsell_0.4.2
## [17] nloptr_1.0.4 nnet_7.3-8 parallel_3.1.1 plyr_1.8.1
## [21] proto_0.3-10 reshape2_1.4 scales_0.2.4 stringr_0.6.2
## [25] tools_3.1.1