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

# Covariates in the detection function
```{marginfigure, echo=answer}
\huge(Solutions)
```
```{marginfigure}
\includegraphics []{images/amakihi1121-300x270.jpg}
Hawaii amakihi *(Hemignathus virens)*
```

We illustrate fitting multiple covariate distance sampling (MCDS) models to point transect data using a bird survey in Hawaii; data on an abundant species, the Hawaii amakihi *(Hemignathus virens)* is used. This practical is based on the case study in \citet[Section 5.3.2]{buckland_distance_2015} which duplicates the analysis presented in \citet{Marques_2007}. This set of data is included in 'Distance for Windows' as one of the Sample Projects: you can open this project (entitled `amakihi.zip`) in the 'Sample projects' directory in the 'My Distance projects' directory residing under 'My Documents'. We describe the analysis of these data using Distance in R \citep{r_language}.

# Objectives of this practical

1. Introduce different types of plots to explore covariates
1. Add covariates to the detection function
1. Plot the detection functions. 

# Importing the data 

Analysis begins by importing the data from a comma-delimited file. This file was created by copying the data from the amakihi Distance project. 

```{r setup_and_read, echo=TRUE, message=FALSE, warning=FALSE}
amakihi <- read.csv(file="https://raw.githubusercontent.com/distanceworkshops/online-course/gh-pages/exercisepdfs/Ch7/datasets/amakihi.csv")
```

\small{file="https://raw.githubusercontent.com/distanceworkshops/online-course/gh-pages/exercisepdfs/Ch7/datasets/amakihi.csv"}

Check that it has been imported correctly.

```{r, echo=TRUE, eval=answer}
head(amakihi, n=3)
```

These data consist of eight columns:

- Study.Area - name of the study area
- Region.Label - survey dates which are used as 'strata'
- Sample.Label - point transect identifier
- Effort - survey effort (1 for all points because they are visited a single time)
- distance - radial distance of detection from observer
- OBS - initials of the observer
- MAS - minutes after sunrise
- HAS - hour after sunrise

Note there is no `Area` field in the dataset.  Detection functions can be fitted to the data, but bird density and abundance cannot be estimated. The latter three columns are the covariates to be considered for possible inclusion into the detection function. 

There a couple of records with missing distances and so can be deleted with the following command:

```{r, echo=TRUE}
amakihi <- amakihi[!is.na(amakihi$distance), ]
```

In this command, 

- records in `amakihi` are selected using the square brackets `[]`
- `amakihi` is a data frame and so selection can be performed on either rows or columns i.e. `[rows, columns]`. In this case, the selection is performed on the rows (because the selection criteria is before the comma) and all columns will be retained
- the rows selected as those where the distances (stored in `amakihi$distance`) are not missing. The function `is.na` selects elements that are missing; the symbol `!` means 'not', and so `!is.na` selects elements that are not missing.

# Exploratory data analysis
```{marginfigure}
__Question:__  Examine the distribution of radial distances of the point transect data of the amakihi.

Basic syntax will be 

`ggplot(amakihi, aes(x=distance)) + geom_histogram(binwidth=1)`

Examine the full dataset, then truncate the data to 82.5m.
```

```{r, echo=FALSE, eval=answer, fig.width=2, fig.height=2, fig.margin=TRUE,results='hold',fig.cap="Two levels of detail examining distribution of detection distances.", results='hide', warning=FALSE}
library(ggplot2)
ggplot(amakihi, aes(x=distance)) + geom_histogram(binwidth=1, colour="black", fill="white")
#ggplot(amakihi[amakihi$distance<82.5,], aes(x=distance)) + 
  geom_histogram(binwidth=2.5, colour="black", fill="white")
ggplot(amakihi[amakihi$distance<82.5,], aes(x=distance)) + 
  geom_histogram(binwidth=8.25, colour="black", fill="white")
```

It is important to gain an understanding of the data prior to fitting detection functions. With this in mind, preliminary analysis of distance sampling data involves:

- assessing the shape of the collected data,
- considering the level of truncation of distances, and
- exploring patterns in potential covariates. 

We begin by assessing the distribution of distances by plotting histograms with different number of bins and different truncation. 

The components of the boxplot are:

- the thick black line indicates the median
- the lower limit of the box is the first quartile (25th percentile) and the upper limit is the third quartile (75th percentile)
- the height of the box is the interquartile range (75th - 25th quartiles)
- the whiskers extend to the most extreme points which are no more than 1.5 times the interquartile range.
- dots indicate 'outliers' if there are any, i.e. points beyond the range of the whiskers.

This format is probably not as useful as a histogram in a distance sampling context but boxplots can be useful to compare the distances for discrete groups in the data. Here we use boxplots to display the distribution of distances recorded by each observer and for each hour after sunrise. Note how the `~` symbol is used to define the groups.

Boxplots of distances by observer: 

```{r, echo=TRUE, eval=answer, fig.width=3.5, fig.height=3.5, fig.margin=TRUE, fig.cap="Detection distance distribution by Observer."}
ggplot(droplevels(amakihi), aes(x=OBS, y=distance)) +
  geom_boxplot() + labs(x="Observer initials", y="Radial distance (m)")
```

Boxplot of distances for each hour after sunrise:

```{r, echo=TRUE, eval=answer, fig.width=3.5, fig.height=3.5, fig.margin=TRUE, fig.cap="Detection distance distribution by hours after sunrise."}
ggplot(amakihi, aes(x=factor(HAS), y=distance)) +
  geom_boxplot() + labs(x="Hours after sunrise", y="Radial distance (m)")
```

For minutes after sunrise (a continuous variable), we create a scatterplot of MAS (on the $x$-axis) against distances (on the $y$-axis).

```{r, echo=TRUE, eval=answer, fig.width=3.5, fig.height=3.5, fig.margin=TRUE, fig.cap="Detection distance distribution by minutes after sunrise."}
ggplot(amakihi, aes(x=MAS, y=distance)) + 
  geom_point(alpha=0.25, size=1.6) + 
  labs(x="Hours after sunrise", y="Radial distance (m)")
```

# Colinearity

Estimating the parameters of a detection function when covariates are involved is complex.  You will recall from multiple linear regression that problems in estimation arise when two covariates in the model are highly correlated.  In the exploratory data analysis, it is useful to look for colinearity in potential covariates.

```{marginfigure}
__Question:__  Compute the correlation of minutes after sunrise and hours after sunrise using the `cor()` function.
```

```{r, echo=FALSE, eval=answer, fig.width=2, fig.height=2, fig.margin=TRUE,results='hold',fig.cap="Diagnostics for the presence of colinearity between HAS and MAS.", warning=FALSE}
correlat <- cor(amakihi$MAS, amakihi$HAS, use = "complete.obs")
ggplot(amakihi, aes(x=MAS, y=HAS)) + geom_point(shape=19, alpha=0.25, size=0.5) +
  labs(x="Minutes after sunrise", y="Hours after sunrise") +
  annotate("text", x=70, y=4, label=paste0("Corr=",round(correlat,3)))

```

To alleviate the potential colinearity difficulty, hours after sunrise could be transformed to a discrete, rather than a continuous variable.

# Adjusting the raw covariates

We would like to treat OBS and HAS as factor variables as in the original analysis; OBS is, by default, treated as a factor variable because it consists of characters rather than numbers. `HAS`, on the other hand, consists of numbers and so by default would be treated as a continuous variable (i.e. non-factor). That is fine if we want the effect of `HAS` to be monotonic (i.e. detectability either increases or decreases as a function of `HAS`). If we want `HAS` to have a non-linear effect on detectability, indicate it is a factor:.  

```{r obsfact, echo=TRUE}
amakihi$HAS <- factor(amakihi$HAS)
```


One final adjustment, and more subtle, is a transformation of the continuous covariate, `MAS`. We are entertaining three possible covariates in our detection function: `OBS`, `HAS` and `MAS`. The first two variables, `OBS` and `HAS`, are both factor variables, and so, essentially, we can think of them as taking on values between 1 and 3 in the case of `OBS`, and 1 to 6 in the case of `HAS`.  However, `MAS` can take on values from -18 (detections before sunrise) to >300 and the disparity in scales of measure between `MAS` and the other candidate covariates can lead to difficulties in the performance of the optimizer fitting the detection functions in R. The solution to the difficulty is to scale `MAS` such that it is on a scale (approx. 1 to 5) comparable with the other covariates.

Scaling the `MAS` measurements accomplishes the desired compaction in the range of the `MAS` covariate without changing the shape of the distribution of `MAS` values. 

```{r MASscale, echo=TRUE}
amakihi$MAS <- scale(amakihi$MAS)
```

Check what this command has done by looking at the range of the adjusted MAS:

```{r, echo=TRUE, eval=answer}
range(amakihi$MAS)
```

# Candidate models

With three potential covariates, there are 8 possible combinations for including them in the detection function:

- No covariates
- OBS
- HAS
- MAS
- OBS + HAS
- OBS + MAS
- HAS + MAS
- OBS + HAS + MAS

Even without considering covariates there are a number of possible key function/adjustment term combinations and if all key function/covariate combinations are considered the number of potential models is large. Note that covariates are not allowed if a uniform key function is chosen and if covariate terms are included, adjustment terms are not allowed. Even with these restrictions, it is not best practice to take a scatter gun approach to detection function model fitting. \citet{buckland_distance_2015} considered 13 combinations of key function/covariates. Here, we look at a subset of these as an illustration of how to incorporate covariates in the detection function.

If it is not already loaded, load the `Distance` package \citep{Miller2019}. 

> We are not using the `convert_units` argument in calls to `ds` because we are not reporting density or abundance estimates.  If we were reporting those estimates, we would be including a unit conversion in our calls to `ds`.

Fit a hazard rate model with no covariates or adjustment terms. By default, line transects are assumed and because our data are point transects, the argument `transect="point"` is specified:

```{r, echo=TRUE, eval=answer, message=FALSE, warning=FALSE}
library(Distance)
conv <- convert_units("meter", NULL, "square kilometer")
hr.model0 <- ds(amakihi, transect="point", key="hr", truncation=82.5,
                adjustment=NULL, convert_units = conv)
```

The fitted model can be investigated using the `summary` function. Make a note of the AIC for this model. 

```{r, echo=TRUE, eval=answer, message=FALSE, warning=FALSE}
summary(hr.model0)
```

```{marginfigure}
__Question:__ Fit a hazard rate model with OBS as a covariate in the detection function and make a note of the AIC. Has the AIC reduced by including a covariate? 
```

```{r, echo=TRUE, eval=answer, message=FALSE, warning=FALSE}
hr.obs <- ds(amakihi, transect="point", key="hr", formula=~OBS, 
             truncation=82.5, convert_units = conv)
print(hr.obs$ddf$criterion)
```
```{marginfigure, echo=answer}
__Answer:__ Yes, AIC of the model with observer covariate is 30 AIC units smaller than the model without this covariate.
```

Fit a hazard rate model with OBS and HAS in the detection function:

```{r, echo=TRUE, eval=answer, message=FALSE, warning=FALSE}
hr.obshas <- ds(amakihi, transect="point", key="hr", 
                formula=~OBS+HAS, truncation=82.5, convert_units = conv)
print(hr.obshas$ddf$criterion)
```

```{marginfigure}
__Question:__ Fit the other candidate models including covariates shown in the list above and decide which model is best in terms of AIC.
```

```{marginfigure, echo=answer}
__Answer__: The model with both observer and hours after sunrise had an AIC score 5 AIC units larger than the model with observer alone; suggesting the inclusion of HAS along with OBS is not a better model than the model with OBS alone.
```

```{r fitmodels, echo=TRUE, eval=answer, message=FALSE, warning=FALSE}
hr.has <- ds(amakihi, transect="point", key="hr", 
             formula=~HAS, truncation=82.5, convert_units = conv)
hr.mas <- ds(amakihi, transect="point", key="hr", 
             formula=~MAS, truncation=82.5, convert_units = conv)
hr.obsmas <- ds(amakihi, transect="point", key="hr", 
                formula=~OBS+MAS, truncation=82.5, convert_units = conv)
hr.hasmas <- ds(amakihi, transect="point", key="hr", 
                formula=~HAS+MAS, truncation=82.5, convert_units = conv)
hr.hasmasobs <- ds(amakihi, transect="point", key="hr",
                   formula=~HAS+MAS+OBS, truncation=82.5, convert_units = conv)
```

A useful function for summarising a candidate model set is `summarize_ds_models`.  The arguments to the function is an enumeration of the candidate model objects.  

```{r summarisemods, echo=TRUE, eval=FALSE}
summarize_ds_models(hr.model0, hr.obs, hr.has, hr.mas,
                    hr.obshas, hr.obsmas, hr.hasmas, hr.hasmasobs)
```

```{r summkable, echo=FALSE, eval=answer}
knitr::kable(summarize_ds_models(hr.model0, hr.obs, hr.has, hr.mas,
                    hr.obshas, hr.obsmas, hr.hasmas, hr.hasmasobs), 
             caption="Candidate model set for Hawaii amakihi covariate analysis.",
             row.names = FALSE, digits=3)
```

# Plotting the detection functions

The detection functions can be investigated using the `plot` function as shown below. A few different plotting options are illustrated.

```{r detnone, echo=TRUE, eval=FALSE}
# Plot simple model
plot(hr.model0, nc=20, main="No covariates", pch=20, pdf=TRUE)
```

```{r detobs, echo=TRUE, eval=FALSE}
# Plot model with OBS
plot(hr.obs, nc=10, main="Model with OBS covariate", 
     cex=0.5, pdf=TRUE, showpoints=FALSE)
add_df_covar_line(hr.obs, data.frame(OBS=c("SGF", "TJS","TKP")),
                  col=c("blue","green","red"), lty=2, pdf=TRUE)
legend("topright", legend=c("SGF", "TJS","TKP"),
       col=c("blue","green","red"), lwd=2)
```

What does the detection function look like for your selected model? 

```{r, echo=FALSE, eval=answer, ref.label="detnone", fig.width=3.5, fig.height=3.5, fig.cap="Detection function fit for model without covariates.", fig.margin=TRUE}
```
```{r, echo=FALSE, eval=answer, ref.label="detobs", fig.width=3.5, fig.height=3.5, fig.cap="Detection function fit for model with observer covariate.", fig.margin=TRUE}
```
```{r, echo=TRUE, eval=answer, fig.width=3.5, fig.height=3.5, message=FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Detection function fit for model with lowest AIC."}
# Fit best model
hr.best <- ds(amakihi, transect="point", key="hr", 
              truncation=82.5, quiet=TRUE,formula=~OBS+MAS)
# Plot model with OBS and MAS
plot(hr.best, nc=10, main="Model with OBS and MAS covariates", 
     pch=".", pdf=TRUE)
```

To see more sophisticated examples of plotting the detection function for the selected model, see the code accompanying \citep{buckland_distance_2015} [Hawaiian Amakihi case study](https://synergy.st-andrews.ac.uk/ds-manda/#hawaiian-amakihi-case-study).
