Package 'mrds'

Description

This package implements mark-recapture distance sampling methods as described in D.L. Borchers, W. Zucchini and Fewster, R.M. (1988), "Mark-recapture models for line transect surveys", Biometrics 54: 1207-1220. and Laake, J.L. (1999) "Distance sampling with independent observers: Reducing bias from heterogeneity by weakening the conditional independence assumption." in Amstrup, G.W., Garner, S.C., Laake, J.L., Manly, B.F.J., McDonald, L.L. and Robertson, D.G. (eds) "Marine mammal survey and assessment methods", Balkema, Rotterdam: 137-148 and Borchers, D.L., Laake, J.L., Southwell, C. and Paxton, C.L.G. "Accommodating unmodelled heterogeneity in double-observer distance sampling surveys". 2006. Biometrics 62:372-378.)

Details

Examples of distance sampling analyses are available at http://examples.distancesampling.org/.

For help with distance sampling and this package, there is a Google Group https://groups.google.com/forum/#!forum/distance-sampling.

Author(s)

Jeff Laake <[email protected]>, David Borchers <[email protected]>, Len Thomas <[email protected]>, David L. Miller <[email protected]>, Jon Bishop <[email protected]>, Felix Petersma <[email protected]>

Description

Add a line or lines to a plot of the detection function which correspond to a a given covariate combination. These can be particularly useful when there is a small number of factor levels or if quantiles of a continuous covariate are specified.

Usage

add.df.covar.line(ddf, data, ndist = 250, pdf = FALSE, breaks = "Sturges", ...)

add_df_covar_line(ddf, data, ndist = 250, pdf = FALSE, breaks = "Sturges", ...)

Arguments

Details

All covariates must be specified in data. Plots can become quite busy when this approach is used. It may be useful to fix some covariates at their median level and plot set values of a covariate of interest. For example setting weather (e.g., Beaufort) to its median and plotting levels of observer, then creating a second plot for a fixed observer with levels of weather.

Arguments to lines are supplied in ... and aesthetics like line type (lty), line width (lwd) and colour (col) are recycled. By default lty is used to distinguish between the lines. It may be useful to add a legend to the plot (lines are plotted in the order of data).

Value

invisibly, the values of detectability over the truncation range.

Author(s)

David L Miller

Examples

## Not run: 
# fit an example model
data(book.tee.data)
egdata <- book.tee.data$book.tee.dataframe
result <- ddf(dsmodel = ~mcds(key = "hn", formula = ~sex),
              data = egdata[egdata$observer==1, ], method = "ds",
              meta.data = list(width = 4))

# make a base plot, showpoints=FALSE makes the plot less busy
plot(result, showpoints=FALSE)

# add lines for sex one at a time
add.df.covar.line(result, data.frame(sex=0), lty=2)
add.df.covar.line(result, data.frame(sex=1), lty=3)

# add a legend
legend(3, 1, c("Average", "sex==0", "sex==1"), lty=1:3)

# alternatively we can add both at once
# fixing line type and varying colour
plot(result, showpoints=FALSE)
add.df.covar.line(result, data.frame(sex=c(0,1)), lty=1,
                  col=c("red", "green"))
# add a legend
legend(3, 1, c("Average", "sex==0", "sex==1"), lty=1,
       col=c("black", "red", "green"))

## End(Not run)

Description

'adj.check.order' checks that the Cosine, Hermite or simple polynomials are of the correct order.

Usage

adj.check.order(adj.series, adj.order, key)

Arguments

Details

Only even functions are allowed as adjustment terms, per p.47 of Buckland et al (2001). If incorrect terms are supplied then an error is throw via stop.

Value

Nothing! Just calls stop if something goes wrong.

Author(s)

David Miller

References

S.T.Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake. 1993. Robust Models. In: Distance Sampling, eds. S.T.Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake. Chapman & Hall.

See Also

adjfct.cos, adjfct.poly, adjfct.herm, detfct, mcds, cds

Description

For internal use only – not to be called by 'mrds' or 'Distance' users directly.

Usage

adj.cos(distance, scaling, adj.order)

Arguments

Value

scalar or vector containing the cosine adjustment term for every value in distance argument

Author(s)

Felix Petersma

Description

For internal use only – not to be called by 'mrds' or 'Distance' users directly.

Usage

adj.herm(distance, scaling, adj.order)

Arguments

Value

scalar or vector containing the Hermite adjustment term for every value in distance argument

Author(s)

Felix Petersma

Description

For internal use only – not to be called by 'mrds' or 'Distance' users directly.

Usage

adj.poly(distance, scaling, adj.order)

Arguments

Value

scalar or vector containing the polynomial adjustment term for every value in distance argument

Author(s)

Felix Petersma

Description

For internal use only – not to be called by 'mrds' or 'Distance' users directly.

Usage

adj.series.grad.cos(
  distance,
  scaling = 1,
  adj.order,
  adj.parm = NULL,
  adj.exp = FALSE
)

Arguments

Value

scalar or vector containing the gradient of the cosine adjustment series for every value in distance argument

Author(s)

Felix Petersma

Description

For internal use only – not to be called by 'mrds' or 'Distance' users directly.

Usage

adj.series.grad.herm(
  distance,
  scaling = 1,
  adj.order,
  adj.parm = NULL,
  adj.exp = FALSE
)

Arguments

Value

scalar or vector containing the gradient of the Hermite adjustment series for every value in distance argument

Author(s)

Felix Petersma

Description

For internal use only – not to be called by 'mrds' or 'Distance' users directly.

Usage

adj.series.grad.poly(
  distance,
  scaling = 1,
  adj.order,
  adj.parm = NULL,
  adj.exp = FALSE
)

Arguments

Value

scalar or vector containing the gradient of the polynomial adjustment series for every value in distance argument

Author(s)

Felix Petersma

Description

Extract the AIC from a fitted detection function.

Usage

## S3 method for class 'ddf'
AIC(object, ..., k = 2)

Arguments

Author(s)

David L Miller

Description

Get the apex for a gamma detection function

Usage

apex.gamma(ddfobj)

Arguments

Value

the distance at which the gamma peaks

Author(s)

Jeff Laake

Description

Assigns default values for argument in list x from argument=value pairs in ... if x$argument doesn't already exist

Usage

assign.default.values(x, ...)

Arguments

Value

x - list with filled values

Author(s)

Jeff Laake

Description

For models with covariates the detection probability for each observation can vary. This function computes an average value for a set of distances to plot an average line to graphically represent the fitted model in plots that compare histograms and the scatter of individual estimated detection probabilities. Averages are calculated over the observed covariate combinations.

Usage

average.line(finebr, obs, model)

Arguments

Value

list with 2 elements

Note

Internal function called from plot functions for ddf objects

Author(s)

Jeff Laake

Description

For models with covariates the detection probability for each observation can vary. This function computes an average value for a set of distances to plot an average line to graphically represent the fitted model in plots that compare histograms and the scatter of individual estimated detection probabilities.

Usage

average.line.cond(finebr, obs, model)

Arguments

Value

list with 2 elements:

Note

Internal function called from plot functions for ddf objects

Author(s)

Jeff Laake

Description

Double platform data collected in a line transect survey of golf tees by 2 observers at St. Andrews. Field sex was actually colour of the golf tee: 0 - green; 1 - yellow. Exposure was either low (0) or high(1) depending on height of tee above the ground. size was the number of tees in an observed cluster.

Format

A list of 4 dataframes, with the list elements named: book.tee.dataframe, book.tee.region, book.tee.samples and book.tee.obs.

book.tee.dataframe is the distance sampling data dataframe. Used in the call to fit the detection function in ddf. Contains the following columns:

object

numeric object id

observer

factor representing observer 1 or 2

detected

numeric 1 if the animal was detected 0 otherwise

distance

numeric value for the distance the animal was detected

size

numeric value for the group size

sex

numeric value for sex of animal

exposure

numeric value for exposure level 0 or 1

book.tee.region: is the region table dataframe. Used to supply the strata areas to the dht function. Contains the following columns:

Region.Label

factor giving the strata labels

Area

numeric value giving the strata areas

book.tee.samples is the samples table dataframe to match the transect ids to the region ids and supply the effort. Used in the dht function. Contains the following columns:

Sample.Label

numeric giving the sample / transect labels

Region.Label

factor giving the strata labels

Effort

numeric value giving the sample / transect lengths

book.tee.obs is the observations table dataframe to match the object ids in the distance data to the transect labels. Used in the dht function. Contains the following columns:

object

numeric value object id

Region.Label

factor giving the strata labels

Sample.Label

numeric giving the sample / transect labels

Description

Find se of average p and N

Usage

calc.se.Np(model, avgp, n, average.p)

Arguments

Author(s)

David L. Miller

Description

Computes cdf values of observed distances from fitted distribution. For a set of observed x it returns the integral of f(x) for the range= (inner, x), where inner is the innermost distance which is observable (either 0 or left if left truncated). In terms of g(x) this is the integral of g(x) over range divided by the integral of g(x) over the entire range of the data (inner, W).

Usage

cdf.ds(model, newdata = NULL)

Arguments

Value

vector of cdf values for each observation

Note

This is an internal function that is not intended to be invoked directly. It is called by qqplot.ddf to compute values for Kolmogorov-Smirnov and Cramer-von Mises tests and the Q-Q plot.

Author(s)

Jeff Laake

See Also

qqplot.ddf

Description

Creates model formula list for conventional distance sampling using values supplied in call to ddf

Usage

cds(
  key = NULL,
  adj.series = NULL,
  adj.order = NULL,
  adj.scale = "width",
  adj.exp = FALSE,
  formula = ~1,
  shape.formula = ~1
)

Arguments

Value

A formula list used to define the detection function model

Author(s)

Jeff Laake; Dave Miller

Description

Simple internal function to check that the optimisation didn't hit bounds. Based on code that used to live in detfct.fit.opt.

Usage

check.bounds(lt, lowerbounds, upperbounds, ddfobj, showit, setlower, setupper)

Arguments

Value

TRUE if parameters are close to the bound, else FALSE

Author(s)

Dave Miller; Jeff Laake

Description

Check that a fitted detection function is monotone non-increasing.

Usage

check.mono(
  df,
  strict = TRUE,
  n.pts = 100,
  tolerance = 1e-08,
  plot = FALSE,
  max.plots = 6
)

Arguments

Details

Evaluates a series of points over the range of the detection function (left to right truncation) then determines:

1. If the detection function is always less than or equal to its value at the left truncation (g(x)<=g(left), or usually g(x)<=g(0)). 2. (Optionally) The detection function is always monotone decreasing (g(x[i])<=g(x[i-1])). This check is only performed when strict=TRUE (the default). 3. The detection function is never less than 0 (g(x)>=0). 4. The detection function is never greater than 1 (g(x)<=1).

For models with covariates in the scale parameter of the detection function is evaluated at all observed covariate combinations.

Currently covariates in the shape parameter are not supported.

Value

TRUE if the detection function is monotone, FALSE if it's not. warnings are issued to warn the user that the function is non-monotonic.

Author(s)

David L. Miller, Felix Petersma

Description

Extract coefficients and provide a summary of parameters and estimates from the output of ddf model objects.

Usage

## S3 method for class 'ds'
coef(object,...)
       ## S3 method for class 'io'
coef(object,...)
       ## S3 method for class 'io.fi'
coef(object,...)
       ## S3 method for class 'trial'
coef(object,...)
       ## S3 method for class 'trial.fi'
coef(object,...)
       ## S3 method for class 'rem'
coef(object,...)
       ## S3 method for class 'rem.fi'
coef(object,...)

Arguments

Value

For coef.ds List of data frames for coefficients (scale and exponent (if hazard))

For all others Data frame containing each coefficient and standard error

Note

These functions are called by the generic function coef for any ddf model object. It can be called directly by the user, but it is typically safest to use coef which calls the appropriate function based on the type of model.

Author(s)

Jeff Laake

Description

Compute individual components of Horvitz-Thompson abundance estimate in covered region for a particular subset of the data depending on value of group = TRUE (do group abundance); FALSE(do individual abundance)

Usage

compute.Nht(pdot, group = TRUE, size = NULL)

Arguments

Value

vector of H-T components for abundance estimate

Note

Internal function called by covered.region.dht

Author(s)

Jeff Laake

Description

Computes H-T abundance within covered region by sample.

Usage

covered.region.dht(obs, samples, group)

Arguments

Value

Nhat.by.sample - dataframe of abundance by sample

Note

Internal function called by dht and related functions

Author(s)

Jeff Laake

Description

This is an internal routine and shouldn't be necessary in normal analyses.

Usage

create.bins(data, cutpoints)

Arguments

Value

argument 'data' with two extra columns 'distbegin' and 'distend'.

Author(s)

David L. Miller

Description

create.command.file

Usage

create.command.file(dsmodel = call(), data, method, meta.data, control)

Arguments

Author(s)

Jonah McArthur

Description

Creates a model.frame for distance detection function fitting. It includes some pre-specified and computed variables with those included in the model specified by user (formula)

Usage

create.model.frame(xmat, scale.formula, meta.data, shape.formula = NULL)

Arguments

Details

The following fields are always included: detected, observer, binned, and optionally distance (unless null), timesdetected (if present in data). If the distance data were binned, include distbegin and distend point fields. If the integration width varies also include int.begin and int.end and include an offset field for an iterative glm, if used. Beyond these fields only fields used in the model formula are included.

Value

model frame for analysis

Note

Internal function and not called by user

Author(s)

Jeff Laake

Description

Creates samples and obs dataframes used to compute abundance and its variance based on a structure of geographic regions and samples within each region. The intent is to generalize this routine to work with other sampling structures.

Usage

create.varstructure(model, region, sample, obs, dht.se)

Arguments

Details

The function performs the following tasks: 1)tests to make sure that region labels are unique, 2) merges sample and region tables into a samples table and issue a warning if not all samples were used, 3) if some regions have no samples or if some values of Area were not valid areas given then issue error and stop, then an error is given and the code stops, 4) creates a unique region/sample label in samples and in obs, 5) merges observations with sample and issues a warning if not all observations were used, 6) sorts regions by its label and merges the values with the predictions from the fitted model based on the object number and limits it to the data that is appropriate for the fitted detection function.

Value

List with 2 elements:

Note

Internal function called by dht

Author(s)

Jeff Laake

Description

Generic function for fitting detection functions for distance sampling with single and double observer configurations. Independent observer, trial and dependent observer (removal) configurations are included. This is a generic function which does little other than to validate the calling arguments and methods and then calls the appropriate method specific function to do the analysis.

Usage

ddf(
  dsmodel = call(),
  mrmodel = call(),
  data,
  method = "ds",
  meta.data = list(),
  control = list(),
  call = NULL
)

Arguments

Details

The fitting code has certain expectations about data. It should be a dataframe with at least the following fields named and defined as follows:

If the data are for clustered objects, the dataframe should also contain a field named size that gives the observed number in the cluster. If the data are for a double observer survey, then there are two records for each observation and each should have the same object number. The code assumes the observations are listed in the same order for each observer such that if the data are subsetted by observer there will be the same number of records in each and each subset will be in the same object order. In addition to these predefined and pre-named fields, the dataframe can have any number and type of fields that are used as covariates in the dsmodel and mrmodel. At present, discrepancies between observations in distance, size and any user-specified covariates cannot be assimilated into the uncertainty of the estimate. The code presumes the values for those fields are the same for both records (observer=1 and observer=2) and it uses the value from observer 1. Thus it makes sense to make the values the same for both records in each pair even when both detect the object or when observer 1 doesn't detect the object the data would have to be taken from observer 2 and would not be consistent.

Five different fitting methods are currently available and these in turn define whether dsmodel and mrmodel need to be defined.

Methods with the suffix ".fi" use the assumption of full independence and do not use the distance sampling portion of the likelihood which is why a dsmodel is not needed. An mrmodel is only needed for double observer surveys and thus is not needed for method ds.

The dsmodel specifies the detection function g(y) for the distance sampling data and the models restrict g(0)=1. For single observer data g(y) is the detection function for the single observer and if it is a double observer survey it is the relative detection function (assuming g(0)=1) of both observers as a team (the unique observations from both observers). In double observer surveys, the detection function is p(y)=p(0)g(y) such that p(0)<1. The detection function g(y) is specified by dsmodel and p(0) estimated from the conditional detection functions (see mrmodel below). The value of dsmodel is specified using a hybrid formula/function notation. The model definition is prefixed with a ~ and the remainder is a function definition with specified arguments. At present there are two different functions, cds and mcds, for conventional distance sampling and multi-covariate distance sampling. Both functions have the same required arguments (key,formula). The first specifies the key function this can be half-normal ("hn"), hazard-rate ("hr"), gamma ("gamma") or uniform ("unif"). The argument formula specifies the formula for the log of the scale parameter of the key function (e.g., the equivalent of the standard deviation in the half-normal). The variable distance should not be included in the formula because the scale is for distance. See Marques, F.F.C. and S.T. Buckland (2004) for more details on the representation of the scale formula. For the hazard rate and gamma functions, an additional shape.formula can be specified for the model of the shape parameter. The default will be ~1. Adjustment terms can be specified by setting adj.series which can have the values: "none", "cos" (cosine), "poly" (polynomials), and "herm" (Hermite polynomials). One must also specify a vector of orders for the adjustment terms (adj.order) and a scaling (adj.scale) which may be "width" or "scale" (for scaling by the scale parameter). Note that the uniform key can only be used with adjustments (usually cosine adjustments for a Fourier-type analysis).

The mrmodel specifies the form of the conditional detection functions (i.e.,probability it is seen by observer j given it was seen by observer 3-j) for each observer (j=1,2) in a double observer survey. The value is specified using the same mix of formula/function notation but in this case the functions are glm and gam. The arguments for the functions are formula and link. At present, only glm is allowed and it is restricted to link=logit. Thus, currently the only form for the conditional detection functions is logistic as expressed in eq 6.32 of Laake and Borchers (2004). In contrast to dsmodel, the argument formula will typically include distance and all other covariates that affect detection probability. For example, mrmodel=~glm(formula=~distance+size+sex) constructs a conditional detection function based on the logistic form with additive factors, distance, size, and sex. As another example, mrmodel=~glm(formula=~distance*size+sex) constructs the same model with an added interaction between distance and size.

The argument meta.data is a list that enables various options about the data to be set. These options include:

point

if TRUE the data are from point counts and FALSE (default) implies line transect data

width

distance specifying half-width of the transect

left

distance specifying inner truncation value

binned

TRUE or FALSE to specify whether distances should be binned for analysis

breaks

if binned=TRUE, this is a required sequence of break points that are used for plotting/gof. They should match distbegin, distend values if bins are fixed

int.range

an integration range for detection probability; either a vector of 2 or matrix with 2 columns

mono

constrain the detection function to be weakly monotonically decreasing (only applicable when there are no covariates in the detection function)

mono.strict

when TRUE constrain the detection function to be strictly monotonically decreasing (again, only applicable when there are no covariates in the detection function)

Using meta.data=list(int.range=c(1,10)) is the same as meta.data=list(left=1,width=10). If meta.data=list(binned=TRUE) is used, the dataframe needs to contain the fields distbegin and distend for each observation which specify the left and right hand end points of the distance interval containing the observation. This is a general data structure that allows the intervals to change rather than being fixed as in the standard distance analysis tools. Typically, if the intervals are changing so is the integration range. For example, assume that distance bins are generated using fixed angular measurements from an aircraft in which the altitude is varying. Because all analyses are truncated (i.e., the last interval does not go to infinity), the transect width (and the left truncation point if there is a blindspot below the aircraft) can potentially change for each observation. The argument int.range can also be entered as a matrix with 2 columns (left and width) and a row for each observation.

The argument control is a list that enables various analysis options to be set. It is not necessary to set any of these for most analyses. They were provided so the user can optionally see intermediate fitting output and to control fitting if the algorithm doesn't converge which happens infrequently. The list values include:

showit

Integer (0-3, default 0) controls the (increasing)amount of information printed during fitting. 0 - none, >=1 - information about refitting and bound changes is printed, >=2 - information about adjustment term fitting is printed, ==3 -per-iteration parameter estimates and log-likelihood printed.

estimate

if FALSE fits model but doesn't estimate predicted probabilities

refit

if TRUE the algorithm will attempt multiple optimizations at different starting values if it doesn't converge

nrefits

number of refitting attempts

initial

a named list of starting values for the dsmodel parameters (e.g. $scale, $shape, $adjustment)

lowerbounds

a vector of lowerbounds for the dsmodel parameters in the order the ds parameters will appear in the par element of the ddf object, i.e. fit.ddf$par where fit.ddf is a fitted ddf model.

upperbounds

a vector of upperbounds for the dsmodel parameters in the order the ds parameters will appear in the par element of the ddf object, i.e. fit.ddf$par where fit.ddf is a fitted ddf model.

limit

if TRUE restrict analysis to observations with detected=1

debug

if TRUE, if fitting fails, return an object with fitting information

nofit

if TRUE don't fit a model, but use the starting values and generate an object based on those values

optimx.method

one (or a vector of) string(s) giving the optimisation method to use. If more than one is supplied, the results from one are used as the starting values for the next. See optimx

optimx.maxit

maximum number of iterations to use in the optimisation.

mono.random.start

By default when monotonicity constraints are enforced, a grid of starting values are tested. Instead random starting values can be used (uniformly distributed between the upper and lower bounds). Set TRUE for random start, FALSE (default) uses the grid method

mono.method

The optimiser method to be used when (strict) monotonicity is enforced. Can be either slsqp or solnp. Default slsqp.

mono.startvals

Controls if the mono.optimiser should find better starting values by first fitting a key function without adjustments, and then use those start values for the key function parameters when fitting the key + adjustment series detection function. Defaults to FALSE

mono.outer.iter

Number of outer iterations to be used by solnp when fitting a monotonic model and solnp is selected. Default 200.

silent

silences warnings within ds fitting method (helpful for running many times without generating many warning/error messages).

optimizer

By default this is set to 'both' for single observer analyses and 'R' for double observer analyses. For single observer analyses where optimizer = 'both', the R optimizer will be used and if present the MCDS optimizer will also be used. The result with the best likelihood value will be selected. To run only a specified optimizer set this value to either 'R' or 'MCDS'. The MCDS optimizer cannot currently be used for detection function fitting with double observer analyses. See mcds_dot_exe for more information.

winebin

Location of the wine binary used to run MCDS.exe. See mcds_dot_exe for more information.

Examples of distance sampling analyses are available at https://examples.distancesampling.org/.

Hints and tips on fitting (particularly optimisation issues) are on the mrds_opt manual page.

Value

model object of class=(method, "ddf")

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

Marques, F.F.C. and S.T. Buckland. 2004. Covariate models for the detection function. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

ddf.ds, ddf.io, ddf.io.fi, ddf.trial, ddf.trial.fi, ddf.rem, ddf.rem.fi, mrds_opt

Examples

# load data
data(book.tee.data)
region <- book.tee.data$book.tee.region
egdata <- book.tee.data$book.tee.dataframe
samples <- book.tee.data$book.tee.samples
obs <- book.tee.data$book.tee.obs

# fit a half-normal detection function
result <- ddf(dsmodel=~mcds(key="hn", formula=~1), data=egdata, method="ds",
              meta.data=list(width=4))

# fit an independent observer model with full independence
result.io.fi <- ddf(mrmodel=~glm(~distance), data=egdata, method="io.fi",
                    meta.data=list(width = 4))

# fit an independent observer model with point independence
result.io <- ddf(dsmodel=~cds(key = "hn"), mrmodel=~glm(~distance),
                 data=egdata, method="io", meta.data=list(width=4))
## Not run: 

# simulated single observer point count data (see ?ptdata.single)
data(ptdata.single)
ptdata.single$distbegin <- (as.numeric(cut(ptdata.single$distance,
                            10*(0:10)))-1)*10
ptdata.single$distend <- (as.numeric(cut(ptdata.single$distance,
                          10*(0:10))))*10
model <- ddf(data=ptdata.single, dsmodel=~cds(key="hn"),
             meta.data=list(point=TRUE,binned=TRUE,breaks=10*(0:10)))

summary(model)

plot(model,main="Single observer binned point data - half normal")

model <- ddf(data=ptdata.single, dsmodel=~cds(key="hr"),
             meta.data=list(point=TRUE, binned=TRUE, breaks=10*(0:10)))

summary(model)

plot(model, main="Single observer binned point data - hazard rate")

dev.new()

# simulated double observer point count data (see ?ptdata.dual)
# setup data
data(ptdata.dual)
ptdata.dual$distbegin <- (as.numeric(cut(ptdata.dual$distance,
                          10*(0:10)))-1)*10
ptdata.dual$distend <- (as.numeric(cut(ptdata.dual$distance,
                        10*(0:10))))*10

model <- ddf(method="io", data=ptdata.dual, dsmodel=~cds(key="hn"),
             mrmodel=~glm(formula=~distance*observer),
             meta.data=list(point=TRUE, binned=TRUE, breaks=10*(0:10)))

summary(model)

plot(model, main="Dual observer binned point data", new=FALSE, pages=1)

model <- ddf(method="io", data=ptdata.dual,
             dsmodel=~cds(key="unif", adj.series="cos", adj.order=1),
             mrmodel=~glm(formula=~distance*observer),
             meta.data=list(point=TRUE, binned=TRUE, breaks=10*(0:10)))

summary(model)

par(mfrow=c(2,3))
plot(model,main="Dual observer binned point data",new=FALSE)


## End(Not run)

Description

Fits a conventional distance sampling (CDS) (likelihood eq 6.6 in Laake and Borchers 2004) or multi-covariate distance sampling (MCDS)(likelihood eq 6.14 in Laake and Borchers 2004) model for the detection function of observed distance data. It only uses key functions and does not incorporate adjustment functions as in CDS/MCDS analysis engines in DISTANCE (Marques and Buckland 2004). Distance can be grouped (binned), ungrouped (unbinned) or mixture of the two. This function is not called directly by the user and is called from ddf,ddf.io, or ddf.trial.

Usage

## S3 method for class 'ds'
ddf(
  dsmodel,
  mrmodel = NULL,
  data,
  method = "ds",
  meta.data = list(),
  control = list(),
  call
)

Arguments

Details

For a complete description of each of the calling arguments, see ddf. The argument model in this function is the same as dsmodel in ddf. The argument dataname is the name of the dataframe specified by the argument data in ddf. The arguments control,meta.data,and method are defined the same as in ddf.

Value

result: a ds model object

Note

If mixture of binned and unbinned distance, width must be set to be >= largest interval endpoint; this could be changed with a more complicated analysis; likewise, if all binned and bins overlap, the above must also hold; if bins don't overlap, width must be one of the interval endpoints; same holds for left truncation Although the mixture analysis works in principle it has not been tested via simulation.

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R. Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

Marques, F.F.C. and S.T. Buckland. 2004. Covariate models for the detection function. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R. Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

flnl, summary.ds, coef.ds, plot.ds,gof.ds

Examples

# ddf.ds is called when ddf is called with method="ds"

data(book.tee.data)
region <- book.tee.data$book.tee.region
egdata <- book.tee.data$book.tee.dataframe
samples <- book.tee.data$book.tee.samples
obs <- book.tee.data$book.tee.obs
result <- ddf(dsmodel = ~mcds(key = "hn", formula = ~1),
              data = egdata[egdata$observer==1, ], method = "ds",
              meta.data = list(width = 4))
summary(result,se=TRUE)
plot(result,main="cds - observer 1")
print(dht(result,region,samples,obs,options=list(varflag=0,group=TRUE),
          se=TRUE))
print(ddf.gof(result))

Description

Generic function that computes chi-square goodness of fit test for detection function models with binned data and Cramer-von Mises and Kolmogorov-Smirnov (if ks=TRUE)tests for exact distance data. By default a Q-Q plot is generated for exact data (and can be suppressed using the qq=FALSE argument).

Usage

ddf.gof(
  model,
  breaks = NULL,
  nc = NULL,
  qq = TRUE,
  nboot = 100,
  ks = FALSE,
  ...
)

Arguments

Details

Formal goodness of fit testing for detection function models using Kolmogorov-Smirnov and Cramer-von Mises tests. Both tests are based on looking at the quantile-quantile plot produced by qqplot.ddf and deviations from the line x=y.

The Kolmogorov-Smirnov test asks the question "what's the largest vertical distance between a point and the y=x line?" It uses this distance as a statistic to test the null hypothesis that the samples (EDF and CDF in our case) are from the same distribution (and hence our model fits well). If the deviation between the y=x line and the points is too large we reject the null hypothesis and say the model doesn't have a good fit.

Rather than looking at the single biggest difference between the y=x line and the points in the Q-Q plot, we might prefer to think about all the differences between line and points, since there may be many smaller differences that we want to take into account rather than looking for one large deviation. Its null hypothesis is the same, but the statistic it uses is the sum of the deviations from each of the point to the line. Note that a bootstrap procedure is required for the Kolmogorov-Smirnov test to ensure that the p-values from the procedure are correct as the we are comparing the cumulative distribution function (CDF) and empirical distribution function (EDF) and we have estimated the parameters of the detection function. The nboot parameter controls the number of bootstraps to use. Set to 0 to avoid computing bootstraps (much faster but with no Kolmogorov-Smirnov results, of course).

One can change the precision of printed values by using the print.ddf.gof method's digits argument.

Value

List of class ddf.gof containing

Author(s)

Jeff Laake

See Also

qqplot.ddf

Description

Mark-Recapture Distance Sampling (MRDS) Analysis of Independent Observer Configuration and Point Independence

Usage

## S3 method for class 'io'
ddf(
  dsmodel,
  mrmodel,
  data,
  method = NULL,
  meta.data = list(),
  control = list(),
  call = ""
)

Arguments

Details

MRDS analysis based on point independence involves two separate and independent analyses of the mark-recapture data and the distance sampling data. For the independent observer configuration, the mark-recapture data are analysed with a call to ddf.io.fi (see likelihood eq 6.8 and 6.16 in Laake and Borchers 2004) to fit conditional distance sampling detection functions to estimate p(0), detection probability at distance zero for the independent observer team based on independence at zero (eq 6.22 in Laake and Borchers 2004). Independently, the distance data, the union of the observations from the independent observers, are used to fit a conventional distance sampling (CDS) (likelihood eq 6.6) or multi-covariate distance sampling (MCDS) (likelihood eq 6.14) model for the detection function, g(y), such that g(0)=1. The detection function for the observer team is then created as p(y)=p(0)*g(y) (eq 6.28 of Laake and Borchers 2004) from which predictions are made. ddf.io is not called directly by the user and is called from ddf with method="io".

For a complete description of each of the calling arguments, see ddf. The argument dataname is the name of the dataframe specified by the argument data in ddf. The arguments dsmodel, mrmodel, control and meta.data are defined the same as in ddf.

Value

result: an io model object which is composed of io.fi and ds model objects

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

ddf.io.fi, ddf.ds,summary.io,coef.io,plot.io, gof.io

Description

Mark-Recapture Analysis of Independent Observer Configuration with Full Independence

Usage

## S3 method for class 'io.fi'
ddf(
  dsmodel = NULL,
  mrmodel,
  data,
  method,
  meta.data = list(),
  control = list(),
  call = ""
)

Arguments

Details

The mark-recapture data derived from an independent observer distance sampling survey can be used to derive conditional detection functions (p_j(y)) for both observers (j=1,2). They are conditional detection functions because detection probability for observer j is based on seeing or not seeing observations made by observer 3-j. Thus, p_1(y) is estimated by p_1|2(y).

If detections by the observers are independent (full independence) then p_1(y)=p_1|2(y),p_2(y)=p_2|1(y) and for the union, full independence means that p(y)=p_1(y) + p_2(y) - p_1(y)*p_2(y) for each distance y. In fitting the detection functions the likelihood given by eq 6.8 and 6.16 in Laake and Borchers (2004) is used. That analysis does not require the usual distance sampling assumption that perpendicular distances are uniformly distributed based on line placement that is random relative to animal distribution. However, that assumption is used in computing predicted detection probability which is averaged based on a uniform distribution (see eq 6.11 of Laake and Borchers 2004).

For a complete description of each of the calling arguments, see ddf. The argument model in this function is the same as mrmodel in ddf. The argument dataname is the name of the dataframe specified by the argument data in ddf. The arguments control,meta.data,and method are defined the same as in ddf.

Value

result: an io.fi model object

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

ddf.io,summary.io.fi,coef.io.fi, plot.io.fi,gof.io.fi,io.glm

Description

Mark-Recapture Distance Sampling (MRDS) Analysis of Removal Observer Configuration and Point Independence

Usage

## S3 method for class 'rem'
ddf(
  dsmodel,
  mrmodel,
  data,
  method = NULL,
  meta.data = list(),
  control = list(),
  call = ""
)

Arguments

Details

MRDS analysis based on point independence involves two separate and independent analyses of the mark-recapture data and the distance sampling data. For the removal observer configuration, the mark-recapture data are analysed with a call to ddf.rem.fi (see Laake and Borchers 2004) to fit conditional distance sampling detection functions to estimate p(0), detection probability at distance zero for the primary observer based on independence at zero (eq 6.22 in Laake and Borchers 2004). Independently, the distance data, the observations from the primary observer, are used to fit a conventional distance sampling (CDS) (likelihood eq 6.6) or multi-covariate distance sampling (MCDS) (likelihood eq 6.14) model for the detection function, g(y), such that g(0)=1. The detection function for the primary observer is then created as p(y)=p(0)*g(y) (eq 6.28 of Laake and Borchers 2004) from which predictions are made. ddf.rem is not called directly by the user and is called from ddf with method="rem".

For a complete description of each of the calling arguments, see ddf. The argument data is the dataframe specified by the argument data in ddf. The arguments dsmodel, mrmodel, control and meta.data are defined the same as in ddf.

Value

result: an rem model object which is composed of rem.fi and ds model objects

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

ddf.rem.fi, ddf.ds

Description

Mark-Recapture Distance Sampling (MRDS) Analysis of Removal Observer Configuration with Full Independence

Usage

## S3 method for class 'rem.fi'
ddf(
  dsmodel = NULL,
  mrmodel,
  data,
  method,
  meta.data = list(),
  control = list(),
  call = ""
)

Arguments

Details

The mark-recapture data derived from an removal observer distance sampling survey can only derive conditional detection functions (p_j(y)) for both observers (j=1) because technically it assumes that detection probability does not vary by occasion (observer in this case). It is a conditional detection function because detection probability for observer 1 is conditional on the observations seen by either of the observers. Thus, p_1(y) is estimated by p_1|2(y).

If detections by the observers are independent (full independence) then p_1(y)=p_1|2(y) and for the union, full independence means that p(y)=p_1(y) + p_2(y) - p_1(y)*p_2(y) for each distance y. In fitting the detection functions the likelihood from Laake and Borchers (2004) are used. That analysis does not require the usual distance sampling assumption that perpendicular distances are uniformly distributed based on line placement that is random relative to animal distribution. However, that assumption is used in computing predicted detection probability which is averaged based on a uniform distribution (see eq 6.11 of Laake and Borchers 2004).

For a complete description of each of the calling arguments, see ddf. The argument model in this function is the same as mrmodel in ddf. The argument dataname is the name of the dataframe specified by the argument data in ddf. The arguments control,meta.data,and method are defined the same as in ddf.

Value

result: an rem.fi model object

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

ddf.io,rem.glm

Description

Mark-Recapture Distance Sampling (MRDS) Analysis of Trial Observer Configuration and Point Independence

Usage

## S3 method for class 'trial'
ddf(
  dsmodel,
  mrmodel,
  data,
  method = NULL,
  meta.data = list(),
  control = list(),
  call = ""
)

Arguments

Details

MRDS analysis based on point independence involves two separate and independent analyses of the mark-recapture data and the distance sampling data. For the trial configuration, the mark-recapture data are analysed with a call to ddf.trial.fi (see likelihood eq 6.12 and 6.17 in Laake and Borchers 2004) to fit a conditional distance sampling detection function for observer 1 based on trials (observations) from observer 2 to estimate p_1(0), detection probability at distance zero for observer 1. Independently, the distance data from observer 1 are used to fit a conventional distance sampling (CDS) (likelihood eq 6.6) or multi-covariate distance sampling (MCDS) (likelihood eq 6.14) model for the detection function, g(y), such that g(0)=1. The detection function for observer 1 is then created as p_1(y)=p_1(0)*g(y) (eq 6.28 of Laake and Borchers 2004) from which predictions are made. ddf.trial is not called directly by the user and is called from ddf with method="trial".

For a complete description of each of the calling arguments, see ddf. The argument dataname is the name of the dataframe specified by the argument data in ddf. The arguments dsmodel, mrmodel, control and meta.data are defined the same as in ddf.

Value

result: a trial model object which is composed of trial.fi and ds model objects

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

ddf.trial.fi, ddf.ds, summary.trial, coef.trial, plot.trial, gof.trial

Description

Mark-Recapture Analysis of Trial Observer Configuration with Full Independence

Usage

## S3 method for class 'trial.fi'
ddf(
  dsmodel = NULL,
  mrmodel,
  data,
  method,
  meta.data = list(),
  control = list(),
  call = ""
)

Arguments

Details

The mark-recapture data derived from a trial observer distance sampling survey can be used to derive a conditional detection function (p_1(y)) for observer 1 based on trials (observations) from observer 2. It is a conditional detection function because detection probability for observer 1 is based on seeing or not seeing observations made by observer 2. Thus, p_1(y) is estimated by p_1|2(y). If detections by the observers are independent (full independence) then p_1(y)=p_1|2(y) for each distance y. In fitting the detection functions the likelihood given by eq 6.12 or 6.17 in Laake and Borchers (2004) is used. That analysis does not require the usual distance sampling assumption that perpendicular distances are uniformly distributed based on line placement that is random relative to animal distribution. However, that assumption is used in computing predicted detection probability which is averaged based on a uniform distribution (see eq 6.13 of Laake and Borchers 2004).

For a complete description of each of the calling arguments, see ddf. The argument model in this function is the same as mrmodel in ddf. The argument dataname is the name of the dataframe specified by the argument data in ddf. The arguments control,meta.data,and method are defined the same as in ddf.

Value

result: a trial.fi model object

Author(s)

Jeff Laake

References

Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete detection at distance zero. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

ddf.trial, summary.trial.fi, coef.trial.fi, plot.trial.fi, gof.trial.fi

Description

Computes delta method variance-covariance matrix of results of any generic function fct that computes a vector of estimates as a function of a set of estimated parameters par.

Usage

DeltaMethod(par, fct, vcov, delta, ...)

Arguments

Details

The delta method (aka propagation of errors is based on Taylor series approximation - see Seber's book on Estimation of Animal Abundance). It uses the first derivative of fct with respect to par. It also uses the variance-covariance matrix of the estimated parameters which is derived in estimating the parameters and is an input argument.

The first argument of fct should be par which is a vector of parameter estimates. It should return a single value (or vector) of estimate(s). The remaining arguments of fct if any can be passed to fct by including them at the end of the call to DeltaMethod as name=value pairs.

Value

a list with values

Note

This is a generic function that can be used in any setting beyond the mrds package. However this is an internal function for mrds and the user does not need to call it explicitly.

Author(s)

Jeff Laake and David L Miller

Description

Creates a series of tables for dual observer data that shows the number missed and detected for each observer within defined distance classes.

Usage

det.tables(model, nc = NULL, breaks = NULL)

Arguments

Value

list object of class "det.tables"

Author(s)

Jeff Laake

Examples

data(book.tee.data)
region <- book.tee.data$book.tee.region
egdata <- book.tee.data$book.tee.dataframe
samples <- book.tee.data$book.tee.samples
obs <- book.tee.data$book.tee.obs
xx <- ddf(mrmodel=~glm(formula=~distance*observer),
          dsmodel=~mcds(key="hn", formula=~sex),
          data=egdata, method="io", meta.data=list(width=4))
tabs <- det.tables(xx, breaks=c(0, 0.5, 1, 2, 3, 4))
par(mfrow=c(2, 2))
plot(tabs, new=FALSE, which=c(1, 2, 5, 6))

Description

Fit detection function to observed distances using the key-adjustment function approach. If adjustment functions are included it will alternate between fitting parameters of key and adjustment functions and then all parameters much like the approach in the CDS and MCDS Distance FORTRAN code. To do so it calls detfct.fit.opt which uses the R optim function which does not allow non-linear constraints so inclusion of adjustments does allow the detection function to be non-monotone.

Usage

detfct.fit(ddfobj, optim.options, bounds, misc.options)

Arguments

Value

fitted detection function model object with the following list structure

Author(s)

Dave Miller; Jeff Laake

Description

Fit detection function to observed distances using the key-adjustment function approach. If adjustment functions are included it will alternate between fitting parameters of key and adjustment functions and then all parameters much like the approach in the CDS and MCDS Distance FORTRAN code. This function is called by the driver function detfct.fit, it then calls the relevant optimisation routine, slsqp, solnp or optimx.

Usage

detfct.fit.opt(ddfobj, optim.options, bounds, misc.options, fitting = "all")

Arguments

Value

fitted detection function model object with the following list structure

Author(s)

Dave Miller; Jeff Laake; Lorenzo Milazzo; Felix Petersma

Description

Compute density and abundance estimates and variances based on Horvitz-Thompson-like estimator.

Usage

dht(
  model,
  region.table,
  sample.table,
  obs.table = NULL,
  subset = NULL,
  se = TRUE,
  options = list()
)

Arguments

Details

Density and abundance within the sampled region is computed based on a Horvitz-Thompson-like estimator for groups and individuals (if a clustered population) and this is extrapolated to the entire survey region based on any defined regional stratification. The variance is based on replicate samples within any regional stratification. For clustered populations, $E(s)$ and its standard error are also output.

Abundance is estimated with a Horvitz-Thompson-like estimator (Huggins 1989, 1991; Borchers et al 1998; Borchers and Burnham 2004). The abundance in the sampled region is simply $1/p_1 + 1/p_2 + ... + 1/p_n$ where $p_i$ is the estimated detection probability for the $i$th detection of $n$ total observations. It is not strictly a Horvitz-Thompson estimator because the $p_i$ are estimated and not known. For animals observed in tight clusters, that estimator gives the abundance of groups (group=TRUE in options) and the abundance of individuals is estimated as $s_1/p_1 + s_2/p_2 + ... + s_n/p_n$, where $s_i$ is the size (e.g., number of animals in the group) of each observation (group=FALSE in options).

Extrapolation and estimation of abundance to the entire survey region is based on either a random sampling design or a stratified random sampling design. Replicate samples (lines) are specified within regional strata region.table, if any. If there is no stratification, region.table should contain only a single record with the Area for the entire survey region. The sample.table is linked to the region.table with the Region.Label. The obs.table is linked to the sample.table with the Sample.Label and Region.Label. Abundance can be restricted to a subset (e.g., for a particular species) of the population by limiting the list the observations in obs.table to those in the desired subset. Alternatively, if Sample.Label and Region.Label are in the data.frame used to fit the model, then a subset argument can be given in place of the obs.table. To use the subset argument but include all of the observations, use subset=1==1 to avoid creating an obs.table.

In extrapolating to the entire survey region it is important that the unit measurements be consistent or converted for consistency. A conversion factor can be specified with the convert.units variable in the options list. The values of Area in region.table, must be made consistent with the units for Effort in sample.table and the units of distance in the data.frame that was analyzed. It is easiest to do if the units of Area is the square of the units of Effort and then it is only necessary to convert the units of distance to the units of Effort. For example, if Effort was entered in kilometres and Area in square kilometres and distance in metres then using options=list(convert.units=0.001) would convert metres to kilometres, density would be expressed in square kilometres which would then be consistent with units for Area. However, they can all be in different units as long as the appropriate composite value for convert.units is chosen. Abundance for a survey region can be expressed as: A*N/a where A is Area for the survey region, N is the abundance in the covered (sampled) region, and a is the area of the sampled region and is in units of Effort * distance. The sampled region a is multiplied by convert.units, so it should be chosen such that the result is in the same units of Area. For example, if Effort was entered in kilometres, Area in hectares (100m x 100m) and distance in metres, then using options=list(convert.units=10) will convert a to units of hectares (100 to convert metres to 100 metres for distance and .1 to convert km to 100m units).

The argument options is a list of variable=value pairs that set options for the analysis. All but two of these have been described above. pdelta should not need to be changed but was included for completeness. It controls the precision of the first derivative calculation for the delta method variance. If the option areas.supplied is TRUE then the covered area is assumed to be supplied in the CoveredArea column of the sample data.frame.

Value

list object of class dht with elements:

The list structure of clusters and individuals are the same:

Uncertainty

If the argument se=TRUE, standard errors for density and abundance is computed. Coefficient of variation and log-normal confidence intervals are constructed using a Satterthwaite approximation for degrees of freedom (Buckland et al. 2001 p. 90). The function dht.se computes the variance and interval estimates.

The variance has two components:

• variation due to uncertainty from estimation of the detection function parameters;

• variation in abundance due to random sample selection;

The first component (model parameter uncertainty) is computed using a delta method estimate of variance (Huggins 1989, 1991, Borchers et al. 1998) in which the first derivatives of the abundance estimator with respect to the parameters in the detection function are computed numerically (see DeltaMethod).

The second component (encounter rate variance) can be computed in one of several ways depending on the form taken for the encounter rate and the estimator used. To begin with there three possible values for varflag to calculate encounter rate:

• 0 uses a binomial variance for the number of observations (equation 13 of Borchers et al. 1998). This estimator is only useful if the sampled region is the survey region and the objects are not clustered; this situation will not occur very often;

• 1 uses the encounter rate $n/L$ (objects observed per unit transect) from Buckland et al. (2001) pg 78-79 (equation 3.78) for line transects (see also Fewster et al, 2009 estimator R2). This variance estimator is not appropriate if size or a derivative of size is used in the detection function;

• 2 is the default and uses the encounter rate estimator $\hat{N}/L$ (estimated abundance per unit transect) suggested by Innes et al (2002) and Marques & Buckland (2004).

In general if any covariates are used in the models, the default varflag=2 is preferable as the estimated abundance will take into account variability due to covariate effects. If the population is clustered the mean group size and standard error is also reported.

For options 1 and 2, it is then possible to choose one of the estimator forms given in Fewster et al (2009) for line transects: "R2", "R3", "R4", "S1", "S2", "O1", "O2" or "O3" by specifying the ervar= option (default "R2"). For points, either the "P2" or "P3" estimator can be selected (>=mrds 2.3.0 default "P2", <= mrds 2.2.9 default "P3"). See varn and Fewster et al (2009) for further details on these estimators.

dht options

Several options are available to control calculations and output:

ci.width

Confidence interval width, expressed as a decimal between 0 and 1 (default 0.95, giving a 95% CI)

pdelta

delta value for computing numerical first derivatives (Default: 0.001)

varflag

0,1,2 (see "Uncertainty") (Default: 2)

convert.units

multiplier for width to convert to units of length (Default: 1)

ervar

encounter rate variance type (see "Uncertainty" and type argument of varn). (Default: "R2" for lines and "P2" for points)

Author(s)

Jeff Laake, David L Miller

References

Borchers, D.L., S.T. Buckland, P.W. Goedhart, E.D. Clarke, and S.L. Hedley. 1998. Horvitz-Thompson estimators for double-platform line transect surveys. Biometrics 54: 1221-1237.

Borchers, D.L. and K.P. Burnham. General formulation for distance sampling pp 10-11 In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

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.

Fewster, R.M., S.T. Buckland, K.P. Burnham, D.L. Borchers, P.E. Jupp, J.L. Laake and L. Thomas. 2009. Estimating the encounter rate variance in distance sampling. Biometrics 65: 225-236.

Huggins, R.M. 1989. On the statistical analysis of capture experiments. Biometrika 76:133-140.

Huggins, R.M. 1991. Some practical aspects of a conditional likelihood approach to capture experiments. Biometrics 47: 725-732.

Innes, S., M.P. Heide-Jorgensen, J.L. Laake, K.L. Laidre, H.J. Cleator, P. Richard, and R.E.A. Stewart. 2002. Surveys of belugas and narwhals in the Canadian High Arctic in 1996. NAMMCO Scientific Publications 4: 169-190.

Marques, F.F.C. and S.T. Buckland. 2004. Covariate models for the detection function. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

See Also

print.dht dht.se

Description

Computes abundance at specified values of parameters for numerical computation of first derivative with respect to parameters in detection function. An internal function called by DeltaMethod which is invoked by dht.se

Usage

dht.deriv(par, model, obs, samples, options = list())

Arguments

Value

vector of abundance estimates at values of parameters specified in par

Note

Internal function; not intended to be called by user

Author(s)

Jeff Laake

See Also

dht, dht.se, DeltaMethod

Description

Computes standard error, cv, and log-normal confidence intervals for abundance and density within each region (if any) and for the total of all the regions. It also produces the correlation matrix for regional and total estimates.

Usage

dht.se(
  model,
  region.table,
  samples,
  obs,
  options,
  numRegions,
  estimate.table,
  Nhat.by.sample
)

Arguments

Details

The variance has two components:

• variation due to uncertainty from estimation of the detection function parameters;

• variation in abundance due to random sample selection;

The first component (model parameter uncertainty) is computed using a delta method estimate of variance (Huggins 1989, 1991, Borchers et al. 1998) in which the first derivatives of the abundance estimator with respect to the parameters in the detection function are computed numerically (see DeltaMethod).

The second component (encounter rate variance) can be computed in one of several ways depending on the form taken for the encounter rate and the estimator used. To begin with there three possible values for varflag to calculate encounter rate:

• 0 uses a binomial variance for the number of observations (equation 13 of Borchers et al. 1998). This estimator is only useful if the sampled region is the survey region and the objects are not clustered; this situation will not occur very often;

• 1 uses the encounter rate $n/L$ (objects observed per unit transect) from Buckland et al. (2001) pg 78-79 (equation 3.78) for line transects (see also Fewster et al, 2009 estimator R2). This variance estimator is not appropriate if size or a derivative of size is used in the detection function;

• 2 is the default and uses the encounter rate estimator $\hat{N}/L$ (estimated abundance per unit transect) suggested by Innes et al (2002) and Marques & Buckland (2004).

In general if any covariates are used in the models, the default varflag=2 is preferable as the estimated abundance will take into account variability due to covariate effects. If the population is clustered the mean group size and standard error is also reported.

For options 1 and 2, it is then possible to choose one of the estimator forms given in Fewster et al (2009). For line transects: "R2", "R3", "R4", "S1", "S2", "O1", "O2" or "O3" can be used by specifying the ervar= option (default "R2"). For points, either the "P2" or "P3" estimator can be selected (>=mrds 2.3.0 default "P2", <= mrds 2.2.9 default "P3"). See varn and Fewster et al (2009) for further details on these estimators.

Exceptions to the above occur if there is only one sample in a stratum. In that case it uses Poisson assumption ($Var(x)=x$) and it assumes a known variance so $z=1.96$ is used for critical value. In all other cases the degrees of freedom for the $t$-distribution assumed for the log(abundance) or log(density) is based on the Satterthwaite approximation (Buckland et al. 2001 pg 90) for the degrees of freedom (df). The df are weighted by the squared cv in combining the two sources of variation because of the assumed log-normal distribution because the components are multiplicative. For combining df for the sampling variance across regions they are weighted by the variance because it is a sum across regions.

A non-zero correlation between regional estimates can occur from using a common detection function across regions. This is reflected in the correlation matrix of the regional and total estimates which is given in the value list. It is only needed if subtotals of regional estimates are needed.

Value

List with 2 elements:

Note

This function is called by dht and it is not expected that the user will call this function directly but it is documented here for completeness and for anyone expanding the code or using this function in their own code.

Author(s)

Jeff Laake

References

see dht

See Also

dht, print.dht

Description

This function has been updated to match distpdf closely, so that it has the same flexibility. Effectively, it gives the gradient of distpdf or detfct, whichever one is specified.

Usage

distpdf.grad(
  distance,
  par.index,
  ddfobj,
  standardize = FALSE,
  width,
  point,
  left = 0,
  pdf.based = TRUE
)

Arguments

Details

Various functions used to specify key and adjustment functions for gradients of detection functions.

So far, only developed for the half-normal, hazard-rate and uniform key functions in combination with cosine, simple polynomial and Hermite polynomial adjustments. It is only called by the gradient-based solver and should not be called by the general user.

distpdf.grad will call either a half-normal, hazard-rate or uniform function with adjustment terms to fit the data better, returning the gradient of detection at that distance w.r.t. the parameters. The adjustments are either cosine, Hermite or simple polynomial.

Value

the gradient of the non-normalised pdf or detection w.r.t. to the parameter with parameter index par.index.

Author(s)

Felix Petersma

Description

Computes values of conditional and unconditional detection functions and probability density functions for for line/point data for single observer or dual observer in any of the 3 configurations (io,trial,rem).

Usage

ds.function(
  model,
  newdata = NULL,
  obs = "All",
  conditional = FALSE,
  pdf = TRUE,
  finebr
)

Arguments

Details

Placeholder – Not functional —-

Value

List containing

Author(s)

Jeff Laake

Description

For a specific set of parameter values, it computes and returns the negative log-likelihood for the distance sampling likelihood for distances that are unbinned, binned and a mixture of both. The function flnl is the function minimized using optim from within ddf.ds.

Usage

flnl(fpar, ddfobj, misc.options, fitting = "all")

Arguments

Details

Most of the computation is in flpt.lnl in which the negative log-likelihood is computed for each observation. flnl is a wrapper that optionally outputs intermediate results and sums the individual log-likelihood values.

flnl is the main routine that manipulates the parameters using getpar to handle fitting of key, adjustment or all of the parameters. It then calls flpt.lnl to do the actual computation of the likelihood. The probability density function for point counts is fr and for line transects is fx. fx=g(x)/mu (where g(x) is the detection function); whereas, f(r)=r*g(r)/mu where mu in both cases is the normalizing constant. Both functions are in source code file for link{detfct} and are called from distpdf and the integral calculations are made with integratepdf.

Value

negative log-likelihood value at the parameter values specified in fpar

Note

These are internal functions used by ddf.ds to fit distance sampling detection functions. It is not intended for the user to invoke these functions but they are documented here for completeness.

Author(s)

Jeff Laake, David L Miller

See Also

flt.var, detfct

Description

The function derives the gradients of the constraint function for all model parameters, in the following order: 1. Scale parameter (if part of key function) 2. Shape parameter (if part of key function) 3. Adjustment parameter 1 4. Adjustment parameter 2 5. Etc.

Usage

flnl.constr.grad.neg(pars, ddfobj, misc.options, fitting = "all")

Arguments

Details

The constraint function itself is formed of a specified number of non-linear constraints, which defaults to 20 and is specified through misc.options$mono.points. The constraint function checks whether the standardised detection function is 1) weakly/strictly monotonic at the points and 2) non-negative at all the points. flnl.constr.grad returns the gradients of those constraints w.r.t. all parameters of the detection function, i.e., 2 times mono.points gradients for every parameter.

This function mostly follows the same structure as flnl.constr in detfct.fit.mono.R.

Value

a matrix of gradients for all constraints (rows) w.r.t to every parameters (columns)

Description

This function derives the gradients of the negative log likelihood function, with respect to all parameters. It is based on the theory presented in Introduction to Distance Sampling (2001) and Distance Sampling: Methods and Applications (2015). It is not meant to be called by users of the mrds and Distance packages directly but rather by the gradient-based solver. This solver is use when our distance sampling model is for single-observer data coming from either line or point transect and only when the detection function contains an adjustment series but no covariates. It is implement for the following key + adjustment series combinations for the detections function: the key function can be half-normal, hazard-rate or uniform, and the adjustment series can be cosine, simple polynomial or Hermite polynomial. Data can be either binned or exact, but a combination of the two has not been implemented yet.

Usage

flnl.grad(pars, ddfobj, misc.options, fitting = "all")

Arguments

Value

The gradients of the negative log-likelihood w.r.t. the parameters

Author(s)

Felix Petersma

Description

Computes hessian to be used for variance-covariance matrix. The hessian is the outer product of the vector of first partials (see pg 62 of Buckland et al 2002).

Usage

flt.var(ddfobj, misc.options)

Arguments

Value

variance-covariance matrix of parameters in the detection function

Note

This is an internal function used by ddf.ds to fit distance sampling detection functions. It is not intended for the user to invoke this function but it is documented here for completeness.

Author(s)

Jeff Laake and David L Miller

References

Buckland et al. 2002

See Also

flnl,flpt.lnl,ddf.ds

Description

Compute value of p(0) using a logit formulation

Usage

g0(beta, z)

Arguments

Value

vector of p(0) values

Author(s)

Jeff Laake

Description

Extracts parameters of a particular type (scale, shape, adjustments or g0 (p(0))) from the vector of parameters in ddfobj. All of the parameters are kept in a single vector for optimization even though they have very different uses. assign.par parses them from the vector based on a known structure and assigns them into ddfobj. getpar extracts the requested types to be extracted from ddfobj.

Usage

getpar(ddfobj, fitting = "all", index = FALSE)

Arguments

Value

index==FALSE, vector of parameters that were requested or index==TRUE, vector of 3 indices for shape, scale, adjustment

Note

Internal functions not intended to be called by user.

Author(s)

Jeff Laake

See Also

assign.par

Description

Compute chi-square goodness-of-fit test for ds models

Usage

gof.ds(model, breaks = NULL, nc = NULL)

Arguments

Value

list with chi-square value, df and p-value

Author(s)

Jeff Laake

See Also

ddf.gof

Description

Computes the integral of distpdf with scale=1 (stdint=TRUE) or specified scale (stdint=FALSE).

Usage

gstdint(
  x,
  ddfobj,
  index = NULL,
  select = NULL,
  width,
  standardize = TRUE,
  point = FALSE,
  stdint = TRUE,
  doeachint = FALSE,
  left = left
)

Arguments

Value

vector of integral values of detection function

Note

This is an internal function that is not intended to be invoked directly.

Author(s)

Jeff Laake and David L Miller

Description

Takes bar heights (height) and cutpoints (breaks), and constructs a line-only histogram from them using the function plot() (if lineonly==FALSE) or lines() (if lineonly==TRUE).

Usage

histline(
  height,
  breaks,
  lineonly = FALSE,
  outline = FALSE,
  ylim = range(height),
  xlab = "x",
  ylab = "y",
  det.plot = FALSE,
  add = FALSE,
  ...
)

Arguments

Value

None

Author(s)

Jeff Laake and David L Miller

Description

Integrates a logistic detection function; a separate function is used because in certain cases the integral can be solved analytically and also because the scale trick used with the half-normal and hazard rate doesn't work with the logistic.

Usage

integratedetfct.logistic(x, scalemodel, width, theta1, integral.numeric, w)

Arguments

Value

vector of integral values

Author(s)

Jeff Laake

Description

Computes integral (analytically) over x from 0 to width of a logistic detection function; For reference see integral #526 in CRC Std Math Table 24th ed

Usage

integratelogistic.analytic(x, models, beta, width)

Arguments

Author(s)

Jeff Laake

Description

Computes integral of pdf of observed distances over x for each observation. The method of computation depends on argument switches set and the type of detection function.

Usage

integratepdf(
  ddfobj,
  select,
  width,
  int.range,
  standardize = TRUE,
  point = FALSE,
  left = 0,
  doeachint = FALSE
)

Arguments

Value

vector of integral values - one for each observation

Author(s)

Jeff Laake & Dave Miller

Description

Gradient of the integral of the detection function, i.e., d beta/d theta in the documentation. This gradient of the integral is the same as the integral of the gradient, thanks to Leibniz integral rule.

Usage

integratepdf.grad(
  par.index,
  ddfobj,
  int.range,
  width,
  standardize = FALSE,
  point = FALSE,
  left = 0,
  pdf.based = TRUE
)

Arguments

Details

For internal use only – not to be called by mrds or Distance users directly.

Author(s)

Felix Petersma

Description

Provides an iterative algorithm for finding the MLEs of detection (capture) probabilities for a two-occasion (double observer) mark-recapture experiment using standard algorithms GLM/GAM and an offset to compensate for conditioning on the set of observations. While the likelihood can be formulated and solved numerically, the use of GLM/GAM provides all of the available tools for fitting, predictions, plotting etc without any further development.

Usage

io.glm(
  datavec,
  fitformula,
  eps = 1e-05,
  iterlimit = 500,
  GAM = FALSE,
  gamplot = TRUE
)

Arguments

Details

Note that currently the code in this function for GAMs has been commented out until the remainder of the mrds package will work with GAMs. This is an internal function that is used as by ddf.io.fi to fit mark-recapture models with 2 occasions. The argument mrmodel is used for fitformula.

Value

list of class("ioglm","glm","lm") or class("ioglm","gam")

Author(s)

Jeff Laake, David Borchers, Charles Paxton

References

Buckland, S.T., J.M. breiwick, K.L. Cattanach, and J.L. Laake. 1993. Estimated population size of the California gray whale. Marine Mammal Science, 9:235-249.

Burnham, K.P., S.T. Buckland, J.L. Laake, D.L. Borchers, T.A. Marques, J.R.B. Bishop, and L. Thomas. 2004. Further topics in distance sampling. pp: 360-363. In: Advanced Distance Sampling, eds. S.T. Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L. Thomas. Oxford University Press.

Description

These functions are used to test whether a logistic detection function is a linear function of distance (is.linear.logistic) or is constant (varies by distance but no other covariates) is.logistic.constant). Based on these tests, the most appropriate manner for integrating the detection function with respect to distance is chosen. The integrals are needed to estimate the average detection probability for a given set of covariates.

Usage

is.linear.logistic(xmat, g0model, zdim, width)

Arguments

Details

If the logit is linear in distance then the integral can be computed analytically. If the logit is constant or only varies by distance then only one integral needs to be computed rather than an integral for each observation.

Value

Logical TRUE if condition holds and FALSE otherwise

Author(s)

Jeff Laake

Description

Determines whether the specified logit model is constant for all observations. If it is constant then only one integral needs to be computed.

Usage

is.logistic.constant(xmat, g0model, width)

Arguments

Value

logical value

Author(s)

Jeff Laake

Description

The key function contains one parameter, the scale. Current implementation assumes that scaled dist is x/scale, not x/width

Usage

keyfct.grad.hn(distance, key.scale)

Arguments

Details

d key / d scale = exp(-y ^ 2 / (2 scale ^ 2)) * (y ^ 2 / scale ^ 3)

Value

vector of derivatives of the half-normal key function w.r.t. the scale parameter

Description

The key function contains two parameters, the scale and the shape, and so the gradient is two-dimensional. Current implementation assumes that scaled dist is x/scale, not x/width

Usage

keyfct.grad.hz(distance, key.scale, key.shape, shape = FALSE)

Arguments

Details

d key / d scale = (shape * exp(-(1/ (x/scale) ^ shape)) / ((x/scale) ^ shape ) * scale) d key / d shape = - ((log(x / scale) * exp(-(1/ (x/scale) ^ shape))) / (x/scale) ^ shape)

When distance = 0, the gradients are also zero. However, the equation below will result in NaN and (-)Inf due to operations such as log(0) or division by zero. We correct for this in line 33.

Value

matrix of derivatives of the hazard rate key function w.r.t. the scale parameter and the shape parameter.

Description

Threshold key function

Usage

keyfct.th1(distance, key.scale, key.shape)

Arguments

Value

vector of probabilities

Description

Threshold key function

Usage

keyfct.th2(distance, key.scale, key.shape)

Arguments

Value

vector of probabilities

Description

The two-part normal detection function of Becker and Christ (2015). Either side of an estimated apex in the distance histogram has a half-normal distribution, with differing scale parameters. Covariates may be included but affect both sides of the function.

Usage

keyfct.tpn(distance, ddfobj)

Arguments

Details

Two-part normal models have 2 important parameters:

• The apex, which estimates the peak in the detection function (where g(x)=1). The log apex is reported in summary results, so taking the exponential of this value should give the peak in the plotted function (see examples).

• The parameter that controls the difference between the sides .dummy_apex_side, which is automatically added to the formula for a two-part normal model. One can add interactions with this variable as normal, but don't need to add the main effect as it will be automatically added.

Value

a vector of probabilities that the observation were detected given they were at the specified distance and assuming that g(mu)=1

Author(s)

Earl F Becker, David L Miller

References

Becker, E. F., & Christ, A. M. (2015). A Unimodal Model for Double Observer Distance Sampling Surveys. PLOS ONE, 10(8), e0136403. doi:10.1371/journal.pone.0136403

Description

These data represent avian point count surveys conducted at 453 point sample survey locations on the 24,000 (approx) live-fire region of Fort Hood in central Texas. Surveys were conducted by independent double observers (2 per survey occasion) and as such we had a maximum of 3 paired survey histories, giving a maximum of 6 sample occasions (see MacKenzie et al. 2006, MacKenzie and Royle 2005, and Laake et al. 2011 for various sample survey design details). At each point, we surveyed for 5 minutes (technically broken into 3 time intervals of 2, 2, and 1 minutes; not used here) and we noted detections by each observer and collected distance to each observation within a set of distance bins (0-25, 25-50, 50-75, 75-100m) of the target species (Black-capped vireo's in this case) for each surveyor. Our primary focus was to use mark-recapture distance sampling methods to estimate density of Black-capped vireo's, and to estimate detection rates for the mark-recapture, distance, and composite model.

Format

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

PointID

Unique identifier for each sample location; locations are the same for both species

VisitNumber

Visit number to the point

Species

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

Distance

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

PairNumber

ID value indicating which observers were paired for that sampling occasion

Observer

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

Detected

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

Date

Date of survey since 15 march 2011

Pred

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

Category

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

Effort

Amount of survey effort at the point

Day

Number of days since 15 March 2011

ObjectID

Unique ID for each paired observations

Details

In addition to detailing the analysis used by Collier et al. (2013, In Review), this example documents the use of mrds for avian point count surveys and shows how density models can be incorporated with occupancy models to develop spatially explicit density surface maps. For those that are interested, for the distance sampling portion of our analysis, we used both conventional distance sampling (cds) and multiple covariate distance sampling (mcds) with uniform and half-normal key functions. For the mark-recapture portion of our analysis, we tended to use covariates for distance (median bin width), observer, and date of survey (days since 15 March 2011).

We combined our mrds density estimates via a Horvitz-Thompson styled estimator with the resource selection function gradient developed in Farrell et al. (2013) and estimated density on an ~3.14ha hexagonal grid across our study area, which provided a density gradient for the Fort Hood military installation. Because there was considerable data manipulation needed for each analysis to structure the data appropriately for use in mrds, rather than wrap each analysis in a single function, we have provided both the Golden-cheeked warbler and Black-capped vireo analyses in their full detail. The primary differences you will see will be changes to model structures and model outputs between the two species.

Author(s)

Bret Collier and Jeff Laake

References

Farrell, S.F., B.A. Collier, K.L. Skow, A.M. Long, A.J. Campomizzi, M.L. Morrison, B. Hays, and R.N. Wilkins. 2013. Using LiDAR-derived structural vegetation characteristics to develop high-resolution, small-scale, species distribution models for conservation planning. Ecosphere 43(3): 42. http://dx.doi.org/10.1890/ES12-000352.1

Laake, J.L., B.A. Collier, M.L. Morrison, and R.N. Wilkins. 2011. Point-based mark recapture distance sampling. Journal of Agricultural, Biological and Environmental Statistics 16: 389-408.

Collier, B.A., S.L. Farrell, K.L. Skow, A. M. Long, A.J. Campomizzi, K.B. Hays, J.L. Laake, M.L. Morrison, and R.N. Wilkins. 2013. Spatially explicit density of endangered avian species in a disturbed landscape. Auk, In Review.

Examples

## Not run: 
data(lfbcvi)
xy=cut(lfbcvi$Pred, c(-0.0001, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1),
  labels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
x=data.frame(lfbcvi, New=xy)

# Note that I scaled the individual covariate of day-helps with
# convergence issues
bird.data <- data.frame(object=x$ObjectID, observer=x$Observer,
                        detected=x$Detected, distance=x$Distance,
                        Region.Label=x$New, Sample.Label=x$PointID,
                        Day=(x$Day/max(x$Day)))

# make observer a factor variable
bird.data$observer=factor(bird.data$observer)

# Jeff Laake suggested this snippet to quickly create distance medians
# which adds bin information to the bird.data dataframe

bird.data$distbegin=0
bird.data$distend=100
bird.data$distend[bird.data$distance==12.5]=25
bird.data$distbegin[bird.data$distance==37.5]=25
bird.data$distend[bird.data$distance==37.5]=50
bird.data$distbegin[bird.data$distance==62.5]=50
bird.data$distend[bird.data$distance==62.5]=75
bird.data$distbegin[bird.data$distance==87.5]=75
bird.data$distend[bird.data$distance==87.5]=100

# Removed all survey points with distance=NA for a survey event;
# hence no observations for use in ddf() but needed later
bird.data=bird.data[complete.cases(bird.data),]

# Manipulations on full dataset for various data.frame creation for
# use in density estimation using dht()

#Samples dataframe
xx=x
x=data.frame(PointID=x$PointID, Species=x$Species,
             Category=x$New, Effort=x$Effort)
x=x[!duplicated(x$PointID),]
point.num=table(x$Category)
samples=data.frame(PointID=x$PointID, Region.Label=x$Category,
                   Effort=x$Effort)
final.samples=data.frame(Sample.Label=samples$PointID,
                         Region.Label=samples$Region.Label,
                         Effort=samples$Effort)

#obs dataframe
obs=data.frame(ObjectID=xx$ObjectID, PointID=xx$PointID)
#used to get Region and Sample assigned to ObjectID
obs=merge(obs, samples, by=c("PointID", "PointID"))
obs=obs[!duplicated(obs$ObjectID),]
obs=data.frame(object=obs$ObjectID, Region.Label=obs$Region.Label,
               Sample.Label=obs$PointID)

region.data=data.frame(Region.Label=c(1, 2, 3,4,5,6,7,8,9, 10),
Area=c(point.num[1]*3.14, point.num[2]*3.14,
       point.num[3]*3.14, point.num[4]*3.14,
       point.num[5]*3.14, point.num[6]*3.14,
       point.num[7]*3.14, point.num[8]*3.14,
       point.num[9]*3.14, point.num[10]*3.14))

# Candidate Models

BV1=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
BV1FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
BV2=ddf(
   dsmodel=~mcds(key="hr",formula=~1),
   mrmodel=~glm(~distance),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
BV3=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance+observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV3FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance+observer),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV4=ddf(
   dsmodel=~mcds(key="hr",formula=~1),
   mrmodel=~glm(~distance+observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV5=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV5FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV6=ddf(
   dsmodel=~mcds(key="hr",formula=~1),
   mrmodel=~glm(~distance*observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV7=ddf(
   dsmodel=~cds(key="hn",formula=~1),
   mrmodel=~glm(~distance*Day),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV7FI=ddf(
   dsmodel=~cds(key="hn",formula=~1),
   mrmodel=~glm(~distance*Day),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV8=ddf(
   dsmodel=~cds(key="hr",formula=~1),
   mrmodel=~glm(~distance*Day),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV9=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer*Day),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV9FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer*Day),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV10=ddf(
   dsmodel=~mcds(key="hr",formula=~1),
   mrmodel=~glm(~distance*observer*Day),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
#BV.DS=ddf(
#    dsmodel=~mcds(key="hn",formula=~1),
#    data=bird.data,
#    method="ds",
#    meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))

#AIC table building code.
AIC = c(BV1$criterion, BV1FI$criterion, BV2$criterion, BV3$criterion,
        BV3FI$criterion, BV4$criterion,  BV5$criterion, BV5FI$criterion,
        BV6$criterion, BV7$criterion, BV7FI$criterion, BV8$criterion,
        BV9$criterion, BV9FI$criterion, BV10$criterion)

#creates a set of row names for me to check my grep() call below
rn = c("BV1", "BV1FI", "BV2", "BV3",  "BV3FI", "BV4", "BV5", "BV5FI",
       "BV6", "BV7", "BV7FI", "BV8", "BV9", "BV9FI", "BV10")

#Number parameters
k = c(length(BV1$par), length(BV1FI$par), length(BV2$par),
      length(BV3$par), length(BV3FI$par), length(BV4$par),
      length(BV5$par),length(BV5FI$par), length(BV6$par),
      length(BV7$par), length(BV7FI$par), length(BV8$par),

#build AIC table
AIC.table=data.frame(AIC = AIC, rn=rn, k=k, dAIC = abs(min(AIC)-AIC) ,
                     likg=exp(-.5*(abs(min(AIC)-AIC))))
#row.names(AIC.table)=grep("BV", ls(), value=TRUE)
AIC.table=AIC.table[with(AIC.table, order(-likg, -dAIC, AIC, k)),]
AIC.table=data.frame(AIC.table, wi=AIC.table$likg/sum(AIC.table$likg))
AIC.table

# Model average N_hat_covered estimates
#  not very clean, but I wanted to show full process, need to use
#  collect.models and model.table here later on
estimate <- c(BV1$Nhat, BV1FI$Nhat, BV2$Nhat, BV3$Nhat, BV3FI$Nhat,
              BV4$Nhat,  BV5$Nhat, BV5FI$Nhat, BV6$Nhat, BV7$Nhat,
              BV7FI$Nhat, BV8$Nhat, BV9$Nhat, BV9FI$Nhat, BV10$Nhat)

AIC.values=AIC

# had to use str() to extract here as Nhat.se is calculated in
# mrds:::summary.io, not in ddf(), so it takes a bit
std.err <- c(summary(BV1)$Nhat.se, summary(BV1FI)$Nhat.se,
             summary(BV2)$Nhat.se, summary(BV3)$Nhat.se,
             summary(BV3FI)$Nhat.se, summary(BV4)$Nhat.se,
             summary(BV5)$Nhat.se, summary(BV5FI)$Nhat.se,
             summary(BV6)$Nhat.se, summary(BV7)$Nhat.se,
             summary(BV7FI)$Nhat.se,summary(BV8)$Nhat.se,
             summary(BV9)$Nhat.se, summary(BV9FI)$Nhat.se,
             summary(BV10)$Nhat.se)

## End(Not run)

## Not run: 
#Not Run
#requires RMark
library(RMark)
#uses model.average structure to model average real abundance estimates for
#covered area of the surveys
  mmi.list=list(estimate=estimate, AIC=AIC.values, se=std.err)
  model.average(mmi.list, revised=TRUE)

#Not Run
#Summary for the top 2 models
 #summary(BV5, se=TRUE)
 #summary(BV5FI, se=TRUE)

#Not Run
#Best Model
 #best.model=AIC.table[1,]

#Not Run
#GOF for models
#ddf.gof(BV5, breaks=c(0, 25, 50, 75, 100))

#Not Run
#Density estimation across occupancy categories
#out.BV=dht(BV5, region.data, final.samples, obs, se=TRUE,
#           options=list(convert.units=.01))

#Plot--Not Run

#Composite Detection Function
#plot(BV5, which=3, showpoints=FALSE, angle=0, density=0, col="black", lwd=3,
# main="Black-capped Vireo",xlab="Distance (m)", las=1, cex.axis=1.25,
# cex.lab=1.25)


## End(Not run)

Description

These data represent avian point count surveys conducted at 453 point sample survey locations on the 24,000 (approx) live-fire region of Fort Hood in central Texas. Surveys were conducted by independent double observers (2 per survey occasion) and as such we had a maximum of 3 paired survey histories, giving a maximum of 6 sample occasions (see MacKenzie et al. 2006, MacKenzie and Royle 2005, and Laake et al. 2011 for various sample survey design details). At each point, we surveyed for 5 minutes (technically broken into 3 time intervals of 2, 2, and 1 minutes; not used here) and we noted detections by each observer and collected distance to each observation within a set of distance bins (0-50, 50-100m; Laake et al. 2011) of the target species (Golden-cheeked warblers in this case) for each surveyor. Our primary focus was to use mark-recapture distance sampling methods to estimate density of Golden-cheeked warblers, and to estimate detection rates for the mark-recapture, distance, and composite model.

Format

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

PointID

Unique identifier for each sample location; locations are the same for both species

VisitNumber

Visit number to the point

Species

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

Distance

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

PairNumber

ID value indicating which observers were paired for that sampling occasion

Observer

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

Detected

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

Date

Date of survey since 15 March 2011, numeric value

Pred

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

Category

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

Effort

Amount of survey effort at the point

Day

Number of days since 15 March 2011, numeric value

ObjectID

Unique ID for each paired observations

Details

In addition to detailing the analysis used by Collier et al. (2013, In Review), this example documents the use of mrds for avian point count surveys and shows how density models can be incorporated with occupancy models to develop spatially explicit density surface maps. For those that are interested, for the distance sampling portion of our analysis, we used both conventional distance sampling (cds) and multiple covariate distance sampling (mcds) with uniform and half-normal key functions. For the mark-recapture portion of our analysis, we tended to use covariates for distance (median bin width), observer, and date of survey (days since 15 March 2011).

We combined our mrds density estimates via a Horvitz-Thompson styled estimator with the resource selection function gradient developed in Farrell et al. (2013) and estimated density on an ~3.14ha hexagonal grid across our study area, which provided a density gradient for Fort Hood. Because there was considerable data manipulation needed for each analysis to structure the data appropriately for use in mrds, rather than wrap each analysis in a single function, we have provided both the Golden-cheeked warbler and Black-capped vireo analyses in their full detail. The primary differences you will see will be changes to model structures and model outputs between the two species.

Author(s)

Bret Collier and Jeff Laake

References

Farrell, S.F., B.A. Collier, K.L. Skow, A.M. Long, A.J. Campomizzi, M.L. Morrison, B. Hays, and R.N. Wilkins. 2013. Using LiDAR-derived structural vegetation characteristics to develop high-resolution, small-scale, species distribution models for conservation planning. Ecosphere 43(3): 42. http://dx.doi.org/10.1890/ES12-000352.1

Laake, J.L., B.A. Collier, M.L. Morrison, and R.N. Wilkins. 2011. Point-based mark recapture distance sampling. Journal of Agricultural, Biological and Environmental Statistics 16: 389-408.

Collier, B.A., S.L. Farrell, K.L. Skow, A.M. Long, A.J. Campomizzi, K.B. Hays, J.L. Laake, M.L. Morrison, and R.N. Wilkins. 2013. Spatially explicit density of endangered avian species in a disturbed landscape. Auk, In Review.

Examples

## Not run: 
data(lfgcwa)
xy <- cut(lfgcwa$Pred, c(-0.0001, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1),
 labels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
x <- data.frame(lfgcwa, New=xy)

# Note that I scaled the individual covariate of day-helps with
# convergence issues
bird.data <- data.frame(object=x$ObjectID, observer=x$Observer,
                        detected=x$Detected, distance=x$Distance,
                        Region.Label=x$New, Sample.Label=x$PointID,
                        Day=(x$Day/max(x$Day)))

# make observer a factor variable
bird.data$observer=factor(bird.data$observer)

# Jeff Laake suggested this snippet to quickly create distance medians
# which adds bin information to the \code{bird.data} dataframe

bird.data$distbegin=0
bird.data$distend=100
bird.data$distend[bird.data$distance==12.5]=50
bird.data$distbegin[bird.data$distance==37.5]=0
bird.data$distend[bird.data$distance==37.5]=50
bird.data$distbegin[bird.data$distance==62.5]=50
bird.data$distend[bird.data$distance==62.5]=100
bird.data$distbegin[bird.data$distance==87.5]=50
bird.data$distend[bird.data$distance==87.5]=100

# Removed all survey points with distance=NA for a survey event;
# hence no observations for use in \code{ddf()} but needed later
bird.data=bird.data[complete.cases(bird.data),]

# Manipulations on full dataset for various data.frame creation
# for use in density estimation using \code{dht()}

# Samples dataframe
xx <- x
x <- data.frame(PointID=x$PointID, Species=x$Species,
                Category=x$New, Effort=x$Effort)
x <- x[!duplicated(x$PointID),]
point.num <- table(x$Category)
samples <- data.frame(PointID=x$PointID, Region.Label=x$Category,
                      Effort=x$Effort)
final.samples=data.frame(Sample.Label=samples$PointID,
                         Region.Label=samples$Region.Label,
                         Effort=samples$Effort)

# obs dataframe
obs <- data.frame(ObjectID=xx$ObjectID, PointID=xx$PointID)
# used to get Region and Sample assigned to ObjectID
obs <- merge(obs, samples, by=c("PointID", "PointID"))
obs <- obs[!duplicated(obs$ObjectID),]
obs <- data.frame(object=obs$ObjectID, Region.Label=obs$Region.Label,
                  Sample.Label=obs$PointID)

#Region.Label dataframe
region.data <- data.frame(Region.Label=c(1,2,3,4,5,6,7,8,9),
                          Area=c(point.num[1]*3.14,
                                 point.num[2]*3.14,
                                 point.num[3]*3.14,
                                 point.num[4]*3.14,
                                 point.num[5]*3.14,
                                 point.num[6]*3.14,
                                 point.num[7]*3.14,
                                 point.num[8]*3.14,
                                 point.num[9]*3.14))

# Candidate Models

GW1=ddf(
   dsmodel=~cds(key="unif", adj.series="cos", adj.order=1,adj.scale="width"),
   mrmodel=~glm(~distance),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW2=ddf(
   dsmodel=~cds(key="unif", adj.series="cos", adj.order=1,adj.scale="width"),
   mrmodel=~glm(~distance+observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW3=ddf(
   dsmodel=~cds(key="unif", adj.series="cos", adj.order=1,adj.scale="width"),
   mrmodel=~glm(~distance*observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW4=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW4FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
GW5=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance+observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW5FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance+observer),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW6=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW6FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW7=ddf(
   dsmodel=~cds(key="hn",formula=~1),
   mrmodel=~glm(~distance*Day),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW7FI=ddf(
   dsmodel=~cds(key="hn",formula=~1),
   mrmodel=~glm(~distance*Day),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW8=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer*Day),
   data=bird.data,
   method="io",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
GW8FI=ddf(
   dsmodel=~mcds(key="hn",formula=~1),
   mrmodel=~glm(~distance*observer*Day),
   data=bird.data,
   method="io.fi",
   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
#GWDS=ddf(
#   dsmodel=~mcds(key="hn",formula=~1),
#   data=bird.data,
#   method="ds",
#   meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))



#### GCWA Summary Metrics

#AIC table building code, not exactly elegant, but I did not
want to add more package dependencies
AIC = c(GW1$criterion, GW2$criterion, GW3$criterion, GW4$criterion,
        GW4FI$criterion, GW5$criterion, GW5FI$criterion,
        GW6$criterion, GW6FI$criterion, GW7$criterion, GW7FI$criterion,
        GW8$criterion, GW8FI$criterion)

#creates a set of row names for me to check my grep() call below
rn <- c("GW1", "GW2", "GW3", "GW4", "GW4FI", "GW5", "GW5FI", "GW6",
        "GW6FI", "GW7","GW7FI", "GW8", "GW8FI")

# number of parameters for each model
k <- c(length(GW1$par), length(GW2$par), length(GW3$par), length(GW4$par),
       length(GW4FI$par), length(GW5$par), length(GW5FI$par),
       length(GW6$par), length(GW6FI$par), length(GW7$par),
       length(GW7FI$par), length(GW8$par), length(GW8FI$par))

# build AIC table and
AIC.table <- data.frame(AIC = AIC, rn=rn, k=k, dAIC = abs(min(AIC)-AIC),
                        likg = exp(-.5*(abs(min(AIC)-AIC))))
# row.names(AIC.table)=grep("GW", ls(), value=TRUE)
AIC.table <- AIC.table[with(AIC.table, order(-likg, -dAIC, AIC, k)),]
AIC.table <- data.frame(AIC.table, wi=AIC.table$likg/sum(AIC.table$likg))
AIC.table

# Model average N_hat_covered estimates
# not very clean, but I wanted to show full process, need to use
# collect.models and model.table here

estimate <- c(GW1$Nhat, GW2$Nhat, GW3$Nhat, GW4$Nhat, GW4FI$Nhat,
              GW5$Nhat, GW5FI$Nhat, GW6$Nhat, GW6FI$Nhat, GW7$Nhat,
              GW7FI$Nhat, GW8$Nhat, GW8FI$Nhat)
AIC.values <- AIC

# Nhat.se is calculated in mrds:::summary.io, not in ddf(), so
# it takes a bit to pull out
std.err <- c(summary(GW1)$Nhat.se, summary(GW2)$Nhat.se,
             summary(GW3)$Nhat.se,summary(GW4)$Nhat.se,
             summary(GW4FI)$Nhat.se, summary(GW5)$Nhat.se,
             summary(GW5FI)$Nhat.se, summary(GW6)$Nhat.se,
             summary(GW6FI)$Nhat.se, summary(GW7)$Nhat.se,
             summary(GW7FI)$Nhat.se,summary(GW8)$Nhat.se,
             summary(GW8FI)$Nhat.se)

## End(Not run)
## Not run: 
#Not Run
#requires RMark
library(RMark)
#uses model.average structure to model average real abundance estimates for
#covered area of the surveys
mmi.list=list(estimate=estimate, AIC=AIC.values, se=std.err)
model.average(mmi.list, revised=TRUE)

#Not Run
#Best Model FI
#best.modelFI=AIC.table[1,]
#best.model
#Best Model PI
#best.modelPI=AIC.table[2,]
#best.modelPI

#Not Run
#summary(GW7FI, se=TRUE)
#summary(GW7, se=TRUE)

#Not Run
#GOF for models
#ddf.gof(GW7, breaks=c(0,50,100))

#Not Run
#Density estimation across occupancy categories
#out.GW=dht(GW7, region.data, final.samples, obs, se=TRUE,
           options=list(convert.units=.01))

#Plots--Not Run
#Composite Detection Function examples
#plot(GW7, which=3, showpoints=FALSE, angle=0, density=0,
#     col="black", lwd=3, main="Golden-cheeked Warbler",
#     xlab="Distance (m)", las=1, cex.axis=1.25, cex.lab=1.25)

#Conditional Detection Function
#dd=expand.grid(distance=0:100,Day=(4:82)/82)
#dmat=model.matrix(~distance*Day,dd)
#dd$p=plogis(model.matrix(~distance*Day,dd)%*%coef(GW7$mr)$estimate)
#dd$Day=dd$Day*82
#with(dd[dd$Day==12,],plot(distance,p,ylim=c(0,1), las=1,
# ylab="Detection probability", xlab="Distance (m)",
#  type="l",lty=1, lwd=3, bty="l", cex.axis=1.5, cex.lab=1.5))
#with(dd[dd$Day==65,],lines(distance,p,lty=2, lwd=3))
#ch=paste(bird.data$detected[bird.data$observer==1],
#         bird.data$detected[bird.data$observer==2],
#         sep="")
#tab=table(ch,cut(82*bird.data$Day[bird.data$observer==1],c(0,45,83)),
# cut(bird.data$distance[bird.data$observer==1],c(0,50,100)))
#tabmat=cbind(colMeans(rbind(tab[3,,1]/colSums(tab[2:3,,1],
#                            tab[3,,1]/colSums(tab[c(1,3),,1])))),
#             colMeans(rbind(tab[3,,2]/colSums(tab[2:3,,2],
#                            tab[3,,2]/colSums(tab[c(1,3),,2])))))
#colnames(tabmat)=c("0-50","51-100")
#points(c(25,75),tabmat[1,],pch=1, cex=1.5)
#points(c(25,75),tabmat[2,],pch=2, cex=1.5)

# Another alternative plot using barplot instead of points
# (this is one in paper)

#ch=paste(bird.data$detected[bird.data$observer==1],
#         bird.data$detected[bird.data$observer==2],
#sep="")
#tab=table(ch,cut(82*bird.data$Day[bird.data$observer==1],c(0,45,83)),
# cut(bird.data$distance[bird.data$observer==1],c(0,50,100)))
#tabmat=cbind(colMeans(rbind(tab[3,,1]/colSums(tab[2:3,,1],
#                            tab[3,,1]/colSums(tab[c(1,3),,1])))),
#colMeans(rbind(tab[3,,2]/colSums(tab[2:3,,2],
#               tab[3,,2]/colSums(tab[c(1,3),,2])))))
#colnames(tabmat)=c("0-50","51-100")
#par(mfrow=c(2, 1), mai=c(1,1,1,1))
#with(dd[dd$Day==12,],
#     plot(distance,p,ylim=c(0,1), las=1,
#          ylab="Detection probability", xlab="",
#          type="l",lty=1, lwd=4, bty="l", cex.axis=1.5, cex.lab=1.5))
#segments(0, 0, .0, tabmat[1,1], lwd=3)
#segments(0, tabmat[1,1], 50, tabmat[1,1], lwd=4)
#segments(50, tabmat[1,1], 50, 0, lwd=4)
#segments(50, tabmat[1,2], 100, tabmat[1,2], lwd=4)
#segments(0, tabmat[1,1], 50, tabmat[1,1], lwd=4)
#segments(100, tabmat[1,2], 100, 0, lwd=4)
#mtext("a",line=-1, at=90)
#with(dd[dd$Day==65,],
#     plot(distance,p,ylim=c(0,1), las=1, ylab="Detection probability",
#          xlab="Distance", type="l",lty=1,
#          lwd=4, bty="l", cex.axis=1.5, cex.lab=1.5))
#segments(0, 0, .0, tabmat[2,1], lwd=4)
#segments(0, tabmat[2,1], 50, tabmat[2,1], lwd=4)
#segments(50, tabmat[2,1], 50, 0, lwd=4)
#segments(50, tabmat[2,2], 50, tabmat[2,1], lwd=4)
#segments(50, tabmat[2,2], 100, tabmat[2,2], lwd=4)
#segments(100, tabmat[2,2], 100, 0, lwd=4)
#mtext("b",line=-1, at=90)

## End(Not run)