Asymptotic_Weight_rept.Rmd

---
title: "Estimating Asymptotic Weight for Lake Trout"
author: "Matt Tyers"
date: "2024-05-31"
output: word_document
---

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

```{r}
library(tidyverse)
library(jagsUI)
library(jagshelper)
```


## Purpose & general outline of methods

Application of the Lester Model to lake trout in Alaska lakes may be substantially improved by direct estimation of asymptotic weight $W_{\infty}$ when possible.  In the computational framework presented by Lester, et al., $W_{\infty}$ is estimated by means of a sequence of relationships: first, the relationship observed between lake area and asymptotic length observed in Canada; then, the aggregate relationship between length and weight observed in Canada.

We have data at a finer resolution: for 84 lakes with documented lake trout presence, in addition to geographical data (lake area and latitude) we have sampling data on 30,000+ individual fish (4089 weights; 32,244 lengths; 1523 ages).  Not only can the relationships between length and weight and between lake area and asymptotic length be estimated for Alaska, they can be estimated specifically for each of the 84 lakes as appropriate, and parameters reported by Lester, et. al. may be incorporated as priors.  Furthermore, in some cases, asymptotic size may be estimated directly without intermediary steps.

Depending on data availability for each lake, $W_{\infty}$ was estimated using one of four methods:

1. $W_{\infty}$ from a weight sample quantile value

2. $L_{\infty}$ from a length ~ age relationship, then $W_{\infty}$ from a weight ~ length relationship

3. $L_{\infty}$ from a length sample quantile value, then $W_{\infty}$ from a weight ~ length relationship

4. $L_{\infty}$ from a $L_{\infty}$ ~ lake area relationship, then $W_{\infty}$ from a weight ~ length relationship

Overall, the relationships between weight and length were strong when data were present, and direct estimation of $L_{\infty}$ from age yielded results that appeared quite reasonable.  Estimates of $L_{\infty}$ appeared to be reasonably estimated from length sample quantiles, and application of the same quantiles to weight samples seemed reasonable due to the strong and monotonic relationship between length and weight.  The relationship between lake area and $L_{\infty}$ was quite weak, however, and estimation of $W_{\infty}$ will be greatly improved with collection of fish-level samples (length, weight, and/or age) from these lakes.

\pagebreak

## Models & model selection

### Weight ~ Length (+ Latitude + Area)

A log-log regression was used to model weight as a function of length, expanding somewhat on the form used by Lester, et al.  Since fish-level data (lengths and weights) were available from several lakes of differing habitat, the length-weight relationships were estimated separately for each lake, with slope and intercept parameters modeled hierarchically as shown below.

$log(W_{j[i]}) = \beta_{0[j]} + \beta_{1[j]}log(L_{j[i]}) + \epsilon_i$

which is equivalent to

$W_{j[i]} = e^{\beta_{0[j]}} \times L^{\beta_{1[j]}} \times e^{\epsilon_i}$

This was extended to allow possible relationships between the lake-level intercept parameters $\beta_{0[j]}$ and slope parameters $\beta_{1[j]}$ and lake-level variables (lake area and latitude), with the selected model expressed below.

$\beta_{0[j]} \sim N(\mu_{\beta_0[j]}, \sigma_{\beta_0}) \text{ for } j \in 1...n_{lakes}$

$\beta_{1[j]} \sim N(\mu_{\beta_1[j]}, \sigma_{\beta_1}) \text{ for } j \in 1...n_{lakes}$

$\mu_{\beta_0[j]} = \gamma_0 + \tau_{0,Area}log(Area_j)  \text{ for } j \in 1...n_{lakes}$

$\mu_{\beta_1[j]} = \gamma_1 +  \tau_{1,Lat}Latitude_j \text{ for } j \in 1...n_{lakes}$

$\epsilon_i \sim N(0,\sigma_{\epsilon})$

Non-informative diffuse Normal priors were used on hyperparameters $\gamma_.$ and $\tau_.$, and appropriately bounded uniform priors were used on all standard deviation parameters $\sigma_.$ on the SD scale.

All possible relationships between intercept and slope parameters ($\beta_{0[j]}$ and $\beta_{1[j]}$, respectively) and lake-level variables (lake area and latitude) were explored using DIC and leave-one-out cross-validation, with the best performing model being that reported above.

For lakes without paired length and weight data, parameters $\beta_{0[j]}$ and $\beta_{0[j]}$ were still modeled from the posterior predictive distribution using lake area and latitude as available.  In the cases in which lake area and latitude were not available, these were imputed for each MCMC sample using draws from Normal distributions with means and standard deviations of the latitude and log area of all lakes.  

Note: Lengths were recentered by subtracting the mean measurement, in order to eliminate collinearity between the intercept and slope parameters and thus aid model convergence.  Raw estimates of the intercept parameters appear quite different from that reported by Lester, et al.; these were then backtransformed for the sake of interpretability and consistency.  Latitude and log area were also recentered by subtracting the respective means.

Intercept and slope parameter estimates are shown below for all lakes with associated weight and length measurements, with log area and latitude of each lake depicted on the X-axes.  Light and heavy vertical bars represent 95% and 50% credible intervals, respectively.  Overlayed trend envelopes represent 95% and 50% credible intervals for the relationships between log area or latitude and regression parameters, as appropriate.

```{r, fig.height=12}
load(file="WL_data_WinfRept.Rdata")

par(mfrow=c(2, 2))
# caterpillar(WL_jags_out$sims.list$b0[,WL_data$whichlakes], 
#             x=WL_data$area[WL_data$whichlakes])
# caterpillar(WL_jags_out$sims.list$b0[,WL_data$whichlakes], 
#             x=WL_data$lat[WL_data$whichlakes])
# caterpillar(WL_jags_out$sims.list$b1[,WL_data$whichlakes], 
#             x=WL_data$area[WL_data$whichlakes])
# caterpillar(WL_jags_out$sims.list$b1[,WL_data$whichlakes], 
#             x=WL_data$lat[WL_data$whichlakes])

lakeswithLWdata <- WL_data$whichlakes[WL_data$whichlakes %in% WL_data$lake[!is.na(WL_data$y)]]
caterpillar(WL_jags_out$sims.list$b0[,lakeswithLWdata], 
            x=WL_data$area[lakeswithLWdata],
            xlab = "log(area) - recentered", main="b0 (intercept)")
envelope(WL_jags_out$sims.list$mu_b0[,lakeswithLWdata], 
            x=WL_data$area[lakeswithLWdata], add=TRUE, dark=.2)
caterpillar(WL_jags_out$sims.list$b0[,lakeswithLWdata], 
            x=WL_data$lat[lakeswithLWdata],
            xlab = "latitude - recentered", main="b0 (intercept)")
caterpillar(WL_jags_out$sims.list$b1[,lakeswithLWdata], 
            x=WL_data$area[lakeswithLWdata],
            xlab = "log(area) - recentered", main="b1 (slope)")
caterpillar(WL_jags_out$sims.list$b1[,lakeswithLWdata], 
            x=WL_data$lat[lakeswithLWdata],
            xlab = "latitude - recentered", main="b1 (slope)")
envelope(WL_jags_out$sims.list$mu_b1[,lakeswithLWdata], 
            x=WL_data$lat[lakeswithLWdata], add=TRUE, dark=.2)
```

\pagebreak

Intercept and slope parameter estimates for all lakes are plotted below, with color denoting data availability.  Estimates presented in green represent lakes with paired lengths and weights; estimates presented in blue represent lakes with area and latitude only; and estimates presented in red represent lakes with no associated data.

Two plots are presented for the intercept parameters: the first being the raw parameters from the model (using recentered length data), and the second being the backtransformed parameters to the original scale.

Dotted horizontal lines represent values reported and used by Lester, et al.  Note that parameter estimates are similar to those used by Lester, et al., but in many cases are substantially different from one another, suggesting variation in weight ~ length relationships between lakes.

```{r, fig.height=11, fig.width=10}
par(mfrow=c(3,1))

## colors are differentiated by whether lat/area data and weight data are present
LatArea_present <- sort(unique(laketrout$LakeNum)) %in% WL_data$whichlakes
weight_present <- table(laketrout$LakeNum,!is.na(laketrout$Weight_g))[,2] > 0
# cols <- 3 - 1*(sort(unique(laketrout$LakeNum)) %in% WL_data$whichlakesc)
cols <- ifelse(weight_present, 3, ifelse(LatArea_present, 4, 2))

caterpillar(WL_jags_out, p="b0", col=cols, ylim=c(-0.3, 0.7),
            xlab="Lake", main="b0 (intercept) - from model")
legend("topright", legend=c("Length & Wt", "Lat & Area only", "no data"),
       lwd=3, col=c(3, 4, 2))
caterpillar(WL_jags_out, p="b0_interp", col=cols, ylim=c(-25, -13),
            xlab="Lake", main="b0 (intercept) - original scale")
legend("topright", legend=c("Length & Wt", "Lat & Area only", "no data"),
       lwd=3, col=c(3, 4, 2))
abline(h=-19.56, lty=3)
caterpillar(WL_jags_out, p="b1", col=cols, ylim=c(2.3, 4.4),
            xlab="Lake", main="b1 (slope)")
legend("topright", legend=c("Length & Wt", "Lat & Area only", "no data"),
       lwd=3, col=c(3, 4, 2))
abline(h=3.2, lty=3)
```

Finally, the modeled relationships between length and weight are shown for each lake, overlayed with available data, and plotted on the log-log scale.  The relationship reported by Lester, et al. is overlayed as a dashed line, and fit envelope color represents the availability of data, using the same color convention as the previous set of plots (green = all data, blue = area and latitude only, red = no data).

```{r, fig.width=15, fig.height=20}
# plot regression bands for each lake, overlay data and lester line
xpredict <- seq(from = min(log(laketrout$ForkLength_mm[!is.na(laketrout$Weight_g)]), na.rm=TRUE),
                  to = max(log(laketrout$ForkLength_mm[!is.na(laketrout$Weight_g)]), na.rm=TRUE),
                  length.out=50)
par(mfrow=c(6, 5))
par(mar=c(4,4,4,1))
lakenames <- levels(as.factor(laketrout$LakeName))
for(ilake in 1:max(laketrout$LakeNum)) {
  logweight_predict <- WL_jags_out$sims.list$b0_interp[,ilake] +
    outer(WL_jags_out$sims.list$b1[,ilake], xpredict)
  plot(NA, xlim=range(xpredict), ylim=range(log(laketrout$Weight_g/1000), na.rm=TRUE),
       xlab="log(Length)", ylab="log(Weight)", main=lakenames[ilake])
  if(all(!is.na(logweight_predict))) envelope(logweight_predict, x=xpredict, add=TRUE, col=cols[ilake])
  points(x = log(laketrout$ForkLength_mm[laketrout$LakeNum==ilake]),
         y = log(laketrout$Weight_g[laketrout$LakeNum==ilake]/1000),
         col=adjustcolor(1,alpha.f=.5))
  abline(a=-19.56, b=3.2, lty=3)
}
```

\pagebreak

### Length ~ Age

The relationship between length and age was modeled using a Von Bertalanffy growth model of the form below, with additive error:

$L_{j[i]} = L_{\infty}[j](1-e^{-k[j](a_{j[i]}-t_0[j])}) + \epsilon_i$

Growth parameters $L_{\infty}$, $k$, and $t_0$ were modeled separately for each lake *j*, but with hierarchical distributions according to the form:

$L_{\infty}[j] \sim N(\mu_{L_{\infty}} , \sigma_{L_{\infty}})$

$k[j] \sim logN(\mu_{k} , \sigma_{k})$

$t_0[j] \sim N(\mu_{t_0} , \sigma_{t_0})$

In this case, the parameter of most inferential interest was asymptotic length $L_{\infty}$.  Hyperparameter $\mu_{L_{\infty}}$ was given a weakly informative prior of $N(700, 100)$, and $\mu_{t_0}$ was mildly constrained by using a prior of $N(0,1)$.  All standard deviation $\sigma_.$ parameters were modeled with uniform priors on the SD scale, bounded as appropriate.

Estimates for all parameters are plotted below for each lake.  Note that parameters were only estimated for lakes that had paired length and age data.

```{r, fig.width=10, fig.height=8}
load(file="LA_data_WinfRept.Rdata")

par(mfrow=c(2,2))
caterpillar(LA_jags_out, p="Linf", ylab="Fork Length (mm)", xlab="Lake")
caterpillar(LA_jags_out, p="k", ylab="/year", xlab="Lake")
caterpillar(LA_jags_out, p="t0", ylab="year", xlab="Lake")
```



Estimated growth curves are plotted below for all lakes with paired length and age data, with length and age data overlayed.  Shaded envelopes correspond to 50% and 95% credible intervals for mean length at a given age.

```{r, fig.height=17}
par(mfrow=c(6,4))
par(mar=c(4,4,4,1))
for(ilake in LA_data$whichlakes) {
  envelope(LA_jags_out$sims.list$yfit[,,ilake],
           ylim=c(0, max(LA_data$y[!is.na(LA_data$x)], na.rm=TRUE)),
           main=lakenames[ilake],
           xlab="Age", ylab="Fork Length (mm)")
  points(x=LA_data$x[LA_data$lake==ilake], y=LA_data$y[LA_data$lake==ilake],
         col=adjustcolor(1, alpha.f=.5))
}
```

For lakes with little or no paired age and length data but sufficient length data, a reasonable estimate of $L_{\infty}$ may be gained from a quantile value of the associated length sample.  To estimate the appropriate quantile value (and with error), posterior distributions of the sample quantile value associated with $L_{\infty}$ were estimated for lakes with sufficient data to estimate both the length ~ age relationship and length sample quantiles (ten lakes with $n_{Age} \geq 10$ and $n_{Length} \geq 100$).  This was calculated for each MCMC sample as the proportion of the length sample less than the value of the given MCMC sample of $L_{\infty}$.  MCMC samples of quantiles were then drawn for the remaining lakes by randomly sampling from the MCMC samples, pooled among all lakes with sufficient paired age-length data.  This can be thought of as was analogous to two-stage sampling (lake then MCMC sample), with all lakes with sufficient data weighted equally. 

Posterior distributions of quantiles associated with $L_{\infty}$ are plotted below, before and after imputation by sampling from lakes with sufficient data.

```{r, fig.height=5, fig.width=10}
# getn <- function(x) sum(!is.na(x))
# laketrout_Winf <- data.frame(n_Weight = tapply(laketrout$Weight_g, laketrout$LakeName, getn),
#                              n_Length = tapply(laketrout$ForkLength_mm, laketrout$LakeName, getn),
#                              n_Age = tapply(laketrout$Age, laketrout$LakeName, getn))

# ## need to redo these criteria
# laketrout_Winf$Wt_quantile <- laketrout_Winf$n_Weight > 200
# laketrout_Winf$Ln_Age <- laketrout_Winf$n_Age > 150
# laketrout_Winf$Ln_quantile <- laketrout_Winf$n_Length > 200

par(mfrow=c(1,1))
cols <- ifelse(laketrout_Winf$n_Age >= 10 & laketrout_Winf$n_Length >= 100, 3, 2)
caterpillar(L_quantile, col=cols, ylim=0:1, xlab="Lake", ylab="Linf quantile")
legend("bottomright", legend=c("n >= 100", "n < 100"), lwd=3, col=c(3, 2))

cols <- ifelse(laketrout_Winf$n_Age >= 10 & laketrout_Winf$n_Length >= 100, 3, 4)
caterpillar(L_quantile_imputed, col=cols, ylim=0:1, xlab="Lake", ylab="Linf quantile")
legend("bottomright", legend=c("n >= 100", "imputed"), lwd=3, col=c(3, 4))
```

### L_inf ~ Area

A relationship was noted between log area and asymptotic length by Lester et al. and incorporated in the computational framework they present.  The same relationship was noted in our dataset, but was not very strong.  Therefore slightly informative priors were taken from values reported by Lester et al.

The model that was used is expressed below, with additive error:

$L_{\infty}[j] = \beta_0 (1 - e^{-\beta_1 * (1 + log(A[j]))}) + \epsilon_j$

```{r}
load(file="LinfArea_data_WinfRept.Rdata")
```


Parameter $\beta_0$ was given a Normal prior with mean equal to the value reported by Lester et al. (957) with a CV (coefficient of variation, equal to SD/mean) of `r 100*LinfArea_data$cv_b0_lester`%.  Parameter $\beta_1$ was similarly given a Normal prior with mean equal to the value reported by Lester et al. (0.14) but with a CV of `r 100*LinfArea_data$cv_b1_lester`%.

It should be noted that asymptotic length $L_{\infty}[j]$ for each lake $j$ was not known.  For the purpose of a data input for this model, quantile values were taken of the length samples associated with each lake, with quantile equal to the posterior median quantile value calculated in the previous section. 

The respective densities of the priors and posteriors of both parameters are plotted below.  Below this, posterior envelopes are plotted for the trend and posterior predictive, respectively, with the lake area and estimated associated asymptotic lengths overlayed.  The curve reported by Lester, et al. is expressed as a dotted line.

```{r, fig.height=8}
par(mfrow=c(2,2))
comparepriors(LinfArea_jags_out, col=c(2,4))
plot(LinfArea_data$x, LinfArea_data$y, log="x",
     xlab="Lake Area (ha)", ylab="L_inf (mm)",
     main="Trend")
envelope(LinfArea_jags_out, p="mu", x=LinfArea_data$x, add=TRUE)
curve(957*(1-exp(-0.14*(1+log(x)))), lty=3, add=TRUE)

plot(LinfArea_data$x, LinfArea_data$y, log="x",
     xlab="Lake Area (ha)", ylab="L_inf (mm)",
     main="Post Predictive")
envelope(LinfArea_jags_out, p="ypp", x=LinfArea_data$x, log="x",
         xlab="Lake Area (ha)", ylab="L_inf (mm)",
         main="Post Predictive", add=TRUE)
# points(LinfArea_data$x, LinfArea_data$y)
curve(957*(1-exp(-0.14*(1+log(x)))), lty=3, add=TRUE)
```


## Multimodel synthesis & method selection

Asymptotic weight $W_{\infty}$ was estimated for all lakes using all four methods described above, depending on data availability.  Following this, an appropriate method was selected for each lake depending on relative sample sizes and estimated variance.

### Estimation of $W_{\infty}$ 

#### 1. $W_{\infty}$ from a weight sample quantile value

For each lake, the vector of MCMC samples corresponding to the quantile value associated with asymptotic size for that lake was applied to the (data) sample of weights from that lake.  

#### 2. $L_{\infty}$ from a length ~ age relationship, then $W_{\infty}$ from a weight ~ length relationship

For each lake in which $L_{\infty}$ could be estimated from a length ~ age relationship, the vectors of MCMC samples of parameter estimates from the weight ~ length relationship were applied to the vector of MCMC samples of $L_{\infty}$.  The intent of this was to "convert" asymptotic length to asymptotic weight, while incorporating the uncertainty in estimating asymptotic length and in estimating the weight ~ length relationship. 

#### 3. $L_{\infty}$ from a length sample quantile value, then $W_{\infty}$ from a weight ~ length relationship

For each lake, the vector of MCMC samples corresponding to the quantile value associated with asymptotic size for that lake was applied to the (data) sample of lengths from that lake.  The vectors of MCMC samples of parameter estimates from the weight ~ length relationship were then applied to this vector, again to "convert" asymptotic length to asymptotic weight while incorporating the uncertainty in estimating asymptotic length and in estimating the weight ~ length relationship.

#### 4. $L_{\infty}$ from a $L_{\infty}$ ~ lake area relationship, then $W_{\infty}$ from a weight-length relationship

Finally, for each lake, the vectors of MCMC samples of parameter estimates from the weight ~ length relationship were applied to the vector of MCMC samples of posterior-predicted $L_{\infty}$ from the $L_{\infty}$ ~ lake area relationship.

\pagebreak

### Method selection

Selection among available methods was motivated by a few guiding principles and assumptions: 

* It would be inappropriate to consider a method that does not meet some associated sample size threshold.  

* A relatively large sample size is required to accurately estimate a population quantile (methods 1 and 3), but direct estimation of asymptotic size (method 2) should be possible using a smaller sample size, particularly since the model incorporates information hierarchically. 

* Methods with fewer steps should be preferred over methods with a greater number of steps: heuristically, it makes sense that methods 1-4 are essentially ranked in order of descending preference.

* Since all methods presented incorporate the uncertainty in estimation due to each step, we can expect that methods incorporating fewer steps will have a smaller associated variance.

* However, since many lakes have much more samples for length than for weight or age, it is possible that methods ranked as less preferable might result in better estimates.

* Method 4 should be avoided unless there is no other viable alternative!

```{r}
load(file="Winf_method_WinfRept.Rdata")
```


From here, model selection proceeded according to a few discrete decision rules.  The minimum sample sizes required to accept methods 1, 2, and 3 were defined as `r n_weight_accept` weights, `r n_age_accept` ages, and `r n_length_accept` lengths, respectively.  Of the methods with minimum sample sizes being met, the method with the minimum variance in estimating $W_{\infty}$ was selected.  If none of methods 1-3 were selected, method 4 was accepted if lake area was known.  

Estimated $W_{\infty}$ from all four methods is plotted below for all lakes, with the selected method highlighted in green and associated sample sizes for the three methods annotated along the X-axis.

```{r, fig.width=12, fig.height=18}
par(mfrow=c(6,5))
par(mar=c(4,4,4,1))
for(i in 1:length(Winf_method)) {
  cols <- rep(2,4)
  cols[Winf_method[i]] <- 3
  caterpillar(Winf_all_all[,i,],
              col = cols,
              xax = c(paste0("n=",laketrout_Winf$n_Weight[i]),
                      paste0("n=",laketrout_Winf$n_Age[i]),
                      paste0("n=",laketrout_Winf$n_Length[i]),
                      ""),
              main = lakenames[i],
              ylab="Est W_inf", xlab="Method")
}
```




## Appendix: Data and data filtering

The raw data consisted of 35,516 measurements from 148 sampling events occurring over 84 lakes  from 1960-2024.  

A few outlying measurements of either length or weight were noted.  

```{r}
morphometry <- read_csv("flat_data/lake_morphometry2.csv", skip=2)
laketrout_all <- read_csv("flat_data/length_weight2.csv", skip=2) %>%
  left_join(morphometry)

laketrout_all %>%
  ggplot(aes(x=ForkLength_mm, y=Weight_g, color=LakeName)) +
  # ggplot(aes(x=ForkLength_mm, y=Weight_g)) +
  geom_point() +
  scale_y_log10() +
  scale_x_log10() +
  theme_bw() +
  theme(legend.position = 'none')
```

A cluster of points in the lower right were best explained by weight having been entered as kilograms rather than grams.  These entries were adjusted according to this assumption.

Obvious outliers and high-leverage points were removed according to the rules: weight > 100,000g, age > 50 years.

Many observations from Paxson Lake showed substantial deviations from the trend.  Residual analysis showed this phenomenon to be related to data from a few specific projects.  It was not clear whether this was an artifact of sampling, or if errors were introduced in data entry or data sorting, etc.  Regardless, the data from these projects were censored from the length-weight analysis.

```{r, fig.height=10}
laketrout_all %>%
  mutate(Weight_g = ifelse(Weight_g < 50, Weight_g*1000, Weight_g)) %>%
  filter(Weight_g < 100000) %>% # is.na(Weight_g) | 
  ggplot(aes(x=ForkLength_mm, y=Weight_g, color=LakeName)) +
  # ggplot(aes(x=ForkLength_mm, y=Weight_g)) +
  geom_point() +
  facet_wrap(facets=~LakeName) +
  scale_y_log10() +
  scale_x_log10() +
  theme_bw() +
  theme(legend.position = 'none')
```

A number of observations still showed large deviations from the trend, likely due to data entry or transcription errors.  These were removed by fitting a global linear regression model, and removing the points with residuals greater than $\pm 4$ times the residual standard deviation, shown below.

```{r}
laketrout1 <- laketrout_all %>%
  mutate(Weight_g = ifelse(Weight_g < 50, Weight_g*1000, Weight_g)) %>%
  filter(is.na(Weight_g) | Weight_g < 100000) %>%
  filter(!ProjectTitle %in% c("Mark-Recapture Event 1 - (September - 2002)",
                              "Mark-Recapture Event 1 - (September - 2003)",
                              "Mark-Recapture Event 1 - (September - 2004)",
                              "Mark-Recapture Event 2 - (May - 2003)"))
  # filter(is.na(ForkLength_mm) | ForkLength_mm )
lm1 <- with(laketrout1, lm(log(Weight_g) ~ log(ForkLength_mm)))
resids1 <- log(laketrout1$Weight_g) - predict(lm1, newdata=laketrout1)
laketrout2 <- filter(laketrout1, is.na(resids1) | abs(resids1) < 4*sd(resids1, na.rm=TRUE))

laketrout <- laketrout2 %>%
  filter(is.na(Age) | Age < 50) %>%
  mutate(LakeNum = as.numeric(as.factor(LakeName)))

# laketrout %>%
#   ggplot(aes(x=ForkLength_mm, y=Weight_g, color=LakeName)) +
#   # ggplot(aes(x=ForkLength_mm, y=Weight_g)) +
#   geom_point() +
#   scale_y_log10() +
#   scale_x_log10(limits=c(100,1200)) +
#   theme_bw() +
#   theme(legend.position = 'none')

outside <- abs(resids1) > 4*sd(resids1, na.rm=TRUE)
plot(laketrout1$ForkLength_mm, laketrout1$Weight_g,
     log="xy", xlim=c(50,1000),
     xlab="Fork Length (mm)", ylab="Weight (g)",
     col=adjustcolor(ifelse(outside, 2, 1), alpha.f=.9),
     pch=ifelse(outside, 16, 1))
legend("topleft", pch=c(1, 16), col=c(1, 2), legend=c("Kept","Removed"))
```

\pagebreak

## Appendix: Model selection for weight ~ length

A log-log regression was used to model weight as a function of length, expanding somewhat on the form used by Lester, et al.  Since fish-level data (lengths and weights) were available from several lakes of differing habitat, the length-weight relationships were fit separately for each lake, with the base model defining slope and intercept parameters as modeled hierarchically, shown below.

$log(W_{j[i]}) = \beta_{0[j]} + \beta_{1[j]}log(L_{j[i]}) + \epsilon_i$

which is equivalent to

$W = e^{\beta_{0[j]}} \times L^{\beta_{1[j]}} \times e^{\epsilon_i}$

in which

$\beta_{0[j]} \sim N(\mu_{\beta_0}, \sigma_{\beta_0}) \text{ for } j \in 1...n_{lakes}$

$\beta_{1[j]} \sim N(\mu_{\beta_1}, \sigma_{\beta_1}) \text{ for } j \in 1...n_{lakes}$

$\epsilon_i \sim N(0,\sigma_{\epsilon})$

This was extended to allow possible relationships between the lake-level intercept parameters $\beta_{0[j]}$ and slope parameters $\beta_{1[j]}$ and both lake area and latitude.  The full model extension is expressed below.

$\beta_{0[j]} \sim N(\mu_{\beta_0[j]}, \sigma_{\beta_0}) \text{ for } j \in 1...n_{lakes}$

$\beta_{1[j]} \sim N(\mu_{\beta_1[j]}, \sigma_{\beta_1}) \text{ for } j \in 1...n_{lakes}$

$\mu_{\beta_0[j]} = \gamma_0 + \tau_{0,Area}log(Area_j) + \tau_{0,Lat}Latitude_j \text{ for } j \in 1...n_{lakes}$

$\mu_{\beta_1[j]} = \gamma_1 + \tau_{1,Area}log(Area_j) + \tau_{1,Lat}Latitude_j \text{ for } j \in 1...n_{lakes}$

$\epsilon_i \sim N(0,\sigma_{\epsilon})$

Models were fit for all possible combinations of relationships, and were compared according to DIC scores, as tabulated below.

```{r, fig.show='hide'}
# loading all the data, but hiding the data exploration plots, etc
source("R/1_laketrout_lwdata.R")

laketrout_lw <- filter(laketrout_lw, !is.na(SurfaceArea_h))

logarea <- log(with(laketrout_lw, tapply(SurfaceArea_h, LakeName, median, na.rm=T)))
lat <- with(laketrout_lw, tapply(Latitude_WGS84, LakeName, median))
lakenames <- names(lat)
```

```{r}
load(file="alloutputs.Rdata")
# alloutputs <- alloutputs1

theDIC <- sapply(alloutputs, function(x) x$DIC)[1:16]
# theDIC
# [1] -3343.813 -3344.528 -3346.021 -3345.798 -3345.114 -3345.044 -3346.590 -3346.775
# [9] -3344.630 -3344.252 -3346.165 -3346.176 -3344.370 -3344.958 -3346.061 -3346.038

# theDIC - min(theDIC)
# [1] 2.9622832 2.2467628 0.7539361 0.9766494 1.6606827 1.7314548 0.1847344 0.0000000
# [9] 2.1452072 2.5227203 0.6104599 0.5986784 2.4050054 1.8175136 0.7138559 0.7370097

# modeldescription[which.min(theDIC)]

controlmat <- expand.grid(0:1,0:1,0:1,0:1)

intmodel <- slopemodel <- rep(NA, nrow(controlmat))
intmodel[controlmat[,1]==0 & controlmat[,2]==0] <- "hierarchical"
intmodel[controlmat[,1]==1 & controlmat[,2]==0] <- "trends with lat"
intmodel[controlmat[,1]==0 & controlmat[,2]==1] <- "trends with area"
intmodel[controlmat[,1]==1 & controlmat[,2]==1] <- "trends with lat & area"
slopemodel[controlmat[,3]==0 & controlmat[,4]==0] <- "hierarchical"
slopemodel[controlmat[,3]==1 & controlmat[,4]==0] <- "trends with lat"
slopemodel[controlmat[,3]==0 & controlmat[,4]==1] <- "trends with area"
slopemodel[controlmat[,3]==1 & controlmat[,4]==1] <- "trends with lat & area"
modeldescription <- paste0("intercept ", intmodel,", ", "slope ", slopemodel)

nparam <- rowSums(controlmat)

# par(mfrow=c(1, 1))
# plot(x=theDIC - min(theDIC), y=rank(theDIC),
#      xlim=c(0, 10))
# text(x=theDIC - min(theDIC), y=rank(theDIC),
#      labels=paste(modeldescription, nparam), pos=4)

# par(mfrow=c(1,1))
# comparecat(alloutputs, p="sig")
# thesigs <- sapply(alloutputs, function(x) jags_df(x, p="sig", exact=T))
# thesigs %>% as.data.frame %>% caterpillar

model_df <- data.frame(Model=modeldescription, delta_DIC = theDIC - min(theDIC),
                       Addl_Param = nparam)
# rownames(model_df) <- modeldescription
highlight <- rep("", nrow(model_df))
bestmodels <- which(model_df$delta_DIC %in% tapply(model_df$delta_DIC, model_df$Addl_Param, min))



###############################################################
bestmodels <- which(theDIC <= min(theDIC)+2)
# bestmodels
## still need to do something equivalent to the .Rmd
###############################################################



highlight[bestmodels] <- " *** "
# model_df$highlight <- highlight

# model_df <- model_df[order(model_df$delta_DIC),]

# knitr::kable(model_df)

par(mfrow=c(1, 1))
plot(x=theDIC - min(theDIC), y=seq_along(theDIC),
     xlim=c(0, 8),
     col=rowSums(controlmat)+1, pch=16,
     ylab="Model", xlab="Delta DIC (smaller is better)",
     ylim=c(16,1))
axis(2, 1:16)
text(x=theDIC - min(theDIC), y=seq_along(theDIC),
     labels=paste0(modeldescription, ", Addl param = ",nparam,highlight), pos=4,
     col=rowSums(controlmat)+1)
```




It is important to compare the differences between inferences in terms of lake-level parameters, among the models considered.  Here the lake-level intercept (b0) and slope (b1) parameters are compared, for the top-performing models (see table below).

```{r}
knitr::kable(model_df[bestmodels,], digits=2)
```

Note that slight differences in inference exist between models, though all differences are very slight.  Dotted lines are overlayed as appropriate, for comparison with the fit reported in Lester et al.

```{r, fig.height=15}
bestoutputs <- alloutputs[bestmodels]

par(mfrow=c(3,1))
comparecat(bestoutputs, p="b0_interp[", col=1:5, main="b0")
abline(h=-19.56, lty=3)
comparecat(bestoutputs, p="b0[", col=1:5, main="b0 (centered)")
comparecat(bestoutputs, p="b1[", col=1:5, main="b1")
abline(h=3.2, lty=3)
```

\pagebreak

Since model fitting was done using recentered data (log length, log area, latitude), some algebraic headache was required to backtransform all parameters back to the data scale, and is not shown here.  Since the net relationship between log length and log weight contains several sources of uncertainty (error in parameter estimation, errors in hierarchical distributions, residual error), error in expressing the length-weight relationship for a new lake is best calculated in-model or directly from the posterior.  However, the formula below should be sufficient for point estimates of regression parameters.

Note: weight is expressed in kilograms, length in millimeters, area in hectares, and latitude in degrees.

```{r, results='asis'}
g0 <- alloutputs[[7]]$q50$b0_int
g1 <- alloutputs[[7]]$q50$b1_int
t0A <- alloutputs[[7]]$q50$b0_area
t1L <- alloutputs[[7]]$q50$b1_lat
meanlogA <- mean(log(laketrout_lw$SurfaceArea_h))
meanlogL <- mean(log(laketrout_lw$ForkLength_mm))
meanlat <- mean(laketrout_lw$Latitude_WGS84)

p1 <- g0 - t0A*meanlogA - g1*meanlogL + t1L*meanlat*meanlogL
p2 <- t0A
p3 <- -t1L*meanlogL
p4 <- g1 - t1L*meanlat
p5 <- t1L

# cat("$$\\hat{log(W)}=", round(p1, 2), "+",
#     round(p2, 3), "\\times log(area) +",
#     round(p3, 3), "\\times lat +",
#     "(", round(p4, 2), "+", round(p5, 3), "\\times lat)\\times log(L)", "$$")

cat("$$\\hat{log(W)}=", round(p1, 2), 
    round(p2, 3), "log(area)+",
    round(p3, 3), "lat +",
    "(", round(p4, 2), round(p5, 3), " lat)log(L)", "$$")
```

The predictive power of the top models was also compared using a cross-validation routine, in which all considered models were re-fit *excluding* the data for each lake in turn, and the model-predicted length-weight relationships were compared to the relationships actually observed at each lake.

This was done two ways: First, the predicted intercept and slope parameters (b0 and b1) were compared to the respective parameters fit from data.  Since inferences for different models *excluding* each lakes' data were compared to  different models *including* each lakes' data, consensus was sought among the ensemble of possible versions of "truth".

```{r}
load(file="lw_loocv.Rdata")
b0_medians <- apply(b0_arr, 2:3, median)
b1_medians <- apply(b1_arr, 2:3, median)

par(mfrow=c(2,1))

comparecat(list(as.data.frame(b0_interp_arr[,,1]),
                as.data.frame(b0_interp_arr[,,2]),
                as.data.frame(b0_interp_arr[,,3]),
                as.data.frame(b0_interp_arr[,,4]),
                as.data.frame(b0_interp_arr[,,5])), col=1:5, 
           main="b0", xlab="lake")
for(i in seq_along(bestmodels)) lines(alloutputs[[bestmodels[i]]]$q50$b0_interp, lty=1)
# legend("topleft", legend = c("from Data", paste("Mod", bestmodels)), 
#        lwd=c(1, rep(3,5)), lty=c(1, rep(1,5)), col=c(1,1:5), pch=c(NA,rep(16,5)))
legend("topleft", legend = c("from Data", "Model prediction"), 
       lwd=c(1, 3), lty=c(1, 1), col=c(1,1), pch=c(NA,16))

comparecat(list(as.data.frame(b1_arr[,,1]),
                as.data.frame(b1_arr[,,2]),
                as.data.frame(b1_arr[,,3]),
                as.data.frame(b1_arr[,,4]),
                as.data.frame(b1_arr[,,5])), col=1:5, 
           main="b1", xlab="lake")
for(i in seq_along(bestmodels)) lines(alloutputs[[bestmodels[i]]]$q50$b1, lty=1)
# legend("topleft", legend = c("from Data", paste("Mod", bestmodels)), 
#        lwd=c(1, rep(3,5)), lty=c(1, rep(1,5)), col=c(1,1:5), pch=c(NA,rep(16,5)))
legend("topleft", legend = c("from Data", "Model prediction"), 
       lwd=c(1, 3), lty=c(1, 1), col=c(1,1), pch=c(NA,16))
```

Root mean squared prediction error is plotted below for predictions from the ensemble of models, as compared to parameters fit from data.  In this case, the "truth" from the ensemble of models is shown (lines and boxplots), but comparison is only made for the top-performing models.

```{r, fig.height=5}
b0_truths <- sapply(alloutputs, function(x) x$q50$b0)
b0_interp_truths <- sapply(alloutputs, function(x) x$q50$b0_interp)
b1_truths <- sapply(alloutputs, function(x) x$q50$b1)

# all models, all truths
b0_rmse <- b0_interp_rmse <- b1_rmse <- matrix(nrow=ncol(b0_truths), ncol=length(bestmodels))
for(itruth in 1:ncol(b0_truths)) {
  for(imodel in seq_along(bestmodels)) {
    b0_rmse[itruth, imodel] <- sqrt(mean((b0_truths[,itruth]-b0_medians[,imodel])^2))
    # b0_interp_rmse[itruth, imodel] <- sqrt(mean((b0_interp_truths[,itruth]-b0_interp_medians[,imodel])^2))
    b1_rmse[itruth, imodel] <- sqrt(mean((b1_truths[,itruth]-b1_medians[,imodel])^2))
  }
}
thecolors <- rcolors(22)
thecolors[bestmodels] <- seq_along(bestmodels)

par(mfrow=c(1,2))
boxplot(b0_rmse, names=bestmodels, xlab="Model used for prediction", main="b0 prediction RMSE")
for(i in 1:ncol(b0_truths)) {
  lines(b0_rmse[i,], col=adjustcolor(thecolors[i], alpha.f=.5), 
        lwd=ifelse(i %in% bestmodels, 2, 1))
}
boxplot(b1_rmse, names=bestmodels, xlab="Model used for prediction", main="b1 prediction RMSE")
for(i in 1:ncol(b0_truths)) {
  lines(b1_rmse[i,], col=adjustcolor(thecolors[i], alpha.f=.5), 
        lwd=ifelse(i %in% bestmodels, 2, 1))
}
# plot(b0_rmse[1,], type="l", ylim=range(b0_rmse))
# plot(b1_rmse[1,], type="l", ylim=range(b1_rmse))
```

Second, the candidate prediction models (that were fit *excluding* the data from each lake in turn) were used to predict log weights from log lengths, using data from the excluded lake.  Model-predicted log weights were compared to actual log-weights for each lake, for each prediction model considered, again in terms of root mean squared prediction error.


```{r, fig.height=5}
lat_all <- with(laketrout_lw, tapply(Latitude_WGS84, LakeName, median))
lakenames_all <- names(lat_all)
pred_r <- pred_sd <- pred_rmse <- matrix(nrow=length(bestmodels), ncol=length(lakenames_all))
for(imodel in seq_along(bestmodels)) {
  for(ilake in seq_along(lakenames_all)) {
    # ypred <- b0_interp_medians[ilake, imodel] +
    #   b1_medians[ilake, imodel]*log(laketrout_lw$ForkLength_mm[laketrout_lw$LakeName==lakenames_all[ilake]])
    ypred <- b0_medians[ilake, imodel] +
      b1_medians[ilake, imodel]*
      (log(laketrout_lw$ForkLength_mm[laketrout_lw$LakeName==lakenames_all[ilake]]) -
         mean(log(laketrout_lw$ForkLength_mm)))
    ytrue <- log(laketrout_lw$Weight_g[laketrout_lw$LakeName==lakenames_all[ilake]]/1000)
    pred_sd[imodel, ilake] <- sd(ypred-ytrue)  # should this be rmse instead??
    pred_rmse[imodel, ilake] <- sqrt(mean(((ypred-ytrue)^2)))
    pred_r[imodel, ilake] <- cor(ypred, ytrue)
  }
}
par(mfrow=c(1,2))
boxplot(t(pred_rmse), names=bestmodels, xlab="Model used for prediction", main="prediction RMSE for all lakes")
apply(pred_rmse, 1, mean) %>% lines(lwd=4)
apply(pred_rmse, 1, mean) %>% plot(main="Mean prediction RMSE error across lakes", type="b", xaxt="n", xlab="Model used for prediction")
axis(side=1, at=seq_along(bestmodels), labels=bestmodels)
```

The results from cross-validation show relatively strong agreement between models, which could perhaps be interpreted as limited benefit from selecting a prediction model that includes lake area and/or latitude.  However, Model 7 (intercept depends on log area, slope depends on latitude) was selected, partly because of the possible inferential benefit in predicting both regression parameters, while still exhibiting strong predictive power on the data level.  Model 7 might be considered a "sweet spot" in predictive power: beyond this number of parameters, the predictive power drops due to overfitting.