Beginner Line-Transect Analysis in Rdistance"

Introduction

This tutorial is meant to be a beginner's guide to analyzing line-transect data using Rdistance. It is assumed that you have some familiarity with Program R, but not necessarily with distance-sampling analysis. This beginning tutorial focuses on input data requirements, fitting a detection function, and estimating abundance (or density). Here, we make use of the example datasets already contained within Rdistance (i.e., line transect surveys of sparrows), so you can complete this tutorial without having any data of your own. 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 input data

Rdistance requires two input datasets. These can be prepared outside of R and read in as data.frames using, for example, read.csv. In the following sections, we make use of the sparrow example datasets already contained within Rdistance.

The first required dataset is a detection data.frame, with a row for each detection, and the following required columns, named as follows:

If the observers recorded sighting distance and sighting angle instead of perpendicular distance (as is often common in line transect surveys), you can use the perpDists function (detailed in Section 3) to calculate the perpendicular distances based on the sighting distances and sighting angles.

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

3: Fit a detection function

After prepping the input data, the first step is to explore your data and fit a detection function. First, load the example dataset of sparrow detections and their distances from the package using data(sparrowDetectionData) and transect characteristices using data(sparrowSiteData). Be sure that you have installed and loaded Rdistance prior to issuing the following commands:

data(sparrowDetectionData)
head(sparrowDetectionData)
data(sparrowSiteData)
head(sparrowSiteData)

Line-transect distance-sampling analysis is done on perpendicular distances (i.e., the distance from each detected group to the transect, not to the observer). We have provided the perpendicular distances (named dist) in the example data, but observers originally recorded sighting distances and sighting angles. Here we use the perpDists function to (re)calculate the perpendicular distances (dist) and remove the sightdist and sightangle columns. See the help documentation for perpDists for details.

sparrowDetectionData$dist<- perpDists(sightDist="sightdist",
                    sightAngle="sightangle", 
                    data = sparrowDetectionData)
sparrowDetectionData <- sparrowDetectionData[,
                    -which(names(sparrowDetectionData) 
                    %in% c("sightdist", "sightangle"))]

Explore the distribution of distances.

hist(sparrowDetectionData$dist, col="grey", main="", 
     xlab="distance (m)")
rug(sparrowDetectionData$dist,quiet = TRUE)
summary(sparrowDetectionData$dist)

Next, fit a detection function (plotted as a red line) using dfuncEstim. For now, we will proceed using the half-normal likelihood as the detection function, but in Section 5 of this tutorial, we demonstrate how to run an automated process that fits multiple detection functions and compares them using AICc. Note that distances greater than 150 m are quite sparse, so here we right-truncate the data, tossing out detections where dist > 150.

dfuncSparrow<- dfuncEstim(formula = dist~1, 
                          detectionData = sparrowDetectionData,
                          siteData = sparrowSiteData, 
                          likelihood = "halfnorm", w.hi = 150)

plot(dfuncSparrow)
dfuncSparrow

The effective strip width (ESW) is the key information from the detection function that will be used to next estimate abundance (or density). The ESW is calculated by integrating under the detection function. A survey with imperfect detection and ESW equal to X effectively covers the same area as a study with perfect detection out to a distance of X. See the help documentation for ESW for details.

4: Estimate abundance given the detection function

Estimating abundance requires the additional information contained in the second required dataset, described earlier, where each row represents one transect.

Next, estimate abundance (or density in this case) using abundEstim. If area=1, then density is given in the squared units of the distance measurements --- in this case, sparrows per square meter. Instead, we set area=10000 in order to convert to sparrows 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), and the detection function fit in each iteration of the bootstrap is plotted as a blue line (if plot.bs=TRUE). Note that, as with all bootstrapping procedures, there may be slight differences in the confidence intervals between runs due to simulation error.

Computing bootstrap confidence intervals for this vignette takes time (>5s), and the masters at CRAN have better things to do than wait for a package to build. In the code below, we supress bootstrap confidence intervals with ci = NULL, then replace confidence endpoints with values from a separate run with R = 500 iterations.

fit <- abundEstim(dfuncSparrow, 
                  detectionData=sparrowDetectionData,
                  siteData=sparrowSiteData,
                  area=10000, 
                  ci=NULL)
# To save vignette build time, we insert values from 
# a separate run with R=500
fit$ci <- c( lowCI=0.6473984, hiCI=1.0167007)
fit

Results of interest (such as the abundance estimate and confidence interval) can be extracted from the resulting object (here called fit).

fit$n.hat
fit$ci

5: Estimate abundance using covariates

For some data, the probability of detection may be dependent upon covariates due to site characteristics such as vegetation type or object characteristics such as animal size. Explore the distribution of distances for sparrows detected in different sagebrush shrub classes.

hist(sparrowDetectionData$dist[sparrowSiteData$shrubclass ==
     "High"], col="grey", main="High Sage",  
     xlab="distance (m)", breaks = 15, xlim = c(0,150))
rug(sparrowDetectionData$dist[sparrowSiteData$shrubclass == 
    "High"],quiet = TRUE)
hist(sparrowDetectionData$dist[sparrowSiteData$shrubclass == 
    "Low"], col="grey", main="Low Sage",  xlab="distance (m)", 
     breaks = 14, xlim = c(0,150))
rug(sparrowDetectionData$dist[sparrowSiteData$shrubclass == 
    "Low"],quiet = TRUE)

We can fit detection functions to these different site or object classes to better estimate abundance. We can specifiy different detection functions for these classes by using the formula argument in dfuncEstim. For the sparrow data, we fit detection functions to observations in high and low sagebrush shrub classes. Like the previous example, to save vignette build time, we supress bootstrap sampling then substitute values from a separate run with R=500.

dfuncShrubClass<- dfuncEstim(formula = dist~shrubclass, 
                          detectionData = sparrowDetectionData,
                          siteData = sparrowSiteData, 
                          likelihood = "halfnorm", w.hi = 150)
dfuncShrubClass <-
structure(list(parameters = c(`(Intercept)` = 4.01969799841488, 
shrubclassHigh = -0.181306851357376), varcovar = structure(c(0.0050270030909592, 
-0.00502975457487661, -0.00502975457487661, 0.00748138480810069
), .Dim = c(2L, 2L)), loglik = 1628.37871595743, convergence = 0L, 
    like.form = "halfnorm", w.lo = 0, w.hi = 150, dist = c(`1` = 16.8, 
    `2` = 12.2, `3` = 24.1, `4` = 3.5, `5` = 69.7, `6` = 10, 
    `7` = 14.9, `8` = 1.7, `9` = 4, `10` = 10, `11` = 9.8, `12` = 67.6, 
    `13` = 73.9, `14` = 74.7, `15` = 21.2, `16` = 12.9, `17` = 23.9, 
    `18` = 22, `19` = 25, `20` = 20.5, `21` = 65, `22` = 21.1, 
    `23` = 24.9, `24` = 31.9, `25` = 3.5, `26` = 70, `27` = 77.1, 
    `28` = 8.6, `29` = 12, `30` = 29.4, `31` = 46.7, `32` = 29.1, 
    `33` = 42.6, `34` = 29, `35` = 29.7, `36` = 22.5, `37` = 42.6, 
    `38` = 35.7, `39` = 27.7, `40` = 45.1, `41` = 5.9, `42` = 103, 
    `43` = 31.5, `44` = 13.5, `45` = 93, `46` = 7.3, `47` = 34.8, 
    `48` = 60.5, `49` = 16.1, `50` = 27.9, `51` = 44.4, `52` = 78, 
    `53` = 5.5, `54` = 75.8, `55` = 48, `56` = 19, `57` = 83.7, 
    `58` = 25.9, `59` = 15.5, `60` = 52.5, `61` = 26.4, `62` = 24, 
    `63` = 66.4, `64` = 8.9, `65` = 32.5, `66` = 26.4, `67` = 64.3, 
    `68` = 115.5, `69` = 5, `70` = 21.2, `71` = 61.9, `72` = 32.5, 
    `73` = 77.3, `74` = 34.5, `75` = 33, `76` = 9.6, `77` = 53.4, 
    `78` = 12.2, `79` = 27.6, `80` = 21.9, `81` = 34.7, `82` = 49.3, 
    `83` = 21.2, `84` = 14.7, `85` = 14, `86` = 16.1, `87` = 36.8, 
    `88` = 0, `89` = 38.4, `90` = 12.2, `91` = 126.6, `92` = 35, 
    `93` = 42.1, `94` = 36.3, `95` = 50, `96` = 46, `97` = 22.9, 
    `98` = 22.9, `99` = 12.7, `100` = 0, `101` = 7.5, `102` = 56.6, 
    `103` = 77.3, `104` = 16.5, `105` = 54.9, `106` = 0, `107` = 28, 
    `108` = 1.4, `109` = 26.8, `110` = 23.5, `111` = 26.8, `112` = 39, 
    `113` = 44, `114` = 57.3, `115` = 20.6, `116` = 32.5, `117` = 10.3, 
    `118` = 115.9, `119` = 50.8, `120` = 92, `121` = 0, `122` = 81, 
    `123` = 52.5, `124` = 19.5, `125` = 62.1, `126` = 100, `127` = 103.9, 
    `128` = 117.1, `129` = 6.4, `130` = 65, `131` = 33.5, `132` = 41, 
    `133` = 31, `135` = 73.9, `136` = 45, `137` = 16.1, `138` = 0, 
    `139` = 17, `140` = 11.6, `141` = 13.4, `142` = 2.6, `143` = 17.5, 
    `144` = 11, `145` = 1.2, `146` = 0, `147` = 0, `148` = 14, 
    `149` = 44.4, `150` = 17.6, `151` = 89.6, `152` = 81.1, `153` = 48.2, 
    `154` = 99.5, `155` = 36.6, `156` = 121.8, `157` = 46.5, 
    `158` = 6.7, `159` = 15.8, `160` = 5.9, `161` = 0, `162` = 67.7, 
    `163` = 7.8, `164` = 47.5, `165` = 52.5, `166` = 26, `167` = 14.6, 
    `168` = 59, `169` = 44.8, `170` = 68.1, `171` = 38, `172` = 2, 
    `173` = 4.6, `174` = 6, `175` = 7.8, `176` = 44.2, `177` = 50.8, 
    `178` = 19.9, `179` = 4.5, `180` = 11.4, `181` = 141.9, `182` = 26.4, 
    `183` = 15.8, `184` = 93, `185` = 40.6, `186` = 101.5, `187` = 24.7, 
    `188` = 23.9, `189` = 14.3, `190` = 34.5, `191` = 91.8, `192` = 109.6, 
    `193` = 53.6, `194` = 32.8, `195` = 6.1, `196` = 1.7, `197` = 40, 
    `198` = 20.5, `199` = 9.1, `200` = 35, `201` = 34.5, `202` = 24.6, 
    `203` = 93.7, `204` = 31.5, `205` = 45.8, `206` = 44.5, `207` = 0, 
    `208` = 7.5, `209` = 63.6, `210` = 10.1, `211` = 0, `212` = 0, 
    `213` = 46, `214` = 72, `215` = 8.5, `216` = 15.5, `217` = 2.2, 
    `218` = 94, `219` = 54, `220` = 14, `221` = 19.1, `222` = 31.2, 
    `223` = 78, `224` = 23, `225` = 110, `226` = 7.1, `227` = 26, 
    `228` = 19, `229` = 0, `230` = 5.2, `231` = 11.8, `232` = 63.9, 
    `233` = 54.4, `234` = 72.1, `235` = 46, `236` = 60.1, `237` = 60.8, 
    `238` = 19.6, `239` = 24.6, `240` = 26, `241` = 38.4, `242` = 144.8, 
    `243` = 25.9, `244` = 17.8, `245` = 32.2, `246` = 31, `247` = 59.1, 
    `248` = 74.6, `249` = 2.1, `250` = 12.5, `251` = 74.3, `252` = 7.2, 
    `254` = 63, `255` = 90.7, `256` = 5.5, `257` = 19.9, `258` = 22.9, 
    `259` = 28, `260` = 87.6, `261` = 114.6, `262` = 61, `263` = 55.4, 
    `264` = 57.5, `265` = 86, `266` = 71.9, `267` = 68.8, `268` = 63, 
    `269` = 39.6, `270` = 8.9, `271` = 71, `272` = 23.8, `273` = 41.5, 
    `274` = 87.6, `275` = 39, `276` = 31, `277` = 76, `278` = 45.2, 
    `279` = 12.2, `280` = 26.2, `281` = 8.5, `282` = 73.9, `283` = 87, 
    `284` = 77.8, `285` = 2.4, `286` = 30.5, `287` = 27.9, `288` = 134, 
    `289` = 30.5, `290` = 17.5, `291` = 4.5, `292` = 55.1, `293` = 19.4, 
    `294` = 8.5, `295` = 21.2, `296` = 79.5, `297` = 21, `298` = 102, 
    `299` = 32.5, `300` = 51.6, `301` = 118, `302` = 115.2, `303` = 44.4, 
    `304` = 117.5, `305` = 39, `306` = 29, `307` = 61, `308` = 42.4, 
    `310` = 21.7, `311` = 12, `312` = 33.4, `313` = 6.9, `314` = 0, 
    `315` = 46, `316` = 0, `317` = 66.5, `318` = 22.5, `319` = 15, 
    `320` = 60.8, `321` = 105, `322` = 0, `323` = 20, `324` = 14.7, 
    `325` = 23.6, `326` = 30.4, `327` = 52, `328` = 60, `329` = 46.8, 
    `330` = 78, `331` = 0, `332` = 24.3, `333` = 14.5, `334` = 9.8, 
    `335` = 46, `336` = 10.9, `337` = 55.3, `338` = 40.7, `339` = 58, 
    `340` = 57.2, `341` = 60.1, `342` = 45, `343` = 0.9, `344` = 28, 
    `345` = 14, `346` = 41.8, `347` = 0, `348` = 5.6, `349` = 14.2, 
    `350` = 35.5, `351` = 14, `352` = 51, `353` = 78.7, `354` = 27.3, 
    `355` = 1, `356` = 3.1), covars = structure(c(1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
    1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), .Dim = c(353L, 
    2L), .Dimnames = list(c("1", "2", "3", "4", "5", "6", "7", 
    "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", 
    "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", 
    "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", 
    "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", 
    "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", 
    "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", 
    "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", 
    "78", "79", "80", "81", "82", "83", "84", "85", "86", "87", 
    "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", 
    "98", "99", "100", "101", "102", "103", "104", "105", "106", 
    "107", "108", "109", "110", "111", "112", "113", "114", "115", 
    "116", "117", "118", "119", "120", "121", "122", "123", "124", 
    "125", "126", "127", "128", "129", "130", "131", "132", "133", 
    "135", "136", "137", "138", "139", "140", "141", "142", "143", 
    "144", "145", "146", "147", "148", "149", "150", "151", "152", 
    "153", "154", "155", "156", "157", "158", "159", "160", "161", 
    "162", "163", "164", "165", "166", "167", "168", "169", "170", 
    "171", "172", "173", "174", "175", "176", "177", "178", "179", 
    "180", "181", "182", "183", "184", "185", "186", "187", "188", 
    "189", "190", "191", "192", "193", "194", "195", "196", "197", 
    "198", "199", "200", "201", "202", "203", "204", "205", "206", 
    "207", "208", "209", "210", "211", "212", "213", "214", "215", 
    "216", "217", "218", "219", "220", "221", "222", "223", "224", 
    "225", "226", "227", "228", "229", "230", "231", "232", "233", 
    "234", "235", "236", "237", "238", "239", "240", "241", "242", 
    "243", "244", "245", "246", "247", "248", "249", "250", "251", 
    "252", "254", "255", "256", "257", "258", "259", "260", "261", 
    "262", "263", "264", "265", "266", "267", "268", "269", "270", 
    "271", "272", "273", "274", "275", "276", "277", "278", "279", 
    "280", "281", "282", "283", "284", "285", "286", "287", "288", 
    "289", "290", "291", "292", "293", "294", "295", "296", "297", 
    "298", "299", "300", "301", "302", "303", "304", "305", "306", 
    "307", "308", "310", "311", "312", "313", "314", "315", "316", 
    "317", "318", "319", "320", "321", "322", "323", "324", "325", 
    "326", "327", "328", "329", "330", "331", "332", "333", "334", 
    "335", "336", "337", "338", "339", "340", "341", "342", "343", 
    "344", "345", "346", "347", "348", "349", "350", "351", "352", 
    "353", "354", "355", "356"), c("(Intercept)", "shrubclassHigh"
    ))), model.frame = structure(list(dist = c(16.8, 12.2, 24.1, 
    3.5, 69.7, 10, 14.9, 1.7, 4, 10, 9.8, 67.6, 73.9, 74.7, 21.2, 
    12.9, 23.9, 22, 25, 20.5, 65, 21.1, 24.9, 31.9, 3.5, 70, 
    77.1, 8.6, 12, 29.4, 46.7, 29.1, 42.6, 29, 29.7, 22.5, 42.6, 
    35.7, 27.7, 45.1, 5.9, 103, 31.5, 13.5, 93, 7.3, 34.8, 60.5, 
    16.1, 27.9, 44.4, 78, 5.5, 75.8, 48, 19, 83.7, 25.9, 15.5, 
    52.5, 26.4, 24, 66.4, 8.9, 32.5, 26.4, 64.3, 115.5, 5, 21.2, 
    61.9, 32.5, 77.3, 34.5, 33, 9.6, 53.4, 12.2, 27.6, 21.9, 
    34.7, 49.3, 21.2, 14.7, 14, 16.1, 36.8, 0, 38.4, 12.2, 126.6, 
    35, 42.1, 36.3, 50, 46, 22.9, 22.9, 12.7, 0, 7.5, 56.6, 77.3, 
    16.5, 54.9, 0, 28, 1.4, 26.8, 23.5, 26.8, 39, 44, 57.3, 20.6, 
    32.5, 10.3, 115.9, 50.8, 92, 0, 81, 52.5, 19.5, 62.1, 100, 
    103.9, 117.1, 6.4, 65, 33.5, 41, 31, 201, 73.9, 45, 16.1, 
    0, 17, 11.6, 13.4, 2.6, 17.5, 11, 1.2, 0, 0, 14, 44.4, 17.6, 
    89.6, 81.1, 48.2, 99.5, 36.6, 121.8, 46.5, 6.7, 15.8, 5.9, 
    0, 67.7, 7.8, 47.5, 52.5, 26, 14.6, 59, 44.8, 68.1, 38, 2, 
    4.6, 6, 7.8, 44.2, 50.8, 19.9, 4.5, 11.4, 141.9, 26.4, 15.8, 
    93, 40.6, 101.5, 24.7, 23.9, 14.3, 34.5, 91.8, 109.6, 53.6, 
    32.8, 6.1, 1.7, 40, 20.5, 9.1, 35, 34.5, 24.6, 93.7, 31.5, 
    45.8, 44.5, 0, 7.5, 63.6, 10.1, 0, 0, 46, 72, 8.5, 15.5, 
    2.2, 94, 54, 14, 19.1, 31.2, 78, 23, 110, 7.1, 26, 19, 0, 
    5.2, 11.8, 63.9, 54.4, 72.1, 46, 60.1, 60.8, 19.6, 24.6, 
    26, 38.4, 144.8, 25.9, 17.8, 32.2, 31, 59.1, 74.6, 2.1, 12.5, 
    74.3, 7.2, 196.6, 63, 90.7, 5.5, 19.9, 22.9, 28, 87.6, 114.6, 
    61, 55.4, 57.5, 86, 71.9, 68.8, 63, 39.6, 8.9, 71, 23.8, 
    41.5, 87.6, 39, 31, 76, 45.2, 12.2, 26.2, 8.5, 73.9, 87, 
    77.8, 2.4, 30.5, 27.9, 134, 30.5, 17.5, 4.5, 55.1, 19.4, 
    8.5, 21.2, 79.5, 21, 102, 32.5, 51.6, 118, 115.2, 44.4, 117.5, 
    39, 29, 61, 42.4, 207, 21.7, 12, 33.4, 6.9, 0, 46, 0, 66.5, 
    22.5, 15, 60.8, 105, 0, 20, 14.7, 23.6, 30.4, 52, 60, 46.8, 
    78, 0, 24.3, 14.5, 9.8, 46, 10.9, 55.3, 40.7, 58, 57.2, 60.1, 
    45, 0.9, 28, 14, 41.8, 0, 5.6, 14.2, 35.5, 14, 51, 78.7, 
    27.3, 1, 3.1), shrubclass = structure(c(2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Low", "High"), class = "factor")), terms = dist ~ 
        shrubclass, row.names = c(NA, 356L), class = "data.frame"), 
    expansions = 0, series = "cosine", call = quote(dfuncEstim(formula = dist ~ 
        shrubclass, detectionData = sparrowDetectionData, siteData = sparrowSiteData, 
        likelihood = "halfnorm", w.hi = 150)), call.x.scl = 0, 
    call.g.x.scl = 1, call.observer = "both", fit = list(par = c(`(Intercept)` = 4.01969799841488, 
    shrubclassHigh = -0.181306851357376), objective = -5686.93412448465, 
        convergence = 0L, iterations = 10L, counts = c(`function` = 14L, 
        gradient = 27L), message = "relative convergence (4)", 
        hessian = structure(c(607.723213611669, 408.573906078495, 
        408.573906078495, 408.350399242343), .Dim = c(2L, 2L))), 
    factor.names = "shrubclass", pointSurvey = FALSE, formula = dist ~ 
        shrubclass, x.scl = 0, g.x.scl = 1), class = "dfunc")
plot(dfuncShrubClass, newdata = data.frame(shrubclass=
                levels(sparrowSiteData$shrubclass)), nbins = 15)
covarFit <- abundEstim(dfuncShrubClass, 
                       detectionData=sparrowDetectionData,
                       siteData=sparrowSiteData,
                       area=10000, 
                       ci=NULL)
# To save vignette build time, insert values from 
# a separate run with R=500
covarFit$ci <- c("2.939477%" = 0.6552051, "97.88301%"= 1.0411813)
covarFit

In this example including the covariate of shrub class changes our estimate from r round(fit$n.hat, 2) sparrows per ha (95% CI = r round(fit$ci[[1]], 2) - r round(fit$ci[[2]], 2)) to r round(covarFit$n.hat, 2) sparrows per ha (95% CI = r round(covarFit$ci[[1]], 2) - r round(covarFit$ci[[2]], 2)). In this case, the density estimate did not change much when we included shrubClass, but this will not always be the case. Covariates can have a large influence on abundance estimates.

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

Alternatively, steps 3 (fitting a detection function) and 4 (estimating abundance) can be automated using the function autoDistSamp. This function attempts to fit multiple detection functions, uses AICc (by default, but see help documentation for AIC.dfunc for other options) to find the 'best' detection function, then proceeds to estimate abundance using that detection function. 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 would return a plot of each detection function. In this example, we fit the three detection function forms combined with 0 or 1 expansion terms, and we don't plot them (plot=FALSE). If autoDistSamp is called from another function, it can be convenient to suppress all output by setting plot.bs=FALSE, plot=FALSE, and showProgress=FALSE.

auto <- autoDistSamp(formula        = dist~1,
                     detectionData  = sparrowDetectionData,
                     siteData       = sparrowSiteData,
                     area           = 10000, 
                     likelihoods    = c("halfnorm","hazrate","uniform"),
                     expansions     = 0:1,
                     series         = "cosine",
                     ci             = 0.95,
                     R              = 500,
                     plot           = FALSE,
                     plot.bs        = FALSE,
                     w.hi           = 150)  
auto<-list(n.hat=1.003543 ,
  ci=list(0.7823177, 1.366216))
## Likelihood   Series  Expans  Converged?  Scale?  AICc
## halfnorm cosine  0   Yes     Ok  3263.443
## halfnorm cosine  1   Yes     Ok  3261.5853
## hazrate      cosine  0   Yes     Ok  3267.6249
## hazrate      cosine  1   Yes     Ok  3263.3098
## uniform      cosine  0   Yes     Ok  3260.7321
## uniform      cosine  1   Yes     Ok  3262.4814
## Computing bootstrap confidence interval on N...
## |===================================================================| 100%
## 101 of 500 iterations did not converge.
## 
## ------------ Abundance Estimate Based on Top-Ranked Detection Function ------------
## Call: dfuncEstim(formula = formula, detectionData = detectionData,
##    siteData = siteData, likelihood = fit.table$like[1], pointSurvey =
##    pointSurvey, w.lo = w.lo, w.hi = w.hi, expansions =
##    fit.table$expansions[1], series = fit.table$series[1])
## Coefficients:
##            Estimate    SE            z         p(>|z|)     
## Threshold  27.7397495  19.537302641  1.419835  1.556557e-01
## Knee        0.0340005   0.005304744  6.409451  1.460443e-10
## 
## Convergence: Success
## Function: UNIFORM  
## Strip: 0 to 150 
## Effective strip width (ESW): 51.34584 
## Probability of detection: 0.3423056 
## Scaling: g(0) = 1
## Log likelihood: 1628.349 
## AICc: 3260.732
## 
## Abundance estimate:  1.003543 ;  95% CI=( 0.7823177 to 1.366216 )
## CI based on 399 of 500 successful bootstrap iterations

You can see that the detection function with the lowest AICc value (and thus selected as the 'best') is the uniform likelihood. It is worth noting that approximately 20% of the bootstrap iterations did not converge, indicating that the uniform detection function is unstable. It might be best in this situation to exclude the uniform from the list of potential distance function forms.

We can also estimate abundances with covariates using the function autoDistSamp. However, automatic selection of covariates does not (yet) occur and none of the likelihoods accept expansions when covariates are present. In this example, we include a covariate in the distance function (shrubclass) and fit three different forms (likelihoods = c("halfnorm", "hazrate", "negexp")) without expansions (expansions = 0).

autoCov <- autoDistSamp(formula     = dist~shrubclass,
                     detectionData  = sparrowDetectionData,
                     siteData       = sparrowSiteData,
                     likelihoods    = c("halfnorm", "hazrate", "negexp"),
                     expansions     = 0,
                     area           = 10000, 
                     ci             = 0.95,
                     R              = 500,
                     plot           = FALSE,
                     plot.bs        = FALSE,
                     w.hi           = 150)  
## Likelihood Series  Expans  Converged?  Scale?  AICc
## halfnorm   cosine  0       Yes         Ok      3260.7917
## hazrate    cosine  0       Yes         Ok      3264.531
## negexp     cosine  0       Yes         Ok      3260.8539
## Computing bootstrap confidence interval on N...
## |===================================================================| 100%
## 
## ------------ Abundance Estimate Based on Top-Ranked Detection Function ------------
## Call: dfuncEstim(formula = formula, detectionData = detectionData,
##    siteData = siteData, likelihood = fit.table$like[1], pointSurvey =
##    pointSurvey, w.lo = w.lo, w.hi = w.hi, expansions =
##    fit.table$expansions[1], series = fit.table$series[1])
## Coefficients:
##                 Estimate    SE          z          p(>|z|)   
## (Intercept)      4.0196980  0.07090136  56.694228  0.00000000
## shrubclassHigh  -0.1813069  0.08649500  -2.096154  0.03606852
## 
## Convergence: Success
## Function: HALFNORM  
## Strip: 0 to 150 
## Average effective strip width (ESW): 62.28349 
## Average probability of detection: 0.4152232 
## Scaling: g(0) = 1
## Log likelihood: 1628.379 
## AICc: 3260.792
## 
## Abundance estimate:  0.8350688 ;  95% CI=( 0.6484067 to 1.036907 )
autoCov<-list(n.hat=0.8350688 ,
  ci=list(0.6484067,1.036907))

Conclusion

Note that the detection function that you select can have a large influence on the resulting abundance estimate. In sections 3 and 4, we fit a half-normal detection function and used that function to estimate sparrow density. Our estimate was r round(fit$n.hat, 2) sparrows per ha (95% CI = r round(fit$ci[[1]], 2) - r round(fit$ci[[2]], 2)). In section 5, we included shrubClass in the model and found that the average density estimate changed little. This will not always be the case as inclusion of covariates can have a substantial effect on estimates. In section 6, we used AICc to determine the best-fitting detection function amoung several we chose to fit and used that function to estimate sparrow density again. Initially, AICc determined that the uniform likelihood fit best and estimated density at r round(auto$n.hat, 2) sparrows per ha (95% CI = r round(auto$ci[[1]], 2) - r round(auto$ci[[2]], 2)). But, the high number of non-convergent bootstap iterations lead to some suspicion of the initial model. We added shrubClass into the model, dropped the uniform likelihood from consideration and AICc determined that the half normal detection function fit better than the hazard rate and negative exponential, but not by much. This latter model, also fitted in section 5, estimated density at r round(autoCov$n.hat, 2) sparrows per ha (95% CI = r round(autoCov$ci[[1]], 2) - r round(autoCov$ci[[2]], 2)). (Note, recall that your estimates may vary slightly from these due to variation inherent in bootstrapping methods). We conclude that the half normal distance function with shrubClass is a reasonable model for these data and that different functional forms for the distance function can have a substantial impact on abundance estimates. The autoDistSamp function can help select a detection function based on AICc, but ultimately the analyst will need to make the final distance function choice based on plots and other diagnostics.



Try the Rdistance package in your browser

Any scripts or data that you put into this service are public.

Rdistance documentation built on May 2, 2019, 3:49 a.m.