```{r setup, include=FALSE, message=FALSE}
knitr::opts_chunk$set(echo=TRUE, message=F, warning=F)
```

The purpose of this exercise is to analyse point transect survey data: it can sometimes be more difficult than line transect data. In the first problem, the data are simulated and so the true density is known. In the second problem, two different data collection methods were used to survey song birds. 

# Objectives

The aim of this practical are to:

1. Practice fitting detection functions to point transect survey data. 
2. Use data from the `dsdata` package. 

# Simulated survey data

Simulated point transect data from 30 points are given in the text file 'IntroDS_5.1.csv'. These data were generated from a half-normal detection function and the true detection function was 79.8 animals/hectare \footnote{1 hectare=0.01km$^2$}. The radial distances were recorded in metres. 

Load the data, check it is OK and plot the radial distances. 

```{r, echo=T, eval=F}
library(Distance)
# Read in data
simptdatafile <- system.file("extdata", "IntroDS_5.1.csv", package = "dsdata")
ptdat <- read.csv(file=simptdatafile, header=TRUE)
conversion.factor <- convert_units("meter", NULL, "hectare")
# What does the data look like
head(ptdat)
hist(ptdat$distance)
```

To fit a point transect detection function, the argument `transect="point"` needs to be specified:

```{r, echo=T, eval=F}
# Fit half-normal point transect detection function, 
ptdat.hn <- ds(data=ptdat, transect="point", key="hn", 
               convert.units=conversion.factor)
```

The `convert.units` argument gives the estimated density in animals per hectare.

The detection function can be plotted as for line transects:

```{r, eval=F}
# Plot detection function
plot(ptdat.hn)
```

To plot the probability density function (pdf), an additional argument is required in the `plot` function:

```{r, eval=F}
# Plot probability density function
plot(ptdat.hn, pdf=TRUE)
```

Experiment with keys other than the half normal (i.e. hazard rate and uniform) to assess whether these data can be satisfactorily analysed using the wrong model: 

+ determine a suitable truncation distance, and 
+ for each key function decide whether any adjustments are needed.

How do the bias and precision compare between models? 

# Wren data (Optional)

A point transect survey of songbirds was conducted at Montrave, Fife, Scotland, in 2004 (Buckland 2006) and for this exercise, the data on winter wrens is used. Several different methods of data collection were used and for this exercise, two point transect methods are used:

1. standard five-minute counts,
2. the 'snapshot' method.

For each method the same 32 point transects were used in 33.2 ha of parkland (Fig. 1) and each point transect was visited twice. Detection distances (recorded in metres) were measured with the aid of a rangefinder.

```{r montrave, echo=FALSE, fig.cap="The study site: the dotted line is a small stream, the short dashed lines are tracks and the thick dashed line is a main road. The 32 points, shown by crosses, are laid out on a systematic grid with 100m separation. The diagonal lines were used for a line transect survey.", out.width = '50%'}
knitr::include_graphics("https://workshops.distancesampling.org/standrews-2019/intro/practicals/figures/Prac_5_Figure_1.png")
```


Data for these methods are available in the `dsdata` package. This is convenient because each data set can be accessed as follows:

```{r, echo=T, eval=F}
# Load library
library(dsdata)
# Extract wren data for method 1
data(wren1)
conversion.factor <- convert_units("meter", NULL, "hectare")
```

You will see that there is an object called `wren1` in your `R` workspace. There is also a data object available in `dsdata` for method 2 (i.e. `wren2`) which can be loaded in the same way. Have a look at the `wren1` data with

```{r, echo=T, eval=F}
head(wren1, n=3)
```

Note the `Effort` field is `2` meaning each point transect was visited twice.  The same applies for `wren2`.

What to do:

1. Select a simple model for exploratory data analysis. Experiment with different truncation distances, $w$, and select a suitable value for each method. Are there any potential problems with any of the data sets? 
2. Try other models and model options. Use plots, AIC values and goodness-of-fit test statistics to determine an adequate model. 
3. Record your estimates of density and corresponding confidence interval for each method. Compare your answers with those of others in the workshop. The conversion units to obtain density in birds per hectare is `convert.units=0.01`. 

# References

Buckland, ST 2006. Point-transect surveys for songbirds: robust methodologies. The Auk 123: 345–57. https://doi.org/10.1642/0004-8038(2006)123[345:PSFSRM]2.0.CO;2.


\newpage

***
**Solution 5. Point transect exercises**

***

# Simulated data

The code for importing and checking these data and fitting various models is shown below. 

```{r, fig.width=4, fig.height=4, message=FALSE}
library(Distance)
# Read in data
simptdatafile <- system.file("extdata", "IntroDS_5.1.csv", package = "dsdata")
ptdat <- read.csv(file=simptdatafile, header=TRUE)

# What does data look like
head(ptdat, n=3)
conversion.factor <- convert_units("meter", NULL, "hectare")
# Fit half-normal detection function, no truncation
ptdat.hn <- ds(data=ptdat, transect="point", key="hn",
               convert.units=conversion.factor)
plot(ptdat.hn, pdf=TRUE)

# Try truncation of 20m based on preliminary fit

# Half normal, no adjustments
ptdat.hn.t20m <- ds(data=ptdat, transect="point", key="hn", truncation=20,
                    convert.units=conversion.factor)
# Hazard rate, no adjustments
ptdat.hr.t20m <- ds(data=ptdat, transect="point", key="hr", truncation=20,
                    convert.units=conversion.factor)
# Uniform, cosine adjustments
ptdat.uf.cos.t20m <- ds(data=ptdat, transect="point", key="unif", 
                        adjustment="cos", truncation=20,
                        convert.units=conversion.factor)
```


```{r, echo=F}
# Harvest results
pt.tab <- data.frame(DetectionFunction=c("Half-nomal","Half-nomal","Hazard rate","Uniform"), Adjustments=c("None","None","None","Cosine"), Truncation=c(34.2,20,20,20), AIC=rep(NA,4), Density=rep(NA,4), D.CV=rep(NA,4), Lower.CI=rep(NA,4), Upper.CI=rep(NA,4))

get.results.f <- function(fit.model) {   return(c(AIC=summary(fit.model$ddf)$aic,
         D=fit.model$dht$individuals$D$Estimate,
         D.CV=fit.model$dht$individuals$D$cv,
         lCL=fit.model$dht$individuals$D$lcl,
         uCL=fit.model$dht$individuals$D$ucl))
}

pt.tab[1,4:8] <- get.results.f(ptdat.hn)
pt.tab[2,4:8] <- get.results.f(ptdat.hn.t20m)
pt.tab[3,4:8] <- get.results.f(ptdat.hr.t20m)
pt.tab[4,4:8] <- get.results.f(ptdat.uf.cos.t20m)

# Print results
pander::pander(pt.tab, caption="Results from simulated point transect data.")
```

```{r, echo=FALSE, eval=FALSE}
# Plot detection functions
par(mfrow=c(2,2))
plot(ptdat.hn, main="Half normal, no truncation")
plot(ptdat.hn.t20m, main="Half normal, truncation 20m")
plot(ptdat.hr.t20m, main="Hazard rate, truncation 20m")
plot(ptdat.uf.cos.t20m, main="Uniform with cosine, truncation 20m")
```

```{r}
# Plot probability density functions
par(mfrow=c(2,2))
plot(ptdat.hn, main="Half normal, no truncation", pdf=TRUE)
plot(ptdat.hn.t20m, main="Half normal, truncation 20m", pdf=TRUE)
plot(ptdat.hr.t20m, main="Hazard rate, truncation 20m", pdf=TRUE)
plot(ptdat.uf.cos.t20m, main="Uniform with cosine, truncation 20m", pdf=TRUE)
```

We see a fair degree of variability between analyses - reliable analysis of point transect data is more difficult than for line transect data. We see greater loss in precision from truncating data relative to line transect sampling, but if we don't truncate data, different models can give widely differing estimates.

# Wren data (Optional)

In the code provided below, each dataset is loaded and then detection functions that we selected are fitted. 

```{r, echo=T, eval=T}
# Load library
library(dsdata)
# Wren data for Method 1 (5 min counts)
data(wren1)
# Wren data for Method 2 (snapshot)
data(wren2)
conversion.factor <- convert_units("meter", NULL, "hectare")
# Method 1
wren1.uf.cos.t110 <- ds(data=wren1, key="unif", adjustment="cos", 
                        transect="point", truncation=110, 
                        convert.units=conversion.factor)
# Method 2
wren2.hr.cos.t110 <- ds(data=wren2, key="hr", adjustment=NULL, 
                        transect="point", truncation=110, 
                        convert.units=conversion.factor)
```

```{r, echo=F}
# Harvest results
n <- 2
wren.tab <- data.frame(Method=1:n, Density=rep(NA,n), 
                       Lower.CI=rep(NA,n), Upper.CI=rep(NA,n))

get.results.f <- function(fit.model) { return(c(D=fit.model$dht$individuals$D$Estimate,
         lCL=fit.model$dht$individuals$D$lcl,
         uCL=fit.model$dht$individuals$D$ucl))
}

wren.tab[1,2:4] <- get.results.f(wren1.uf.cos.t110)
wren.tab[2,2:4] <- get.results.f(wren2.hr.cos.t110)

# Print results
pander::pander(wren.tab, caption="Winter wren density estimates from 5 minute counts and snapshot moment.")
```

```{r, fig.height=4}
# Plot detection functions
par(mfrow=c(1,2))
plot(wren1.uf.cos.t110, main="5 minute count")
plot(wren2.hr.cos.t110, main="Snapshot moment")
```

As the detection distance histograms indicate, winter wren showed evidence of observer avoidance, more than other species in the Montrave study.  We show the detection function graph rather than the PDF to emphasise the evasive movement aspect of the data.  If you conduct the goodness of fit test, using `gof_ds()`, you will find that the models still suitably fit the data.

