BCMR_output.Rmd

---
title: "Beaver Creek Mark-Recapture Output"
author: "Matt Tyers"
date: "2024-09-18"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, fig.width = 10, fig.height = 7, dpi=300,
                      warning=FALSE, message=FALSE)
```

```{r fig.show='hide'}
# source("R/1_BCMR_data.R")
source("R/3_BCMR_estimates.R")

# # # doing this with the object that will be used to tally!
# # bc_cap2_recaps$sample[bc_cap2_recaps$Tag == 515] <- "float"
# # 
# # 
# # n1 <- table(bc_cap1$sample)    # this should work
# # n2 <- table(bc_cap2$sample)    # this should work
# # # m2 <- table(recaps$sample1)    # wrong, there are Recap $Tag's in bc_cap2
# # m2 <- table(bc_cap2_recaps$sample)  # this will work now
# n1_Stratum1 <- table(bc_cap1$Stratum1)    # this should work
# n2_Stratum1 <- table(bc_cap2$Stratum1)    # this should work
# m2_Stratum1 <- table(bc_cap2_recaps$Stratum1)  # this will work now
# 
# n1_Stratum2 <- table(bc_cap1$Stratum2)    # this should work
# n2_Stratum2 <- table(bc_cap2$Stratum2)    # this should work
# m2_Stratum2 <- table(bc_cap2_recaps$Stratum2)  # this will work now


## defining a function to illustrate differences in distribution, with ks test
ksplot <- function(x1, x2, legend=c("x1","x2"), main="", col=c(1,1), lty=c(1,1), xlab="") {
  d1 <- density(x1, na.rm=TRUE)
  d2 <- density(x2, na.rm=TRUE)
  plot(d1, main=main, col=col[1], lty=lty[1],
       xlim=range(d1$x, d2$x), ylim=range(d1$y, d2$y), xlab=xlab)
  lines(density(x2, na.rm=TRUE), col=col[2], lty=lty[2])
  legend("topright", lty=lty, col=col,
         legend=paste0(legend, " (n=",c(sum(!is.na(x1)), sum(!is.na(x2))),")"))

  ksks <- suppressWarnings(ks.test(x1, x2))
  plot(ecdf(x1), main=c(main,
                        paste0("D=", signif(ksks$statistic, digits=3), ", ",
                        "pval=", signif(ksks$p.value, digits=3), collapse=NULL)),
       col=col[1], lty=lty[1], xlab=xlab)
  plot(ecdf(x2), col=col[2], lty=lty[2], add=TRUE)
  legend("bottomright", lty=lty, col=col,
         legend=paste0(legend, " (n=",c(sum(!is.na(x1)), sum(!is.na(x2))),")"))
}
```

## Evaluation of Assumptions

### Length Selectivity

Length selectivity was assessed using Kolmogorov-Smirnov (KS) tests, using the lengths of Arctic grayling captured in each event tested against the lengths of fish that were captured in both events.  KS tests showed strong evidence of length selectivity in both events when considered without stratification.

```{r}
par(family="serif")
par(mfrow=c(2,2))

ksplot(x1=bc_cap1$Length, x2=bc_cap1_recaps$Length,
       main="Event 1", legend=c("All","Recaps"), col=c(1,2),
       xlab="Fork Length (mm)")
ksplot(x1=bc_cap2$Length, x2=bc_cap2_recaps$Length,
       main="Event 2", legend=c("All","Recaps"), col=c(1,4),
       xlab="Fork Length (mm)")
```

However, after splitting the data by geographic stratum (Beaver Creek vs. Nome Creek), there is no longer any evidence of length selectivity for Beaver Creek, and only evidence of length selectivity in Nome Creek in the second event (KS test with Event 1 data).  Results are plotted below.

It should be noted that one tagged fish (tag # 515) was tagged in Beaver Creek and recaptured in Nome Creek.  Since this fish was tagged 
`r round(abs(42.80986 - bc_cap1_recaps$upstream[bc_cap1_recaps$Tag==515]), 2)`
river kilometers (rkm) from the confluence and recaptured 
`r round(abs(42.80986 - bc_cap2_recaps$upstream[bc_cap2_recaps$Tag==515]), 2)`
rkm from the confluence, it was treated as being part of the Beaver Creek throughout the study for the sake of consistency in stratification.  Since only one tagged fish was captured in another spatial stratum, it was deemed sufficient to reassign this fish and make the assumption of negligible movement between strata, rather than use a partially-stratified (Darroch) estimator. 

```{r}
par(family="serif")
par(mfrow=c(2,2))

par(mfrow=c(2,2))
# ks.test(bc_cap1$Length[bc_cap1$Stratum1=="Beaver"], bc_cap1_recaps$Length[bc_cap1_recaps$Stratum1=="Beaver"])
# D = 0.10478, p-value = 0.7664
ksplot(bc_cap1$Length[bc_cap1$Stratum1=="Beaver"], 
       bc_cap1_recaps$Length[bc_cap1_recaps$Stratum1=="Beaver"],
       main="Event 1 - Beaver", legend=c("All","Recaps"), col=c(1,2),
       xlab="Fork Length (mm)")

# ks.test(bc_cap2$Length[bc_cap2$Stratum1=="Beaver"], bc_cap2_recaps$Length[bc_cap2_recaps$Stratum1=="Beaver"])
# D = 0.067667, p-value = 0.9869
ksplot(bc_cap2$Length[bc_cap2$Stratum1=="Beaver"], 
       bc_cap2_recaps$Length[bc_cap2_recaps$Stratum1=="Beaver"],
       main="Event 2 - Beaver", legend=c("All","Recaps"), col=c(1,4),
       xlab="Fork Length (mm)")

par(mfrow=c(2,2))
# ks.test(bc_cap1$Length[bc_cap1$Stratum1=="Nome"], bc_cap1_recaps$Length[bc_cap1_recaps$Stratum1=="Nome"])
# D = 0.16239, p-value = 0.03623
ksplot(bc_cap1$Length[bc_cap1$Stratum1=="Nome"], 
       bc_cap1_recaps$Length[bc_cap1_recaps$Stratum1=="Nome"],
       main="Event 1 - Nome", legend=c("All","Recaps"), col=c(1,2),
       xlab="Fork Length (mm)")

# ks.test(bc_cap2$Length[bc_cap2$Stratum1=="Nome"], bc_cap2_recaps$Length[bc_cap2_recaps$Stratum1=="Nome"])
# D = 0.14644, p-value = 0.09537
ksplot(bc_cap2$Length[bc_cap2$Stratum1=="Nome"], 
       bc_cap2_recaps$Length[bc_cap2_recaps$Stratum1=="Nome"],
       main="Event 2 - Nome", legend=c("All","Recaps"), col=c(1,4),
       xlab="Fork Length (mm)")
```

It was also of interest to estimate abundance and length composition for the stretch of river comparable to the 2000 study, discussed further.  There was no evidence of length selectivity in this stratum, for either capture event.

```{r}
par(mfrow=c(2,2))
# ks.test(bc_cap1$Length[bc_cap1$Stratum2=="2000 Study"], bc_cap1_recaps$Length[bc_cap1_recaps$Stratum2=="2000 Study"])
# D = 0.15153, p-value = 0.1002
ksplot(bc_cap1$Length[bc_cap1$Stratum2=="2000 Study"], bc_cap1_recaps$Length[bc_cap1_recaps$Stratum2=="2000 Study"],
       main="Event 1 - 2000 Study", legend=c("All","Recaps"), col=c(1,2),
       xlab="Fork Length (mm)")

# ks.test(bc_cap2$Length[bc_cap2$Stratum2=="2000 Study"], bc_cap2_recaps$Length[bc_cap2_recaps$Stratum2=="2000 Study"])
# D = 0.14634, p-value = 0.1611
ksplot(bc_cap2$Length[bc_cap2$Stratum2=="2000 Study"], bc_cap2_recaps$Length[bc_cap2_recaps$Stratum2=="2000 Study"],
       main="Event 2 - 2000 Study", legend=c("All","Recaps"), col=c(1,4),
       xlab="Fork Length (mm)")
```

\pagebreak

### Spatial Selectivity

#### KS tests

Spatial selectivity was similarly evaluated using KS tests, in these cases testing the upriver position (distance in rkm from the lowest point in the study area) of the full samples in each event versus the subset of recaptured fish.  Results are plotted below.  Similarly, testing showed strong evidence of spatial selectivity when pooling the Beaver Creek sample together with the Nome Creek sample; however, there was no evidence of spatial selectivity when the two samples were considered separately.

```{r}
## -------------- spatial selectivity? ---------------- ##
par(family="serif")
par(mfrow=c(2,2))

# ks.test(bc_cap1$upstream, bc_cap1_recaps$upstream)
# D = 0.22825, p-value = 8.778e-06
ksplot(bc_cap1$upstream, bc_cap1_recaps$upstream,
       main="Event 1", legend=c("All","Recaps"), col=c(1,2),
       xlab="Upstream Position (rkm)")

# ks.test(bc_cap2$upstream, bc_cap2_recaps$upstream)
# D = 0.49899, p-value < 2.2e-16
ksplot(bc_cap2$upstream, bc_cap2_recaps$upstream,
       main="Event 2", legend=c("All","Recaps"), col=c(1,4),
       xlab="Upstream Position (rkm)")
```

```{r}
## -------------- spatial selectivity? ---------------- ##
par(family="serif")
par(mfrow=c(2,2))

## actually I don't think we need to worry about this one
### actually yes we do!!  This suggests much more fish in the lower river
### might even be float vs hike!   - IT IS, KEEP THIS STRATIFICATION!!!
par(mfrow=c(2,2))
# ks.test(bc_cap1$upstream[bc_cap1$Stratum1=="Beaver"], bc_cap1_recaps$upstream[bc_cap1_recaps$Stratum1=="Beaver"])
# D = 0.11782, p-value = 0.6284
ksplot(bc_cap1$upstream[bc_cap1$Stratum1=="Beaver"], bc_cap1_recaps$upstream[bc_cap1_recaps$Stratum1=="Beaver"],
       main="Event 1 - Beaver", legend=c("All","Recaps"), col=c(1,2),
       xlab="Upstream Position (rkm)")

# ks.test(bc_cap2$upstream[bc_cap2$Stratum1=="Beaver"], bc_cap2_recaps$upstream[bc_cap2_recaps$Stratum1=="Beaver"])
# D = 0.12232, p-value = 0.5175
ksplot(bc_cap2$upstream[bc_cap2$Stratum1=="Beaver"], bc_cap2_recaps$upstream[bc_cap2_recaps$Stratum1=="Beaver"],
       main="Event 2 - Beaver", legend=c("All","Recaps"), col=c(1,4),
       xlab="Upstream Position (rkm)")

par(mfrow=c(2,2))
# ks.test(bc_cap1$upstream[bc_cap1$Stratum1=="Nome"], bc_cap1_recaps$upstream[bc_cap1_recaps$Stratum1=="Nome"])
# D = 0.1425, p-value = 0.09111
ksplot(bc_cap1$upstream[bc_cap1$Stratum1=="Nome"], bc_cap1_recaps$upstream[bc_cap1_recaps$Stratum1=="Nome"],
       main="Event 1 - Nome", legend=c("All","Recaps"), col=c(1,2),
       xlab="Upstream Position (rkm)")

# ks.test(bc_cap2$upstream[bc_cap2$Stratum1=="Nome"], bc_cap2_recaps$upstream[bc_cap2_recaps$Stratum1=="Nome"])
# D = 0.11333, p-value = 0.3215
ksplot(bc_cap2$upstream[bc_cap2$Stratum1=="Nome"], bc_cap2_recaps$upstream[bc_cap2_recaps$Stratum1=="Nome"],
       main="Event 2 - Nome", legend=c("All","Recaps"), col=c(1,4),
       xlab="Upstream Position (rkm)")
```

#### $\chi^2$ tests

Consistency tests ($\chi^2$), while more sensitive to detecting selectivity, confirm these results, with output as shown below.  It should be noted that spatial strata for Beaver Creek are broken by the confluence with Trail Creek, and the spatial strata for Nome Creek are broken by the lower and upper endpoints of the stretch comparable to the 2000 study, discussed further.

**Full Study Area (all data pooled) - No tests satisfied**

```{r}
library(recapr)  # for automated consistency tests

# with(subset(bc_all, event=="mark"),
     # table(cut(Site, breaks=c(0,15,23,28,32,35,42))))

chisq_strat <- rep(NA, nrow(bc_all))
chisq_strat[bc_all$seg==7 | bc_all$seg==8] <- 1
chisq_strat[bc_all$seg==5 | bc_all$seg==6| bc_all$seg==4] <- 2
chisq_strat[bc_all$Stratum2=="Lower Nome"] <- 3
chisq_strat[bc_all$Stratum2=="2000 Study"] <- 4
chisq_strat[bc_all$Stratum2=="Upper Nome"] <- 5

bc_all$chisq_strat2000 <- ifelse(bc_all$upstream <= 57, 1,
                          ifelse(bc_all$upstream > 62, 3, 2))

# table(chisq_strat)
n1 <- table(chisq_strat[bc_all$event=="mark"])
n2 <- table(chisq_strat[bc_all$event=="recap"])

recap_tags <- bc_cap1_recaps$Tag
m2strata1 <- chisq_strat[bc_all$Tag %in% recap_tags & bc_all$event=="mark"][order(bc_all$Tag[bc_all$Tag %in% recap_tags & bc_all$event=="mark"])]
m2strata2 <- chisq_strat[bc_all$Tag %in% recap_tags & bc_all$event=="recap"][order(bc_all$Tag[bc_all$Tag %in% recap_tags & bc_all$event=="recap"])]

# consistencytest(n1=n1, n2=n2, m2strata1=m2strata1, m2strata2 = m2strata2)

with(bc_all, {
  n1 <- table(chisq_strat[event=="mark"])
  n2 <- table(chisq_strat[event=="recap"])

  m2strata1 <- chisq_strat[Tag %in% recap_tags & event=="mark"][order(Tag[Tag %in% recap_tags & event=="mark"])]
  m2strata2 <- chisq_strat[Tag %in% recap_tags & event=="recap"][order(Tag[Tag %in% recap_tags & event=="recap"])]

  consistencytest(n1=n1, n2=n2, m2strata1=m2strata1, m2strata2 = m2strata2)
})  # tests are not satisfied
```

**Beaver Creek - Tests 2 and 3 satisfied**

```{r}
bc_all$chisq_strat <- chisq_strat
with(subset(bc_all, Stratum1=="Beaver"), {
  n1 <- table(chisq_strat[event=="mark"])
  n2 <- table(chisq_strat[event=="recap"])

  m2strata1 <- chisq_strat[Tag %in% recap_tags & event=="mark"][order(Tag[Tag %in% recap_tags & event=="mark"])]
  m2strata2 <- chisq_strat[Tag %in% recap_tags & event=="recap"][order(Tag[Tag %in% recap_tags & event=="recap"])]

  consistencytest(n1=n1, n2=n2, m2strata1=m2strata1, m2strata2 = m2strata2)
}) # tests 2 and 3 are satisfied

```

**Nome Creek - Test 3 satisfied**

```{r}
with(subset(bc_all, Stratum1=="Nome"), {
  n1 <- table(chisq_strat[event=="mark"])
  n2 <- table(chisq_strat[event=="recap"])

  m2strata1 <- chisq_strat[Tag %in% recap_tags & event=="mark"][order(Tag[Tag %in% recap_tags & event=="mark"])] %>% as.factor %>% as.numeric
  m2strata2 <- chisq_strat[Tag %in% recap_tags & event=="recap"][order(Tag[Tag %in% recap_tags & event=="recap"])]  %>% as.factor %>% as.numeric

  consistencytest(n1=n1, n2=n2, m2strata1=m2strata1, m2strata2 = m2strata2)
}) # test 3 is satisfied
```

#### 2000 Study Area - KS tests

It was also of interest to estimate abundance and length composition for the stretch of river comparable to the 2000 study, discussed further.  There was no evidence of spatial selectivity in this stratum for either capture event.

```{r}
par(mfrow=c(2,2))
par(family="serif")
# ks.test(bc_cap1$upstream[bc_cap1$Stratum2=="2000 Study"], bc_cap1_recaps$upstream[bc_cap1_recaps$Stratum2=="2000 Study"])
# D = 0.13215, p-value = 0.205
ksplot(bc_cap1$upstream[bc_cap1$Stratum2=="2000 Study"], bc_cap1_recaps$upstream[bc_cap1_recaps$Stratum2=="2000 Study"],
       main="Event 1 - 2000 Study", legend=c("All","Recaps"), col=c(1,2),
       xlab="Upstream Position (rkm)")

# ks.test(bc_cap2$upstream[bc_cap2$Stratum2=="2000 Study"], bc_cap2_recaps$upstream[bc_cap2_recaps$Stratum2=="2000 Study"])
# D = 0.072168, p-value = 0.919
ksplot(bc_cap2$upstream[bc_cap2$Stratum2=="2000 Study"], bc_cap2_recaps$upstream[bc_cap2_recaps$Stratum2=="2000 Study"],
       main="Event 2 - 2000 Study", legend=c("All","Recaps"), col=c(1,4),
       xlab="Upstream Position (rkm)")
```

#### 2000 Study Area - $\chi^2$ tests

**Test 2 satisfied**

```{r}
with(subset(bc_all, Stratum2=="2000 Study"), {
  n1 <- table(chisq_strat2000[event=="mark"])
  n2 <- table(chisq_strat2000[event=="recap"])

  m2strata1 <- chisq_strat2000[Tag %in% recap_tags & event=="mark"][order(Tag[Tag %in% recap_tags & event=="mark"])] %>% as.factor %>% as.numeric
  m2strata2 <- chisq_strat2000[Tag %in% recap_tags & event=="recap"][order(Tag[Tag %in% recap_tags & event=="recap"])]  %>% as.factor %>% as.numeric

  consistencytest(n1=n1, n2=n2, m2strata1=m2strata1, m2strata2 = m2strata2)
}) # test 3 is satisfied
```

### Growth Recruitment

```{r}
## ---------------- growth recruitment? ----------------- ##
length1 <- bc_cap1_recaps$Length[order(bc_cap1_recaps$Tag)]
length2 <- bc_cap2_recaps_justtags$Length[order(bc_cap2_recaps_justtags$Tag)]

diffs <- length2 - length1
diffs <- diffs[abs(diffs) < 80]  # what happens when we remove that big outlier
```

With the relatively short time period between mark and recapture events, growth recruitment was not anticipated.  Quantitative assessment of growth recruitment confirms this: after culling one outlying measurement, the mean difference in recorded length was `r round(mean(diffs), 1)` mm, as compared to a standard deviation of `r round(sd(diffs), 1)` mm.  This yielded a paired t-test p-value of `r round(t.test(diffs)$p.value, 2)` when testing the null hypothesis of zero growth, giving no evidence of appreciable growth recruitment.

### Immigration & Emigration

With the high degree of summer site fidelity observed in the telemetry portion of this study, it was not anticipated that a meaningful degree of immigration or emigration would occur in the time period between the mark and recapture events.

```{r}
up1 <- bc_cap1_recaps$upstream[order(bc_cap1_recaps$Tag)]
up2 <- bc_cap2_recaps_justtags$upstream[order(bc_cap2_recaps_justtags$Tag)]
diffs <- up2-up1

# data.frame(tag = bc_cap1_recaps$Tag[order(bc_cap1_recaps$Tag)],
#            upstream1 = up1, upstream2 = up2, upstream_diff = diffs) %>%
#   write.csv("R_output/upstream_diffs.csv")
```

Movement of marked fish may be approximately inferred from movement observed in recaptured fish. Movement is illustrated below, with upriver position (distance in rkm from the lowest point in the study area) expressed on the respective x-axes, with movement between events expressed on the y-axes.  As is evident from both figures, the majority of recaptured fish did not travel any appreciable distance, with only `r round(100*mean(abs(diffs)>.5))`% of fish observed to travel more than 0.5 rkm. Additionally, those that did move did not travel far in comparison to the size of the study area, and relatively few fish were observed near the endpoints of the study area.

```{r fig.height=10}
# plot(diffs)
# sd(diffs)
# mean(abs(diffs))
# median(abs(diffs))
# mean(abs(diffs)>.5)
# par(mfrow=c(1,1))
# plot(ecdf(abs(diffs)), xlim=c(0,5))

par(mfrow=c(2,1))
par(family="serif")

plot(x=up1, y=diffs, asp=1,
     xlab="Upriver Position - Mark Event (rkm)",
     main="", ylab="Movement (rkm)")
plot(x=up2, y=diffs, asp=1,
     xlab="Upriver Position - Recapture Event (rkm)",
     main="", ylab="Movement (rkm)")  # hardly any fish near endpoints, it's probably just fine

```

```{r}

nboot <- 10000
count_table <- matrix(nrow=nboot, ncol=4)
for(i in 1:nboot) {
  upstream_boot <- bc_cap1$upstream + sample(diffs, size=nrow(bc_cap1), replace=TRUE)
  # table(cut(upstream_boot,
  #           c(-100, 0, 42.81, 82.52, 200)))
  #           # levels=c("(-100,0]",   " (0,42.8]", "(42.8,82.5]",  "(82.5,200])")))
  count_table[i,1] <- sum(upstream_boot <= 0)
  count_table[i,2] <- length(upstream_boot %s_inside(]% c(0, 42.81))
  count_table[i,3] <- length(upstream_boot %s_inside(]% c(42.81, 82.5))
  count_table[i,4] <- sum(upstream_boot > 82.5)
}
prop_table <- count_table/nrow(bc_cap1)
# colMeans(prop_table)
# we might lose 1% of tagged fish out the bottom, and the top is negligible


```

Since fairly equal numbers of recaptured fish were observed to move up- and downstream, it seems reasonable to assume that an approximately equal degree of movement occurs in and out of the study area at the endpoints and at the confluences of tributary streams, and between the spatial strata (e.g. Beaver vs Nome Creeks).  It is therefore reasonable to assume that any change in abundance in the population of inference (i.e. immigration or emigration) is negligible.  However, the potential does exist that marked fish may have left the study area between events, thus imparting some small bias.  This was investigated by approximating the movement of tagged fish by resampling the differences in observed position with replacement 10,000 times and adding them to the upriver positions observed in the mark event.  In this simulation, a negligible mean proportion (`r round(100*colMeans(prop_table)[4], 2)`%) of tagged fish had an ending positions above the top end of the study area, and a small mean proportion (`r round(100*colMeans(prop_table)[1], 2)`%) had an ending position below the bottom end.  It may be worth accounting for the bias this would impart, but it is likely that the bias would be small.

## Abundance Estimation

Abundance was estimated for each geographic stratum *t* (e.g. Beaver Creek vs. Nome Creek) using a Bailey estimator shown below:

$$\hat{N}_t = \frac{n_{1,t}(n_{2,t}+1)}{m_{2,t}+1}$$

and 

$$\hat{V}(\hat{N}_t) =  \frac{n_{1,t}^2(n_{2,t}+1)(n_{2,t}-m_{2,t})}{(m_{2,t}+1)^2(m_{2,t}+2)}$$

The abundances and associated variances estimated for each stratum were summed for the purpose of combined estimates:

$$\hat{N} = \hat{N}_{Beaver} + \hat{N}_{Nome}$$

and

$$\hat{V}(\hat{N}) = \hat{V}(\hat{N}_{Beaver}) + \hat{V}(\hat{N}_{Nome})$$

It should be noted that during the recapture event, one fish got away before its tag number was recorded, and one fish had no tag but did have a fin clip, indicating capture.  Both were recorded as recaptures.

```{r}
######### ----------------- abundance estimates ----------------- #########

# library(recapr)
# 
# nhat_Stratum1 <- NBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1) %>% 
#   unname %>% round %>% formatC(format="d", big.mark=",")  # estimated abundance
# 
# se_nhat_Stratum1 <- seBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1) %>% 
#   unname %>% round  %>% formatC(format="d", big.mark=",") # SE of estimated abundance
# 
# 
# nhat_combined_Stratum1 <- Nstrat(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1, 
#                                  estimator = "Bailey") %>% 
#   unname %>% round %>% formatC(format="d", big.mark=",")   # estimated abundance
# se_nhat_combined_Stratum1 <- sestrat(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1, 
#                                      estimator = "Bailey") %>% 
#   unname %>% round %>% formatC(format="d", big.mark=",")  # SE of estimated abundance
# 
# 
# nhat_Stratum2 <- NBailey(n1=n1_Stratum2, n2=n2_Stratum2, m2=m2_Stratum2)[c(2,3,1,4)] %>% 
#   # unname %>% 
#   round %>% formatC(format="d", big.mark=",")  # estimated abundance
# 
# se_nhat_Stratum2 <- seBailey(n1=n1_Stratum2, n2=n2_Stratum2, m2=m2_Stratum2)[c(2,3,1,4)] %>% 
#   # unname %>% 
#   round  %>% formatC(format="d", big.mark=",") # SE of estimated abundance
# 
# 
# nhat_combined_Stratum2 <- Nstrat(n1=n1_Stratum2, n2=n2_Stratum2, m2=m2_Stratum2, 
#                                  estimator = "Bailey") %>% 
#   unname %>% round %>% formatC(format="d", big.mark=",")   # estimated abundance
# se_nhat_combined_Stratum2 <- sestrat(n1=n1_Stratum2, n2=n2_Stratum2, m2=m2_Stratum2, 
#                                      estimator = "Bailey") %>% 
#   unname %>% round %>% formatC(format="d", big.mark=",")  # SE of estimated abundance
```

### Results - Beaver & Nome Creeks 

The abundance was estimated as `r unname(nhat_Stratum1["Beaver"])` and `r unname(nhat_Stratum1["Nome"])` for Beaver Creek and Nome Creek, respectively (SE `r unname(se_nhat_Stratum1["Beaver"])` and `r unname(se_nhat_Stratum1["Nome"])`, respectively), and `r nhat_combined_Stratum1` (SE `r se_nhat_combined_Stratum1`) for the full study area.

### Results - Finer stratification

Also of interest was the reach corresponding to the 2000 study, defined as the stretch of Nome Creek between the airstrip and the confluence with Moose Creek.  To that end, the Nome Creek stratum was further subdivided into Lower, Middle (2000 Study), and Upper strata.  In addition to the one fish (Tag # 515) observed to travel between Nome Creek and Beaver Creek, three additional fish were observed to move between strata of Nome Creek, as tabulated below:

```{r}
# problemtable1
```

```{r}
problemtable2 
```

All three fish were located near the stratum boundary during one event and substantially further from the boundary in the other event.  In all cases, recaptured fish were assigned to the spatial stratum reflecting the largest movement from the boundary.

If the movement among strata is accepted and assumed negligible, the finer stratification of Nome Creek yields the following abundance estimates.  

```{r}
thenames <- c("Beaver", "Lower Nome", "2000 Study", "Upper Nome")
cbind(n1=formatC(unname(n1_Stratum2[thenames]), format="d", big.mark=","),
           n2=formatC(unname(n2_Stratum2[thenames]), format="d", big.mark=","),
           m2=unname(m2_Stratum2[thenames]), 
           Nhat = nhat_Stratum2[thenames], 
           SE_Nhat = se_nhat_Stratum2[thenames]) %>% knitr::kable()
```

## Length Distribution

```{r, results='asis'}
######### ----------------- ASL estimates (length bins) ----------------- #########

# lengthbreaks <- c(250, 270, 300, 350, 400, 450) # c(250, 270, 300, 400)
# 
# bc_all_justonce <- rbind(bc_cap1,
#                          subset(bc_cap2, Tag=="-"))

lengthbreakcount <- ifelse(length(lengthbreaks)-1 <= 9, 
                           c("one", "two", "three", "four", "five", "six", "seven", "eight", "nine")[length(lengthbreaks)-1],
                           length(lengthbreaks)-1)
lengthbinwords <- paste(lengthbreaks[-length(lengthbreaks)],
                        lengthbreaks[-1]-1, sep="-")
lengthbinphrase <- paste(paste(lengthbinwords[-length(lengthbinwords)], collapse=", "),
      lengthbinwords[length(lengthbinwords)], sep=", and ")

boilerplate <- capture.output(
  ASL_boilerplate(stratum = as.numeric(as.factor(bc_all$Stratum1)),
                # length = bc_all$Length,
                age = cut(bc_all$Length, lengthbreaks, right=FALSE),
                Nhat = as.numeric(NBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1)),
                se_Nhat = as.numeric(seBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1))))

# replacing "age" with "length"
boilerplate_fix <- gsub(pattern=" age", replacement=" length", x=boilerplate)

# splitting the boilerplate text: up to references, and references onward
boilerplate_uptoref <- boilerplate_fix[1:(which(boilerplate_fix=="## References")-1)]
boilerplate_refon <- boilerplate_fix[(which(boilerplate_fix=="## References")):length(boilerplate_fix)]

# printing the boilerplate up to references
cat(boilerplate_uptoref, sep="\n")
```

The length composition of Arctic grayling was described by considering lengths in 
`r lengthbreakcount` categories: `r lengthbinphrase` mm FL.

The length composition was tabulated and stratified estimators were calculated as appropriate for the following cases.  For cases in which length measurements from both events (i.e. marking and recapture) were used, recaptured fish (that is, present in both the marking and recapture events) were only incorporated once.

\pagebreak

### Full study area (stratified) - just using Event 1 data

```{r, results='asis'}

digitvec <- c(0,3,3,0,0)

# just using first event sample
# ASL_table(stratum = as.numeric(as.factor(bc_cap1$Stratum1)),
#           # length = bc_cap1$Length,
#           age = cut(bc_cap1$Length, lengthbreaks, right=FALSE),
#           Nhat = as.numeric(NBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1)),
#           se_Nhat = as.numeric(seBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1))) %>%
#   knitr::kable(digits=digitvec)
ASL_firstevent %>% knitr::kable(digits=digitvec)
```

### Beaver Creek only - using both events

```{r, results='asis'}

# pooling both events, but separating by capture stratum
# with(subset(bc_all, Stratum1 == "Beaver"),
#      ASL_table(
#        # length = Length,
#        age = cut(Length, lengthbreaks, right=FALSE),
#        Nhat = as.numeric(NBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1))[1],
#        se_Nhat = as.numeric(seBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1))[1])) %>%
#   knitr::kable(digits=digitvec)
ASL_Beaver %>% knitr::kable(digits=digitvec)
```

### Nome Creek only - just using Event 1 data

```{r, results='asis'}

# with(subset(bc_all, Stratum1 == "Nome"),
#      ASL_table(
#        # length = Length,
#        age = cut(Length, lengthbreaks, right=FALSE),
#        Nhat = as.numeric(NBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1))[2],
#        se_Nhat = as.numeric(seBailey(n1=n1_Stratum1, n2=n2_Stratum1, m2=m2_Stratum1))[2])) %>%
#   knitr::kable(digits=digitvec)
ASL_Nome %>% knitr::kable(digits=digitvec)
```

### 2000 Study Area only - using both events

```{r, results='asis'}

# with(subset(bc_all, Stratum2 == "2000 Study"),
#      ASL_table(
#        # length = Length,
#        age = cut(Length, lengthbreaks, right=FALSE),
#        Nhat = as.numeric(NBailey(n1=n1_Stratum2, n2=n2_Stratum2, m2=m2_Stratum2))[1],
#        se_Nhat = as.numeric(seBailey(n1=n1_Stratum2, n2=n2_Stratum2, m2=m2_Stratum2))[1])) %>%
#   knitr::kable(digits=digitvec)
ASL_2000 %>% knitr::kable(digits=digitvec)
```

```{r, results='asis'}
# printing references
cat(boilerplate_refon, sep="\n")
```