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?
# Clear out working memory
rm(list = ls())
# set global option for output when producing html file
opts_chunk$set(comment=NA, fig.width=6, fig.height=6)
load libraries
library(lme4)
library(doBy)
library(gplots)
library(mosaic)
library(AICcmodavg)
library(MuMIn)
Function for getting SE
SE<-function(x){sqrt(var(x)/length(x))}
Read in data:
mtg<-read.csv("../Data/MTG_Sight_Alaska.csv")
Look at number of observations for each animal
table(mtg$ID)
LG003 LG004 LG005 LG007 LG008 LG009 LG010 LG011 LG015 LG016 LG017 LG020
1 1 1 1 1 1 1 1 1 1 1 1
LG021 LG022 LG023 LG024 LG025 LG026 LG027 LG029 LG030 LG036 LG037 LG038
1 1 1 2 1 1 2 2 1 2 2 2
LG039 LG040 LG041 LG045 LG049 LG050 LG051 LG052 LG053 LG054 LG055 LG056
1 3 2 2 1 1 2 2 2 1 1 1
LG057 LG058 LG059 LG060 LG061 LG062 LG063 LG064 LG066 LG068 LG069 LG070
1 1 1 2 3 2 1 2 1 1 2 1
LG071 LG074 LG075 LG076 LG077 LG078 LG080 LG081 LG083 LG084 LG086 LG087
1 1 1 2 2 1 1 1 2 2 2 2
LG088 LG089 LG090 LG091 LG092 LG093 LG094 LG095 LG096 LG097 LG098 LG099
2 2 1 1 1 1 1 1 1 1 1 1
LG100 LG101 LG102 LG103 LG104 LG105 LG106 LG108 LG109 LG110 LG111 LG112
1 1 1 1 1 1 1 1 1 1 1 1
LG113 LG114 LG115 LG116 LG117 LG118 LG119 LG120
1 1 1 1 1 1 1 1
Distribution of time since capture
plot(density(mtg$Day.Since.Capture, from=0), xlab="Days since capture", ylab="Density", main="")
Since so few animals were observed more than once, fit a glm (assuming independence), and then use bootstrap for standard errors if necessary (e.g., if statistically significant when independence is assumed).
fit.mtg<-glm(Seen.~Day.Since.Capture, family=binomial(), data=mtg)
summary(fit.mtg)
Call:
glm(formula = Seen. ~ Day.Since.Capture, family = binomial(),
data = mtg)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.19 -1.18 -1.11 1.18 1.25
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.018938 0.252848 0.07 0.94
Day.Since.Capture -0.000416 0.001044 -0.40 0.69
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 164.89 on 118 degrees of freedom
Residual deviance: 164.73 on 117 degrees of freedom
AIC: 168.7
Number of Fisher Scoring iterations: 3
Could also consider controlling for area if pdetect depeneded on area….
fit.mtg2<-glm(Seen.~Day.Since.Capture+Area2, family=binomial(), data=mtg)
summary(fit.mtg2)
Call:
glm(formula = Seen. ~ Day.Since.Capture + Area2, family = binomial(),
data = mtg)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.344 -1.139 -0.961 1.202 1.430
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.56e+01 1.46e+03 0.01 0.99
Day.Since.Capture -6.21e-04 1.09e-03 -0.57 0.57
Area2East Berners -1.59e+01 1.46e+03 -0.01 0.99
Area2Lions Head -1.52e+01 1.46e+03 -0.01 0.99
Area2Sinclair -1.56e+01 1.46e+03 -0.01 0.99
Area2Villard -1.59e+01 1.46e+03 -0.01 0.99
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 164.89 on 118 degrees of freedom
Residual deviance: 161.41 on 113 degrees of freedom
AIC: 173.4
Number of Fisher Scoring iterations: 14
Days since capture does not appear to be significant. What about “years since capture” (since the number of days since capture is bi-modal)
mtg$Itime<-I(mtg$Day.Since.Capture > 200)
fit.mtg3<-glmer(Seen.~Itime + (1|ID), family=binomial(), nAGQ=10, data=mtg)
summary(fit.mtg3)
Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
Family: binomial ( logit )
Formula: Seen. ~ Itime + (1 | ID)
Data: mtg
AIC BIC logLik deviance df.resid
170.6 178.9 -82.3 164.6 116
Scaled residuals:
Min 1Q Median 3Q Max
-1.013 -1.013 -0.913 0.987 1.095
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 1.73e-14 1.32e-07
Number of obs: 119, groups: ID, 92
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.0267 0.2310 0.12 0.91
ItimeTRUE -0.2090 0.3808 -0.55 0.58
Correlation of Fixed Effects:
(Intr)
ItimeTRUE -0.607
Test and confidence interval for a difference in proportions
prop.test(Seen.~Itime, data=mtg)
2-sample test for equality of proportions with continuity
correction
data: t(table_from_formula)
X-squared = 0.129, df = 1, p-value = 0.7195
alternative hypothesis: two.sided
95 percent confidence interval:
-0.1515 0.2558
sample estimates:
prop 1 prop 2
0.5455 0.4933
with(data=mtg,prop.table(table(Seen., Itime), margin=2))
Itime
Seen. FALSE TRUE
0 0.4933 0.5455
1 0.5067 0.4545
Read in data:
mtgwa<-read.csv("../Data/WAgoats.csv")
Fit full model, with interaction between yr.s.cap & method, then model without the interaction
fit.glmer1<-glmer(observed~grpsize+vegmidpt + terrainob+ method + yr.s.cap + method*yr.s.cap + (1|id), family=binomial(), data=mtgwa, 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 ~ grpsize + vegmidpt + terrainob + method + yr.s.cap +
method * yr.s.cap + (1 | id)
Data: mtgwa
AIC BIC logLik deviance df.resid
97.4 119.1 -40.7 81.4 103
Scaled residuals:
Min 1Q Median 3Q Max
-3.087 0.016 0.264 0.489 1.051
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.228 0.477
Number of obs: 111, groups: id, 52
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.6060 1.2516 1.28 0.199
grpsize 0.2508 0.1295 1.94 0.053 .
vegmidpt -0.0113 0.0109 -1.03 0.302
terrainobY -0.5268 0.6647 -0.79 0.428
methodHeli -0.1484 1.4895 -0.10 0.921
yr.s.cap -0.4628 0.7389 -0.63 0.531
methodHeli:yr.s.cap 0.3297 1.0136 0.32 0.745
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) grpsiz vgmdpt trrnbY mthdHl yr.s.c
grpsize -0.231
vegmidpt -0.071 0.090
terrainobY -0.438 -0.132 -0.161
methodHeli -0.601 -0.108 -0.221 0.167
yr.s.cap -0.766 -0.034 -0.087 0.131 0.621
mthdHl:yr.. 0.540 0.148 0.145 -0.146 -0.889 -0.743
Drop interaction
fit.glmer2<-glmer(observed~grpsize+vegmidpt + method + terrainob+ yr.s.cap + (1|id), family=binomial(), data=mtgwa, 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 ~ grpsize + vegmidpt + method + terrainob + yr.s.cap +
(1 | id)
Data: mtgwa
AIC BIC logLik deviance df.resid
95.5 114.5 -40.8 81.5 104
Scaled residuals:
Min 1Q Median 3Q Max
-3.253 0.017 0.273 0.507 1.105
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.189 0.435
Number of obs: 111, groups: id, 52
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.3949 1.0233 1.36 0.173
grpsize 0.2452 0.1265 1.94 0.052 .
vegmidpt -0.0118 0.0107 -1.10 0.270
methodHeli 0.2835 0.6767 0.42 0.675
terrainobY -0.4988 0.6522 -0.76 0.444
yr.s.cap -0.2871 0.4864 -0.59 0.555
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) grpsiz vgmdpt mthdHl trrnbY
grpsize -0.370
vegmidpt -0.179 0.068
methodHeli -0.300 0.050 -0.204
terrainobY -0.434 -0.119 -0.142 0.089
yr.s.cap -0.632 0.111 0.029 -0.163 0.022
Drop other non-significant variables (control only for groupsize) to test for robustness. First, include interaction
fit.glmer3<-glmer(observed~grpsize+ method + yr.s.cap + method*yr.s.cap+(1|id), family=binomial(), data=mtgwa, 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 ~ grpsize + method + yr.s.cap + method * yr.s.cap +
(1 | id)
Data: mtgwa
AIC BIC logLik deviance df.resid
95.6 111.8 -41.8 83.6 105
Scaled residuals:
Min 1Q Median 3Q Max
-3.189 0.013 0.291 0.521 0.805
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.108 0.329
Number of obs: 111, groups: id, 52
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.005 1.057 0.95 0.342
grpsize 0.261 0.125 2.10 0.036 *
methodHeli -0.285 1.399 -0.20 0.839
yr.s.cap -0.483 0.692 -0.70 0.485
methodHeli:yr.s.cap 0.398 0.955 0.42 0.677
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) grpsiz mthdHl yr.s.c
grpsize -0.309
methodHeli -0.641 -0.078
yr.s.cap -0.815 -0.032 0.608
mthdHl:yr.. 0.552 0.141 -0.888 -0.726
Now, drop interaction
fit.glmer3b<-glmer(observed~grpsize+ method + yr.s.cap +(1|id), family=binomial(), data=mtgwa, nAGQ=10 )
summary(fit.glmer3b)
Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
Family: binomial ( logit )
Formula: observed ~ grpsize + method + yr.s.cap + (1 | id)
Data: mtgwa
AIC BIC logLik deviance df.resid
93.7 107.3 -41.9 83.7 106
Scaled residuals:
Min 1Q Median 3Q Max
-3.433 0.015 0.287 0.536 0.741
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.095 0.308
Number of obs: 111, groups: id, 52
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.768 0.854 0.90 0.368
grpsize 0.255 0.123 2.08 0.038 *
methodHeli 0.236 0.641 0.37 0.713
yr.s.cap -0.277 0.469 -0.59 0.554
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) grpsiz mthdHl
grpsize -0.477
methodHeli -0.372 0.106
yr.s.cap -0.709 0.104 -0.157
Now, method and yr.s.cap as only predictors (along with groupsize)
fit.glmer4<-glmer(observed~grpsize+ method + (1|id), family=binomial(), data=mtgwa, 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 ~ grpsize + method + (1 | id)
Data: mtgwa
AIC BIC logLik deviance df.resid
92.1 102.9 -42.1 84.1 107
Scaled residuals:
Min 1Q Median 3Q Max
-3.1374 0.0131 0.2939 0.5505 0.6739
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.126 0.355
Number of obs: 111, groups: id, 52
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.413 0.606 0.68 0.496
grpsize 0.265 0.125 2.12 0.034 *
methodHeli 0.183 0.638 0.29 0.774
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) grpsiz
grpsize -0.573
methodHeli -0.690 0.126
fit.glmer5<-glmer(observed~grpsize+ yr.s.cap + (1|id), family=binomial(), data=mtgwa, 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 ~ grpsize + yr.s.cap + (1 | id)
Data: mtgwa
AIC BIC logLik deviance df.resid
91.9 102.7 -41.9 83.9 107
Scaled residuals:
Min 1Q Median 3Q Max
-3.280 0.015 0.302 0.561 0.729
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.0255 0.16
Number of obs: 111, groups: id, 52
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.887 0.787 1.13 0.260
grpsize 0.251 0.121 2.08 0.037 *
yr.s.cap -0.250 0.459 -0.54 0.586
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) grpsiz
grpsize -0.465
yr.s.cap -0.840 0.122
Yr.s.cap is not even close to significant in any of the models. Construct plots to look at effect size.
AK mountain goats: basic proportions and confidence intervals for proportion seen vs. time since cap
mtg.ak<-summaryBy(Seen.~Itime, FUN=c(mean, SE),data=mtg)
Model based predictions for Alaskan mtgs
xmat<-matrix(rbind(c(1,0), c(1,1)), 2,2)
mt.ak.prd<- matrix(NA,2,3)
mt.ak.prd[,1]<-plogis(xmat%*%fixef(fit.mtg3))
se.prd<-diag(xmat%*%vcov(fit.mtg3)%*% t(xmat))
mt.ak.prd[,3]<- plogis(xmat%*%fixef(fit.mtg3)+1.96*se.prd)
mt.ak.prd[,2]<- plogis(xmat%*%fixef(fit.mtg3)-1.96*se.prd)
Now for WA mtgs. Use fit.glmer5, with grpsize and yr.s.cap (lowest AIC of models containing yr.s.cap) Get predictions for median group size (which varied a bit by yr.s.cap, but not that much)
tapply(mtgwa$grpsize, mtgwa$yr.s.cap, mean)
0 1 2 3
7.000 9.065 8.821 46.000
tapply(mtgwa$grpsize, mtgwa$yr.s.cap, median)
0 1 2 3
6 4 3 46
mgs<-median(mtgwa$grpsize)
xmat<-matrix(rbind(c(1,mgs,0), c(1,mgs,1), c(1,mgs,2)), 3,3)
mt.wa.prd<- matrix(NA,3,3)
mt.wa.prd[,1]<-plogis(xmat%*%fixef(fit.glmer5))
se.prd.wa<-diag(xmat%*%vcov(fit.glmer5)%*% t(xmat))
mt.wa.prd[,3]<- plogis(xmat%*%fixef(fit.glmer5)+1.96*se.prd.wa)
mt.wa.prd[,2]<- plogis(xmat%*%fixef(fit.glmer5)-1.96*se.prd.wa)
Now, get proportion observed as function of time since cap
mtgwa$Iobs<-I(mtgwa$observed=="Y")
mtg.wa.ob<-summaryBy(Iobs~yr.s.cap,FUN=c(mean, SE), data=mtgwa)
# Now, plot with helicopter/ground as well for washington data, using glmer3 (with interaction) between method & yr since cap
xmat<-matrix(rbind(c(1,mgs,0,0,0),c(1,mgs,0,1,0),c(1,mgs,0,2,0),
c(1,mgs,1,0,1),c(1,mgs,1,1,1),c(1,mgs,1,2,1)), 6,5)
mt.wa.prd<- matrix(NA,6,3)
mt.wa.prd[,1]<-plogis(xmat%*%fixef(fit.glmer3))
se.prd.wa<-diag(xmat%*%vcov(fit.glmer3)%*% t(xmat))
mt.wa.prd[,3]<- plogis(xmat%*%fixef(fit.glmer3)+1.96*se.prd.wa)
mt.wa.prd[,2]<- plogis(xmat%*%fixef(fit.glmer3)-1.96*se.prd.wa)
par(mfrow=c(1,1), bty="L", cex.lab=1.8, cex=1, cex.axis=1.8, mar=c(6,6,4.1,2.1) , xpd=T, oma=c(0,1,0,0))
plotCI(c(1,2), mt.ak.prd[,1], ui=mt.ak.prd[,3], li=mt.ak.prd[,2], xlab="Years since capture", ylab="Probability of detection", ylim=c(0,1), xlim=c(0,2), gap=0, pch=17, cex=2,xaxt="n")
axis(side=1, c(0,1,2))
plotCI(c(c(0,1,2)-.02, c(0,1,2)+.02), mt.wa.prd[,1], ui=mt.wa.prd[,3], li=mt.wa.prd[,2], xlab="", ylab="",add=T, ylim=c(0,1), xlim=c(0,2), gap=0, pch=c(rep(16,3), rep(1,3)),cex=2, xaxt="n")
legend("bottomright", pch=c(17,16,1), cex=2,lty=1, title="", legend=c("Alaska Helicopter/Helicopter", "Washington Ground/Helicopter", "Washington Helicopter/Helicopter"), bty="n")
For WA, where there were no missing values, also calculate model-averaged coefficients to get at effect sizes (use all models except those with interactions between method and years since capture; these models had larger AIC's anyway)
Cands <- list(fit.glmer2, fit.glmer3b, fit.glmer4, fit.glmer5)
model selection
aictab(cand.set = Cands)
Warning:
Model names have been supplied automatically in the table
Model selection based on AICc :
K AICc Delta_AICc AICcWt Cum.Wt LL
Mod4 4 92.27 0.00 0.42 0.42 -41.94
Mod3 4 92.48 0.21 0.38 0.80 -42.05
Mod2 5 94.32 2.06 0.15 0.95 -41.87
Mod1 7 96.59 4.32 0.05 1.00 -40.75
model average for difference between period 1 vs period 4
summary(mod.avg(fit.glmer2, fit.glmer3b, fit.glmer4, fit.glmer5))
Error: error in evaluating the argument 'object' in selecting a method for function 'summary': Error: could not find function "mod.avg"
Effect sizes:
c(exp(-0.161), exp(-0.161-1.96*0.46395), exp(-0.161+1.96*0.46395))
[1] 0.8513 0.3429 2.1135
c(exp(0.1188361),exp(0.1188361-1.96*0.4996),exp(0.1188361+1.96*0.4996))
[1] 1.126 0.423 2.998
# 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] RMark_2.1.8 MuMIn_1.10.5 AICcmodavg_2.00 gplots_2.14.1
[5] knitr_1.6 mosaic_0.9.1-3 lattice_0.20-29 dplyr_0.2
[9] car_2.0-21 doBy_4.5-10 MASS_7.3-33 survival_2.37-7
[13] nlme_3.1-117 ggplot2_1.0.0 lme4_1.1-7 Rcpp_0.11.2
[17] Matrix_1.1-4
loaded via a namespace (and not attached):
[1] assertthat_0.1 bitops_1.0-6 caTools_1.17
[4] coda_0.16-1 colorspace_1.2-4 digest_0.6.4
[7] evaluate_0.5.5 expm_0.99-1.1 formatR_1.0
[10] gdata_2.13.3 ggdendro_0.1-14 grid_3.1.1
[13] gridExtra_0.9.1 gtable_0.1.2 gtools_3.4.1
[16] htmltools_0.2.4 KernSmooth_2.23-12 labeling_0.3
[19] markdown_0.7.4 matrixcalc_1.0-3 mime_0.1.2
[22] minqa_1.2.3 mosaicData_0.9.1 msm_1.4
[25] munsell_0.4.2 mvtnorm_1.0-0 nloptr_1.0.4
[28] nnet_7.3-8 parallel_3.1.1 plyr_1.8.1
[31] proto_0.3-10 reshape2_1.4 rmarkdown_0.2.49
[34] scales_0.2.4 snowfall_1.84-6 stringr_0.6.2
[37] tools_3.1.1 yaml_2.1.13