```{r, echo=FALSE, warning=FALSE, message=FALSE}
library(knitr)
answer <- TRUE
```


# More complex analyses
```{marginfigure, echo=answer}
\huge(Solutions)
```
```{marginfigure}
\includegraphics[]{images/robin-300x256.jpg}
European robin _(Erithacus rubecula)_; one of the species in the Montrave study \citep{Buckland_2006}.

\includegraphics[]{images/montrave-map.jpg}
Aerial view of Montrave study area.  White diagonal lines represent transects walked for data analysed here.
```

This practical is based on the Montrave songbird case study in \citet[Section 5.2.2.3]{buckland_distance_2015}, with computer code under [Montrave songbird case study](https://synergy.st-andrews.ac.uk/ds-manda/#montrave-songbird-case-study). Both point and line transect surveys were conducted and here we use the data from the __line transect__ data, although the issues (and solutions) will be similar. 

These data are provided in a 'flat file' format (i.e. it contains all the necessary columns to estimate a detection function, density and abundance). While both formats are equally valid, the 'flat file' approach has a particular idiosyncrasy which we exploit here to introduce more functions and data manipulation.  

Several species of birds were identified but not all species were detected on all transects. If a simple data selection is performed to select records for a particular species, then not all of the transects will be included in the resulting data (because that species may not have been seen). This doesn't matter if we are only interested in fitting detection functions, but will matter if we wish to estimate density and abundance because the effort will be too low since some of the transects are missing. To correct for this, some data frame manipulation is required. There is generally more than one way to do something in R \citep{Miller2019} - for an alternative way see the computer code 'Montrave song bird case study' associated with \citet{buckland_distance_2015}, as well as Section \ref{functions} below. 

# Objectives of the practical

1. Data frame selection and manipulation
1. Extracting estimates from `dht` object
1. Customising detection function plots
1. Improve re-usability of code with functions

# Importing the data

The data is in a 'flat file' format and contains the following columns:

- Region.Label - name of study
- Area - size of study region (hectares)
- repeats - number of visits to transect
- Sample.Label - line transect identifier
- Effort - length of transect (km)
- distance - perpendicular distance (m)
- species - species of bird (c=chaffinch, g=great tit, r=robin and w=wren)
- visit - on which visit bird was detected.

Use the following command to import the data from the website associated with \citet{buckland_distance_2015} and then use the `head` command to examine it. 

```{r, echo=TRUE}
birds <- read.csv(file="https://raw.githubusercontent.com/distanceworkshops/online-course/gh-pages/exercisepdfs/Ch7/datasets/montrave-line.csv")
```

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

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


```{marginfigure}
__Question:__ Explore the data. How many transects are there? 
```
  
```{r, echo=TRUE, eval=answer}
length(unique(birds$Sample.Label))
```

For now, save the transect labels to a new object as we will use them later on:

```{r, echo=TRUE}
tran.lab <- unique(birds$Sample.Label)
```

The `table` command is a quick way to determine how many detections there are of each species: 

```{r, echo=TRUE, eval=answer}
table(birds$species)
```

As a hint of things to come, create a two-way table showing the number of detections by transect and by species.  If there are zeroes in this table, it will create a challenge.

```{r twoway, echo=TRUE, eval=FALSE}
with(birds, table(species, Sample.Label))
```

```{r kabletwoway, echo=FALSE, eval=TRUE}
knitr::kable(with(birds, table(species, Sample.Label)))
```

Each of the line transects was visited twice which is not taken into account at present. However, it is straightforward to do so: 

```{r, echo=TRUE}
birds$Effort <- birds$Effort * birds$repeats
```

# Getting the effort right for the robin data {#manip}

For the purposes of this practical, we are interested in estimating the density of robins and so we select only these records: 

```{r, echo=TRUE}
robins <- birds[birds$species=="r", ]
```

```{marginfigure}
__Question:__ On how many transects were robins detected?
```

```{r, echo=TRUE, eval=answer}
length(unique(robins$Sample.Label))
```

If we were to use the `robins` data as it is at present to estimate density, then density would be __incorrect__ because the search effort associated with three transects is missing. Resolution of this problem is to use additional arguments to the `ds()` function.  These arguments carry information about the survey design that cannot be deduced from the flatfile `montrave-lines.csv`.

## Additional arguments

Rather than proceeding with the `robins` data set, which will produce erroneous results, let's back up and work with the entire multiple species data set.

We will create two data frames to contain information about Buckland's survey.  These data frames will be passed to `ds()` via the arguments `region_table` and `sample_table`.  The contents of the data frames are described in the help for the `ds()` function

> - region_table	data.frame with two columns:
>   - Region.Label label for the region
>   - Area area of the region
>   - region_table has one row for each stratum. If there is no stratification then region_table has one entry with Area corresponding to the total survey area. If Area is omitted density estimates only are produced.
> - sample_table	data.frame mapping the regions to the samples (i.e. transects). There are three columns:
>   - Sample.Label label for the sample
>   - Region.Label label for the region that the sample belongs to.
>   - Effort the effort expended in that sample (e.g. transect length).

The code in the next chunk performs the magic we need.
```{r buildtables}
bird.regiontab <- data.frame(Region.Label=as.factor(c(1,2)), Area=c(33.2,33.2))
bird.sampletab <- data.frame(Region.Label=as.factor(rep(c(1,2), each=19)),
                            Sample.Label=rep(1:19, times=2),
                            Effort=c(0.208, 0.401, 0.401, 0.299, 0.350,
                                     0.401, 0.393, 0.405, 0.385, 0.204,
                                     0.039, 0.047, 0.204, 0.271, 0.236,
                                     0.189, 0.177, 0.200, 0.020))
```

## Dissecting the code

The `region_table` data frame is the most simple.  We have indicated there are two strata (the two visits; estimation is not required for each visit, but it does represent how the sampling was carried out).  Because the strata simply repeat visits to the same transects, the size of the strata is identical and equal to the size of the study area.

The `sample_table` data frame needs to know the "names" of each stratum (here simply "1" and "2") and the names of the transects (numbers 1-19) and the stratum to which each transect belongs.  Finally the length of each transect is required.  When completed, the data frame contains 38 rows (the 19 transects duplicated once).  The 19 transects belonging to stratum 1 (first visit) and again the same 19 transects belonging to stratum 2 (second visit).

When these data frames are provided to `ds()` the integrity of the survey design will be preserved regardless of the data filtering that might take place.

# Analysis

Before we fit any models, have a quick look at the histogram of distances:

```{r, echo=TRUE, eval=answer, message=FALSE, fig.height=3.5, fig.width=3.5, fig.margin=TRUE, fig.cap="Perpendicular distances of robins in Montrave study."}
hist(robins$distance, breaks=20)
```

Consistent with \citet{buckland_distance_2015}, three detection functions are fitted using the `ds()` function in the R package `Distance` \citep{Miller2019} using the truncation distance used by Prof Buckland:

```{r, echo=TRUE, message=FALSE, warning=FALSE}
library(Distance)
robins$Region.Label <- robins$visit
robin.hn.herm <- ds(robins, truncation=95, transect="line",
                    region_table = bird.regiontab, sample_table = bird.sampletab,
                    key="hn", adjustment="herm", convert_units=0.1)
robin.uni.cos <- ds(robins, truncation=95, transect="line", 
                    region_table = bird.regiontab, sample_table = bird.sampletab,
                    key="unif", adjustment="cos", convert_units=0.1)
robin.haz.simp <- ds(robins, truncation=95, transect="line", 
                     region_table = bird.regiontab, sample_table = bird.sampletab,
                     key="hr", adjustment="poly", convert_units=0.1)
```

```{marginfigure}
__Question:__ What is the preferred model for the robin data?
```

```{r modsel, echo=TRUE, eval=FALSE}
summarize_ds_models(robin.hn.herm, robin.uni.cos, robin.haz.simp)
```

```{r modselprint, echo=FALSE, eval=answer}
kable(summarize_ds_models(robin.hn.herm, robin.uni.cos, robin.haz.simp)[,c(1,4,5,6,7)], 
      digits = c(rep(x = 3, 4)), row.names = FALSE,
      caption="Model selection for robin data from Montrave line transect survey.")
```

```{marginfigure, echo=answer}
__Note:__  All three detection function fit the data (based upon the C-vM test of exact distances).  The estimated detection probability is very similar for all models, and the $\Delta$AIC values of all models is $<$ 1.  Hence all models will give very similar estimates of density.
```

# Examining the `dht` object

The fitted model object (e.g. `robin.uni.cos`) is made up of two parts; the detection function in the `ddf` part and the estimates in the `dht` part. In this section, we look at the `dht` part. 

To list the elements that are contained in `dht`, use the `names` function:

```{r, echo=TRUE, eval=answer}
names(robin.uni.cos$dht)
```

Detections were of individual birds and so group size was not included in these data - if it had been included (in a column called `size`), then as well as `individuals` there would have been elements `clusters` and `Expected.S`. 

The estimates stored in the `individuals` object can be listed in a similar manner:

```{r, echo=TRUE, eval=answer}
names(robin.uni.cos$dht$individuals)
```

To collect together the density estimates (and estimates of precision) from all the fitted models, we can use the following command:

```{r, echo=TRUE}
model.results <- rbind(robin.uni.cos$dht$individuals$D[3,],
                       robin.haz.simp$dht$individuals$D[3,], 
                       robin.hn.herm$dht$individuals$D[3,])
```

```{marginfigure}
__Question:__ Examine the three sets of density estimates to see if the previous suggestion (that the density estimates are similar) is confirmed.
```

```{r, echo=answer, eval=FALSE}
model.results
```

```{r dentable, echo=FALSE, eval=answer}
knitr::kable(model.results, row.names = FALSE, digits=4,
      caption = "Density estimates for Montrave robins under three fitted detection functions.")
```

# Goodness of fit

Here we look at goodness of fit test with unequal bin intervals and just consider one of the fitted models. First we specify the required bin intervals.

```{r, echo=TRUE, message=FALSE}
robin.brks <- c(0, 12.5, 22.5, 32.5, 42.5, 52.5, 62.5, 77.5, 95.0)
```

Perform the tests using both exact distance data for the Cramer-von Mises test and specified breakpoints for $\chi^2$ test for the uniform-cosine model that had the (slightly) smallest AIC score.

```{r, echo=TRUE, eval=answer, fig.width=3.5, fig.height=3.5, fig.margin=TRUE}
testout <- gof_ds(robin.uni.cos, breaks=robin.brks, 
                  chisq = TRUE, main="QQ plot unif-cos for robins")
nice <- data.frame(statistic=c(testout$chisquare$chi1$chisq,
                               testout$dsgof$CvM$W),
                   pvalue=c(testout$chisquare$chi1$p,
                            testout$dsgof$CvM$p))
row.names(nice) <- c("Chisquare", "CvM test")
knitr::kable(nice, row.names = TRUE, digits=3,
             caption="Goodness of fit tests for unicos model.")
```

```{marginfigure, echo=answer}
__Note:__ The detections fall close to the diagonal line of the qq plot, suggesting an adquate fit for the uniform cosine model.  The _p-value_ of the Cramer-von Mises test (at bottom of printout) confirms this.  Similarly the _p-value_ for the $\chi^2$ test also suggests an adequate fit.
```

# Customising the detection function plot

The `plot` function provides a basic plot of the fitted detection function overlaid onto the scaled distribution of distances:

```{r, echo=TRUE, eval=answer, fig.width=3.5, fig.height=3.5, fig.margin=TRUE}
plot(robin.uni.cos)
```

# Advanced:  modularising R code to work with multiple species {#functions}

When analysing a multi-species survey, it is likely that the invesigator will want to analyse all (or at least many) of the species encountered during the survey.  This will necessitate some repetitive calculation fitting multiple detection functions for each species.

To facilitate the repetitive nature of such analyses, it is useful to take advantage of the programmatic nature of the R language to create _functions_ that can be called repeatedly with arguments to accommodate changing circumstances.  The code below demonstrates such a modular approach with `fit.hn.uni.haz()` defined to aide in the repeated analyses.

`fit.hn.uni.haz()` fits three candidate models to a dataset provided as the first argument.  The second argument is the truncation distance.  The final argument determines whether the `summarize_ds_models()` table is printed.

```{r fitfn, echo=TRUE, eval=TRUE, tidy=FALSE}
fit.hn.uni.haz <- function(data, trunc, print=TRUE) {
  #  Purpose:  fit three key functions to transect data, 
  #            perform model selection and
  #            print model selection table
  #  Input: data to analyse, truncation distance, print flag
  #  Output: fitted model object (class `dsmodel`)
  #  Rexstad August 2018
  hn.herm <- ds(data, trun=trunc, key="hn", adj="herm",
                convert_units = .1)
  uni.cos <- ds(data, trun=trunc, key="unif", adj="cos", 
                convert_units = .1)
  haz.simp <- ds(data, trun=trunc, key="hr", adj="poly", 
                 convert_units =.1)
  mods <- summarize_ds_models(hn.herm, uni.cos, haz.simp, output="plain")
  if(print) print(knitr::kable(mods))
  names(mods) <- c("mod","key","form","fit","pa","sepa","daic")
  if(mods[1,1]=="hn.herm") {
    result <- hn.herm
  } else {
    if(mods[1,1]=="uni.cos") {
    result <- uni.cos
  } else {
    result <- haz.simp
  }
  }  
  return(result)
}
```

The function is used in the calling code below.  Note the `for` loop that iterates through three of the four species detected in the Montrave survey (great tit not analysed because there were few detections).  Note finally, that even though the model with the smallest AIC for the winter wren was the hazard rate model, this was not the model used for inference by \citet{Buckland_2006} as the shape (perfect detectability to 70m) was determined to be biologically implausible.  This underscores one of the challenges associated with "assembly line" analysis of multiple species; biological insight needs to be maintained in the model selection process.

```{r usefunct, echo=TRUE, eval=answer, fig.width=3.5, fig.height=3.5, fig.margin=TRUE, results='hold', message=FALSE, tidy=FALSE, warning=FALSE}
for(thisspecies in c("r", "c", "w")) {
  best.model <- fit.hn.uni.haz(birds[birds$species==thisspecies, ],
                               trunc= 100, print=FALSE)
  plot(best.model, 
       main=paste("Montrave lines, species ", thisspecies,
       "\nD-hat=", round(best.model$dht$individuals$D$Estimate,4), 
       "SE=", round(best.model$dht$individuals$D$se, 4),
       "\n",best.model$ddf$name.message))
}
```

