```{r setup, include=FALSE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
answer <- TRUE
```

# Using survey data

```{marginfigure, echo=answer}
\huge(Solutions)
```
```{marginfigure}
\includegraphics [] {images/monte-vista-nwr02.jpg}

The Monte Vista National Wildlife Refuge.
```
We use field data to fit different detection function models and estimate density and abundance. The data were collected during line transect surveys of duck nests in Monte Vista National Wildlife Refuge, Colorado, USA in 1967 and 1968. Twenty transects, each 25.75km in length were walked 5 times over the two years.  Total transect length of 128.75km (25.75 $\times$ 5) and a distance out to 2.4m was searched.  Consult \citet{anderson_posphala} for a description of the survey.  Distances of detected nests have been provided in a 'csv' text file in a basic format required by 'Distance'. The columns in the the file are:

- Study.Area - name of the study region (Monte Vista NWR)
- Region.Label - identifier of regions or strata (in this case there is only one region and it is set to 'Default')
- Area - size of the stratum
- Sample.Label - line transect identifier
- Effort - length of each transect
- distance - perpendicular distances (m).

The distances allow different key functions/adjustments to be fitted in the detection function model and, by including the transect lengths and area of the region, density and abundance can be estimated.

# Objectives of the practical

1. Import a text file
1. Understand the structure of a data frame
1. Fit different key functions/adjustments in the detection function model using `ds`
1. Examine the results of an analysis, i.e. `ddf` and `dht` components of a `dsmodel` object  

# Importing the data

The file containing the duck nest survey data is located on the online workshop website.  Either download the file following [this link](https://workshops.distancesampling.org/online-course/exercisepdfs/Ch7/datasets/ducks-area-effort.csv) or the content of the data file can be read directly from the web into an R \citep{r_language} data frame named `nests` via the following command 

```{r, echo=TRUE}
nests <- read.csv(file="https://workshops.distancesampling.org/online-course/exercisepdfs/Ch7/datasets/ducks-area-effort.csv",
                  header=TRUE)
```

The URL of the file location is quite long, we repeat it here so it is more legible.

\scriptsize{https://workshops.distancesampling.org/online-course/exercisepdfs/Ch7/datasets/ducks-area-effort.csv}

\normalsize

This command is made up of several components:

- `read.csv` is a function to read data files of type 'csv' (comma-separated values),
- the function has two arguments specified; `file` specifies the name of the data file and `header=TRUE` specifies that the first row of the data file contains the names of the data columns. 
- the `<-` symbol has assigned the data set to an object called `nests`. Note that there is now an object called `nests` listed on the 'Environment' tab.

To check that the data file has been read into `R` correctly, use the `head` and `tail` 'functions' to look at the top and bottom rows of the data, respectively. To look at the first few rows of `nests` type the following command.

```{r,echo=TRUE, eval=answer}
head(nests, n=2)
```

The `head` function as used above displays the first 6 records of the named object. The argument `n` controls the number of rows to display. To display the _last_ 2 records in the data, type the command:

```{r,echo=TRUE, eval=answer}
tail(nests, n=2)
```

The object `nests` is a dataframe object made up of rows and columns. Use the function `dim` to find out the dimensions of the data set (i.e. the total number of rows and columns):

```{r, echo=TRUE, eval=FALSE}
dim(nests)
```

```{marginfigure, echo=!answer}
__Question:__  How many rows and columns are there in the dataframe `nests`?
```

```{r, echo=FALSE, eval=answer}
dim(nests)
```

Another way to look at a data frame is to move to the 'Environment tab' in R-Studio and click on the rectangle (with the grid); this opens a new tab showing the data.

# Summarising the perpendicular distances

To access an individual column within a data frame use the `$` symbol, for example to summarise the distances:

```{r, echo=TRUE, eval=answer}
summary(nests$distance)
```

Similarly to plot the histogram of distances, the command is:

```{r, echo=TRUE, eval=answer, fig.margin=TRUE, fig.width=3.5, fig.height=3.5}
hist(nests$distance, xlab="Distance (m)")
```

# Fitting different models

To use the `ds` function, load the `Distance` package \citep{distance_pkg}.

The function `ds` requires a data frame to have a column called `distance`, we  specify the name of the data frame as follows:

```{r, echo=TRUE, eval=TRUE}
library(Distance)
conversion <- convert_units("meter", "kilometer", "square kilometer")
nest.model1 <- ds(nests, key="hn", adjustment=NULL, convert_units = conversion)
```

The `convert_units` argument ensures that the correct units are specified - in this example, distances are in metres, lengths in km and the area in km$^2$.  Think of this argument as a divider used to transform units of transect effort into units of perpendicular distance (e.g., 1km / 0.001 = 1000m).  The `convert_units` function performs these calculation for you if you provide the units in which perpendicular distances, transects and study area size are recorded in the data.

This call to `ds` fits a half-normal key function with no adjustment terms. Summarise the fitted model:

```{r, echo=TRUE, eval=answer}
summary(nest.model1$ddf)
```

Plot the detection function with the histogram having 12 bins:

```{r, echo=TRUE, eval=answer, message=FALSE, fig.margin=TRUE, fig.width=3.5, fig.height=3.5}
plot(nest.model1, nc=12)
```

To fit different detection functions, change the `key` and `adjustment` arguments. For example to fit a half-normal key function with cosine adjustment terms, use the command:

```{r, echo=TRUE, eval=answer, message=FALSE}
nest.model2 <- ds(nests, key="hn", adjustment="cos", convert_units = conversion)
```

By default, AIC selection will be used to fit adjustment terms of up to order 5. 

```{marginfigure}
__Question:__ Have any adjustment terms been selected?  
```

```{r, echo=FALSE, eval=answer, message=FALSE}
summary(nest.model2$ddf)
```

```{marginfigure, echo=answer}
__Answer:__ No adjustment terms have been included in the preferred model, because there is only a single row in the list of `Detection function parameters`.
```

To fit a hazard rate key function with Hermite polynomial adjustment terms, then use the command:

```{r, echo=TRUE, eval=answer, message=FALSE, warning=FALSE}
nest.model3 <- ds(nests, key="hr", adjustment="herm", convert_units = conversion)
```

```{r, echo=FALSE, message=FALSE, warning=FALSE, eval=FALSE}
# used for the final question, executed under all situations
hr.model <- ds(nests, key="hr", adjustment="herm", convert_units = conversion)
```

```{r, echo=TRUE, eval=answer, message=FALSE}
summary(nest.model3$ddf)
```

Use the `help` command to find out what other key functions and adjustment terms are available.

```{marginfigure}
__Question:__ Is there much difference in the probability of detecting a nest given it is within the maximum detection distance between the three models fitted to the duck nest data?
```

```{marginfigure, echo=answer}
__Answer:__ No, there is barely any difference between the estimated probability of detection.  For the first two models, $\hat{P}_a$ = 0.869 and for the hazard rate model, $\hat{P}_a$ = 0.889.
```

# The `ds` object

The objects created with `ds` (e.g. `nest.model1`) are made up of two parts. We can list them using the `names` function as below:

```{r, echo=TRUE, eval=answer}
names(nest.model1)
```

The detection function information is in the `ddf` part and the density and abundance estimates are stored in the `dht` part. To access each part, then the `$` can be used (as with columns in a data frame). For example to see what information is stored in the `ddf` part, we can use the `names` function again:

```{r, echo=TRUE, eval=answer}
names(nest.model1$ddf)
```

# Goodness of fit

Before making inference from the detection function we have fitted, we should evaluate the model.  First assessment is goodness of fit, accomplished using the function `gof_ds`:

```{r, echo=TRUE, eval=answer, fig.margin=TRUE, fig.width=3.5, fig.height=3.5}
gof_ds(nest.model1)
```

Calling the function `gof_ds` with the default arguments and exact distance data, a _q-q_ plot is produced along with the unweighted Cramer-von Mises goodness of fit test.  

```{marginfigure}
__Question:__ Interpret the q-q plot and CvM test results for the duck nest data.
```

```{marginfigure, echo=answer}
__Answer:__ The half-normal detection function fits the duck nest data very well.  All points of the q-q plot fall on the diagonal line and the p-value associated with the CvM test statisic is very large (>> 0.05) indicating a good fit of the model to the data.
```

# Estimating density and abundance

So far, we have concentrated on the detection function but, with more information such as transect lengths and the area of the region, we can estimate density and abundance.  The second component of a `ds` object, contains this additional information.  This information can be viewed with:

```{r, echo=TRUE, eval=answer}
str(nest.model1$dht$individuals, max=1)
```

This `dht` object contains considerable information.  However, focus upon three tables generated by the `summary()` function: summary, abundance and density.  Dig more deeply into the content of these tables.

## Summary information

This provides information about the survey:

- size of study area,
- area covered by sampling effort
- length of all transects
- number of detections
- number of transects
- encounter rate (ER) number of detections per unit transect length and its associated variability

```{r echo=TRUE, eval=FALSE}
nest.model1$dht$individuals$summary
```

```{r fig.margin=FALSE, eval=answer, echo=FALSE}
knitr::kable(nest.model1$dht$individuals$summary, digits=3)
```

## Abundance estimates

Estimated density multiplied by the size of the study area.

```{r echo=TRUE, eval=FALSE}
nest.model1$dht$individuals$N
```

```{r fig.margin=FALSE, eval=answer, echo=FALSE}
knitr::kable(nest.model1$dht$individuals$N, digits=2)
```

## Density estimates

Density estimated using the formula

$$
\hat{D} = \frac{n}{a \hat{P}_a}
$$

where n (number of nests)=`r nest.model1$dht$individuals$summary$n`, a (covered area)=`r nest.model1$dht$individuals$summary$CoveredArea` and $\hat{P_a}$ (probability of detection)=`r round(nest.model1$dht$individuals$average.p,4)`

```{r echo=TRUE, eval=FALSE}
nest.model1$dht$individuals$D
```

```{r fig.margin=FALSE, eval=answer, echo=FALSE}
knitr::kable(nest.model1$dht$individuals$D, digits=2)
```

```{marginfigure}
__Question:__ Compute (by hand) the density estimate resulting from the estimated probability of detection arising from the hazard rate detection function: $\hat{P}_a$= 0.8891.
```

```{marginfigure, echo=answer}
__Answer:__ $\hat{D} = \frac{534}{12.36 \cdot 0.8891} =$ 48.59 nests$\cdot \text{km}^{-2}$ compared to 49.70 nests$\cdot \text{km}^{-2}$ (a difference of 2.3%).  In other words, the difference in estimated density is quite small; particularly when uncertainty is taken into account.
```

<!-- the 2015 book is not cited in the narrative, but I want it in the bibliography -->
\nocite{buckland_distance_2015}


