Black-capped vireo 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-25, 25-50, 50-75, 75-100m) of the target species (Black-capped vireo's in this case) for each surveyor. Our primary focus was to use mark-recapture distance sampling methods to estimate density of Black-capped vireo's, 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
Predicted occupancy value for that survey hexagon based on Farrell et al. (2013)
Region.Label categorization, see mrds help file for details on data structure
Amount of survey effort at the point
Number of days since 15 March 2011
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 the Fort Hood military installation. 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(lfbcvi)
xy=cut(lfbcvi$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(lfbcvi, 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 bird.data dataframe
bird.data$distbegin=0
bird.data$distend=100
bird.data$distend[bird.data$distance==12.5]=25
bird.data$distbegin[bird.data$distance==37.5]=25
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]=75
bird.data$distbegin[bird.data$distance==87.5]=75
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 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 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.data=data.frame(Region.Label=c(1, 2, 3,4,5,6,7,8,9, 10),
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, point.num[10]*3.14))
# Candidate Models
BV1=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)))
BV1FI=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)))
BV2=ddf(
dsmodel=~mcds(key="hr",formula=~1),
mrmodel=~glm(~distance),
data=bird.data,
method="io",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
BV3=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)))
BV3FI=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)))
BV4=ddf(
dsmodel=~mcds(key="hr",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)))
BV5=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)))
BV5FI=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)))
BV6=ddf(
dsmodel=~mcds(key="hr",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)))
BV7=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)))
BV7FI=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)))
BV8=ddf(
dsmodel=~cds(key="hr",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)))
BV9=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)))
BV9FI=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)))
BV10=ddf(
dsmodel=~mcds(key="hr",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)))
#BV.DS=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)))
#AIC table building code.
AIC = c(BV1$criterion, BV1FI$criterion, BV2$criterion, BV3$criterion,
BV3FI$criterion, BV4$criterion, BV5$criterion, BV5FI$criterion,
BV6$criterion, BV7$criterion, BV7FI$criterion, BV8$criterion,
BV9$criterion, BV9FI$criterion, BV10$criterion)
#creates a set of row names for me to check my grep() call below
rn = c("BV1", "BV1FI", "BV2", "BV3", "BV3FI", "BV4", "BV5", "BV5FI",
"BV6", "BV7", "BV7FI", "BV8", "BV9", "BV9FI", "BV10")
#Number parameters
k = c(length(BV1$par), length(BV1FI$par), length(BV2$par),
length(BV3$par), length(BV3FI$par), length(BV4$par),
length(BV5$par),length(BV5FI$par), length(BV6$par),
length(BV7$par), length(BV7FI$par), length(BV8$par),
#build AIC table
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("BV", 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 later on
estimate <- c(BV1$Nhat, BV1FI$Nhat, BV2$Nhat, BV3$Nhat, BV3FI$Nhat,
BV4$Nhat, BV5$Nhat, BV5FI$Nhat, BV6$Nhat, BV7$Nhat,
BV7FI$Nhat, BV8$Nhat, BV9$Nhat, BV9FI$Nhat, BV10$Nhat)
AIC.values=AIC
# had to use str() to extract here as Nhat.se is calculated in
# mrds:::summary.io, not in ddf(), so it takes a bit
std.err <- c(summary(BV1)$Nhat.se, summary(BV1FI)$Nhat.se,
summary(BV2)$Nhat.se, summary(BV3)$Nhat.se,
summary(BV3FI)$Nhat.se, summary(BV4)$Nhat.se,
summary(BV5)$Nhat.se, summary(BV5FI)$Nhat.se,
summary(BV6)$Nhat.se, summary(BV7)$Nhat.se,
summary(BV7FI)$Nhat.se,summary(BV8)$Nhat.se,
summary(BV9)$Nhat.se, summary(BV9FI)$Nhat.se,
summary(BV10)$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
#Summary for the top 2 models
#summary(BV5, se=TRUE)
#summary(BV5FI, se=TRUE)
#Not Run
#Best Model
#best.model=AIC.table[1,]
#Not Run
#GOF for models
#ddf.gof(BV5, breaks=c(0, 25, 50, 75, 100))
#Not Run
#Density estimation across occupancy categories
#out.BV=dht(BV5, region.data, final.samples, obs, se=TRUE,
# options=list(convert.units=.01))
#Plot--Not Run
#Composite Detection Function
#plot(BV5, which=3, showpoints=FALSE, angle=0, density=0, col="black", lwd=3,
# main="Black-capped Vireo",xlab="Distance (m)", las=1, cex.axis=1.25,
# cex.lab=1.25)
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.