mrds (version 2.2.4)

lfgcwa: Golden-cheeked warbler mark-recapture distance sampling analysis

Description

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.

Format

The format is a data frame with the following covariate metrics.

VisitNumber

Visit number to the point

Species

Species designation, either Golden-cheeked warbler (GW) or Black-capped Vireo (BV)

Distance

Distance measure, which is either NA (representing no detection), or the median of the binned detection distances

PairNumber

ID value indicating which observers were paired for that sampling occasion

Observer

Observer ID, either primary(1), or secondary (2)

Detected

Detection of a bird, either 1 = detected, or 0 = not detected

Date

Date of survey since 15 March 2011, numeric value

Pred

Predicted occupancy value for that survey hexagon based on Farrell et al. (2013)

Category

Region.Label categorization, see R package mrds help file for details on data structure

Effort

Amount of survey effort at the point

Day

Number of days since 15 March 2011, numeric value

ObjectID

Unique ID for each paired observations

Details

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.

References

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.

Examples

# NOT RUN {
# }
# 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)
# }
# 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)
# }
# NOT RUN {
# }