Beginner Point-Transect Analysis in Rdistance"
In Rdistance: Distance-Sampling Analyses for Density and Abundance Estimation

Introduction

This tutorial is a beginner's guide to doing point transect distance-sampling analysis using Rdistance. Topics covered include input data requirements, fitting a detection function, estimating abundance (or density), and selecting the best fit detection function using AICc. We use the internal datasets thrasherDetectionData and thrasherSiteData (point transect surveys of brown thrashers). This tutorial is current as of version r packageVersion("Rdistance") of Rdistance.

1: Install and load Rdistance

If you haven't already done so, install the latest version of Rdistance. In the R console, issue install.packages("Rdistance"). After the package is installed, it can be loaded into the current session as follows:

require(Rdistance)

2: Read in the input data

For this tutorial, we use two datasets collected by J. Carlisle on brown thrashers in central Wyoming that are included with Rdistance.

The first dataset, thrasherDetectionData, is a detection data.frame with one row for each detected object. Columns in the data frame are:

siteID = Factor, the site (point) and transect surveyed. Levels are five character codes like 'TTXPP' where TT is transect number and PP is the point within the transect.
groupsize = Numeric, the number of individuals within the detected group.
dist = Numeric, the radial distance from the point to the detected group. Obtain access to the example dataset of thrasher detections and observed distances (thrasherDetectionData) using the following commands:

data("thrasherDetectionData")
head(thrasherDetectionData)

The second required dataset, thrasherSiteData, is a transect data.frame, with one row for each transect surveyed, and the following required columns:

siteID = Factor, the site (point) and transect surveyed.
... = Any additional transect-level covariate columns (these will not be used in this tutorial).

Load the example dataset of thrasher transects (thrasherSiteData) using the following commands:

data("thrasherSiteData")
head(thrasherSiteData)

3: Fit a detection function

Once the data are imported, the first step is to fit a detection function. Before we do so, explore the distribution of the distances:

hist(thrasherDetectionData$dist, n=40, col="grey", main="", xlab="distance (m)")
summary(thrasherDetectionData$dist)

Next, we fit a detection function using dfuncEstim to the radial distances collected from the point transects and plot it. We specify point transects using option PointSurvey = TRUE in the call to dfuncEstim and specify the the half-normal distance function using option likelihood = "halfnorm". In section 5, we demonstrate an automated process to fit multiple detection functions and compare them using AICc.

dfunc <- dfuncEstim(formula = dist ~ 1,
                    detectionData = thrasherDetectionData, 
                    pointSurvey = TRUE, 
                    likelihood = "halfnorm")

plot(dfunc)
dfunc

The effective detection radius (EDR) is the essential information from the detection function that will be used to estimate abundance in section 4. The EDR is calculated by integrating the detection function to compute area under the detection function. See the help documentation for EDR for details.

4: Estimate abundance given the detection function

Estimating abundance requires the additional information contained in the the thrasher site dataset, described in section 2, where each row represents one transect. Load the example dataset of surveyed thrasher transects from the package.

We estimate abundance (or density in this case) using abundEstim. If area = 1, density is given in the squared units of the distance measurements --- in this case, thrashers per square meter. If we set area = 10000, we convert to thrashers per hectare (1 ha == 10,000 m^2^). The equation used to calculate the abundance estimate is detailed in the help documentation for abundEstim.

Confidence intervals for abundance are calculated using a bias-corrected bootstrapping method (see abundEstim). Note that, as with all bootstrapping procedures, there may be slight differences in the confidence intervals between runs. Increasing the number of bootstrap iterations (R = 100 used here for brevity) may be necessary to stabilize CI estimates.

# Estimate Abundance - Density; fatalities per m2
fit <- abundEstim(dfunc          = dfunc,
                   detectionData = thrasherDetectionData,
                   siteData      = thrasherSiteData,
                   area           = 10000,       # density per hectare
                   R              = 100, 
                   ci             = 0.95)

fit

The abundance estimate can be extracted from the fit object.

fit$n.hat

The confidence interval (in this case 95%) can be extracted from the fit object.

fit$ci

5: Use AICc to select a detection function and estimate abundance

Fitting several detection functions, choosing the best fitting, and estimating abundance (sections 3 and 4) can be automated using the function autoDistSamp. The function attempts to fit multiple detection functions, uses AICc (by default, but see help documentation for autoDistSamp under criterion for other options) to find the 'best' detection function, then proceeds to estimate abundance using the best fit detection function (the distance function with lowest AICc). By default, autoDistSamp tries a large subset of Rdistance's built-in detection functions, but you can control exactly which detection functions are attempted (see help documentation for autoDistSamp). Specifying plot=TRUE produces a plot of each detection function as it is estimated. Specifying, plot.bs=TRUE plots the selected distance function each iteration of the bootstrap procedure. In this example, we fit the half-normal, hazard rate, exponential, and uniform likelihoods with no expansion terms, we do not plot all fitted functions (plot=FALSE), but we plot the best distance function fitted during each bootstrap iteration.

# Automated Fit - fit several models, choose the best model based on AIC
autoDS <- autoDistSamp(formula       = thrasherDetectionData$dist ~ 1,
                     detectionData = thrasherDetectionData,
                     siteData      = thrasherSiteData,
                     pointSurvey   = TRUE,
                     expansions    = c(0),
                     likelihoods   = c("halfnorm", "hazrate", "negexp", "uniform"),
                     plot          = FALSE,
                     area          = 10000,
                     R             = 100,
                     ci            = 0.95,
                     plot.bs       = FALSE)

autoDS

The detection function with the lowest AICc value (and thus selected as the 'best') is the hazard rate likelihood with 0 cosine expansion terms.

Conclusion

In sections 3 and 4, we fitted a half-normal detection function and used that function to estimate thrasher density. Our estimate was r round(fit$n.hat, 2) thrashers per ha (95% CI = r round(fit$ci[[1]], 2) to r round(fit$ci[[2]], 2)). In section 5, we used AICc to estimate a better fitting detection function and used it to estimate thrasher density. The thrasher density estimated by the better-fitting model was r round(autoDS$n.hat, 2) thrashers per ha (95% CI = r round(autoDS$ci[[1]], 2) to r round(autoDS$ci[[2]], 2)). (Note, CI estimates may vary slightly from these due to minor 'simulation slop' inherent in bootstrapping methods).