Golden-cheeked warbler mark-recapture distance sampling analysis
These data represent avian point count surveys conducted at 453 point sample survey locations on the 24,000 (approx) live-fire region of Fort Hood in central Texas. Surveys were conducted by independent double observers (2 per survey occasion) and as such we had a maximum of 3 paired survey histories, giving a maximum of 6 sample occasions (see MacKenzie et al. 2006, MacKenzie and Royle 2005, and Laake et al. 2011 for various sample survey design details). At each point, we surveyed for 5 minutes (technically broken into 3 time intervals of 2, 2, and 1 minutes; not used here) and we noted detections by each observer and collected distance to each observation within a set of distance bins (0-50, 50-100m; Laake et al. 2011) of the target species (Golden-cheeked warblers in this case) for each surveyor. Our primary focus was to use mark-recapture distance sampling methods to estimate density of Golden-cheeked warblers, and to estimate detection rates for the mark-recapture, distance, and composite model.
The format is a data frame with the following covariate metrics.
Unique identifier for each sample location; locations are the same for both species
Visit number to the point
Species designation, either Golden-cheeked warbler (GW) or Black-capped Vireo (BV)
Distance measure, which is either NA (representing no detection), or the median of the binned detection distances
ID value indicating which observers were paired for that sampling occasion
Observer ID, either primary(1), or secondary (2)
Detection of a bird, either 1 = detected, or 0 = not detected
Date of survey since 15 March 2011, numeric value
Predicted occupancy value for that survey hexagon based on Farrell et al. (2013)
Region.Label categorization, see R package mrds help file for details on data structure
Amount of survey effort at the point
Number of days since 15 March 2011, numeric value
Unique ID for each paired observations
In addition to detailing the analysis used by Collier et al. (2013, In Review), this example documents the use of mrds for avian point count surveys and shows how density models
can be incorporated with occupancy models to develop spatially explicit density surface maps. For those that are interested, for the distance sampling portion of our analysis, we
used both conventional distance sampling (cds) and multiple covariate distance sampling (mcds) with uniform and half-normal key functions. For the mark-recapture portion of
our analysis, we tended to use covariates for distance (median bin width), observer, and date of survey (days since 15 March 2011).
We combined our mrds density estimates via a Horvitz-Thompson styled estimator with the resource selection function gradient developed
in Farrell et al. (2013) and estimated density on an ~3.14ha hexagonal grid across our study area, which provided a density gradient for Fort Hood. Because there
was considerable data manipulation needed for each analysis to structure the data appropriately for use in mrds, rather than wrap each analysis in a single function, we have provided
both the Golden-cheeked warbler and Black-capped vireo analyses in their full detail. The primary differences you will see will be changes to
model structures and model outputs between the two species.
Bret Collier and Jeff Laake
Farrell, S.F., B.A. Collier, K.L. Skow, A.M. Long, A.J. Campomizzi, M.L. Morrison, B. Hays, and R.N. Wilkins. 2013. Using LiDAR-derived structural vegetation characteristics to develop high-resolution, small-scale, species distribution models for conservation planning. Ecosphere 43(3): 42. http://dx.doi.org/10.1890/ES12-000352.1
Laake, J.L., B.A. Collier, M.L. Morrison, and R.N. Wilkins. 2011. Point-based mark recapture distance sampling. Journal of Agricultural, Biological and Environmental Statistics 16: 389-408.
Collier, B.A., S.L. Farrell, K.L. Skow, A.M. Long, A.J. Campomizzi, K.B. Hays, J.L. Laake, M.L. Morrison, and R.N. Wilkins. 2013. Spatially explicit density of endangered avian species in a disturbed landscape. Auk, In Review.
## Not run:
data(lfgcwa)
xy <- cut(lfgcwa$Pred, c(-0.0001, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1),
labels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
x <- data.frame(lfgcwa, New=xy)
# Note that I scaled the individual covariate of day-helps with
# convergence issues
bird.data <- data.frame(object=x$ObjectID, observer=x$Observer,
detected=x$Detected, distance=x$Distance,
Region.Label=x$New, Sample.Label=x$PointID,
Day=(x$Day/max(x$Day)))
# make observer a factor variable
bird.data$observer=factor(bird.data$observer)
# Jeff Laake suggested this snippet to quickly create distance medians
# which adds bin information to the \code{bird.data} dataframe
bird.data$distbegin=0
bird.data$distend=100
bird.data$distend[bird.data$distance==12.5]=50
bird.data$distbegin[bird.data$distance==37.5]=0
bird.data$distend[bird.data$distance==37.5]=50
bird.data$distbegin[bird.data$distance==62.5]=50
bird.data$distend[bird.data$distance==62.5]=100
bird.data$distbegin[bird.data$distance==87.5]=50
bird.data$distend[bird.data$distance==87.5]=100
# Removed all survey points with distance=NA for a survey event;
# hence no observations for use in \code{ddf()} but needed later
bird.data=bird.data[complete.cases(bird.data),]
# Manipulations on full dataset for various data.frame creation
# for use in density estimation using \code{dht()}
# Samples dataframe
xx <- x
x <- data.frame(PointID=x$PointID, Species=x$Species,
Category=x$New, Effort=x$Effort)
x <- x[!duplicated(x$PointID),]
point.num <- table(x$Category)
samples <- data.frame(PointID=x$PointID, Region.Label=x$Category,
Effort=x$Effort)
final.samples=data.frame(Sample.Label=samples$PointID,
Region.Label=samples$Region.Label,
Effort=samples$Effort)
# obs dataframe
obs <- data.frame(ObjectID=xx$ObjectID, PointID=xx$PointID)
# used to get Region and Sample assigned to ObjectID
obs <- merge(obs, samples, by=c("PointID", "PointID"))
obs <- obs[!duplicated(obs$ObjectID),]
obs <- data.frame(object=obs$ObjectID, Region.Label=obs$Region.Label,
Sample.Label=obs$PointID)
#Region.Label dataframe
region.data <- data.frame(Region.Label=c(1,2,3,4,5,6,7,8,9),
Area=c(point.num[1]*3.14,
point.num[2]*3.14,
point.num[3]*3.14,
point.num[4]*3.14,
point.num[5]*3.14,
point.num[6]*3.14,
point.num[7]*3.14,
point.num[8]*3.14,
point.num[9]*3.14))
# Candidate Models
GW1=ddf(
dsmodel=~cds(key="unif", adj.series="cos", adj.order=1,adj.scale="width"),
mrmodel=~glm(~distance),
data=bird.data,
method="io",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW2=ddf(
dsmodel=~cds(key="unif", adj.series="cos", adj.order=1,adj.scale="width"),
mrmodel=~glm(~distance+observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW3=ddf(
dsmodel=~cds(key="unif", adj.series="cos", adj.order=1,adj.scale="width"),
mrmodel=~glm(~distance*observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW4=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance),
data=bird.data,
method="io",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW4FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW5=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance+observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW5FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance+observer),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW6=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW6FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW7=ddf(
dsmodel=~cds(key="hn",formula=~1),
mrmodel=~glm(~distance*Day),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW7FI=ddf(
dsmodel=~cds(key="hn",formula=~1),
mrmodel=~glm(~distance*Day),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW8=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer*Day),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW8FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer*Day),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
#GWDS=ddf(
# dsmodel=~mcds(key="hn",formula=~1),
# data=bird.data,
# method="ds",
# meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
#### GCWA Summary Metrics
#AIC table building code, not exactly elegant, but I did not
want to add more package dependencies
AIC = c(GW1$criterion, GW2$criterion, GW3$criterion, GW4$criterion,
GW4FI$criterion, GW5$criterion, GW5FI$criterion,
GW6$criterion, GW6FI$criterion, GW7$criterion, GW7FI$criterion,
GW8$criterion, GW8FI$criterion)
#creates a set of row names for me to check my grep() call below
rn <- c("GW1", "GW2", "GW3", "GW4", "GW4FI", "GW5", "GW5FI", "GW6",
"GW6FI", "GW7","GW7FI", "GW8", "GW8FI")
# number of parameters for each model
k <- c(length(GW1$par), length(GW2$par), length(GW3$par), length(GW4$par),
length(GW4FI$par), length(GW5$par), length(GW5FI$par),
length(GW6$par), length(GW6FI$par), length(GW7$par),
length(GW7FI$par), length(GW8$par), length(GW8FI$par))
# build AIC table and
AIC.table <- data.frame(AIC = AIC, rn=rn, k=k, dAIC = abs(min(AIC)-AIC),
likg = exp(-.5*(abs(min(AIC)-AIC))))
# row.names(AIC.table)=grep("GW", ls(), value=TRUE)
AIC.table <- AIC.table[with(AIC.table, order(-likg, -dAIC, AIC, k)),]
AIC.table <- data.frame(AIC.table, wi=AIC.table$likg/sum(AIC.table$likg))
AIC.table
# Model average N_hat_covered estimates
# not very clean, but I wanted to show full process, need to use
# collect.models and model.table here
estimate <- c(GW1$Nhat, GW2$Nhat, GW3$Nhat, GW4$Nhat, GW4FI$Nhat,
GW5$Nhat, GW5FI$Nhat, GW6$Nhat, GW6FI$Nhat, GW7$Nhat,
GW7FI$Nhat, GW8$Nhat, GW8FI$Nhat)
AIC.values <- AIC
# Nhat.se is calculated in mrds:::summary.io, not in ddf(), so
# it takes a bit to pull out
std.err <- c(summary(GW1)$Nhat.se, summary(GW2)$Nhat.se,
summary(GW3)$Nhat.se,summary(GW4)$Nhat.se,
summary(GW4FI)$Nhat.se, summary(GW5)$Nhat.se,
summary(GW5FI)$Nhat.se, summary(GW6)$Nhat.se,
summary(GW6FI)$Nhat.se, summary(GW7)$Nhat.se,
summary(GW7FI)$Nhat.se,summary(GW8)$Nhat.se,
summary(GW8FI)$Nhat.se)
## End(Not run)
## Not run:
#Not Run
#requires RMark
library(RMark)
#uses model.average structure to model average real abundance estimates for
#covered area of the surveys
mmi.list=list(estimate=estimate, AIC=AIC.values, se=std.err)
model.average(mmi.list, revised=TRUE)
#Not Run
#Best Model FI
#best.modelFI=AIC.table[1,]
#best.model
#Best Model PI
#best.modelPI=AIC.table[2,]
#best.modelPI
#Not Run
#summary(GW7FI, se=TRUE)
#summary(GW7, se=TRUE)
#Not Run
#GOF for models
#ddf.gof(GW7, breaks=c(0,50,100))
#Not Run
#Density estimation across occupancy categories
#out.GW=dht(GW7, region.data, final.samples, obs, se=TRUE,
options=list(convert.units=.01))
#Plots--Not Run
#Composite Detection Function examples
#plot(GW7, which=3, showpoints=FALSE, angle=0, density=0,
# col="black", lwd=3, main="Golden-cheeked Warbler",
# xlab="Distance (m)", las=1, cex.axis=1.25, cex.lab=1.25)
#Conditional Detection Function
#dd=expand.grid(distance=0:100,Day=(4:82)/82)
#dmat=model.matrix(~distance*Day,dd)
#dd$p=plogis(model.matrix(~distance*Day,dd)%*%coef(GW7$mr)$estimate)
#dd$Day=dd$Day*82
#with(dd[dd$Day==12,],plot(distance,p,ylim=c(0,1), las=1,
# ylab="Detection probability", xlab="Distance (m)",
# type="l",lty=1, lwd=3, bty="l", cex.axis=1.5, cex.lab=1.5))
#with(dd[dd$Day==65,],lines(distance,p,lty=2, lwd=3))
#ch=paste(bird.data$detected[bird.data$observer==1],
# bird.data$detected[bird.data$observer==2],
# sep="")
#tab=table(ch,cut(82*bird.data$Day[bird.data$observer==1],c(0,45,83)),
# cut(bird.data$distance[bird.data$observer==1],c(0,50,100)))
#tabmat=cbind(colMeans(rbind(tab[3,,1]/colSums(tab[2:3,,1],
# tab[3,,1]/colSums(tab[c(1,3),,1])))),
# colMeans(rbind(tab[3,,2]/colSums(tab[2:3,,2],
# tab[3,,2]/colSums(tab[c(1,3),,2])))))
#colnames(tabmat)=c("0-50","51-100")
#points(c(25,75),tabmat[1,],pch=1, cex=1.5)
#points(c(25,75),tabmat[2,],pch=2, cex=1.5)
# Another alternative plot using barplot instead of points
# (this is one in paper)
#ch=paste(bird.data$detected[bird.data$observer==1],
# bird.data$detected[bird.data$observer==2],
#sep="")
#tab=table(ch,cut(82*bird.data$Day[bird.data$observer==1],c(0,45,83)),
# cut(bird.data$distance[bird.data$observer==1],c(0,50,100)))
#tabmat=cbind(colMeans(rbind(tab[3,,1]/colSums(tab[2:3,,1],
# tab[3,,1]/colSums(tab[c(1,3),,1])))),
#colMeans(rbind(tab[3,,2]/colSums(tab[2:3,,2],
# tab[3,,2]/colSums(tab[c(1,3),,2])))))
#colnames(tabmat)=c("0-50","51-100")
#par(mfrow=c(2, 1), mai=c(1,1,1,1))
#with(dd[dd$Day==12,],
# plot(distance,p,ylim=c(0,1), las=1,
# ylab="Detection probability", xlab="",
# type="l",lty=1, lwd=4, bty="l", cex.axis=1.5, cex.lab=1.5))
#segments(0, 0, .0, tabmat[1,1], lwd=3)
#segments(0, tabmat[1,1], 50, tabmat[1,1], lwd=4)
#segments(50, tabmat[1,1], 50, 0, lwd=4)
#segments(50, tabmat[1,2], 100, tabmat[1,2], lwd=4)
#segments(0, tabmat[1,1], 50, tabmat[1,1], lwd=4)
#segments(100, tabmat[1,2], 100, 0, lwd=4)
#mtext("a",line=-1, at=90)
#with(dd[dd$Day==65,],
# plot(distance,p,ylim=c(0,1), las=1, ylab="Detection probability",
# xlab="Distance", type="l",lty=1,
# lwd=4, bty="l", cex.axis=1.5, cex.lab=1.5))
#segments(0, 0, .0, tabmat[2,1], lwd=4)
#segments(0, tabmat[2,1], 50, tabmat[2,1], lwd=4)
#segments(50, tabmat[2,1], 50, 0, lwd=4)
#segments(50, tabmat[2,2], 50, tabmat[2,1], lwd=4)
#segments(50, tabmat[2,2], 100, tabmat[2,2], lwd=4)
#segments(100, tabmat[2,2], 100, 0, lwd=4)
#mtext("b",line=-1, at=90)
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.