In this exercise we look at mark-recapture distance sampling (MRDS) models. The first part of this exercise involves analysis of a survey of a known number of golf tees. This is intended mainly to familiarise you with the double-platform data structure and analysis features in the R function `mrds` (Laake \textit{et al}. 2019). The second part of the practical involves analysis of the pack-ice seal survey data of Borchers \textit{et al}. (2006) and Southwell \textit{et al}. (2007).

To help understand the terminology using in MRDS and the output produced by `mrds`, there is a guide available on the workshop website called 'Interpreting MRDS output: making sense of all the numbers'. 

# Aims

The aims of this practical are to learn how to model 

1. trial and independent-observer configuration

2. full and point independence assumptions,

3. include covariates in the detection function(s) and 

4. select between competing models. 

# 1. Golf tee survey

## Golf tee data

These data come from a survey of golf tees which conducted by statistics students at the University of St Andrews. The data were collected along transect lines, 210 metres in total. A distance of 4 metres out from the centre line was searched and, for the purposes of this exercise, we assume that this comprised the total study area, which was divided into two strata. There were 250 clusters of tees in total and 760 individual tees in total. 

The population was independently surveyed by two observer teams. The following data were recorded for each detected group: perpendicular distance, cluster size, observer (team 1 or 2), 'sex' (males are yellow and females are green and golf tees occur in single-sex clusters) and 'exposure'. Exposure was a subjective judgment of whether the cluster was substantially obscured by grass (exposure=0) or not (exposure=1). The lengths of grass varied along the transect line and the grass was slightly more yellow along one part of the line compared to the rest. 

The golf tee dataset is provided as part of the `mrds` package.

Open R and load the `mrds` package and golf tee dataset (called `book.tee.data`). The elements required for an MRDS analysis are contained within the object dataset. These data are in a hierarchical structure (rather than in a 'flat file' format) so that there are separate elements for observations, samples and regions. In the code below, each of these tables is extracted to avoid typing long names.   

```{r preliminaries, message=FALSE, comment=NA, echo=TRUE}
# Load libraries
library(knitr)
library(mrds)

# Access the golf tee data
data(book.tee.data)

# Investigate the structure of the dataset
str(book.tee.data)

# Extract the list elements from the dataset into easy-to-access objects
detections <- book.tee.data$book.tee.dataframe # detection information
region <- book.tee.data$book.tee.region # region info
samples <- book.tee.data$book.tee.samples # transect info
obs <- book.tee.data$book.tee.obs # links detections to transects and regions
```

Let's have a close look at the columns in the `detections` data because it has a particular structure.

```{r headdet, echo=T}
# Check detections
head(detections)
```
The structure of the detection is as follows:

+ each detected object (in this case the object was a group or cluster of golf tees) is given a unique number in the `object` column, 

+ each `object` occurs twice - once for observer 1 and once for observer 2, 

+ the `detected` column indicates whether the object was seen (`detected=1`) or not seen (`detected=0`) by the observer,

+ perpendicular distance is in the `distance` column and cluster size is in the `size` column (the same default names as for the `ds` function). 

To ensure that the variables `sex` and `exposure` are treated correctly, define them as factor variables.

```{r, echo=T}
# Define sex and exposure as factor variables 
detections$sex <- as.factor(detections$sex)
detections$exposure <- as.factor(detections$exposure)
```

## Golf tee survey analyses

### Estimation of \textit{p}(0): distance only

We will start by analysing these data assuming that Observer 2 was generating trials for Observer 1 but not vice versa, i.e. trial configuration where Observer 1 is the primary and Observer 2 is the tracker. (The data could also be analysed in independent observer configuration - you are welcome to try this for yourself). We begin by assuming full independence (i.e. detections between observers are independent at all distances): this requires only a mark-recapture (MR) model and, to start with, perpendicular distance will be included as the only covariate. 

Remember that `?` or `help` can be used to find out more about any of the functions used -- e.g., `?ddf` will tell you more about the ddf function. 

```{r, fit-nocovar, message=FALSE, comment=NA, fig.height=4, fig.cap="Goodness of fit (FI-trial) to golftee data."}
# Fit trial configuration with full independence model
fi.mr.dist <- ddf(method='trial.fi', mrmodel=~glm(link='logit',formula=~distance),
                  data=detections, meta.data=list(width=4))
```

#### Examining `mrds` output

Having fitted the model, we can create tables summarizing the detection data. In the commands below, the tables are created using the `det.tables` function and saved to `detection.tables`. 

```{r}
# Create a set of tables summarizing the double observer data 
detection.tables <- det.tables(fi.mr.dist)
# Print these detection tables
detection.tables
```

The information in detection summary tables can be plotted, but, in the interest of space, only one (out of six possible plots) is shown below. 

```{r dettabplot, echo=T, fig.height=4}
# Plot detection information, change number to see other plots
plot(detection.tables, which=1)
```
The plot numbers are:

1. Histograms of distances for detections by either, or both, observers. The shaded regions show the number for observer 1. 

2. Histograms of distances for detections by either, or both, observers. The shaded regions show the number for observer 2. 

3.	Histograms of distances for duplicates (detected by both observers).

4. Histogram of distances for detections by either, or both, observers. Not shown for trial configuration. 

5. Histograms of distances for observer 2. The shaded regions indicate the number of duplicates - for example, the shaded region is the number of clusters in each distance bin that were detected by Observer 1 given that they were also detected by Observer 2 (the "|" symbol in the plot legend means "given that"). 

6. Histograms of distances for observer 1. The shaded regions indicate the number of duplicates as for plot 5. Not shown for trial configuration. 

Note that if an independent observer configuration had been chosen, all plots would be available. 

A summary of the detection function model is available using the `summary` function. The Q-Q plot has the same interpretation as a Q-Q plot in a conventional, single platform analysis. 

```{r fi.summary, echo=T, fig.height=4}
# Produce a summary of the fitted detection function object
summary(fi.mr.dist)

# Produce goodness of fit statistics and a qq plot
gof.result <- ddf.gof(fi.mr.dist, 
                      main="Full independence, trial configuration\ngoodness of fit Golf tee data")

# Extract chi-square statistics for reporting
chi.distance <- gof.result$chisquare$chi1$chisq
chi.markrecap <- gof.result$chisquare$chi2$chisq
chi.total <- gof.result$chisquare$pooled.chi
```

Abbreviated $\chi^2$ goodness-of-fit assessment shows the $\chi^2$ contribution from the distance sampling model to be `r round(chi.distance,1)` and the $\chi^2$ contribution from the mark-recapture model to be `r round(chi.markrecap,1)`.  The combination of these elements produces a total $\chi^2$ of `r round(chi.total$chisq,1)` with `r chi.total$df` degrees of freedom, resulting in a $p$-value of `r round(chi.total$p,3)`

The (two) detection functions can be plotted:

```{r plotdf, echo=T, fig.height=3.5}
# Divide the plot region into 2 columns
par(mfrow=c(1,2))
# Plot detection functions
plot(fi.mr.dist)
```

The plot headed 

+ "Observer=1 detections" shows a histogram of Observer 1 detections with the estimated Observer 1 detection function overlaid on it and adjusted for \textit{p}(0). The dots show the estimated detection probability for all Observer 1 detections.

+ "Conditional detection probability" shows the proportion of Obs 2's detections that were detected by Obs 1 (also see the detection tables). The fitted line is the estimated detection probability function for Obs 1 (given detection by Obs 2) - this is the MR model. Dots are estimated detection probabilities for each Obs 1 detection.

There is some evidence of unmodelled heterogeneity in that the fitted line in the left-hand plot declines more slowly than the histogram as the distance increases. 

#### Estimating abundance

Abundance is estimated using the `dht` function. In this function, we need to supply information about the transects and survey regions.  

```{r, abund-from-dist}
# Calculate density estimates using the dht function
tee.abund <- dht(model=fi.mr.dist, region.table=region, sample.table=samples, 
                 obs.table=obs)

# Print out results in a nice format
pander::pander(tee.abund$individuals$summary, digits=2, 
      caption="Survey summary statistics for golftees")
pander::pander(tee.abund$individuals$N, digits=2, 
      caption="Abundance estimates for golftee population with two strata")
```

The estimated abundance is 593 (recall that the true abundance is 760) and so this is somewhat negatively biased. The 95\% confidence interval does not include the true value. 

### Estimation of \textit{p}(0): distance and other explanatory variables

How about including the other covariates, `size`, `sex` and `exposure`, in the MR model? Which MR model would you use? Don't spend too long on this - just try a couple of models. In the command below, `distance` and `sex` are included in the detection function - remember `sex` was defined as a factor earlier on.

In the code below, all possible models (excluding interaction terms) are fitted. 

```{r}
# Full independence model
# Set up list with possible models
mr.formula <- c("~distance","~distance+size","~distance+sex","~distance+exposure","~distance+size+sex","~distance+size+exposure","~distance+sex+exposure","~distance+size+sex+exposure")
num.mr.models <- length(mr.formula)
# Create dataframe to store results
fi.results <- data.frame(MRmodel=mr.formula, AIC=rep(NA,num.mr.models))
# Loop through all MR models
for (i in 1:num.mr.models) {
  fi.model  <- ddf(method='trial.fi', 
                   mrmodel=~glm(link='logit',formula=as.formula(mr.formula[i])),
                  data=detections, meta.data=list(width=4))
  fi.results$AIC[i] <- summary(fi.model)$aic
}
# Calculate delta AIC
fi.results$deltaAIC <- fi.results$AIC - min(fi.results$AIC)
# Order by delta AIC
fi.results <- fi.results[order(fi.results$deltaAIC), ]
# Print results in pretty way
pander::pander(fi.results)
```

We see that the preferred model contains `distance + sex + exposure` and so let's check the goodness-of-fit statistics and detection function plots. 

```{r, }
# Fit chosen model
fi.mr.dist.sex.exp  <- ddf(method='trial.fi', mrmodel=~glm(link='logit',formula=~distance+sex+exposure),
                           data=detections, meta.data=list(width=4))

# Check goodness-of-fit 
ddf.gof(fi.mr.dist.sex.exp, main="FI trial mode\nMR=dist+sex+exp")

# Plot detection functions
par(mfrow=c(1,2))
plot(fi.mr.dist.sex.exp)

# Get abundance estimates 
tee.abund.fi <- dht(model=fi.mr.dist.sex.exp, region.table=region, sample.table=samples, obs.table=obs)

# Print results
tee.abund.fi
```

This model incorporates the effect of more variables causing the heterogeneity. The estimated abundance is 651 which is less biased than the previous estimate and the 95\% confidence intervals (494, 858) contain the true value. 

The model is a reasonable fit to the data (i.e. non-significant $\chi^2$ and Cramer von Mises tests). This model has a lower AIC (405.7) than the model with only distance (452) and so is to be preferred. 

### Point independence

A less restrictive assumption than full independence is point independence, which assumes that detections are only independent on the transect centre line i.e. at perpendicular distance zero. 

Let's start by seeing if a simple point independence model is better than a simple full independence one. This requires that a distance sampling (DS) model is specified as well a MR model. Here we try a half-normal key function for the DS model. 

```{r pit-nocovar, echo=T, fig.height=3.5}
# Fit trial configuration with point independence model
pi.mr.dist <- ddf(method='trial', 
                  mrmodel=~glm(link='logit', formula=~distance),
                  dsmodel=~cds(key='hn'), 
                  data=detections, meta.data=list(width=4))

# Summary pf the model 
summary(pi.mr.dist)

# Produce goodness of fit statistics and a qq plot
gof.results <- ddf.gof(pi.mr.dist, main="Point independence, trial configuration\n goodness of fit Golftee data")
```

Compare the results with the corresponding full independence model. Which has the lower AIC? Which has an estimate closer to known true abundance. 

The AIC for this point independence model is 452.02 which is marginally smaller than the first full independence model that was fitted and hence is to be preferred. 

```{r,}
# Get abundance estimates 
tee.abund.pi <- dht(model=pi.mr.dist, region.table=region, sample.table=samples, obs.table=obs)

# Print results
tee.abund.pi
```

This results in an estimated abundance of 687. Can we do better if more covariates are included in the DS model? 

#### Covariates in the DS model

To include covariates in the DS detection function, we need to specify an MCDS model as follows:

```{r pit-mod1, echo=T, fig.height=3.5}
# Fit the PI-trial model - DS sex and MR distance 
pi.mr.dist.ds.sex <- ddf(method='trial', 
                         mrmodel=~glm(link='logit',formula=~distance),
                         dsmodel=~mcds(key='hn',formula=~sex), 
                         data=detections, meta.data=list(width=4))
```

Use the `summary` function to check the AIC and decide if you are going to include any additional covariates in the detection function. 

Now try a point independence model that has the preferred MR model from your full independence analyses. Which has the lower AIC and bias? 


```{r}
# Point independence model
# Include covariates in DS model
# Use selected MR model
# DS models
ds.formula <- c("~size","~sex","~exposure","~size+sex","~size+exposure","~sex+exposure",
                "~size+sex+exposure")
num.ds.models <- length(ds.formula)
# Create dataframe to store results
pi.results <- data.frame(DSmodel=ds.formula, AIC=rep(NA,num.ds.models))
# Loop through ds models - use selected MR model from earlier
for (i in 1:num.ds.models) {
  pi.model <- ddf(method='trial', mrmodel=~glm(link='logit',formula=~distance+sex+exposure),
                  dsmodel=~mcds(key='hn',formula=as.formula(ds.formula[i])), 
                  data=detections, meta.data=list(width=4))
  pi.results$AIC[i] <- summary(pi.model)$AIC
}
# Calculate delta AIC
pi.results$deltaAIC <- pi.results$AIC - min(pi.results$AIC)
# Order by delta AIC
pi.results <- pi.results[order(pi.results$deltaAIC), ]
# Print results in pretty way
pander::pander(pi.results)
```

This indicates that `sex` should be included in the DS model. We do this and check the goodness of fit and obtain abundance.  

```{r}
# Fit chosen model
pi.ds.sex <- ddf(method='trial', mrmodel=~glm(link='logit',formula=~distance+sex+exposure),
                dsmodel=~mcds(key='hn',formula=~sex), data=detections,  
                meta.data=list(width=4))
summary(pi.ds.sex)

# Check goodness-of-fit 
ddf.gof(pi.ds.sex, main="PI trial configutation\nGolfTee DS model sex")

# Get abundance estimates 
tee.abund.pi <- dht(model=pi.ds.sex, region.table=region, sample.table=samples, obs.table=obs)
tee.abund.pi
```

This model estimated an abundance of 695, which is closest to the true value of all the models - it is still less than the true value indicating, perhaps, some unmodelled heterogeneity on the trackline (or perhaps just bad luck - remember this was only one survey). 

Was this complex modelling worthwhile? In this case, the estimated $p(0)$ for the best model was 0.97 (which is very close to 1). If we ran a conventional distance sampling analysis, pooling the data from the two observers, we should get a very robust estimate of true abundance. 

# 2. Crabeater seal survey

## Crabeater seal data 

This analysis is described in Borchers \textit{et al}. (2006) and Southwell \textit{et al}. (2007). These data come from a helicopter survey of crabeater seals conducted by the Australian Antarctic Division within the pack-ice seals programme. The helicopter could only operate within a relatively short distance from the ice-breaker ship which acted as its base. The ice-breaker could only go where the pack ice was thin enough and so the aerial transects could not be located at random. This means that design-based estimation was not a valid option and so, in the published analysis, abundance was estimated using density surface modelling. For the purposes of this exercise, we concentrate on detection function estimation and create an artificial region as a device to produce abundance estimates. 

There were four independent observers in the helicopter, two on each side (front and back). The front observers were considered to be one 'team' and the back observers were considered to be the other 'team'. Various environmental factors were recorded. In addition to perpendicular distance and cluster size, the following explanatory variables are available:

* side - the side of the helicopter from which seal were seen (L and R)
* exp - the experience (in survey hours) of the observer
* fatigue - the number of minutes the observer had been on duty on the current flight
* gscat - group size category (1, 2 and greater than or equal to 3)
* vis - visibility category (Poor, Good and Excellent)
* glare - whether there was glare (Yes or No)
* ssmi - a measure of ice cover
* altitude - the height of the aircraft in metres
* obsname - unique identifier of observer

The data from the survey has been saved in a `.csv` file. This file is read into R using `read.csv`.

```{r, crabdat, comment=NA, echo=T, message=F}
crabseal <- read.csv("crabbieMRDS.csv")

```

## Crabeater seal analyses

The observer teams acted independently and so an 'independent observer' (IO) configuration can be specified. The code below fits simple models (i.e. distance only) with the full independence assumption and the point independence assumption. For each model, make a note of the estimated values for $p(0)$ for each observer and the observers combined. Check goodness-of-fit and plot the detection function. 

#### Full independence

Fit an IO configuration assuming full independence: 

```{r}
# IO configuration - full independence
# MR model - distance only
# Truncation 700m
fi.mr.dist <- ddf(method="io.fi", mrmodel=~glm(link="logit", formula=~distance),
                 data=crabseal, meta.data=list(width=700))
summary(fi.mr.dist)
# Goodness of fit
ddf.gof(fi.mr.dist, main="Crabeater - FI")
# Detection functions
par(mfrow=c(2,2))
plot(fi.mr.dist, which=c(1,2,3,4))
plot(fi.mr.dist, which=c(5,6))

```

#### Point independence

Next fit IO configuration assuming full independence: specify a half-normal key function for the DS model: again only include perpendicular distance in the MR model.

```{r crabdet, echo=T}
# IO configuration - point independence
# MR model - distance only
# DS model - half normal detection function, no additional covars
# Truncation at 700m
pi.mr.dist.ds.hn <- ddf(method="io", dsmodel=~cds(key="hn"),
                 mrmodel=~glm(link="logit", formula=~distance),
                 data=crabseal, meta.data=list(width=700))
summary(pi.mr.dist.ds.hn)

# Goodness of fit
ddf.gof(pi.mr.dist.ds.hn, main="Crabeater - PI")
# Detection functions
par(mfrow=c(2,2))
plot(pi.mr.dist.ds.hn, which=c(1,2,3,4))
plot(pi.mr.dist.ds.hn, which=c(5,6))
```

Which of these models do you prefer? 

The FI model had an AIC of 25681 but is not adequate - clearly seen in the goodness-of-fit plot and statistics and the top two detection function plots. Nevertheless, the conditional detection function plots are pretty good (bottom plots) - the model seems adequate for modelling the conditional probability of one observer detecting a group, given the other detected it, but is not adequate for the unconditional probability of detecting a group (which is what we want to estimate). This implies there is something about detected groups that tends to make them more detectable than other groups (i.e. some unmodelled heterogeneity). Treating the conditional detection functions as if they applied to undetected groups will lead to positive bias in estimating detection probability and negative bias in estimating abundance.

The PI model had an AIC of 25325, hence the PI model is preferred on that account and also the goodness-of-fit statistics do not give cause for alarm. 

If you have time, try a few models with covariates - warning, this is a large dataset and so fitting models and obtaining model summaries can take a long time with complicated models.

The models given below are those that were used in Borchers \textit{et al}. (2006) and Southwell \textit{et al}. (2007). We do not run them here because they take a long time to fit but are provided for comparison with your models. 

```{r}
# Full independence like papers
#mr.model <- "~distance + observer + gscat + fatigue + obsname + exp + ssmi + side + vis + distance:obsname + distance:observer"
#fi.paper <- ddf(method="io.fi", mrmodel=~glm(link="logit", formula=as.formula(mr.model)),
#               data=crabseal, meta.data=list(width=700))

# Point independence like papers
#pi.paper <- ddf(method="io.fi", dsmodel=~mcds(key="hn", formula=~side+vis),
#                mrmodel=~glm(link="logit", formula=as.formula(mr.model)),
#                data=crabseal, meta.data=list(width=700))
```

#### Estimating abundance

Following model criticism and selection, the abundance can be estimated - below we use the simple point independence model. The estimates of abundance for the study area are arbitrary because inference of the study was restricted to the covered region, hence, the estimates of abundance here are artificial. Nevertheless, we illustrate the method to estimate abundance. We require tables of the region, transects and detections - these can easily be created from the data using the `checkdata` function in the `Distance` package (the `:::` is shorthand for accessing a function without loading a package). Using these tables, Horvitz-Thompson-like estimators can be applied to produce estimates of $\hat{N}$. The use of `convert.units=0.001` adjusts the units of perpendicular distance measurement (m) to units of transect effort (km). 

```{r, crabsummary, echo=T}
# Create tables for estimating abundance 
# Selecting observer==1 ensures that observations in the obs.table are unique 
tables <- Distance:::checkdata(crabseal[crabseal$observer==1, ])

# Estimate abundance in covered region
pi.abund <- dht(model=pi.mr.dist.ds.hn,
                region=tables$region.table, 
                sample=tables$sample.table, obs=tables$obs.table,
                se=TRUE, options=list(convert.units=0.001))

# Pretty tables of data summary
pander::pander(pi.abund$individuals$summary, digits=3,
      caption="Summary information from crabeater seal aerial survey.")
```

```{r, crabestimates}
# Pretty tables of estimates of individual abundance
pander::pander(pi.abund$individual$N, digits=3,
      caption="Crabeater seal abundance estimates for study area of arbitrary size.")
```

### Crabeater seals with MCDS (optional)

We can also analyse the crabeater seals data as if it were single platform data (i.e. ignoring that $p(0)$ is less than 1). These data are available in `crabbieMCDS.csv`.

This short exercise guides you through the import of these data into R and fits a simple half-normal detection function examining the possible improvement of the model by incorporating _side of plane_ and _visibility_ covariates.

```{r, message=FALSE, comment=NA}
# Load Distance for MCDS 
library(Distance)
# Read in data
crab.covar <- read.csv("crabbieMCDS.csv")
# Check data imported OK
head(crab.covar, n=3)
```

After checking that the data have been read into `R` appropriately, we are ready to fit a detection function.  

As before, _side of plane_ and _visibility_ are assigned characters and so we need to tell R to treat them as factors. 

```{r, comment=NA}
# Define factor variables
crab.covar$side <- as.factor(crab.covar$side)
crab.covar$vis <- as.factor(crab.covar$vis)
```

With two potential explanatory variables, there are a number of possible models. We start by fitting a detection function with _side of plane_ as a covariate using a half-normal key function. 

```{r, comment=NA}
# Fit HN key function with side of plane
ds.side <- ds(crab.covar, key="hn", formula=~side, truncation=700)
```

We would now like to assess the fit of this function to our data. Two visual assessments are provided by the panels below: histogram and fitted function on the left and Q-Q plot on the right.

```{r, fig.cap="Histogram and fitted half-normal detection function on left.  Q-Q plot of detection function and data on right."}
# Divide plot region
par(mfrow = c(1, 2))
# Create a title for the plot
plot.title <- "Two sets of points\none for each 'side' of plane"
# Plot model
plot(ds.side, pch=19, cex=0.5, main=plot.title)
# Plot qq plot
gof.result <- ds.gof(ds.side, lwd = 2, lty = 1, pch = ".", cex = 0.5)
# Extract gof statistics
message <- paste("CVM GOF p-value=", round(gof.result$dsgof$CvM$p, 4))
# Add gof stats to plot
text(0.6, 0.2, message, cex=0.5)
```

The code below fits the model without any covariates. 

```{r, comment=NA}
# Fit HN key function with no covars and no adjustments
ds.nocov <- ds(crab.covar, key="hn", adjustment=NULL, truncation=700)
```

AIC score for model without covariates is `r options(scipen=9)` `r round(ds.nocov$ddf$criterion)` and AIC score for model with _side_ as a covariate is  `r round(ds.side$ddf$criterion)` so the model without _side_ as a covariate is very slightly preferred.  

We could also fit further detection functions and contrast the resulting models:

* with _visibility_ only
* with _side of plane_ and _visibility_ (excluding an interaction). 

Out of the four possible models which is to be preferred? 

We could go on to produce abundance estimates from our preferred model using the `dht` function if we had provided information about the size of the crabeater seal study area.

In the code below, we fit all possible models and go on examine the selected model, based on AIC. 

```{r}
# Create list of formula
mcds.formula <- c("~1", "~side","~vis","~side+vis")
num.mcds.models <- length(mcds.formula)
# Set up table to save results
crab.results <- data.frame(MCDSmodel=mcds.formula, AIC=rep(NA,num.mcds.models))
# Loop through models
for (i in 1:num.mcds.models) {
  ds.model <- ds(crab.covar, key="hn", formula=as.formula(mcds.formula[i]), adjustment=NULL,
  truncation=700)
  crab.results$AIC[i] <- summary(ds.model$ddf)$aic
}
# Calculate delta AIC
crab.results$deltaAIC <- crab.results$AIC - min(crab.results$AIC)
# Order by delta AIC
crab.results <- crab.results[order(crab.results$deltaAIC), ]
# Print results in pretty way
pander::pander(crab.results)

# Fit selected model
crab.mcds <- ds(crab.covar, key="hn", formula=~side+vis, adjustment=NULL,
                truncation=700)
summary(crab.mcds$ddf)
plot(crab.mcds)
```

Rather than fitting an MRDS model, as above, would an MCDS analyses have been adequate? 

We can see that, from fitting the MRDS models earlier, the estimate of $p(0)$ for both platforms combined is 0.988 hence, the conventional distance sampling assumption that $p(0)=1$ is very nearly satisfied. Therefore, in this example the MCDS model seems adequate. 

## References

Borchers DL, Laake JL, Southwell C and Paxton CGM (2006) Accommodating unmodeled heterogeneity in double-observer distance sampling surveys. Biometrics 62: 371-378

Laake JL, Borchers DL, Thomas L, Miller DL and Bishop JRB (2019) mrds: Mark-Recapture Distance Sampling. R package version 2.2.1.

Southwell C, Borchers DL, Paxton CGM, Burt ML and de la Mare W (2007) Estimation of detection probability in aerial surveys of Antarctic pack-ice seals. Journal of Agricultural, Biological and Environmental Statistics 12:41-54

Thomas L, Buckland ST, Rexstad EA, Laake JL, Strindberg S, Hedley SL, Bishop JRB, Marques TA, and Burnham KP (2010)  Distance software: design and analysis of distance sampling surveys for estimating population size. Journal of Applied Ecology 47: 5-14.  DOI: 10.1111/j.1365-2664.2009.01737.x