Package 'Rdistance'

Description

Rdistance contains functions and associated routines to analyze distance-sampling data collected on point or line transects. Some of Rdistance's features include:

• Accommodation of both point and line transect analyses in one routine (dfuncEstim()).

• Regression-like formula for inclusion of distance function covariates (dfuncEstim()).

• Automatic bootstrap confidence intervals (abundEstim()).

• Parallel processing of bootstrap iterations (parallel argument of abundEstim()).

• Availability of both study-area and site-level abundance estimates (help("predict.dfunc")).

• Rigorous physical measurement requirements and automated conversion when necessary (%#%, setUnits()).

• Classic parametric distance functions (halfnorm.like(), hazrate.like(), negexp.like()), and expansion functions (cosine.expansion(), hermite.expansion(), simple.expansion()).

• Mixture distance functions for non-standard shapes and thresholds (oneStep.like(), triangle.like(), and huber.like()).

• Automated distance function fitting and selection autoDistSamp().

• print, plot, predict, coef, and summary methods for distance function objects and abundance classes.

Background:

Distance-sampling is a popular method for abundance estimation in ecology. Line transect surveys are conducted by traversing randomly placed transects in a study area with the objective of sighting animals and estimating density or abundance. Data collected during line transect surveys consists of target sightings, either of individuals or groups, and off-transect distances to the original location of the target. When targets are sighted in groups, data include the number of individuals in the group.

Point transect surveys are similar except that observers stop one or more times along the transect to observe targets. Point transects are popular avian survey methods where detections are often auditory cues. Point transects are also for studies using automated auditory detectors or trail cameras. Point transect data consists of radial distances from the observer to the target.

The defining feature of distance sampling is the tendency for probability of detection to decline as off-transect or radial distances increase. Targets far from the observer are generally harder to detect than those closer. In most line transect studies, targets on the transect (off-transect distance = 0) are assumed to be sighted with 100% probability. This assumption allows researchers to estimate the proportion of missed targets and in turn adjust the number of sighted targets for missed detections. Some studies utilize two observers searching the same areas and are able to estimate the proportion of individuals missed on the transect line and thereby eliminate the assumption that all individuals on the line have been observed.

Purpose:

The author's aims are to provide an easy-to-use, rigorous, and flexible analysis option in R for distance-sampling data. The authors believe that beginning users need easy-to-use and easy-to-understand software, while advanced users require greater flexibility and customization, and their aim is to meet the demands of both groups.

Data sets:

Rdistance contains the following example data sets:

• Line-transect sampling of Brewers sparrows in central Wyoming (sparrowDf()).

• Point-transect sampling of Sage Thrashers in central Wyoming (thrasherDf()).

References

Buckland, S.T., Anderson, D.R., Burnham, K.P. and Laake, J.L. 1993. Distance Sampling: Estimating Abundance of Biological Populations. Chapman and Hall, London.

Author(s)

Main author and maintainer: Trent McDonald [email protected]

Coauthors: Ryan Nielson, Jason Carlisle, and Aidan McDonald

Contributors: Ben Augustine, James Griswald, Joel Reynolds, Pham Quang, Earl Becker, Aaron Christ, Brook Russelland, Patrick McKann, Lacey Jeroue, Abigail Hoffman, Michael Kleinsasser, and Ried Olson

See Also

Useful links:

• https://mcdonalddatasciences.com/Rdistance.html

• Report bugs at https://github.com/tmcd82070/Rdistance/issues

Description

Helper functions for assigning physical measurement units in Rdistance. All are convenience wrappers for ⁠units::⁠units::set_units().

Usage

x %#% u

x %acre% .

x %cm% .

x %ft% .

x %ha% .

x %inches% .

x %km% .

x %km^2% .

x %m% .

x %m^2% .

x %mi% .

x %mi^2% .

x %yd% .

dropUnits(x)

setUnits(x, u)

Arguments

Details

The fixed unit assignment operators are designed to behave like unary operators (i.e., 1 argument); but, R does not allow user defined unary operators. Technically, these fixed unit assignment operators are instances of R's user-defined infix operator, and as such they require two arguments. Their syntax must be ⁠x %<units>% <something>⁠; but, the second argument is ignored and '.' is suggested. See Examples.

Value

For ⁠%#%⁠ and setUnits, argument x with units u attached.

For all the fixed unit assignment operators (i.e., %units%), argument x with the respective units assigned.

For dropUnits, argument x with no units. If input x has no units, x is returned unchanged.

Examples

2 %#% "m" 
setUnits(2,"km")
x <- 2 %#% "km^2"
10 %#% units(x)
2 %#% "km^2" %#% "acres" # Convert km^2 to acres
x %#% "acres"            # Same
x %#% NULL    # Drop units
dropUnits(x)  # Same

# %#%'s precedence is below "^" but above "+" and "*"
# The following fails: 
# 2 %#% "m" / (2 %#% "ha") %#% "in/acre"
# The following succeeds:
(2 %#% "m" / (2 %#% "ha")) %#% "in/acre"
1 %m%. ^ 2      # [m]
(1 %m%.) ^ 2    # [m^2]


# For fixed unit assignment, 2nd argument does not matter
# All of the following are equivalent
2 %m%.
2 %m% x
2 %m% 3
2 %m% NULL
2 %m% NA

# Conversion:
x <- 10 %#% "ft"
x %m%.

Description

Estimate abundance (or density) from an estimated detection function and supplemental information on observed group sizes, transect lengths, area surveyed, etc. Computes confidence intervals on abundance (or density) using a the bias corrected bootstrap method.

Usage

abundEstim(
  object,
  area = NULL,
  propUnitSurveyed = 1,
  ci = 0.95,
  R = 500,
  plot.bs = FALSE,
  showProgress = TRUE,
  parallel = TRUE
)

Arguments

Details

The abundance estimate for line-transect surveys (if no covariates are included in the detection function and both sides of the transect are observed) is

$N =\frac{n(A)}{2(ESW)(L)}$

where n is total number of sighted individuals (i.e., sum(groupSizes(dfunc))), L is the total length of surveyed transect (i.e., sum(effort(dfunc))), and ESW is effective strip width computed from the estimated distance function (i.e., ESW(dfunc)). If only one side of transects were observed, the "2" in the denominator is not present (or, replaced with a "1").

The abundance estimate for point transect surveys (if no covariates are included) is

$N =\frac{n(A)}{\pi(ESR^2)(P)}$

where n is total number of sighted individuals (i.e., sum(groupSizes(dfunc))), P is the total number of surveyed points (i.e., sum(effort(dfunc))), and ESR is effective search radius computed from the estimated distance function (i.e., ESR(dfunc)).

This routine, abundEstim, estimates abundance on the entire study area. Site-specific density estimates are computed by predict(x, type = "density"), which returns a tibble containing density and abundance on the area surveyed by every transect.

Value

An Rdistance 'abundance estimate' object, which is a list of class c("abund", "dfunc"), containing all the components of a "dfunc" object (see dfuncEstim()), plus the following:

Bootstrap Confidence Intervals:

Rdistance's nested data frames (produced by RdistDf()) contain all information required to estimate bootstrap CIs. To compute bootstrap CIs, Rdistance resamples, with replacement, the rows of the ⁠$data⁠ component contained in Rdistance fitted models. Rdistance assumes each row of ⁠$data⁠ contains information on one transect. The ⁠$data⁠ component also contains information on which observations inform the detection function, which observations should be counted as detected targets, and which transects count toward transect length. After resampling rows of ⁠$data⁠, Rdistance refits the distance function using non-missing distances, recomputes the detected number of targets using non-missing group sizes on transects with non-missing length, and re-computes total transect length from transects with non-missing lengths. By default, R = 500 bootstrap iterations are performed, after which bias corrected confidence intervals are computed (Manly, 1997, section 3.4).

The distance function is not re-selected during bootstrap resampling. The model of the input object is re-fitted every iteration.

During bootstrap iterations, the distance function can fail. An iteration can fail for a two reasons: (1) no detections on the iteration, and (2) a bad configuration of distances that push the distance function's parameters to their limits. When an iteration fails, Rdistance skips the iteration and effectively ignores the failed iterations. If the proportion of failed iterations is small (less than 20% by default), the resulting abundance confidence interval is probably valid and no warning is issued. If the proportion of non-convergent iterations is not small (exceeds 20% by default), a warning is issued. The warning can be modified by re-setting option "Rdistance_maxBSFailPropForWarning" to the acceptable proportion of failures.. Setting options(Rdistance_masBSFailPropForWarning = 1.0) will turn suppress the warning. Setting options(Rdistance_masBSFailPropForWarning = 0.0) will warn if any iteration failed. Results (density and effective sampling distance) from all successful iterations are contained in the non-NA rows of data frame 'B' in the output object.

Missing Transect Lengths

Transect lengths can be missing in the RdistDf object. Missing length transects are equivalent to 0 [m] transects and do not count toward total surveyed units, nor do group sizes on these transects count toward total detected individuals. Use NA-length transects to include their associated distances during distance function estimation, but not when estimating abundance. For example, this allows estimation of abundance on one study area using off-transect distances from another. This allows sightability to be estimated using two or more similar targets (e.g., two similar species), but abundance to be estimated separate for each target type. Include NA-length transects by including the "extra" distance observations in the detection data frame, with valid site IDs, but set the length of those site IDs to NA in the site data frame.

Point Transect Lengths

Point transects do not have a physical measurement for length. The "length" of point transects is the number of points on the transect. Point transects can contain only one point. Rdistance treats transects of points as independent and bootstrap resamples them to estimate variance. The number of points on each point transect must exist in the RdistDf and cannot have physical measurement units (it is a count, not a distance).

References

Manly, B.F.J. (1997) Randomization, bootstrap, and Monte-Carlo methods in biology, London: Chapman and Hall.

Buckland, S.T., D.R. Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. (2001) Introduction to distance sampling: estimating abundance of biological populations. Oxford University Press, Oxford, UK.

See Also

dfuncEstim(), autoDistSamp(), predict.dfunc() with 'type = "density"'.

Examples

# Load example sparrow data (line transect survey type)
# sparrowDf <- RdistDf(sparrowSiteData, sparrowDetectionData)
data(sparrowDf)

# Fit half-normal detection function
dfunc <- sparrowDf |> 
  dfuncEstim(formula=dist ~ groupsize(groupsize)
           , likelihood="halfnorm"
           , w.hi=150 %m%.
  )

# Estimate abundance - Convenient for programming 
abundDf <- estimateN(dfunc
                   , area = 4105 %km^2%.
           )

# Same - Nicer output 
# Set ci=0.95 (or another value) to estimate bootstrap CI's 
fit <- abundEstim(dfunc
                , area = 4105 %km^2%.
                , ci = NULL
                )
summary(fit)

Description

Computes AICc, AIC, or BIC for estimated distance functions.

Usage

## S3 method for class 'dfunc'
AIC(object, ..., criterion = "AICc")

Arguments

Details

Regular Akaike's information criterion (https://en.wikipedia.org/wiki/Akaike_information_criterion) ($AIC$) is

$AIC = LL + 2p,$

where $LL$ is the maximized value of the log likelihood (the minimized value of the negative log likelihood) and $p$ is the number of coefficients estimated in the detection function. For dfunc objects, $AIC$ = obj$loglik + 2*length(coef(obj)).

A correction for small sample size, $AIC_c$, is

$AIC_c = LL + 2p + \frac{2p(p+1)}{n-p-1},$

where $n$ is sample size or number of detected groups for distance analyses. By default, this function computes $AIC_c$. $AIC_c$ converges quickly to $AIC$ as $n$ increases.

The Bayesian Information Criterion (BIC) is

$BIC = LL + log(n)p,$

.

Value

A scalar, the requested fit statistic for object.

References

Burnham, K. P., and D. R. Anderson, 2002. Model selection and multi-model inference: A practical information-theoretic approach, Second ed. Springer-Verlag. ISBN 0-387-95364-7.

McQuarrie, A. D. R., and Tsai, C.-L., 1998. Regression and time series model selection. World Scientific. ISBN 981023242X

See Also

coef(), dfuncEstim()

Examples

data(sparrowDf)
dfunc <- sparrowDf |> dfuncEstim(dist~1)
  
# Fit statistics
AIC(dfunc)  # AICc
AIC(dfunc, criterion="AIC")  # AIC
AIC(dfunc, criterion="BIC")  # BIC

Description

Perform automated likelihood, expansion, and series selection for a classic distance sampling analysis. Estimate abundance using the best fitting likelihood, expansion, and series.

Usage

autoDistSamp(
  data,
  formula,
  likelihoods = c("halfnorm", "hazrate", "negexp"),
  w.lo = setUnits(0, "m"),
  w.hi = NULL,
  expansions = 0:3,
  series = c("cosine"),
  x.scl = w.lo,
  g.x.scl = 1,
  warn = TRUE,
  outputUnits = NULL,
  area = NULL,
  propUnitSurveyed = 1,
  ci = 0.95,
  R = 500,
  plot.bs = FALSE,
  showProgress = TRUE,
  parallel = TRUE,
  plot = TRUE,
  criterion = "AICc"
)

Arguments

Details

During distance function selection, all combinations of likelihoods, series, and number of expansions is fitted. For example, if likelihoods has 3 elements, series has 2 elements, and expansions has 4 elements, this routine fits a total of 3 (likelihoods) * 2 (series) * 4 (expansions) = 24 models. Default parameters fit 9 detection functions, i.e., all combinations of "halfnorm", "hazrate", and "negexp" likelihoods and 0 through 3 expansions. Other combinations are specified through values of likelihoods, series, and expansions.

Suppress all intermediate output using plot.bs=FALSE, showProgress=FALSE, and plot=FALSE.

The returned abundance estimate object contains an additional component, the fitting table (a list of models fitted and criterion values) in component ⁠$fitTable⁠.

Value

An Rdistance 'abundance estimate' object, which is a list of class c("abund", "dfunc"), containing all the components of a "dfunc" object (see dfuncEstim()), plus the following:

See Also

dfuncEstim(), abundEstim().

Examples

# Load example sparrow data (line transect survey type)
data(sparrowDf)

autoDistSamp(data = sparrowDf
           , formula = dist ~ groupsize(groupsize)
           , likelihoods = c("halfnorm","negexp")
           , expansions = 0
           , plot = FALSE
           , ci = NULL
           , area = 1 %ha%.
)

## Not run: 
autoDistSamp(data = sparrowDf
    , formula = dist ~ 1 + groupsize(groupsize)
    , ci = 0.95
    , area = 1 %ha%.
)     

## End(Not run)

Description

Calculate bias-corrected confidence intervals for bootstrap data using methods in Manly textbook.

Usage

bcCI(x.bs, x, ci = 0.95)

Arguments

Value

A named vector containing the lower and upper endpoints of the bias-corrected bootstrap confidence interval.

Description

Performs bootstrap resampling iterations, either in parallel across CPU cores or in serial on a single core.

Usage

bootstrap(
  object,
  area,
  propUnitSurveyed,
  R,
  plot.bs,
  plotCovValues,
  showProgress,
  parallel,
  cores
)

Arguments

Value

A data frame containing density, abundance, and other relevant statistics for every bootstrap iteration. Number of rows is R. If the model from one iteration failed for any reason (e.g., non-convergence), the entire row except the ID column is missing.

See Also

abundEstim(); oneBsIter()

Description

Computes spline, specifically b-spline, expansion terms that modify the shape of distance likelihood functions.

Usage

bspline.expansion(x, expansions)

Arguments

Value

A 3D array of size nrow(x) X ncol(x) X expansions. The 'pages' (3rd dimension) of this array are the cosine expansions of x. i.e., page 1 is the first expansion term of x, page 2 is the second expansion term of x, etc.

See Also

dfuncEstim(), cosine.expansion(), hermite.expansion(), simple.expansion().

Examples

x <- matrix(seq(0, 1, length = 200), ncol = 1)
spl.expn <- bspline.expansion(x, 5)
plot(range(x), range(spl.expn), type="n")
matlines(x, spl.expn[,1,1:5], col=rainbow(5), lty = 1)

Description

Check that number of integration intervals is odd and sufficiently large. Plus, compute and store the Simpson's Composite rule coefficients.

Usage

checkNEvalPts(nEvalPts)

Arguments

Value

The first element of nEvalPts is returned if it is acceptable. If nEvalPts is not acceptable, an error is thrown. Simpson's composite coefficients are store in options()

Description

Check for the presence of physical measurement units on key columns of an RdistDf data frame.

Usage

checkUnits(ml)

Arguments

Value

The input ml list, with units of various quantities converted to common units. If a check fails, for example, a quantity does not have units, an error is thrown.

Description

Extract distance model coefficients from an estimated detection function object.

Usage

## S3 method for class 'dfunc'
coef(object, ...)

Arguments

Value

The estimated coefficient vector for the detection function. Length and interpretation of values vary depending on the form of the detection function and expansion terms.

See Also

AIC(), dfuncEstim()

Examples

data(sparrowDfuncObserver) # pre-estimated dfunc

# Same as sparrowDfuncObserver$par 
coef(sparrowDfuncObserver) 

## Not run: 
data(sparrowDf)
dfunc <- sparrowDf |> dfuncEstim(dist~bare + observer,
                      w.hi = 150 %m%.)
coef(dfunc)

## End(Not run)

Description

Add ANSI color to a string using the crayon package, if the R environment accepts color. This function is needed because some output cannot be colorized. Color determination is made by crayon::has_color().

Rdistance results often include units, e.g., "25 [m]". Only colorize the numbers, not the units. Everything between the first set of square brackets (⁠[...]⁠) is NOT colorized. Subsequent sets of brackets (e.g., "25 [m] + 30 [ft]") ARE colorized (i.e., "[ft]" is color).

Usage

colorize(STR, col = NULL, bg = NULL)

Arguments

Value

If color is not allowed in the terminal, the input string is returned unperturbed. If color is allowed, the input string is returned with color and background ANSI code surrounding the initial part of the string from character 1 to the character before the [ in the first pair of ⁠[]⁠.

See Also

crayon::style

Description

Computes the cosine expansion terms that modify the shape of distance likelihood functions.

Usage

cosine.expansion(x, expansions)

Arguments

Details

The cosine expansion used here is:

• First term:

  $h_1(x)=\cos(2\pi x),$

• Second term:

  $h_2(x)=\cos(4\pi x),$

• Third term:

  $h_3(x)=\cos(6\pi x),$

• Fourth term:

  $h_4(x)=\cos(8\pi x),$

• Fifth term:

  $h_5(x)=\cos(10\pi x),$

The maximum number of expansion terms is 5.

The cosine expansion frequency in Rdistance is 2*pi. Each term is two pi more than the previous. The sine expansion frequency in Rdistance is pi. Consequently, the sine and cosine expansions fit different models.

Value

A 3D array of size nrow(x) X ncol(x) X expansions. The 'pages' (3rd dimension) of this array are the cosine expansions of x. i.e., page 1 is the first expansion term of x, page 2 is the second expansion term of x, etc.

See Also

dfuncEstim(), hermite.expansion(), simple.expansion()

Examples

x <- matrix(seq(0, 1, length = 200), ncol = 1)
cos.expn <- cosine.expansion(x, 5)
plot(range(x), range(cos.expn), type="n")
matlines(x, cos.expn[,1,1:5], col=rainbow(5), lty = 1)

Description

Fits a detection function to off-transect distances collected by multiple observers.

Usage

dE.multi(
  data,
  formula,
  likelihood = "halfnorm",
  w.lo = setUnits(0, "m"),
  w.hi = NULL,
  expansions = 0,
  series = "cosine",
  x.scl = setUnits(0, "m"),
  g.x.scl = 1,
  warn = TRUE,
  outputUnits = NULL
)

Arguments

Value

An object of class 'dfunc' with the following components:

• par: The vector of estimated parameter values. Length of this vector is the sum of the following:

  1. The number of columns of the design matrix. This equals the number of covariates in the distance function plus one for the intercept, assuming an intercept is included.

  2. The number of constant parameters in the distance function. Constant parameters are those not related to covariates. For example, the exponent 'k' parameter for hazard rate likelihood, or the mixing fraction 'p' for the oneStep likelihood. This can be zero.

  3. The number of expansion functions called for. This equals the input expansions.

• loglik: The maximized value of the log likelihood.

• convergence: The convergence code. This code is returned by the optimizing routine (e.g., optim or nlminb). Values other than 0 indicate suspect convergence.

• message: If maximization did not converge (convergence != 0), this is the reason given by the optimizing routine.

• varcovar: The variance-covariance matrix for coefficients of the distance function, either estimated by the inverse of the fit's Hessian or by bootstrapping. If the likelihood is smooth (i.e., those listed by ⁠Rdistance:::differentiableLikelihoods())⁠, Rdistance initially estimates the variance-covariance matrix using the second derivative of the log likelihood surface at the final estimates, where second derivatives are estimated by numeric differentiation (in routine secondDeriv(). The variance-covariance matrix is re-set to NULL if the Hessian is not positive-definite. If bootstrap resampling has been performed (using abundEstim()), the variance-covariance matrix is re-estimated using the bootstrap values of parameters and automatically reset. Error estimates derived from bootstrapping are generally preferable to the asymptotic estimates, hence the automatic re-set.

• limits: A list containing the lower and upper limits of parameters.

• evaluations: The number of likelihood evaluations performed by the optimizer.

• mf: An R 'model frame' containing the detections (within the strip or circle) used in the fit, covariates specified in the formula, and groupsizes. Column 'dist' contains the observed distances. The intercept, if included in the model, is not included as a column in this model frame. (Test whether an intercept is included using attr(terms(return$mf), "intercept")). Column offset(...) contains group sizes associated with the values of dist. Name of the group size column is "offset(...)", not "groupsize(...)", so that group sizes can be treated offsets in other R routines. The mf component is a proper model.frame and contains both terms and contrasts attributes. This model frame contains only non-missing distances between w.lo and w.hi.

• data: The original nested data frame subset to information required to complete distance estimation. This data frame contains information on replication (i.e., rows are sites and are re-sampled during bootstrapping), missing distances, missing transect lengths, and distances outside the observation strip (below w.lo or above w.hi).

• formula: The distance function's formula.

• dataName: Name of the original nested data frame.

• likelihood: The name of the likelihood fitted to observation distances.

• w.lo: Left-truncation value used during the fit.

• w.hi: Right-truncation value used during the fit.

• expansions: The number of expansion terms.

• series: The type of expansion used during estimation. This is only relevant if expansions > 0.

• x.scl: The distance at which the function has been scaled to some value. This is the x at which g(x) = g.x.scl.

• g.x.scl: The height of the distance function at distance x.scl.

• outputUnits: A list of type symbolic_units containing the physical measurement units used during estimation.

• asymptoticSE: A logical scalar indication whether the variance-covariance matrix in component varcovar is asymptotic (TRUE; estimated from the Hessian) or bootstrap (FALSE; estimated by bootstrap resampling).

• optimizer: The optimizing routine used.

• call: The original function call.

• nCovars: The number of exogenous covariates fitted in the distance function. Does not include the intercept.

• LhoodType: The type of likelihood fitted. Currently, only 'parametric' types are fitted.

Description

Fits a detection function to off-transect distances collected by a single observer.

Usage

dE.single(
  data,
  formula,
  likelihood = "halfnorm",
  w.lo = setUnits(0, "m"),
  w.hi = NULL,
  expansions = 0,
  series = "cosine",
  x.scl = w.lo,
  g.x.scl = 1,
  warn = TRUE,
  outputUnits = NULL,
  asymptoticSE = TRUE
)

Arguments

Details

Optimization and estimation controls can be modified using options(). See RdistanceControls().

Value

An object of class 'dfunc' with the following components:

• par: The vector of estimated parameter values. Length of this vector is the sum of the following:

  1. The number of columns of the design matrix. This equals the number of covariates in the distance function plus one for the intercept, assuming an intercept is included.

  2. The number of constant parameters in the distance function. Constant parameters are those not related to covariates. For example, the exponent 'k' parameter for hazard rate likelihood, or the mixing fraction 'p' for the oneStep likelihood. This can be zero.

  3. The number of expansion functions called for. This equals the input expansions.

• loglik: The maximized value of the log likelihood.

• convergence: The convergence code. This code is returned by the optimizing routine (e.g., optim or nlminb). Values other than 0 indicate suspect convergence.

• message: If maximization did not converge (convergence != 0), this is the reason given by the optimizing routine.

• varcovar: The variance-covariance matrix for coefficients of the distance function, either estimated by the inverse of the fit's Hessian or by bootstrapping. If the likelihood is smooth (i.e., those listed by ⁠Rdistance:::differentiableLikelihoods())⁠, Rdistance initially estimates the variance-covariance matrix using the second derivative of the log likelihood surface at the final estimates, where second derivatives are estimated by numeric differentiation (in routine secondDeriv(). The variance-covariance matrix is re-set to NULL if the Hessian is not positive-definite. If bootstrap resampling has been performed (using abundEstim()), the variance-covariance matrix is re-estimated using the bootstrap values of parameters and automatically reset. Error estimates derived from bootstrapping are generally preferable to the asymptotic estimates, hence the automatic re-set.

• limits: A list containing the lower and upper limits of parameters.

• evaluations: The number of likelihood evaluations performed by the optimizer.

• mf: An R 'model frame' containing the detections (within the strip or circle) used in the fit, covariates specified in the formula, and groupsizes. Column 'dist' contains the observed distances. The intercept, if included in the model, is not included as a column in this model frame. (Test whether an intercept is included using attr(terms(return$mf), "intercept")). Column offset(...) contains group sizes associated with the values of dist. Name of the group size column is "offset(...)", not "groupsize(...)", so that group sizes can be treated offsets in other R routines. The mf component is a proper model.frame and contains both terms and contrasts attributes. This model frame contains only non-missing distances between w.lo and w.hi.

• data: The original nested data frame subset to information required to complete distance estimation. This data frame contains information on replication (i.e., rows are sites and are re-sampled during bootstrapping), missing distances, missing transect lengths, and distances outside the observation strip (below w.lo or above w.hi).

• formula: The distance function's formula.

• dataName: Name of the original nested data frame.

• likelihood: The name of the likelihood fitted to observation distances.

• w.lo: Left-truncation value used during the fit.

• w.hi: Right-truncation value used during the fit.

• expansions: The number of expansion terms.

• series: The type of expansion used during estimation. This is only relevant if expansions > 0.

• x.scl: The distance at which the function has been scaled to some value. This is the x at which g(x) = g.x.scl.

• g.x.scl: The height of the distance function at distance x.scl.

• outputUnits: A list of type symbolic_units containing the physical measurement units used during estimation.

• asymptoticSE: A logical scalar indication whether the variance-covariance matrix in component varcovar is asymptotic (TRUE; estimated from the Hessian) or bootstrap (FALSE; estimated by bootstrap resampling).

• optimizer: The optimizing routine used.

• call: The original function call.

• nCovars: The number of exogenous covariates fitted in the distance function. Does not include the intercept.

• LhoodType: The type of likelihood fitted. Currently, only 'parametric' types are fitted.

Group Sizes

To specify non-unity group sizes, use groupsize() on the RHS of formula. When group sizes are not all 1, they must appear in a column of the 'detections' list-column of data. For example, d ~ habitat + groupsize(number) specifies distances in column d, one covariate named habitat, and that column number contains the number of individuals associated with each detection. If group sizes are not specified, all group sizes are assumed to be 1.

Contrasts

Factor contrasts in Rdistance are specified the same way as in lm or glm. By default, Rdistance uses contrasts in getOption("contrasts"). To change contrasts, use a statement like options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly")). Or, to set contrasts for a specific factor in the input data frame, use contrasts(df$A) <- "contr.sum" or similar. See contrasts() or the contrasts.arg of model.matrix().

Transect types

Rdistance accommodates two kinds of transects: continuous and point. Detections can occur at any point on continuous transects. Rdistance calls these 'line-transects' even though routes are not necessarily a straight line. On point transects, detections occur at a series of stops (points). Rdisance calls these point-transects. Transects are the basic sampling unit in both cases. Rdistance assumes each row of data contains information from one transect. See RdistDf() for more details.

Measurement Units

As of Rdistance version 3.0.0, measurement units are require on all physical distances. Requiring units ensures that internal calculations and results (e.g., ESW and abundance) are correct and that output units are clear. Physical distances are required on off-transect distances, radial distances, truncation distances (w.lo, unless it is zero; and w.hi, unless it is NULL), scale locations (x.scl, unless it is zero), line-transect lengths, and study area size. All units are 1-dimensional except those on study area, which are 2-dimensional.

Physical measurement units can vary. For example, off-transect distances can be meters ("m"), w.hi can be inches ("in"), and w.lo can be kilometers ("km"). Internally, all distances are converted to the units specified by outputUnits (or the units of input distances if outputUnits is NULL), and all output is reported in units of outputUnits. Valid conversions must exist between units or an error is thrown (e.g., meters cannot convert into hectares).

Measurement units can be assigned using one of Rdistance's unit helper routines (see help(unitHelpers)), Rdistance's setUnits() function, or units::set_units() See units::valid_udunits() for a list of valid symbolic units.

If measurements are truly unit-less, or measurement units are unknown, set options(Rdist_requireUnits = FALSE). This suppresses all unit checks and conversions. Users are on their own here and must make sure all inputs are scaled correctly so that internal computations are correct and output units are known.

References

Buckland, S.T., D.R. Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. (2001) Introduction to distance sampling: estimating abundance of biological populations. Oxford University Press, Oxford, UK.

See Also

abundEstim(), autoDistSamp(). Likelihood-specific help files (e.g., halfnorm.like()).

Examples

# Load example sparrow data (line transect survey type)
data(sparrowDf)

dfunc <- dfuncEstim(data = sparrowDf
                  , formula = dist ~ 1)
dfunc
plot(dfunc)

Description

Fits a detection function using maximum likelihood.

Usage

dfuncEstim(data, ...)

Arguments

Details

Optimization and estimation controls can be modified using options(). See RdistanceControls().

Value

An object of class 'dfunc' with the following components:

• par: The vector of estimated parameter values. Length of this vector is the sum of the following:

  1. The number of columns of the design matrix. This equals the number of covariates in the distance function plus one for the intercept, assuming an intercept is included.

  2. The number of constant parameters in the distance function. Constant parameters are those not related to covariates. For example, the exponent 'k' parameter for hazard rate likelihood, or the mixing fraction 'p' for the oneStep likelihood. This can be zero.

  3. The number of expansion functions called for. This equals the input expansions.

• loglik: The maximized value of the log likelihood.

• convergence: The convergence code. This code is returned by the optimizing routine (e.g., optim or nlminb). Values other than 0 indicate suspect convergence.

• message: If maximization did not converge (convergence != 0), this is the reason given by the optimizing routine.

• varcovar: The variance-covariance matrix for coefficients of the distance function, either estimated by the inverse of the fit's Hessian or by bootstrapping. If the likelihood is smooth (i.e., those listed by ⁠Rdistance:::differentiableLikelihoods())⁠, Rdistance initially estimates the variance-covariance matrix using the second derivative of the log likelihood surface at the final estimates, where second derivatives are estimated by numeric differentiation (in routine secondDeriv(). The variance-covariance matrix is re-set to NULL if the Hessian is not positive-definite. If bootstrap resampling has been performed (using abundEstim()), the variance-covariance matrix is re-estimated using the bootstrap values of parameters and automatically reset. Error estimates derived from bootstrapping are generally preferable to the asymptotic estimates, hence the automatic re-set.

• limits: A list containing the lower and upper limits of parameters.

• evaluations: The number of likelihood evaluations performed by the optimizer.

• mf: An R 'model frame' containing the detections (within the strip or circle) used in the fit, covariates specified in the formula, and groupsizes. Column 'dist' contains the observed distances. The intercept, if included in the model, is not included as a column in this model frame. (Test whether an intercept is included using attr(terms(return$mf), "intercept")). Column offset(...) contains group sizes associated with the values of dist. Name of the group size column is "offset(...)", not "groupsize(...)", so that group sizes can be treated offsets in other R routines. The mf component is a proper model.frame and contains both terms and contrasts attributes. This model frame contains only non-missing distances between w.lo and w.hi.

• data: The original nested data frame subset to information required to complete distance estimation. This data frame contains information on replication (i.e., rows are sites and are re-sampled during bootstrapping), missing distances, missing transect lengths, and distances outside the observation strip (below w.lo or above w.hi).

• formula: The distance function's formula.

• dataName: Name of the original nested data frame.

• likelihood: The name of the likelihood fitted to observation distances.

• w.lo: Left-truncation value used during the fit.

• w.hi: Right-truncation value used during the fit.

• expansions: The number of expansion terms.

• series: The type of expansion used during estimation. This is only relevant if expansions > 0.

• x.scl: The distance at which the function has been scaled to some value. This is the x at which g(x) = g.x.scl.

• g.x.scl: The height of the distance function at distance x.scl.

• outputUnits: A list of type symbolic_units containing the physical measurement units used during estimation.

• asymptoticSE: A logical scalar indication whether the variance-covariance matrix in component varcovar is asymptotic (TRUE; estimated from the Hessian) or bootstrap (FALSE; estimated by bootstrap resampling).

• optimizer: The optimizing routine used.

• call: The original function call.

• nCovars: The number of exogenous covariates fitted in the distance function. Does not include the intercept.

• LhoodType: The type of likelihood fitted. Currently, only 'parametric' types are fitted.

Group Sizes

To specify non-unity group sizes, use groupsize() on the RHS of formula. When group sizes are not all 1, they must appear in a column of the 'detections' list-column of data. For example, d ~ habitat + groupsize(number) specifies distances in column d, one covariate named habitat, and that column number contains the number of individuals associated with each detection. If group sizes are not specified, all group sizes are assumed to be 1.

Contrasts

Factor contrasts in Rdistance are specified the same way as in lm or glm. By default, Rdistance uses contrasts in getOption("contrasts"). To change contrasts, use a statement like options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly")). Or, to set contrasts for a specific factor in the input data frame, use contrasts(df$A) <- "contr.sum" or similar. See contrasts() or the contrasts.arg of model.matrix().

Measurement Units

As of Rdistance version 3.0.0, measurement units are require on all physical distances. Requiring units ensures that internal calculations and results (e.g., ESW and abundance) are correct and that output units are clear. Physical distances are required on off-transect distances, radial distances, truncation distances (w.lo, unless it is zero; and w.hi, unless it is NULL), scale locations (x.scl, unless it is zero), line-transect lengths, and study area size. All units are 1-dimensional except those on study area, which are 2-dimensional.

Physical measurement units can vary. For example, off-transect distances can be meters ("m"), w.hi can be inches ("in"), and w.lo can be kilometers ("km"). Internally, all distances are converted to the units specified by outputUnits (or the units of input distances if outputUnits is NULL), and all output is reported in units of outputUnits. Valid conversions must exist between units or an error is thrown (e.g., meters cannot convert into hectares).

Measurement units can be assigned using one of Rdistance's unit helper routines (see help(unitHelpers)), Rdistance's setUnits() function, or units::set_units() See units::valid_udunits() for a list of valid symbolic units.

If measurements are truly unit-less, or measurement units are unknown, set options(Rdist_requireUnits = FALSE). This suppresses all unit checks and conversions. Users are on their own here and must make sure all inputs are scaled correctly so that internal computations are correct and output units are known.

References

Buckland, S.T., D.R. Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. (2001) Introduction to distance sampling: estimating abundance of biological populations. Oxford University Press, Oxford, UK.

See Also

abundEstim(), autoDistSamp(). Likelihood-specific help files (e.g., halfnorm.like()).

Examples

# Sparrow line transect example
data(sparrowDetectionData)
data(sparrowSiteData)

sparrowDf <- RdistDf(sparrowSiteData, sparrowDetectionData)

dfunc <- dfuncEstim(sparrowDf, 
                    formula = dist ~ 1
                  )
summary(dfunc)

  
data(sparrowDfuncObserver) # pre-estimated object
## Not run:                  
# Commands that produced 'sparrowDfuncObserver'
sparrowDfuncObserver <- sparrowDf |> 
         dfuncEstim( 
           formula = dist ~ observer
         )

## End(Not run)     
sparrowDfuncObserver
summary(sparrowDfuncObserver)
plot(sparrowDfuncObserver)
plot(sparrowDfuncObserver
   , newdata = data.frame(observer = c("obs1", "obs2", "obs3"
                                     , "obs4", "obs5")))

Description

Utility function to produce error messages suitable for stop

Usage

dfuncEstimErrMessage(txt, attri)

Arguments

Value

A string

Description

Return a character vector of likelihoods which are differentiable and hence the second derivative method can be used to estimate variance-covariance. Any likelihoods not in this list must use bootstrapping. This vector is polled in a few Rdistance routines, notably parseModel and varcovarEstim.

Usage

differentiableLikelihoods()

Value

A character vector of differentiable likelihoods.

Description

Extract the observation distances (i.e., responses for an Rdistance model) from an Rdistance model frame.

Usage

distances(ml, na.rm = TRUE, ...)

Arguments

Value

A vector containing observation distances contained in the Rdistance model frame.

Examples

data(sparrowDf)
sparrowModel <- parseModel( sparrowDf, dist ~ observer )
stats::model.response(sparrowModel$mf)
distances(sparrowModel) # same, but future-proof

Description

Computes Effective Detection Radius (EDR) for estimated detection functions on point transects. See ESW() is for line transects.

Usage

EDR(object, newdata = NULL)

Arguments

Details

Effective Detection Radius is the integral under the detection function times distance.

Value

If newdata is present, the returned value is a vector of effective sampling distances associated with covariate values in newdata. Length of return in this case is the number of rows in newdata. If newdata is NULL, the returned value is a vector of effective sampling distances associated with covariate values in object. Length of return in this case is the number of detected groups. The returned vector has measurement units, i.e., object$outputUnits.

See Also

dfuncEstim(), ESW(), effectiveDistance()

Examples

# Load example thrasher data (point transect survey type)
data(thrasherDf)

# Fit half-normal detection function
dfunc <- thrasherDf |> dfuncEstim(formula=dist~bare)

# Compute effective detection radius (EDR)
EDR(dfunc) # vector length 192
effectiveDistance(dfunc) # same
EDR(dfunc, newdata = data.frame(bare=30)) # vector length 1

Description

Computes Effective Strip Width (ESW) for line-transect detection functions, or the analogous Effective Detection Radius (EDR) for point-transect detection functions.

Usage

effectiveDistance(object, newdata = NULL)

Arguments

Details

Serves as a wrapper for ESW() and EDR().

Effective distances are areas under scaled distance functions (i.e., area under g(x)). Areas are exact for functions whose integral is known (e.g., negexp). Numeric integration is used for all others.

Value

If newdata is present, the returned value is a vector of effective sampling distances associated with covariate values in newdata. Length of return in this case is the number of rows in newdata. If newdata is NULL, the returned value is a vector of effective sampling distances associated with covariate values in object. Length of return in this case is the number of detected groups. The returned vector has measurement units, i.e., object$outputUnits.

See Also

dfuncEstim(), ESW(), EDR(), integrateNumeric(), integrateNegexpLines()

Description

Extract effort information from an Rdistance data frame. Effort is length of line-transects or number of points on point-transects.

Usage

effort(x, ...)

Arguments

Value

A vector containing effort. If line-transects, return is length of transects, with units. If point-transects, return is number of points (integers, no units). Vector length is number of transects. If input is not an RdistDf or estimated distance function, return is NULL.

Examples

data(sparrowDf)
effort(sparrowDf)
fit <- dfuncEstim(sparrowDf, dist ~ 1)
effort(fit)

Description

Constructs a string stating what is "unknown" that is suitable for use in warning and error functions.

Usage

errDataUnk(txt, attri)

Arguments

Value

A descriptive string, suitable for warning or error.

Description

Estimate abundance from an Rdistance fitted model. This function is called internally by abundEstim. Most users will call abundEstim to estimate abundance unless they are running simulations or bootstrapping.

Usage

estimateN(object, area = NULL, propUnitSurveyed = 1)

Arguments

Details

The abundance estimate for line-transect surveys (if no covariates are included in the detection function and both sides of the transect are observed) is

$N =\frac{n(A)}{2(ESW)(L)}$

where n is total number of sighted individuals (i.e., sum(groupSizes(dfunc))), L is the total length of surveyed transect (i.e., sum(effort(dfunc))), and ESW is effective strip width computed from the estimated distance function (i.e., ESW(dfunc)). If only one side of transects were observed, the "2" in the denominator is not present (or, replaced with a "1").

The abundance estimate for point transect surveys (if no covariates are included) is

$N =\frac{n(A)}{\pi(ESR^2)(P)}$

where n is total number of sighted individuals (i.e., sum(groupSizes(dfunc))), P is the total number of surveyed points (i.e., sum(effort(dfunc))), and ESR is effective search radius computed from the estimated distance function (i.e., ESR(dfunc)).

This routine, abundEstim, estimates abundance on the entire study area. Site-specific density estimates are computed by predict(x, type = "density"), which returns a tibble containing density and abundance on the area surveyed by every transect.

Value

A list containing the following components:

For line-transects that do not involve covariates, object$density is object$n.seen / (2 * propUnitSurveyed * object$w * object$pDetection * object$surveyedUnits)

See Also

dfuncEstim(), abundEstim()

Description

Returns effective strip width (ESW) for line-transect detection functions. See EDR() is for point transects.

Usage

ESW(object, newdata = NULL)

Arguments

Details

ESW is area under a scaled distance function between its left-truncation limit (obj$w.lo) and its right-truncation limit (obj$w.hi).

If detection does not decline with distance, the detection function is flat (horizontal), and area under the detection function is $g(0)(w.hi - w.lo)$. If, in this case, $g(0) = 1$, effective sampling distance is the half-width of the surveys, $(w.hi - w.lo)$

Value

If newdata is present, the returned value is a vector of effective sampling distances associated with covariate values in newdata. Length of return in this case is the number of rows in newdata. If newdata is NULL, the returned value is a vector of effective sampling distances associated with covariate values in object. Length of return in this case is the number of detected groups. The returned vector has measurement units, i.e., object$outputUnits.

See Also

dfuncEstim(), EDR(), effectiveDistance()

Examples

data(sparrowDfuncObserver)

ESW(sparrowDfuncObserver) # vector length 356 = number of groups
ESW(sparrowDfuncObserver, newdata = data.frame(observer = 
   c("obs2", "obs4"))) # vector length 2

Description

The domain to which expansion factors are applied varies by likelihood. When the domain varies, it depends on parameter values. For example, oneStep expansion factors are applied between 0 and the likelihood's first parameter ($\theta$), which varies by covariates. This function computes the likelihood-specific expansion domains

Usage

expandW(ml, params, k)

Arguments

Value

A plain vector of length k containing expansion domains, with units.

Description

Compute "expansion" terms that modify the shape of a base distance function.

Usage

expansionTerms(a, d, series, nexp, w)

Arguments

Details

Expansion terms modify the base likelihood function and are used to incorporate "wiggle". The modified distance function is, key * expTerms where key is a vector of values in the base likelihood function (e.g., halfnorm.like()$L.unscaled) and expTerms is the matrix returned by this routine. In equation form,

$f(x_i|\beta,a_1,a_2,\dots,a_m) = f(x_i|\beta)(1 + \sum_{k=1}^{m} a_k h_{k}(x_i/w)).$

, where $m$ = the the number of expansions (nexp), $h_{j}(x)$ are expansion terms for distance $x$, and $a_1, a_2, \dots, a_m$ are the (estimated) expansion term coefficients.

Value

If nexp equals 0, 1 is returned. If nexp is greater than 0, a matrix of size $n$X$k$ containing expansion terms, where $n$ = length(d) and $k$ = nrow(a). The expansion series associated with row $j$ of a are in column $j$ of the return. i.e., element ($i$,$j$) of the return is

$1 + \sum_{k=1}^{m} a_{jk} h_{k}(x_{i}/w).$

Examples

a1 <- c(log(40), 0.5, -.5)
a2 <- c(log(40), 0.25, -.5)
dists <- seq(0, 100, length = 100) %m%.
w = 100 %m%.

expTerms1 <- expansionTerms(a1, dists, "cosine", 2, w)
expTerms2 <- expansionTerms(a2, dists, "cosine", 2, w)
plot(dists, expTerms2, ylim = c(0,2.5))
points(dists, expTerms1, pch = 16)

# Same as above
a <- rbind(a1, a2)
w <- rep(w, nrow(a))
expTerms <- expansionTerms(a, dists, "cosine", 2, w)
matlines(dists, expTerms, lwd=2, col=c("red", "blue"), lty=1)

# Showing key and expansions
key <- halfnorm.like(log(40), dists, matrix(1,length(dists),1))$L.unscaled
plot(dists, key, type = "l", col = "blue", ylim=c(0,1.5))
lines(dists, key * expTerms1, col = "red")
lines(dists, key * expTerms2, col = "purple")

Description

Evaluate the gamma distance function for sighting distances, potentially including covariates and expansion terms

Usage

Gamma.like(a, dist, covars, w.hi = NULL)

Arguments

Details

The Rdistance implementation of a Gamma distance function follows Becker and Quang (2009). Rdistance's Gamma distance function is

$f(d|\alpha, \sigma) = \left(\frac{d}{m}\right)^{\alpha - 1}e^{-(d-m)/\sigma},$

where $\alpha$ is the shape parameter, $\sigma$ is the scale parameter, and $m = (\alpha-1)\sigma$. $m$ is the mode of the Gamma function, and in Rdistance it's scaled to have a maximum of 1.0 at $m$. The scale parameter is a function of the shape parameter and sighting covariates, i.e.,

$\sigma = k [exp(x'\beta)],$

where $x$ is a vector of covariate values associated with distance $d$ (i.e., a row of covars), $\beta$ is a vector of the first $q$ (=ncol(covars)) values of the first argument of the function (a), and $k$ is a function of the shape parameter, i.e.,

$k = \frac{1}{\Gamma(\alpha)} \left(\frac{a - 1}{e^1} \right)^{a - 1}.$

The shape parameter $\alpha$ is the $q+1$-st value in the function's first argument and is constrained to be strictly greater than 1.0.

See Examples for use of GammaReparam() to compute $\alpha$ and $\sigma$ from fitted object coefficients.

Value

A list containing the following two components:

• L.unscaled: A matrix of size nXk (n = length dist; k = number of cases = nrow(a)) containing likelihood values evaluated at distances in dist. Each row is associated with a single distance, and each column is associated with a single case (row of a). Values in this matrix are the distance function $g(d)$ which generally have $g(0) = 1$. These values are "unscaled" likelihood values; they must be scaled (divided by) with the area under $g(x)$ between w.lo and w.hi to form proper likelihood values.

• params: A nXbXk array of the likelihood's (canonical) parameters in link space (i.e., on log scale; b = number of canonical parameters in the likelihood; k = number of cases). Rows correspond to distances in dist. Columns correspond to parameters (columns of a), and pages correspond to cases (rows of a).

References

Becker, E. F., and P. X. Quang, 2009. A Gamma-Shaped Detection Function for Line-Transect Surveys with Mark-Recapture and Covariate Data. Journal of Agricultural, Biological, and Environmental Statistics 14(2):207-223.

See Also

dfuncEstim(), abundEstim(), other ⁠<likelihood>.like⁠ functions

Examples

x <- seq(0, 100, length=100)
covars <- matrix(1,100,1)

# Plots showing changes in scale
plot(x, Gamma.like(c(log(20),2.5), x, covars)$L.unscaled, type="l", col="red")
lines(x, Gamma.like(c(log(40),2.5), x, covars)$L.unscaled, col="blue")

# Plots showing changes in shape
plot(x, Gamma.like(c(log(20),1.5), x, covars)$L.unscaled, type="l", col="red")
lines(x, Gamma.like(c(log(20),2.5), x, covars)$L.unscaled, col="blue")
lines(x, Gamma.like(c(log(20),4.5), x, covars)$L.unscaled, col="green")

# Roll-your-own plot, showing "re-parameterization":
# Assume fitted object coefficients are c(log(20), 4.5)
fit <- list(par = c(log(20), 4.5))

# The distance function is then,
gammaPar <- GammaReparam( scl = exp(fit$par[1])
                        , shp = fit$par[2] ) # returns scl=k*exp(x'B)
scl <- gammaPar$scl
shp <- gammaPar$shp
m <- (shp - 1) * scl
g <- (x / m)^(shp - 1) * exp(-(x - m) / scl) # distance function
lines(x, g, lwd = 3, lty = 2, col="green3")

Description

Compute starting values and limits for the Gamma distance function.

Usage

Gamma.start.limits(ml)

Arguments

Value

A list containing the following components

The length of each vector in the return is: ⁠(Num expansions) + 1 + 1*(like %in% c("hazrate")) + (Num Covars)⁠.

Description

Compute mode (i.e., maximum) of Rdistance's version of the gamma distribution.

Usage

GammaModes(params)

Arguments

Value

A vector of the locations of the gamma modes associated with each row in the params matrix.

Description

Transform Rdistance's version of the Gamma distribution parameters, which is that of Becker and Quan, into the version for use in R::dgamma() and elsewhere.

Usage

GammaReparam(shp, scl)

Arguments

Value

A list with components $shp and $scl, which are the re-parameterized versions of the input parameters suitable for us in R::dgamma().

See Also

'Details' section of Gamma.like() for Rdistance's Gamma distribution

Examples

# Rdistance Gamma parameters
Rd.scl <- 50  # must be >0
Rd.shp <- 1.5 # must be >1

# dgamma parameters
dgParams <- GammaReparam(Rd.shp, Rd.scl)
dgParams

# Gamma distribution with (Rd.scl, Rd.shp) from 0 to 100
curve(dgamma(x, shape=dgParams$shp, scale = dgParams$scl)
           , from = 0
           , to = 100)
 
# Rdistance's version: same curve but scaled so maximum = 1
x <- seq(0, 100, length = 200) 
scl <- dgParams$scl
shp <- dgParams$shp
m <- (shp - 1) * scl
g <- (x / m)^(shp - 1) * exp(-(x - m) / scl) # distance function
plot(x, g, type = "l")

Description

Set the number of cores for parallel operations. Also convert integer 'parallel' to TRUE-FALSE for ease.

Usage

getNCores(parallel)

Arguments

Details

Input parallel <= 0 is converted to 1. Input parallel > maxCores is converted to maxCores. If input parallel is numeric, is first converted to integer by rounding down.

Value

A list with components $parallel and $cores. $parallel is logical, T for parallel operations, F o.w. $cores is the number of cores to use during parallel operations. $cores is integer in the range 1, 2, ..., max(Available cores).

Description

Extract the group size information from an Rdistance model frame.

Usage

groupSizes(ml, ...)

Arguments

Value

A vector containing group sizes contained in the Rdistance model frame or fitted object.

Examples

data("sparrowDf")
sparrowModel <- parseModel( sparrowDf, dist ~ observer )
stats::model.offset(sparrowModel$mf)
groupSizes(sparrowModel)  # same, but future-proof

sparrowModel <- parseModel( sparrowDf
                 , dist ~ observer + groupsize(groupsize) )
groupSizes(sparrowModel)

Description

Estimate distance function scaling factor , g(0) or g(x), for a specified distance function.

Usage

gxEstim(fit)

Arguments

Details

This routine scales sightability such that g(x.scl) = g.x.scl, where g() is the sightability function. Specification of x.scl and g.x.scl covers several estimation cases:

1. g(0) = 1 : (the default) Inputs are x.scl = 0, g.x.scl = 1. If w.lo > 0, x.scl will be set to w.lo so technically this case is g(w.low) = 1.

2. User supplied probability at specified distance: Inputs are x.scl = a number greater than or equal to w.lo, g.x.scl = a number between 0 and 1. This case covers situations where sightability on the transect (distance 0) is not perfect. This case assumes researchers have an independent estimate of sightability at distance x.scl off the transect. For example, researchers could be using multiple observers to estimate that sightability at distance x.scl is g.x.scl.

3. Maximum sightability specified: Inputs are x.scl="max", g.x.scl = a number between 0 and 1. In this case, g() is scaled such that its maximum value is g.x.scl. This routine computes the distance at which g() is maximum, sets g()'s height there to g.x.scl, and returns x.max where x.max is the distance at which g is maximized. This case covers the common aerial survey situation where maximum sightability is slightly off the transect, but the distance at which the maximum occurs is unknown.

4. Double observer system: Inputs are x.scl="max", g.x.scl = \<a data frame\>. In this case, g(x) = h, where x is the distance that maximizes g and h is the height of g() at x computed from the double observer data frame (see below for structure of the double observer data frame).

5. Distance of independence specified, height computed from double observer system: Inputs are x.scl = a number greater than or equal to w.lo g.x.scl = a data frame. In this case, g(x.scl) = h, where h is computed from the double observer data frame (see below for structure of the double observer data frame).

When x.scl, g.x.scl, or observer are NULL, the routine will look for and use ⁠$call.x.scl⁠, or ⁠$call.g.x.scl⁠, or ⁠$call.observer⁠ components of the fit object for whichever of these three parameters is missing. Later, different values can be specified in a direct call to F.gx.estim without having to re-estimate the distance function. Because of this feature, the default values in dfuncEstim are x.scl = 0 and g.x.scl = 1 and observer = "both".

Value

A list comprised of the following components:

Structure of the double observer data frame

When g.x.scl is a data frame, it is assumed to contain the components ⁠$obsby.1⁠ and ⁠$obsby.2⁠ (no flexibility on names). Each row in the data frame contains data from one sighted target. The ⁠$obsby.1⁠ and ⁠$obsby.2⁠ components are TRUE/FALSE (logical) vectors indicating whether observer 1 (obsby.1) or observer 2 (obsby.2) spotted the target.

See Also

dfuncEstim()

Examples

data(sparrowDf)
fit <- dfuncEstim(sparrowDf, dist ~ groupsize(groupsize))
gxEstim(fit)
  
fit <- dfuncEstim(sparrowDf, dist ~ groupsize(groupsize)
                , x.scl = 50 %m%.
                , g.x.scl = 0.75)
gxEstim(fit)
plot(fit)
abline(h=0.75)
abline(v=50%m%.)

Description

Evaluate the half-normal distance function, for sighting distances, potentially including covariates and expansion terms

Usage

halfnorm.like(a, dist, covars, w.hi = NULL)

Arguments

Details

The half-normal distance function is

$f(d|\sigma) = \exp(-\frac{d^2}{2\sigma^2})$

where $\sigma = exp(x'a)$, $x$ is a vector of covariate values associated with distance $d$ (i.e., a row of covars), and $a$ is a vector of the first $q$ (=ncol(covars)) values in argument a.

Some authors parameterize the halfnorm without the "2" in the denominator of the exponent. Rdistance includes "2" in this denominator to make quantiles of the half normal agree with the standard normal. This means that half-normal coefficients in Rdistance (i.e., $\sigma = exp(x'a)$) can be interpreted as normal standard errors. Approximately 95% of distances should occur between 0 and 2$\sigma$.

Value

A list containing the following two components:

• L.unscaled: A matrix of size nXk (n = length dist; k = number of cases = nrow(a)) containing likelihood values evaluated at distances in dist. Each row is associated with a single distance, and each column is associated with a single case (row of a). Values in this matrix are the distance function $g(d)$ which generally have $g(0) = 1$. These values are "unscaled" likelihood values; they must be scaled (divided by) with the area under $g(x)$ between w.lo and w.hi to form proper likelihood values.

• params: A nXbXk array of the likelihood's (canonical) parameters in link space (i.e., on log scale; b = number of canonical parameters in the likelihood; k = number of cases). Rows correspond to distances in dist. Columns correspond to parameters (columns of a), and pages correspond to cases (rows of a).

See Also

dfuncEstim(), abundEstim(), other ⁠<likelihood>.like⁠ functions

Examples

d <- seq(0, 100, length=100)
covs <- matrix(1,length(d),1)
halfnorm.like(log(20), d, covs)

plot(d, halfnorm.like(log(20), d, covs)$L.unscaled, type="l", col="red")
lines(d, halfnorm.like(log(40), d, covs)$L.unscaled, col="blue")

# Evaluate 3 functions at once using matrix of coefficients:
# sigma ~ 20, 30, 40
coefs <- matrix(log(c(7.39,7.33, 4.48,44.80, 2.72,216.54))
               , byrow = TRUE
               , ncol=2) # (3 coef vectors)X(2 covars)
covs <- matrix(c(rep(1,length(d))
               , rep(0.5,length(d)))
               , nrow = length(d)) # 100 X 2
L <- halfnorm.like( coefs, d, covs ) 
L$L.unscaled # 100 X (3 coef vectors)
L$params     # 100 X (3 coef vectors); ~ log(c(20,30,40))
matplot(d, L$L.unscaled, type="l")

Description

Compute starting values and limits for the half normal distance function.

Usage

halfnorm.start.limits(ml)

Arguments

Value

A list containing the following components

The length of each vector in the return is: ⁠(Num expansions) + 1 + 1*(like %in% c("hazrate")) + (Num Covars)⁠.

Description

Computes the hazard rate distance function.

Usage

hazrate.like(a, dist, covars, w.hi = NULL)

Arguments

Details

The hazard rate likelihood is

$f(x|\sigma,k) = 1 - \exp(-(x/\sigma)^{-k})$

where $\sigma$ determines location (i.e., distance at which the function equals 1 - exp(-1) = 0.632), and $k$ determines slope of the function at $\sigma$ (i.e., larger k equals steeper slope at $\sigma$). For distance analysis, the valid range for both $\sigma$ and k is $\geq 0$.

Value

A list containing the following two components:

• L.unscaled: A matrix of size nXk (n = length dist; k = number of cases = nrow(a)) containing likelihood values evaluated at distances in dist. Each row is associated with a single distance, and each column is associated with a single case (row of a). Values in this matrix are the distance function $g(d)$ which generally have $g(0) = 1$. These values are "unscaled" likelihood values; they must be scaled (divided by) with the area under $g(x)$ between w.lo and w.hi to form proper likelihood values.

• params: A nXbXk array of the likelihood's (canonical) parameters in link space (i.e., on log scale; b = number of canonical parameters in the likelihood; k = number of cases). Rows correspond to distances in dist. Columns correspond to parameters (columns of a), and pages correspond to cases (rows of a).

See Also

dfuncEstim(), abundEstim(), other ⁠<likelihood>.like⁠ functions

Examples

d <- seq(0, 100, length=100)
covs <- matrix(1,length(d),1)
hazrate.like(c(log(20), 5), d, covs)

# Changing location parameter
plot(d, hazrate.like(c(log(20), 5), d, covs)$L.unscaled, type="l", col="red")
lines(d, hazrate.like(c(log(40), 5), d, covs)$L.unscaled, col="blue")
abline(h = 1 - exp(-1), lty = 2)
abline(v = c(20,40), lty = 2)

# Changing slope parameter
plot(d, hazrate.like(c(log(50), 20), d, covs)$L.unscaled, type="l", col="red")
lines(d, hazrate.like(c(log(50), 2), d, covs)$L.unscaled, col="blue")
abline(h = 1 - exp(-1), lty = 2)
abline(v = 50, lty = 2)

Description

Compute starting values and limits for the hazard rate distance function.

Usage

hazrate.start.limits(ml)

Arguments

Value

A list containing the following components

The length of each vector in the return is: ⁠(Num expansions) + 1 + 1*(like %in% c("hazrate")) + (Num Covars)⁠.

Description

Computes Hermite expansion terms for use in distance analysis. The Hermite (and other expansions) allow "wiggle" in estimated distance functions.

Usage

hermite.expansion(x, expansions)

Arguments

Details

There are, in general, several expansions that can be called Hermite. Let $w = 4x - 2$. Rdistance's Hermite expansions are:

• First term:

  $h_1(w) = w + 2,$

• Second term:

  $h_2(w) = w^2 - 4,$

• Third term:

  $h_3(w) = w^3 - 3w + 2,$

• Fourth term:

  $h_4(w) = w^4 - 6w^2 + 8,$

The maximum number of expansion terms computed is 4.

Value

A 3D array of size nrow(x) X ncol(x) X expansions. The 'pages' (3rd dimension) of this array are the cosine expansions of x. i.e., page 1 is the first expansion term of x, page 2 is the second expansion term of x, etc.

See Also

dfuncEstim() , cosine.expansion() , sine.expansion() , simple.expansion().

Examples

x <- matrix(seq(0, 1, length = 200), ncol = 1)
herm.expn <- hermite.expansion(x, 4)
plot(range(x), range(herm.expn), type="n")
matlines(x, herm.expn[,1,1:4], col=rainbow(4), lty = 1)

Description

Call R native function 'nlminb' to perform optimization.

Usage

HookeJeeves(ml, strt.lims)

Arguments

Value

A list with following named components:

• par = parameters

• loglik = objective function value at minimum

• convergence = 0 for yes, other for no

• iterations = number of iterations

• evaluations = function evaluations

• message = a convergence message

• varcovar = a variance covariance matrix of parameters

• limits = low and high limits

Description

Computes the cumulative function of the huber likelihood. The cumulative function is proportional to the huber cumulative distribution function, differing only by appropriate scaling constant.

Usage

huber.cumFunc(x, t1, t2, p, w)

Arguments

Value

A vector of values from the huber cumulative function. The huber cumulative function is

$F(x|\theta_1,\theta_2,p) = \int_0^x f(y|\theta_1,\theta_2,p) dy,$

where $f(y|\theta_1,\theta_2,p)$ is Rdistance's huber likelihood. The only difference between this cumulative function, and the cumulative distribution function is the scaling constant. That is, the maximum of the cumulative function is greater than 1 while the maximum cumulative distribution function is exactly 1.

See Also

huber.like()

Examples

d <- -10:210

# Cumulative function
fd <- huber.cumFunc(d, 125, 25, .05, 200) 
plot(d, fd, type="l")

# Cumulative distribution function
Fd <- fd / huber.cumFunc(200, 125, 25, .05, 200)

Description

Computes the Huber likelihood for use as a distance function. The Huber likelihood is a mixture of an inverted Huber loss and a uniform distribution. The distance function is quadratic from the lowest distance to its first parameter, linear from the first parameter to the sum of its first and second parameter, and constant after that.

Usage

huber.like(a, dist, covars, w.hi = NULL)

Arguments

Details

The 'huber' likelihood is an inverted version of the Huber loss function, mixed with a uniform at the upper end, and contains three canonical parameters. The 'huber' distance function is a negative quadratic between 0 and its first parameter $\theta_1$, linear between $\theta_1$ and $\theta_1 + \theta_2$, and constant after $\theta_1 + \theta_2$. Specifically, the 'huber' likelihood is,

$f(d|\theta_1,\theta_2, p) = \left(1 - \frac{(1-p)h(d)}{m}\right) \, I(0 \leq d \leq \Theta) \ + \ p \, I(\Theta < d \leq w)$

where $\Theta = \theta_1 + \theta_2$, $w$ = w.hi - w.lo, $m = \theta_1(\Theta - 0.5\theta_1)$ and h(d) is Huber loss between 0 and $\Theta$, i.e.,

$h(d) = 0.5d^2\, I(0 \leq d \leq \theta_1) \ + \ \theta_1(d - 0.5\theta_1)\, I(\theta_1 < d \leq \Theta).$

The first parameter, $\theta_1$, is related to covariates, i.e., $log(\theta_1) = \beta_0 + \beta_1x_1 + \ldots + \beta_qx_q$. $\theta_2$ and $p$ are constant across covariate values.

Value

A list containing the following two components:

• L.unscaled: A matrix of size nXk (n = length dist; k = number of cases = nrow(a)) containing likelihood values evaluated at distances in dist. Each row is associated with a single distance, and each column is associated with a single case (row of a). Values in this matrix are the distance function $g(d)$ which generally have $g(0) = 1$. These values are "unscaled" likelihood values; they must be scaled (divided by) with the area under $g(x)$ between w.lo and w.hi to form proper likelihood values.

• params: A nXbXk array of the likelihood's (canonical) parameters in link space (i.e., on log scale; b = number of canonical parameters in the likelihood; k = number of cases). Rows correspond to distances in dist. Columns correspond to parameters (columns of a), and pages correspond to cases (rows of a).

See Also

See Rdistance tutorials for a method to generate random observations from the Huber likelihood.

Examples

t1 <- c(65,80,120)
t2 <- 60
p  <- 0.05
a <- matrix(c(log(t1)
       , rep(t2,3)
       , rep(p,3))
       , 3,3)
d <- seq(0, 200, length=201)
X <- matrix(1,length(d),1)
y <- huber.like(a, d, X)

# Plot showing covariate effects 
plot(range(d), range(y$L.unscaled)
  , type = "n", xlab = "d", ylab = "Huber(d)")
for(i in 1:3){
  # Distance functions
  lines(d
    , y$L.unscaled[,i]
    , col = i
    , lwd= 2)
  # Quadradic to linear transitions
  points(exp(a[i,1])
    , y$L.unscaled[(t1[i]-0.1) < d & d < (t1[i]+0.01),i]
    , pch = 16
    , col = i )
  # Linear to constant transition
  Theta <- exp(a[i,1]) + a[i,2]
  points(Theta
    , y$L.unscaled[(Theta-0.1) < d & d < (Theta+0.1),i]
    , pch = 15
    , col = i )
}

Description

Computes starting values and limits for the 'huber' distance function.

Usage

huber.start.limits(ml)

Arguments

Value

A list containing the following components

The length of each vector in the return is: ⁠(Num expansions) + 1 + 1*(like %in% c("hazrate")) + (Num Covars)⁠.

Description

Compute break points in a oneStep likelihood and insert them into a sequence of distances. The idea is to insert a point just left and just right of the breaks so that they plot as vertical lines.

Usage

insertOneStepBreaks(obj, newData, xseq)

Arguments

Value

A vector like xseq, but with the break points inserted.

Description

Integrates under distances functions using exact integrals when possible. If exact integrals are not known, numerical integration is used.

Usage

integrateDfuncs(object, ml)

Arguments

Details

Let K be the integral under distance function g(x) (i.e., the output from this routine). In distance analysis, the observation likelihood being evaluated for maximization is the density, f(x) = g(x)/K. K is a key quantity in distance analysis and is called the "effective sampling distance".

Value

A vector of areas under the distance functions represented in object. If object is a distance function and newdata is specified, the returned vector's length is nrow(newdata). If object is a distance function and newdata is NULL, returned vector's length is length(distances(object)). If object is a matrix, return's length is nrow(object).

Note

Users will not normally call this function. It is called internally by nLL() and effectiveDistance().

Examples

# Faking a model frame
ml <- list( likelihood = "halfnorm"
          , expansions = 0
          , w.lo = 0 %m% .
          , w.hi = 100 %m% .
          , Units = "m"
          )
attr(ml, "transType") <- "line"

parms <- matrix(75, nrow = 1) 
integrateDfuncs(parms, ml)

# check: Normal, 0 to 100, sd = 75, scaled to mode = 1
(pnorm(q = 100, mean = 0, sd = 75) - 0.5) * sqrt(2*pi)*75

Description

Compute integral of the Gamma distance function for line-transect surveys.

Usage

integrateGammaLines(
  object,
  newdata = NULL,
  w.lo = NULL,
  w.hi = NULL,
  Units = NULL
)

Arguments

Details

#' Returned integrals are

$\int_0^{w} \left(\frac{x}{m}\right)^{\alpha -1} e^{-(x - m)/\sigma_i}dx,$

where $w = w.hi - w.lo$, $\sigma_i$ is the i-th estimated scale parameter for the Gamma distance function, and $m$ is the mode of Gamma (i.e., $(\alpha - 1)\sigma_i$. Rdistance computes the integral using R's base function pgamma(), which for all intents and purposes is exact. See also Gamma.like().

Value

A vector of areas under the distance functions represented in object. If object is a distance function and newdata is specified, the returned vector's length is nrow(newdata). If object is a distance function and newdata is NULL, returned vector's length is length(distances(object)). If object is a matrix, return's length is nrow(object).

Note

Users will not normally call this function. It is called internally by nLL() and effectiveDistance().

See Also

integrateNumeric(); integrateNegexpLines(); integrateOneStepLines()

Examples

# Fake distance function object w/ minimum inputs for integration
d <- rep(1,4) %m%. # Only units needed, not values
obs <- factor(rep(c("obs1", "obs2"), 2))
beta <- c(4.0, -0.5, 1.5) # {'Intercept', b_1, shape}
w.hi <- 125
w.lo <- 20
ml <- list(
    mf = model.frame(d ~ obs)
  , par = beta 
  , likelihood = "Gamma"
  , w.lo = w.lo %#% "m"
  , w.hi = w.hi %#% "m"
  , expansions = 0
)
class(ml) <- "dfunc"
integrateGammaLines(ml)

# Check: Integral of Gamma density from 0 to w.hi-w.lo
b <- exp(c(beta[1], beta[1] + beta[2]))
B <- Rdistance::GammaReparam(shp = beta[3], scl = b)
m <- (B$shp - 1)*B$scl
f.at.m <- dgamma(m, shape = B$shp, scale = B$scl)
intgral <- pgamma(q = w.hi - w.lo, shape = B$shp, scale = B$scl) / f.at.m
intgral

Description

Compute integral of the half-normal distance function for line-transect surveys.

Usage

integrateHalfnormLines(
  object,
  newdata = NULL,
  w.lo = NULL,
  w.hi = NULL,
  Units = NULL
)

Arguments

Details

Returned integrals are

$\int_0^{w} e^{-x^2/2\sigma_i^2} dx = \sqrt{2\pi}\sigma_i (\Phi(w) - 0.5),$

where $w = w.hi - w.lo$, $\sigma_i$ is the estimated half-normal distance function parameter for the i-th observed distance, and $\Phi$ is the standard normal cumulative probability function. Rdistance uses R's base function pnorm(), which for all intents and purposes is exact.

Value

A vector of areas under the distance functions represented in object. If object is a distance function and newdata is specified, the returned vector's length is nrow(newdata). If object is a distance function and newdata is NULL, returned vector's length is length(distances(object)). If object is a matrix, return's length is nrow(object).

Note

Users will not normally call this function. It is called internally by nLL() and effectiveDistance().

See Also

integrateNumeric(); integrateNegexpLines(); integrateOneStepLines()

Examples

# Fake distance function object w/ minimum inputs for integration
d <- rep(1,4) %m%. # Only units needed, not values
obs <- factor(rep(c("obs1", "obs2"), 2))
beta <- c(3.5, -0.5)
w.hi <- 125
w.lo <- 20
ml <- list(
    mf = model.frame(d ~ obs)
  , par = beta 
  , likelihood = "halfnorm"
  , w.lo = w.lo %#% "m"
  , w.hi = w.hi %#% "m"
  , expansions = 0
)
class(ml) <- "dfunc"
integrateHalfnormLines(ml)

# Check: Integral of exp(-x^2/(2*s^2)) from 0 to w.hi-w.lo
b <- exp(c(beta[1], beta[1] + beta[2]))
intgral <- (pnorm(w.hi, mean=w.lo, sd = b) - 0.5) * sqrt(2*pi)*b
intgral

Description

Compute integral of the half-normal distance function for point surveys.

Usage

integrateHalfnormPoints(
  object,
  newdata = NULL,
  w.lo = NULL,
  w.hi = NULL,
  Units = NULL
)

Arguments