```{r echo=FALSE}
library(knitr)
```

In this example we use a data set of simulated minke whales in the Southern Ocean to examine data collected in two strata.

```{r whalesetup, comment=NA}
library(Distance)
whales <- read.csv("minke.csv")
head(whales)
```

In this dataset, the `Region.Label` takes on two values, `South` and `North` relating to the two strata in the data set.

###  Model selection: stratum-specific detection functions or not?

For full geographic stratification we make two separate calls to `ds()` behaving as if the two strata have nothing in common; which is what the full stratification analysis presumes.

```{r full-geo-strat, message=FALSE, warning=FALSE}
whale.trunc <- 1.5
whale.full.strat1 <- ds(whales[whales$Region.Label=="South",], truncation=whale.trunc, 
                        key="hr",  adjustment=NULL)
whale.full.strat2 <- ds(whales[whales$Region.Label=="North",], truncation=whale.trunc,
                        key="hr",  adjustment=NULL)
```
The model selection metric for the full stratification analysis is the sum of AIC for the two distinct analyses.
```{r full-results}
full.aic <- summary(whale.full.strat1)$ds$aic + summary(whale.full.strat2)$ds$aic
```
AIC scores for the two strata analysed separately: `r summary(whale.full.strat1)$ds$aic` +  `r summary(whale.full.strat2)$ds$aic` = `r full.aic`

Contrast `r exp(whale.full.strat1$ddf$par)` with (3.182, 0.5547) and `r exp(whale.full.strat2$ddf$par)` with (2.770, 0.9706).  The first pair of estimates were produced with this current R analysis whereas the second set of each pair was produced by Program Distance 6.2 (numerical results may be slightly different on your computer).

#### Using stratum as a covariate

Here we manufacture a new variable in our dataset `stratum` based upon the `Region.Label`.  The new variable is then used in a formula as a discrete covariate.
```{r covariatefit, message=FALSE, warning=FALSE}
whales$stratum <- ifelse(whales$Region.Label=="North", "N", "S")
whale.strat.covariate <- ds(whales, truncation=whale.trunc, quiet=TRUE,
                  formula = ~as.factor(stratum),
                  key="hr",  adjustment=NULL)
```
AIC score for this model with stratum as a covariate is `r summary(whale.strat.covariate)$ds$aic`.

```{r pooledfit, message=FALSE, warning=FALSE}
whale.pooledf0 <- ds(whales, truncation=whale.trunc,
                        key="hr",  adjustment=NULL)
```
 AIC score for this model pooling sightings from both strata into a single detection function is `r summary(whale.pooledf0)$ds$aic`.

##### Model selection results
Our model selection table for this stratified survey design

Model | Num. parameters |AIC
------|--------|--------------
Full geographic stratification | 4 | `r full.aic`
Detection function shared between strata | 2 | `r summary(whale.pooledf0)$ds$aic`
Stratum as covariate | 3 | `r summary(whale.strat.covariate)$ds$aic`


This shows that the pooled analysis (AIC=`r summary(whale.pooledf0)$ds$aic`) is not preferred to the full geographic stratification analysis (AIC=`r full.aic`).  The model with the smaller AIC is preferable.  However if we introduce `stratum` as a covariate, this forms a halfway house between the extremes, with an added parameter causing the two detection functions to share the same basic shape, but detectability falls off more slowly in one stratum compared to the other (see following figure) and the lowest of the three AIC scores is `r summary(whale.strat.covariate)$ds$aic`.

The detection function that falls off most rapidly is the detection function for the southern stratum, nearer the Antarctic coast where observation conditions were understandable poorer.

```{r covariate-fig, fig.cap="Detection function and qq-plot goodness of fit for minke whales.", fig.width=9}
par(mfrow=c(1,2))
plot(whale.strat.covariate, main="Minke whales, \ndetection function uses stratum as covariate")
covar.fit <- ddf.gof(whale.strat.covariate$ddf)
message <- paste("Cramer von-Mises W=", round(covar.fit$dsgof$CvM$W,3), 
                 "\nP=", round(covar.fit$dsgof$CvM$p,3))
text(0.6, 0.1, message, cex=0.8)
par(mfrow=c(1,1))
```
What remains is to examine the estimated abundance produced by the three models.

### Stratum specific abundance estimates
#### Full geographic stratification
Estimates for the two strata analysed individually.

```{r indivabundfull, echo=FALSE, results="asis", warning=FALSE}
full.1 <- summary(whale.full.strat1)
full.2 <- summary(whale.full.strat2)
strat.spec <- rbind(full.2$dht$individual$N, full.1$dht$individual$N)
kable(strat.spec, format="markdown")
```

#### Pooled detection function
Estimates of group abundance when a detection function is fitted to data pooled across strata.  These results are actually **incorrect** because effort was not equally allocated between the strata.  The southern stratum is much smaller (84000km^2^) than the northern (630000km^2^).  But the southern stratum is more desirable habitat for the minke whales because it is closer to the ice edge in Antarctica.  The southern stratum had much greater survey effort per area than the northern stratum.  This is not represented in the *pooled* analysis. 

```{r indivabundpooled, echo=FALSE,  results="asis"}
pooled <- summary(whale.pooledf0)
kable(pooled$dht$individual$N, format="markdown")
```

#### Stratum as covariate in detection function

Group abundance estimates when strata is a covariate in the detection function (this model was preferred in the model selection exercise).

```{r indivabundcovar, echo=FALSE, results="asis"}
strat <- summary(whale.strat.covariate)
kable(strat$dht$individual$N, format="markdown")
```

