growth_curve.Rmd

---
title: "fit_growth_curve"
author: "Kosmas Hench"
date: "5/24/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(DT)
```

```{r, warning = FALSE, message = FALSE}
# code from https://stackoverflow.com/questions/14190883/fitting-a-curve-to-specific-data
# function from http://www.pisces-conservation.com/growthhelp/index.html?logistic_curve.htm
library(tidyverse)
library(tidymodels)
```


```{r}
data <- tibble(x = c(225, 225, 225, 225, 240, 240, 240, 240, 
                     245, 245, 245, 245, 185, 185, 185, 185, 190, 190, 190, 190, 243, 
                     243, 243, 243, 258, 258, 258, 258, 215, 215, 215, 215, 256, 255, 
                     256, 256, 147, 146, 147, 147, 244, 244, 244, 244, 175, 175, 175, 
                     175, 243, 242, 243, 243, 262, 262, 262, 262, 156, 156, 156, 156, 
                     268, 268, 268, 268, 153, 153, 153, 153, 134, 134, 134, 134, 128, 
                     128, 128, 128, 252, 251, 252, 252, 168, 168, 168, 168, 186, 185, 
                     186, 186, 261, 261, 261, 261, 209, 209, 209, 209, 133, 133, 133, 
                     133, 188, 188, 188, 188, 193, 193, 193, 193, 181, 181, 181, 181, 
                     220, 219, 220, 220, 162, 162, 162, 162, 220, 220, 220, 220, 267, 
                     267, 267, 267, 164, 164, 164, 164, 191, 191, 191, 191, 247, 247, 
                     247, 247, 166, 165, 166, 166, 143, 142, 143, 143, 179, 178, 179, 
                     179, 213, 213, 213, 213, 93, 93, 93, 93, 243, 242, 243, 243, 
                     244, 244, 244, 244, 89, 89, 89, 89, 120, 120, 120, 120, 165, 
                     165, 165, 165, 196, 196, 196, 196, 157, 157, 157, 157, 127, 127, 
                     127, 127, 175, 174, 175, 175, 141, 141, 141, 141, 141, 140, 141, 
                     141, 167, 167, 167, 167, 139, 139, 139, 139, 172, 172, 172, 172, 
                     117, 116, 117, 117, 129, 129, 129, 129, 172, 172, 172, 172, 97, 
                     97, 97, 97, 156, 156, 156, 156, 100, 100, 100, 100, 198, 198, 
                     198, 198, 188, 188, 188, 188, 200, 200, 200, 200, 192, 192, 192, 
                     192, 203, 202, 203, 203, 186, 185, 186, 186, 180, 180, 180, 180, 
                     134, 133, 134, 134, 201, 201, 201, 201, 130, 130, 130, 130, 146, 
                     146, 146, 146, 107, 107, 107, 107, 171, 171, 171, 171, 221, 221, 
                     221, 221, 134, 134, 134, 134, 186, 185, 186, 186, 119, 119, 119, 
                     119, 149, 149, 149, 149, 195, 195, 195, 195, 224, 224, 224, 224, 
                     148, 148, 148, 148, 204, 204, 204, 204, 214, 214, 214, 214, 224, 
                     224, 224, 224, 224, 224, 224, 224, 187, 187, 187, 187, 204, 204, 
                     204, 204, 126, 126, 126, 126, 131, 130, 131, 131, 171, 171, 171, 
                     171, 199, 199, 199, 199, 147, 147, 147, 147, 122, 122, 122, 122, 
                     116, 116, 116, 116, 155, 155, 155, 155, 138, 138, 138, 138, 157, 
                     157, 157, 157, 88, 88, 88, 88, 90, 90, 90, 90, 110, 109, 110, 
                     110, 70, 70, 70, 70, 147, 146, 147, 147, 119, 119, 119, 119, 
                     156, 155, 156, 156, 76, 76, 76, 76, 132, 132, 132, 132, 113, 
                     113, 113, 113, 135, 135, 135, 135, 77, 77, 77, 77, 50, 49, 50, 
                     50, 101, 100, 101, 101, 129, 129, 129, 129, 86, 86, 86, 86, 114, 
                     114, 114, 114, 91, 91, 91, 91, 79, 79, 79, 79, 58, 57, 58, 58, 
                     145, 144, 145, 145, 128, 127, 128, 128, 132, 132, 132, 132),
               y = c(163, 163, 163, 163, 160, 160, 160, 160, 154, 154, 154, 
                     154, 151, 151, 151, 151, 149, 149, 149, 149, 144, 144, 144, 
                     144, 137, 137, 137, 137, 136, 136, 136, 136, 134, 134, 134, 
                     134, 128, 128, 128, 128, 127, 127, 127, 127, 126, 126, 126, 
                     126, 125, 125, 125, 125, 123, 123, 123, 123, 122, 122, 122, 
                     122, 122, 122, 122, 122, 122, 122, 122, 122, 122, 122, 122, 
                     122, 121, 121, 121, 121, 121, 121, 121, 121, 118, 118, 118, 
                     118, 117, 117, 117, 117, 116, 116, 116, 116, 115, 115, 115, 
                     115, 115, 115, 115, 115, 112, 112, 112, 112, 111, 111, 111, 
                     111, 111, 111, 111, 111, 111, 111, 111, 111, 110, 110, 110, 
                     110, 110, 110, 110, 110, 110, 110, 110, 110, 108, 108, 108, 
                     108, 107, 107, 107, 107, 106, 106, 106, 106, 104, 104, 104, 
                     104, 103, 103, 103, 103, 103, 103, 103, 103, 102, 102, 102, 
                     102, 101, 101, 101, 101, 101, 101, 101, 101, 100, 100, 100, 
                     100, 98, 98, 98, 98, 93, 93, 93, 93, 93, 93, 93, 93, 93, 
                     93, 93, 93, 91, 91, 91, 91, 90, 90, 90, 90, 90, 90, 90, 90, 
                     86, 86, 86, 86, 85, 85, 85, 85, 85, 85, 85, 85, 83, 83, 83, 
                     83, 75, 75, 75, 75, 73, 73, 73, 73, 70, 70, 70, 70, 62, 62, 
                     62, 62, 61, 61, 61, 61, 59, 59, 59, 59, 55, 55, 55, 55, 109, 
                     109, 109, 109, 104, 104, 104, 104, 104, 104, 104, 104, 103, 
                     103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 103, 102, 
                     102, 102, 102, 99, 99, 99, 99, 98, 98, 98, 98, 95, 95, 95, 
                     95, 94, 94, 94, 94, 94, 94, 94, 94, 93, 93, 93, 93, 92, 92, 
                     92, 92, 87, 87, 87, 87, 86, 86, 86, 86, 86, 86, 86, 86, 85, 
                     85, 85, 85, 84, 84, 84, 84, 83, 83, 83, 83, 81, 81, 81, 81, 
                     79, 79, 79, 79, 79, 79, 79, 79, 77, 77, 77, 77, 76, 76, 76, 
                     76, 75, 75, 75, 75, 67, 67, 67, 67, 67, 67, 67, 67, 65, 65, 
                     65, 65, 65, 65, 65, 65, 64, 64, 64, 64, 64, 64, 64, 64, 64, 
                     64, 64, 64, 62, 62, 62, 62, 61, 61, 61, 61, 61, 61, 61, 61, 
                     58, 58, 58, 58, 58, 58, 58, 58, 56, 56, 56, 56, 55, 55, 55, 
                     55, 52, 52, 52, 52, 50, 50, 50, 50, 50, 50, 50, 50, 49, 49, 
                     49, 49, 48, 48, 48, 48, 47, 47, 47, 47, 46, 46, 46, 46, 44, 
                     44, 44, 44, 41, 41, 41, 41, 40, 40, 40, 40, 38, 38, 38, 38, 
                     35, 35, 35, 35, 35, 35, 35, 35, 34, 34, 34, 34, 31, 31, 31, 
                     31, 28, 28, 28, 28, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 
                     27, 27, 25, 25, 25, 25),
               sex = rep(c('m','f'), each = 240))
```

```{r}
(p1 <- ggplot(data, aes(x, y,
                     fill = factor(sex), 
                     color = factor(sex), 
                     group = sex)) + 
  geom_point(shape = 21, color = 'black') +
  labs( x = 'Age (days)', y = 'Size (mm)' )+
  scale_color_brewer('Sex', palette = 'Set1') +
  scale_fill_brewer('Sex', palette = 'Set1') +
  scale_linetype(guide = FALSE)+
  theme_minimal() +
  theme(legend.position = c(.95,.1)))
```

```{r, warning = FALSE, message = FALSE}
lmax <- 160
```

**simple linear model**

$$l_{t} = at + b$$

```{r}
model_lm <- function(tib, ...){
   lm(y ~ x, 
      data = tib)
}
```


**von Bertalanffy growth**

$$l_{t} = L_{\infty}( 1 - e^{-K(t -t_{0})})$$

- $l$ is length (or some other measure of size, `y` in our data set)
- $t$ is age/ time (`x` in our data set)
1. $L_{\infty}$ termed *$L$ infinity* in fisheries science, is the asymptotic length at which growth is zero (`L`/`lmax`)
2. $K$ is the growth rate (represented by `k` in the equation and set by `b_k`)
3. $t_{0}$ is the initial size of the organism (`t0`/`b_t0`)

```{r, warning = FALSE, message = FALSE}
b_k <- .005
b_t0 <- 0

model_bert <- function(tib, L = lmax){
   nls(y ~ I( L * (1 - exp(-k*(x-t0)))), 
      data = tib,
      start = list( L = lmax, k = b_k, t0 = b_t0))
}
```

**Logistic groth**

$$l_{t} = \frac{L_{\infty}}{1 + e^{-k(t-I)}})$$

- $l$ is length (or some other measure of size, `y` in our data set)
- $t$ is age/ time (`x` in our data set)
1. $L_{\infty}$ termed *$L$ infinity* in fisheries science, is the asymptotic length at which growth is zero (`L`/`lmax`)
2. $k$ is the growth rate (represented by `k` in the equation and set by `l_k`)
3. $I$ is the age at the inflection point (`i`/`l_i`)

```{r, warning = FALSE, message = FALSE}
l_k <- .01
l_i <- 170


model_logistic <- function(tib, L = lmax, ...){
  nls(y ~ I( L / (1 + exp(-k*(x-i)))) , 
      data = tib, 
      start = list( L = L, k = l_k, i = l_i),
      control = list(minFactor = .001, maxiter = 500))
}
```


```{r,fig.asp=.5}
p1 + 
  stat_smooth(method = lm, data = data %>% mutate(mod = 'Linear'), 
              se = FALSE,
              formula = as.formula(model_lm(tib = data %>% filter(sex == 'm'))),
              aes(linetype = 'linear')) +
  stat_smooth(method = nls, data = data %>% mutate(mod = 'Logistic'), 
              se = FALSE,
              formula = as.formula(model_logistic(tib = data %>% filter(sex == 'm'))),
              aes(linetype = 'logistic'),
              method.args = list(start = list( L = lmax, k = l_k, i = l_i)))+
  stat_smooth(method = nls, data = data %>% mutate(mod = 'von Bertalanffy'), 
              se = FALSE,
              formula = as.formula(model_bert(tib = data %>% filter(sex == 'm'))),
              aes(linetype = 'von Bertalanffy'),
              method.args = list(start = list( L = lmax, k = b_k, t0 = b_t0)))+
  facet_grid(.~mod)
```


```{r, warning = FALSE, message = FALSE}
fit_model <- function(data, model_fun,model_name, sex_querry, ...){

  A <- data %>%
  filter(sex == sex_querry) %>% # subset males
  model_fun(., ...) %>% AIC()
    
  data %>%
  filter(sex == sex_querry) %>% # subset males
  model_fun(., ...) %>%
  tidy() %>%
  mutate(model = model_name,
           sex = sex_querry,
         aic = A)
}
```


```{r}
tibble(model_fun = rep(c(model_lm,model_bert,model_logistic), each = 2),
       model_name = rep(c('Linear','von Bertalanffy','Logistic'), each = 2),
       sex_querry = rep(c('m','f'),3)
       ) %>%
  purrr::pmap(fit_model, data = data) %>%
  bind_rows() %>% 
  knitr::kable(align = 'crrrrrcr')
```