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?
Steps:
rm(list = ls())
# Load libraries
library(RMark)
library(knitr)
library(lme4)
library(ggplot2)
library(doBy)
# set global option for output
opts_chunk$set(comment=NA, fig.width=6, fig.height=6)
MSH<-import.chdata("../Data/MSH.inp", header=F, field.names=c("ch", "AdF", "BA"))
head(MSH)
ch AdF BA
533 149.550 COWEEMAN 11...... 1 0
535 150.110 COWEEMAN 111110.. 1 0
505 150.390 COWEEMAN 11110111 1 0
537 150.450 COWEEMAN 0000.... 0 1
536 150.030 COWEEMAN 10...... 1 0
525 150.260 COWEEMAN 11100000 1 0
There are 8 capture histories (some of which are missing), create a “long” data set from this file.
cdat.MSH<-NULL
for(i in 1:nrow(MSH)){
temp<-MSH[i,]
first.cap=NA
for(j in 1:8){
temp2<-as.numeric(substr(temp$ch,j,j))
if(is.na(temp2)==F){
if(is.na(first.cap)==T){first.cap<-(j+1)%/%2}#first.cap<-j note, first captures were all in week 1
cdat.MSH<-rbind(cdat.MSH, data.frame(id=rownames(temp), ID=i, year=(j+1)%/%2, week=(j+1)%%2+1, observed=temp2, first.cap=first.cap, IAdF=temp$AdF))}
}
}
# Time since capture (assumes week=1 or 2 does not matter...)
cdat.MSH$time.cap<-cdat.MSH$year-cdat.MSH$first.cap
head(cdat.MSH, 20)
id ID year week observed first.cap IAdF time.cap
1 533 149.550 COWEEMAN 1 1 1 1 1 1 0
2 533 149.550 COWEEMAN 1 1 2 1 1 1 0
3 535 150.110 COWEEMAN 2 1 1 1 1 1 0
4 535 150.110 COWEEMAN 2 1 2 1 1 1 0
5 535 150.110 COWEEMAN 2 2 1 1 1 1 1
6 535 150.110 COWEEMAN 2 2 2 1 1 1 1
7 535 150.110 COWEEMAN 2 3 1 1 1 1 2
8 535 150.110 COWEEMAN 2 3 2 0 1 1 2
9 505 150.390 COWEEMAN 3 1 1 1 1 1 0
10 505 150.390 COWEEMAN 3 1 2 1 1 1 0
11 505 150.390 COWEEMAN 3 2 1 1 1 1 1
12 505 150.390 COWEEMAN 3 2 2 1 1 1 1
13 505 150.390 COWEEMAN 3 3 1 0 1 1 2
14 505 150.390 COWEEMAN 3 3 2 1 1 1 2
15 505 150.390 COWEEMAN 3 4 1 1 1 1 3
16 505 150.390 COWEEMAN 3 4 2 1 1 1 3
17 537 150.450 COWEEMAN 4 1 1 0 1 0 0
18 537 150.450 COWEEMAN 4 1 2 0 1 0 0
19 537 150.450 COWEEMAN 4 2 1 0 1 0 1
20 537 150.450 COWEEMAN 4 2 2 0 1 0 1
Lets just plot the proportion seen versus year since capture.
par(bty="L")
plot(0:3,prop.table(table(cdat.MSH$observed, cdat.MSH$time.cap), 2)[2,], ylim=c(0,1), xlab="Years since capture", ylab="Proportion observed")
No apparent relationship, but to be sure, lets fit a few mixed binary regression models.
First, look at a model with Year (as a factor), Week, ID (random)+time since capture + AdF
fit.temp1a<-glmer(observed~time.cap+as.factor(week)+as.factor(year)+(1|id)+IAdF, family=binomial(), data=cdat.MSH)
summary(fit.temp1a)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: observed ~ time.cap + as.factor(week) + as.factor(year) + (1 |
id) + IAdF
Data: cdat.MSH
AIC BIC logLik deviance df.resid
782.9 817.6 -383.5 766.9 558
Scaled residuals:
Min 1Q Median 3Q Max
-1.696 -0.895 0.596 0.875 1.361
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.507 0.712
Number of obs: 566, groups: id, 142
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1291 0.3220 -0.40 0.689
time.cap 0.0505 0.1274 0.40 0.692
as.factor(week)2 -0.0964 0.1791 -0.54 0.590
as.factor(year)2 -0.1211 0.2972 -0.41 0.684
as.factor(year)3 -0.0568 0.3173 -0.18 0.858
as.factor(year)4 -0.4900 0.3513 -1.40 0.163
IAdF1 0.4502 0.2662 1.69 0.091 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tim.cp as.fctr(w)2 as.fctr(y)2 as.()3 as.()4
time.cap 0.197
as.fctr(w)2 -0.276 0.001
as.fctr(y)2 -0.537 -0.259 0.001
as.fctr(y)3 -0.573 -0.483 0.000 0.637
as.fctr(y)4 -0.545 -0.570 0.004 0.609 0.724
IAdF1 -0.640 -0.159 -0.005 0.026 0.081 0.076
Not much there, but perhaps a significant sex effect (AdF). Also, there is quite a bit of variation from individual to individual. Lets try modeling year and week as random
fit.temp1b<-glmer(observed~time.cap+(1 | year/week) +(1|id)+IAdF, family=binomial(), data=cdat.MSH)
summary(fit.temp1b)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: observed ~ time.cap + (1 | year/week) + (1 | id) + IAdF
Data: cdat.MSH
AIC BIC logLik deviance df.resid
782.6 808.7 -385.3 770.6 560
Scaled residuals:
Min 1Q Median 3Q Max
-1.495 -0.880 0.593 0.854 1.488
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.505 0.71
week:year (Intercept) 0.000 0.00
year (Intercept) 0.000 0.00
Number of obs: 566, groups: id, 142; week:year, 8; year, 4
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.3185 0.2334 -1.36 0.172
time.cap -0.0393 0.1024 -0.38 0.701
IAdF1 0.4748 0.2641 1.80 0.072 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tim.cp
time.cap -0.185
IAdF1 -0.826 -0.137
Lets try dropping the year and week effects, which have small variance components.
Lets also consider a random slope effect for time since capture. I.e., are there some individuals that
happen to respond behaviorally more than others?
fit.temp2<-glmer(observed~time.cap+(1+time.cap|id)+IAdF, family=binomial(), data=cdat.MSH)
summary(fit.temp2)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: observed ~ time.cap + (1 + time.cap | id) + IAdF
Data: cdat.MSH
AIC BIC logLik deviance df.resid
782.6 808.7 -385.3 770.6 560
Scaled residuals:
Min 1Q Median 3Q Max
-1.478 -0.885 0.615 0.856 1.486
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 0.535645 0.732
time.cap 0.000574 0.024 -1.00
Number of obs: 566, groups: id, 142
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.3188 0.2344 -1.36 0.174
time.cap -0.0372 0.1026 -0.36 0.717
IAdF1 0.4740 0.2649 1.79 0.074 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tim.cp
time.cap -0.189
IAdF1 -0.824 -0.141
Conculsion regarding time.cap is robust to how we look at the data. Lets construct a plot to look at individual-level variation in predicted probabilities of detection. Appears to be significant variation among individuals in their propensity to be detected.
cdat.MSH$phat<-predict(fit.temp2, type="resp")
qplot(time.cap, phat, data=cdat.MSH, geom="line", group=id, colour=IAdF, xlab="Years since capture", ylab="Predicted probability of being observed")
# Look at distribution of p_observed across individuals
pobs<-summaryBy(IAdF+observed~id, data=cdat.MSH)
pobs[,2]<-pobs[,2]-1
pobs[,2]<-as.character(pobs[,2])
ggplot(pobs, aes(observed.mean))+geom_histogram(aes(y=..density..), binwidth=0.01)+facet_grid(IAdF.mean~.)
Warning: position_stack requires constant width: output may be incorrect
Warning: position_stack requires constant width: output may be incorrect
Nook<-import.chdata("../Data/Nook.inp", header=F, field.names=c("ch", "Female", "Male"))
head(Nook)
ch Female Male
149.060 110001111100.... 1 0
149.150 1100............ 1 0
149.530 11.............. 1 0
149.570 110001111100.... 1 0
149.820 01100111111010.. 1 0
149.850 01111110111110.. 1 0
There are 8 capture histories (some of which are missing), create a “long” data set from this file.
cdat.NS<-NULL
for(i in 1:nrow(Nook)){
temp<-Nook[i,]
first.cap=NA
for(j in 1:16){
temp2<-as.numeric(substr(temp$ch,j,j))
if(is.na(temp2)==F){
if(is.na(first.cap)==T){first.cap<-(j+1)%/%2}#first.cap<-j note, first captures were all in week 1
cdat.NS<-rbind(cdat.NS, data.frame(id=rownames(temp), ID=i, year=(j+1)%/%2, week=(j+1)%%2+1, observed=temp2, first.cap=first.cap, Female=temp$Female))}
}
}
# Time since capture (assumes week=1 or 2 does not matter...)
cdat.NS$time.cap<-cdat.NS$year-cdat.NS$first.cap
cdat.NS$Female<-as.factor(cdat.NS$Female)
head(cdat.NS, 20)
id ID year week observed first.cap Female time.cap
1 149.060 1 1 1 1 1 1 0
2 149.060 1 1 2 1 1 1 0
3 149.060 1 2 1 0 1 1 1
4 149.060 1 2 2 0 1 1 1
5 149.060 1 3 1 0 1 1 2
6 149.060 1 3 2 1 1 1 2
7 149.060 1 4 1 1 1 1 3
8 149.060 1 4 2 1 1 1 3
9 149.060 1 5 1 1 1 1 4
10 149.060 1 5 2 1 1 1 4
11 149.060 1 6 1 0 1 1 5
12 149.060 1 6 2 0 1 1 5
13 149.150 2 1 1 1 1 1 0
14 149.150 2 1 2 1 1 1 0
15 149.150 2 2 1 0 1 1 1
16 149.150 2 2 2 0 1 1 1
17 149.530 3 1 1 1 1 1 0
18 149.530 3 1 2 1 1 1 0
19 149.570 4 1 1 1 1 1 0
20 149.570 4 1 2 1 1 1 0
Lets just plot the proportion seen versus year since capture.
par(bty="L")
plot(0:7,prop.table(table(cdat.NS$observed, cdat.NS$time.cap), 2)[2,], ylim=c(0,1), xlab="Years since capture", ylab="Proportion observed")
Nothing too strong there…
First, look at a model with Year (as a factor), Week, ID (random)+time since capture + sex
fit.temp1a<-glmer(observed~time.cap+as.factor(week)+as.factor(year)+(1|id)+Female, family=binomial(), data=cdat.NS)
Warning: Model failed to converge with max|grad| = 0.00455583 (tol =
0.001, component 11)
summary(fit.temp1a)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: observed ~ time.cap + as.factor(week) + as.factor(year) + (1 |
id) + Female
Data: cdat.NS
AIC BIC logLik deviance df.resid
1210.4 1268.3 -593.2 1186.4 908
Scaled residuals:
Min 1Q Median 3Q Max
-1.664 -0.840 -0.448 0.885 2.760
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.306 0.553
Number of obs: 920, groups: id, 121
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3392 0.3083 -4.34 1.4e-05 ***
time.cap 0.0766 0.0575 1.33 0.1826
as.factor(week)2 -0.2528 0.1422 -1.78 0.0753 .
as.factor(year)2 0.7576 0.2587 2.93 0.0034 **
as.factor(year)3 0.0273 0.2859 0.10 0.9238
as.factor(year)4 0.5166 0.3133 1.65 0.0991 .
as.factor(year)5 0.6659 0.3460 1.92 0.0543 .
as.factor(year)6 -0.2579 0.3860 -0.67 0.5040
as.factor(year)7 0.1233 0.3709 0.33 0.7396
as.factor(year)8 -0.6050 0.4169 -1.45 0.1467
Female1 1.0949 0.2593 4.22 2.4e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tim.cp as.fctr(w)2 as.fctr(y)2 as.()3 as.()4 as.()5
time.cap 0.288
as.fctr(w)2 -0.210 -0.005
as.fctr(y)2 -0.546 -0.233 -0.015
as.fctr(y)3 -0.509 -0.387 -0.001 0.530
as.fctr(y)4 -0.508 -0.469 -0.008 0.515 0.548
as.fctr(y)5 -0.473 -0.573 -0.010 0.497 0.552 0.581
as.fctr(y)6 -0.488 -0.663 0.003 0.475 0.552 0.585 0.628
as.fctr(y)7 -0.507 -0.661 -0.001 0.489 0.562 0.594 0.636
as.fctr(y)8 -0.460 -0.645 0.005 0.447 0.521 0.554 0.598
Female1 -0.758 -0.313 -0.017 0.178 0.171 0.201 0.188
as.()6 as.()7 as.()8
time.cap
as.fctr(w)2
as.fctr(y)2
as.fctr(y)3
as.fctr(y)4
as.fctr(y)5
as.fctr(y)6
as.fctr(y)7 0.674
as.fctr(y)8 0.635 0.653
Female1 0.248 0.250 0.225
Pretty big year-to-year variation, females also more likely to be seen. Week also marginallky significant. Lets turn year and week into a random effects.
fit.temp1a<-glmer(observed~time.cap+(1|year/week)+(1|id)+Female, family=binomial(), data=cdat.NS)
summary(fit.temp1a)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: observed ~ time.cap + (1 | year/week) + (1 | id) + Female
Data: cdat.NS
AIC BIC logLik deviance df.resid
1211.0 1240.0 -599.5 1199.0 914
Scaled residuals:
Min 1Q Median 3Q Max
-1.737 -0.834 -0.467 0.864 2.804
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.3217 0.567
week:year (Intercept) 0.2207 0.470
year (Intercept) 0.0225 0.150
Number of obs: 920, groups: id, 121; week:year, 16; year, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.2896 0.2674 -4.82 1.4e-06 ***
time.cap 0.0475 0.0501 0.95 0.34
Female1 1.1593 0.2577 4.50 6.8e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tim.cp
time.cap -0.220
Female1 -0.681 -0.272
Now, consider a random slope for time since capture along with the other randome effects Big difference between males and females, and significant animal-to-animal variation. No strong signature wrt time since capture.
fit.temp2<-glmer(observed~time.cap+(1+time.cap|id)+Female+(1|year/week), family=binomial(), data=cdat.NS)
summary(fit.temp2)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula:
observed ~ time.cap + (1 + time.cap | id) + Female + (1 | year/week)
Data: cdat.NS
AIC BIC logLik deviance df.resid
1214.9 1253.5 -599.4 1198.9 912
Scaled residuals:
Min 1Q Median 3Q Max
-1.754 -0.830 -0.463 0.865 2.799
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 0.263788 0.5136
time.cap 0.000499 0.0223 1.00
week:year (Intercept) 0.222075 0.4712
year (Intercept) 0.021419 0.1464
Number of obs: 920, groups: id, 121; week:year, 16; year, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.2755 0.2669 -4.78 1.8e-06 ***
time.cap 0.0462 0.0507 0.91 0.36
Female1 1.1454 0.2565 4.47 8.0e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tim.cp
time.cap -0.226
Female1 -0.686 -0.249
Lets create a plot to look at individual-level variation in predicted probabilities of detectoin. Appears to be significant variation among individuals in their propensity to be detected.
b.id<-ranef(fit.temp2)$id
names(b.id)<-c("Int", "beta.time.cap")
b.id$id<-rownames(b.id)
beta.hat<-fixef(fit.temp2)
cdat.NS.new<-merge(cdat.NS, b.id, by="id", all=T)
cdat.NS.new$phat<-plogis(cdat.NS.new$Int+beta.hat[1]+cdat.NS.new$time.cap*(cdat.NS.new$beta.time.cap+beta.hat[2]))
# Similar plots
qplot(time.cap, phat, data=cdat.NS.new, geom=c("line", "point"), group=id, colour=Female, xlab="Years since capture", ylab="Predicted probability of being observed")
# Look at distribution of p_observed across individuals
pobs<-summaryBy(Female+observed~id, data=cdat.NS)
pobs[,2]<-pobs[,2]-1
pobs[,2]<-as.character(pobs[,2])
ggplot(pobs, aes(observed.mean))+geom_histogram(aes(y=..density..), binwidth=0.01)+facet_grid(Female.mean~.)
Warning: position_stack requires constant width: output may be incorrect
Warning: position_stack requires constant width: output may be incorrect
Now, combine everything into 1 plot
#x11(width=4, height=4)
par(mfrow=c(2,2), bty="L", cex.lab=1.8, cex.axis=1.8, mar=c(5,6,4.1,2.1) , xpd=T, oma=c(0,1,0,0))
with(cdat.MSH,plot(time.cap, phat, type="n", xlab="Years since capture", ylab="Probability of detection"))
uid<-unique(cdat.MSH$id)
luid<-length(uid)
cols<-c("black", "gray")
lwds<-c(2,1)
for(i in 1:luid){
tempdat<-cdat.MSH[cdat.MSH$i==uid[i],]
with(tempdat,lines(time.cap, phat, col=cols[as.numeric(IAdF)], lwd=lwds[as.numeric(IAdF)]))
}
legend("top", inset=-3, col=c("black", "gray"), lty=1, lwd=c(2,1),bty="n", legend=c("No", "Yes"), title="Adult Female", ncol=2, cex=1.4)
mtext(outer=F, line=2, side=3, cex=1.4, "A)", adj=0)
# Look at distribution of p_observed across individuals
pobs.MSH<-with(cdat.MSH, tapply(observed,id,mean))
hist(pobs.MSH, xlab="", ylab="Distribution", main="")
mtext(outer=F, line=2, side=3, cex=1.4, "B) Proportion of times observed", adj=0)
with(cdat.NS.new,plot(time.cap, phat, type="n", xlab="Years since capture", ylab="Probability of Detection"))
uid<-unique(cdat.NS.new$id)
luid<-length(uid)
cols<-c("black", "gray")
lwds<-c(2,1)
for(i in 1:luid){
tempdat<-cdat.NS.new[cdat.NS.new$i==uid[i],]
with(tempdat,lines(time.cap, phat,lwd=lwds[as.numeric(Female)], col=cols[as.numeric(Female)]))
}
mtext(outer=F, line=2, side=3, cex=1.4, "C)", adj=0)
# Look at distribution of p_observed across individuals
pobs.NS<-with(cdat.NS.new, tapply(observed,id,mean))
hist(pobs.NS, xlab="", ylab="Distribution", main="")
mtext(outer=F, line=2, side=3, cex=1.4, "D) Proportion of times observed", adj=0)
# 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