##Mortality estimates of the Chen & Watanabe (1989) model###
###Function to be used##############
############################################################
##########################################################
######################################################## 
MChen <- function(x, kmin, kmax, Linfmin, Linfmax, t0min, t0max, tmmin, tmmax, n, tm){
  k <- runif(n, min = kmin, max = kmax)
  t0 <- runif(n, min = t0min, max = t0max)
  Linf <- runif(n, min = Linfmin, max = Linfmax)
    a0 <- 1-exp(-k*(tm-t0))
  a1 <- k*exp(-k*(tm-t0))
  a2 <- -0.5*k^2*(exp(-k*(tm-t0)))
  if(x <= tm){
    M <- k/(1-exp(-k*(x-t0)))
  }
  if(x > tm){
    M <- (k/(a0)+(a1*(x- tm))+a2*(x- tm)^2)
  }
  M <- subset(M, M > 0) #Only positive results are viable
  return(M)
 }

 
 ####Input parameters####
kmin <- 0.2   #minimum coefficient growth from von Bertalanffy model
kmax <- 0.7 #Maximum growth coefficient from von Bertalanffy model
Linfmin <- 23.5 #Minimum Asymptotic length coefficient from von Bertalanffy model
Linfmax <- 38 ##Maximum Asymptotic length coefficient from von Bertalanffy model
t0min <- -0.8 ##Minimum t0 parameter from von Bertalanffy model
t0max <- -0.4  #Maximum t0 parameter from von Bertalanffy model
n <- 1000  #Number of samples


#####For Early maturation (relative age = 1.5 years)######
tm1.5 <- 1.5  ##Age at maturity (first scenario)
 
age <- 0:10 ###vector of relative ages contemplated


#####Loop to estimate M at-age
M1.5 <- list()
for(i in 1:11){
M1.5[[i]] <- MChen(age[i], kmin = kmin, kmax = kmax, Linfmin = Linfmin, Linfmax = Linfmax, t0min = t0min, t0max = t0max, n = n, tm = tm1.5)
M1.5
} 
 
 
##############Create a data frames with M estimations at-age
list1.5 <- list()
for(i in 1:11){
list1.5[[i]] <- data.frame(rep(age[i], length(M1.5[[i]])), M1.5[[i]])
list1.5
}

######merge all the estimates into a single table#####
predictData1.5 <- do.call(rbind.data.frame, list1.5)
names(predictData1.5) <- c("Age", "M")





#####For Late maturation (relative age = 5 years)######
tm5 <- 5 #Age at maturity (second scenario)

#####Loop to estimate M at-age
M5 <- list()
for(i in 1:11){
M5[[i]] <- MChen(age[i], kmin = kmin, kmax = kmax, Linfmin = Linfmin, Linfmax = Linfmax, t0min = t0min, t0max = t0max, n = n, tm = tm5)
M5
} 
 
 
 ##############Create a data frames with M estimations at-age
list5 <- list()
for(i in 1:11){
list5[[i]] <- data.frame(rep(age[i], length(M5[[i]])), M5[[i]])
list5
}

######merge all the estimates into a single table#####
predictData5 <- do.call(rbind.data.frame, list5)
names(predictData5) <- c("Age", "M")



######merge both scenarios into a single table#####
predictData1.5$Stage <- rep("Early Maturation", length(predictData1.5$M))
predictData5$Stage <- rep("Late Maturation", length(predictData5$M))
dato <- rbind(predictData1.5, predictData5)


require(ggplot2) ###package needed

##create Figure 2 of manuscript
p <- ggplot(dato, aes(x = Age, y = M, group = Age)) + facet_wrap(~Stage) + geom_boxplot(fill = "gray")
p + theme_bw(15) + xlab("Relative age (years)") + ylab(expression("Natural Mortality"~~(yr^-1))) +
scale_y_continuous(breaks = seq(0, 3, 0.2), labels = seq(0, 3, 0.2)) +
scale_x_continuous(breaks = seq(0, 10, 1), labels = seq(0, 10, 1))