5.1_L50.Rmd

---
title: "L50"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_depth: 3
    toc_float:
      collapsed: FALSE
    code_folding: show
    number_sections: TRUE
---

Done extremely quickly without much thought. just exploratory

    
# SET-UP
```{r  markdown, echo=F}
library(knitr)
opts_chunk$set(echo = T, collapse = T, fig.align = "center", fig.width = 12, fig.height = 6)
options(width = 140)
```    

## settings
```{r settings, message = F}
source('0.0_settings.R')
```

## load all data
```{r data, message = F}
load(paste0(dir.rdat, "bio.Rdata"))                                  # from 1.0 read
```

# CALCULATE
```{r calculates, message = F}
table(bio$matur_class,bio$sex)

### subset
bio.l50 <- bio[bio$month %in% 6:7,c('year','nafo','length','agef','matur_class')]
names(bio.l50)[4] <- 'age'

# remove NAs
bio.l50 <- bio.l50[!is.na(bio.l50$matur_class) & !is.na(bio.l50$length) & !is.na(bio.l50$age),]

# remove absurd ages (see 3.0_caa)
bio.l50 <- bio.l50[bio.l50$age<18,]      
bio.l50 <- ddply(bio.l50,c('age'),transform,outlier=outlier(length,coef=3))  
bio.l50[bio.l50$age==0 & bio.l50$length>300,'outlier'] <- TRUE             
table(bio.l50$matur_class,bio.l50$outlier)

bio.l50 <- bio.l50[bio.l50$outlier==FALSE,]
bio.l50$outlier <- NULL

### transform
# cohorts
bio.l50$cohort <- with(bio.l50,year-age)
#mature vs immature
bio.l50$mat <- ifelse(bio.l50$matur_class=='immature',0,1)

### check
table(bio.l50$mat,bio.l50$cohort)

### l50
get.l50 <- function(x){
    cf <- coef(x)
    l <- (log(0.5/(1-0.5))-cf[1])/cf[2]
    names(l) <- 'l50'
    return(l)
}

l50 <- ddply(bio.l50[bio.l50$cohort>=min(bio.l50$year),],c('cohort'),function(x){
    x <- x[!is.na(x$length) &!is.na(x$mat),]
    tab <- c(mat=nrow(x[x$mat==1,]),immat=nrow(x[x$mat==0,]),tot=nrow(x))
    mod <- suppressWarnings(try(glm(mat~length,data=x,family = binomial(link=logit)), silent=TRUE))
    if ('try-error' %in% class(mod)){ 
        ret <- c(l50=NA,'(Intercept)'=NA,length=NA)
    }else{
        ret <- c(get.l50(mod),coef(mod))
        # plot(mat~length,data=x)
        # l <- seq(min(x$length),max(x$length),1)
        # m <-predict(mod,list(length=l),type='response')
        # lines(l, m,col='red',lwd=3)
        # abline(h=0.5,col='grey',lty=3)
        # abline(v=ret[1],col='grey',lty=3)
    }
    return(c(ret,tab))
})

```
## PLOT

Something is off. too low. 
Test an approach like mat~length+.... doy/gear/zone
```{r l50 plot, message = F,fig.width=12}
ggplot(l50,aes(x=cohort,y=l50))+
  geom_point()+
  geom_text(aes(label=tot,y=Inf),vjust=3,size=2)+
  geom_text(aes(label=immat,y=Inf),vjust=6,size=2)+
  geom_segment(aes(x=1973,xend=2013,y=250,yend=250),col='darkred')+
  geom_segment(aes(x=2014,xend=2020,y=263,yend=263),col='darkred')+
  scale_y_continuous(limits=c(0,350),expand=c(0,0))
```