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

For this exercise, we use data from a survey of Antarctic minke whales. The study region was divided into two strata (North and South) the two strata were surveyed by different vessels at the same time. The whales tend to be found in high densities against the ice edge, where they feed, and so the density in the southern stratum is typically higher than in the northern stratum (Figure 1). This is the primary reason for using a stratified survey design. It is also the reason for covering the southern stratum more intensely: in the southern stratum the transect length per unit area is more than 2.5 times that of the northern stratum.

![](https://workshops.distancesampling.org/standrews-2019/intro/practicals/figures/Prac_7_Figure_1.png){width=48%}![An example of the sort of survey design used (left) and a typical whale density gradient (right). The irregular bottom border is the ice-edge. The 'stepped' black line defines the boundary between the strata; dotted lines are transects and dots are detections.](https://workshops.distancesampling.org/standrews-2019/intro/practicals/figures/Prac_7_Figure_2.png){width=42%}

Figure 1. An example of the sort of survey design used and a typical whale density gradient. The irregular bottom border is the ice-edge. The 'stepped' black line defines the boundary between the strata; dotted lines are transects and dots are detections.

# Objectives

The objectives of this exercise are to:

1. Create subsets of the data
2. Decide whether to fit separate detection functions or a pooled detection function 
3. Specify different stratification options using the `dht2` function. 

# Getting started

Begin by reading in the data and in the commands below the data is renamed to avoid overwriting the original dataset. Distances are in kilometers and a truncation distance of 1.5km is specified and used in the following detection function fitting. Perpendicular distances, transect lengths and study area size are all measured in kilometers; hence `convert.units` argument to `ds` is 1 and has been omitted.  To keep things simple, a hazard rate detection function with no adjustments is used for all detection functions. 

```{r, echo=T, eval=F}
library(Distance)
data(minke)
# Rename data
whale <- minke
# Check data
head(whale)
# Specify truncation distance
whale.trunc <- 1.5
```

You will see that these data contain a column called 'Region.Label': this contains values 'North' or 'South'. 

# Full geographical stratification

First, we want to fit encounter rate and detection function separately in each strata. This is easily performed by splitting the data by region and using  `ds` on each subset. The commands below do this for the southern region (note, there are alternative ways to select a subset of data). 

```{r, echo=T, eval=F}
# Create dataset for South 
whale.S <- whale[whale$Region.Label=="South", ]
# Fit df to south
whale.df.S.strat <- ds(whale.S, key="hr", adjustment=NULL, 
                       truncation=whale.trunc)
# Print results
summary(whale.df.S.strat)
```

Make a note of the AIC. Perform a similar commands to obtain estimates for the northern region. What is the total AIC? 

Also make a note of the abundance in each region. What is the total abundance in the study region? 

# Fitting a pooled detection function

We want to compare the total AIC found previously with the AIC from fitting a detection function to all data combined. This is easy to obtain:

```{r, echo=T, eval=F}
# Fit to pooled data
whale.df.all <- ds(whale, truncation=whale.trunc, key="hr", adjustment=NULL)

summary(whale.df.all)
```

Given the AIC value for the detection function from the pooled data, would you fit a separate detection function in each strata or not? 

# Stratification options using `dht2`

The command `summary(whale.df.all)` will provide the abundance estimates for each region and the total and for this simple example, this is sufficient. However, if we want to consider different stratification options, then the `dht2` function is useful.  

After fitting a detection function, the `dht2` function, allows abundance estimates to be computed over some specified regions. In the command below, the pooled detection function is used to obtain estimates in each strata and over all (like the summary function previously used). 

```{r, echo=T, eval=F}
dht2(ddf=whale.df.all, flatfile=whales, strat_formula=~Region.Label,
     stratification="geographical")
```

The arguments are:

+ `ddf` 
  - the detection function (fitted by `ds`)
+ `flatfile` 
  - the data object containing all the necessary information \footnote{Data  is referred to as being in a 'flatfile' format if it contains information on region, transects and observations. An alternative is to use a hierarchical structure and have region, transect and observation information in separate data files with links between them to ensure that  transects are mapped to the relevant region and observations to the relevant transect.}
+ `strat_formula=~Region.Label` 
  - formula (hence the `~`) giving the stratification structure
+ `stratification="geographical"` 
  - in this example, we specify that each strata (specified in `strat_formula`) represents a geographical region. 

Make a note the total abundance in the study region. 

What happens if we were to ignore the regions and treat the data as though it came from one large study region?  This can be done by changing the stratification formula, as shown below.

```{r, echo=T, eval=F}
dht2(ddf=whale.df.all, flatfile=whales, strat_formula=~1, 
     stratification="geographical")
```

Has this changed the abundance estimate? 

\newpage

***
**Solution 7. Analysis of stratified data**

***

Reading in the data:

```{r, message=FALSE}
library(Distance)
# Load data
data(minke)
# Rename data
whale <- minke
# Check data
head(whale, n=3)
# Specify truncation distance
whale.trunc <- 1.5
```

Fit detection function and encounter rate separately in each strata. 

```{r, echo=T, eval=T, message=FALSE}
## Fit to each region separately - full geographical stratification
# Create dataset for South
whale.S <- whale[whale$Region.Label=="South", ]
# Fit detection function to South data
whale.df.S.strat <- ds(whale.S, truncation=whale.trunc, key="hr", 
                       adjustment=NULL)
summary(whale.df.S.strat)
# Combine selection and df fitting for North
whale.df.N.strat <- ds(whale[whale$Region.Label=="North", ],
                       truncation=whale.trunc, key="hr", adjustment=NULL)
summary(whale.df.N.strat)
```

Next we fitted a pooled detection function. 

```{r, echo=T, eval=T, message=FALSE}
## Fit to pooled data
whale.df.all <- ds(whale, truncation=whale.trunc, key="hr", adjustment=NULL)
summary(whale.df.all)
```

```{r, echo=F, eval=T}
# Save AIC values
aic.all <- summary(whale.df.all$ddf)$aic
aic.S <- summary(whale.df.S.strat$ddf)$aic
aic.N <- summary(whale.df.N.strat$ddf)$aic
aic.SN <- aic.S + aic.N

```

The AIC value for the detection function in the South was `r format(aic.S, digits=4)` and the AIC for the North was `r format(aic.N, digits=4)`. This gives a total AIC of `r format(aic.SN, digits=4)`. The AIC value for the pooled detection function was `r format(aic.all, digits=4)`. Because `r format(aic.all, digits=4)` is greater than `r format(aic.SN, digits=4)`, estimation of separate detection functions in each stratum is preferable. 

```{r, echo=F}
# Plot detection functions
par(mfrow=c(2,2))
plot(whale.df.S.strat, main="South")
plot(whale.df.N.strat, main="North")
plot(whale.df.all, main="All")
```

# Stratification options

```{r, echo=F}
# Harvest abundance estimates
est.full <- data.frame(Label=c("North","South","Total"),Estimate=rep(NA,3))
est.full[1,2] <- whale.df.N.strat$dht$individuals$N$Estimate
est.full[2,2] <- whale.df.S.strat$dht$individuals$N$Estimate
est.full[3,2] <- est.full[1,2] + est.full[2,2]

est.dfpool <- whale.df.all$dht$individuals$N[ ,c(1,2)]
```

In the full geographical stratification, both encounter rate and detection function were estimated separately for each region (or strata). This resulted in the following abundances:

```{r, echo=F}
pander::pander(est.full, caption="Abundance estimates using full geographical stratification.")
```

Next, the distances were combined to fit a pooled detection function but encounter rate was obtained for each region. This resulted in the following abundances:

```{r, echo=F}
pander::pander(est.dfpool, caption="Abundance estimates calculating encounter rate by strata and a pooled detection function.")
```

An equivalent result could be produced using the `dht2` function, which does not require

```{r, echo=T, eval=T}
# Geographical stratification with pooled detection function 
newdht <- dht2(ddf=whale.df.all, flatfile=whale, strat_formula=~Region.Label,
               stratification="geographical")
print(newdht, report="abundance")
```

If geographic stratification were ignored, the abundance estimate of would be 18,293 whales. This estimate is substantially larger than the estimates above. The reason is that the survey design was geographically stratified with less survey effort in the north and more in the south and this is being ignored in the unstratified analysis.  

<!--
What is not included in these data are sizes of the observed minke whale groups (we didn't want to clutter the analyses with these details). However, there is also a story in geographic variation in group sizes: group sizes tend to be higher in the southern stratum than in the northern stratum. Ignoring the stratification will mean that transects from both strata are treated as if they are representative of the whole survey region: this results in a positively biased estimate for the region as a whole. The estimate of group size from an unstratified analysis would result in a positively biased estimate of group size for the North stratum and a negatively biased estimate of group size for the South stratum. When it is applied to both strata, the result is a positively biased estimate of whale abundance because the North stratum is much larger and contains roughly twice as many whales as the south stratum.

The important message from this exercise is __do not perform analyses without taking the survey design into account__!

-->