In einarhjorleifsson/datrasdoodle: Does not matter.

SR model {#models}

... Something about the model

Stock recruitment relationship

library(tidyverse)
library(modelr)
library(broom)

rby <-
  fishvice::read_separ("/net/hafkaldi/export/u2/reikn/hoski/Mackerel/HCRSimulations/2periods/",
                       "HockeyAcf00",
                       run = "HockeyEstAcf",
                       mName = "sep")$rby %>% 
  select(year, rec = r, ssb) %>% 
  filter(year < 2016) %>% 
  as_tibble() %>% 
  mutate(ssb = ssb / 1e3,
         rec = rec / 1e3)

# segreg
segreg <- function(df, a, b) {
  log(ifelse(df$ssb >= b, a * b, a * df$ssb))
}

initial <- function(df) {
    c(a = log(stats::median(df$rec/df$ssb, na.rm = TRUE)), 
      b = log(stats::median(df$ssb)),
      cv = 0)
  }

nllik <- function(param, data, model) {

    pred <- segreg(data, a = exp(param[1]), b = exp(param[2]))
    -1 * sum(stats::dnorm(log(data$rec), pred, exp(param[3]), log = TRUE))

}

fit <- stats::nlminb(initial(rby), objective = nllik, data = rby, model = segreg)


#segreg <- function(x) {ifelse(x >= b, a * b, a * x)}
#p <- 
#  ggplot(rby, aes(x = ssb, y = rec)) +
#  geom_point() +
#  expand_limits(x = 0, y = 0)
#b <- min(rby$ssb)
#a <- median(rby$rec)/b
#p + stat_function(fun = segreg)

x <- 
  expand.grid(b = seq(0, 5, by = 0.01),
              a = seq(1, 15, by = 0.01),
              ssb = rby$ssb) %>% 
  left_join(rby) %>% 
  mutate(pY = ifelse(ssb >= b, a * b, a * ssb),
         res = log(rec) - log(pY),
         res2 = res^2) %>% 
  group_by(a, b) %>% 
  summarise(sse = sum(res2),
            nllik = -sum(dnorm(log(rec), log(pY), sd(res), log = TRUE)))
library(viridis)
x3 <-
  x %>% 
  ungroup(a, b) %>% 
  filter(nllik == min(nllik))
sbreak <- exp(fit$par[[2]])
slope <- exp(fit$par[[1]])
x %>% 
  mutate(nllik = ifelse(nllik > 30, NA, nllik)) %>% 
  ggplot(aes(b, a, fill = nllik)) +
  geom_raster() +
  scale_fill_viridis(option = "B", direction = -1) +
  geom_vline(xintercept = sbreak) +
  geom_hline(yintercept = slope) +
  geom_point(data = x3, aes(a, b), colour = "blue")

# Following also leads to a false minimum
initial <- function(df) {
    c(a = log(stats::median(df$rec) / min(df$ssb)), 
      b = log(min(df$ssb)),
      cv = 0)
}
fit <- stats::nlminb(initial(rby), objective = nllik, data = rby, model = segreg)
sbreak <- exp(fit$par[[2]])
slope <- exp(fit$par[[1]])
x %>% 
  mutate(nllik = ifelse(nllik > 30, NA, nllik)) %>% 
  ggplot(aes(b, a, fill = nllik)) +
  geom_raster() +
  scale_fill_viridis(option = "B", direction = -1) +
  geom_vline(xintercept = sbreak) +
  geom_hline(yintercept = slope) +
  geom_point(data = x3, aes(a, b), colour = "blue")

# lets try husky
hockey <- function(df, rmax, sbreak) {
  log(ifelse(df$ssb >= rmax, rmax, rmax / sbreak * df$ssb))
}
initial <- function(df) {
    c(sbreak = log(min(df$ssb)),
      rmax = log(median(df$rec)/min(df$ssb)), 
      cv = 0)
}
fit <- stats::nlminb(initial(rby), objective = nllik, data = rby, model = hockey)
exp(fit$par)
sbreak <- exp(fit$par[[1]])
rmax   <- exp(fit$par[[2]])
x %>% 
  mutate(nllik = ifelse(nllik > 30, NA, nllik)) %>% 
  ggplot(aes(b, a, fill = nllik)) +
  geom_raster() +
  scale_fill_viridis(option = "B", direction = -1) +
  geom_vline(xintercept = sbreak) +
  geom_hline(yintercept = rmax) +
  geom_point(data = x3, aes(a, b), colour = "blue")
# problem??: the message is "false convergence (8)"

# -------------------------------------------------------------------------
# Try a smooth damper
llik <-  function (param, data, logpar = FALSE) 
{
    if (logpar) {
        pred <- Segreg(list(a = exp(param[1]), b = exp(param[2])), 
            data$ssb)
        sum(dnorm(log(data$rec), pred, exp(param[3]), log = TRUE))
    }
    else {
        pred <- Segreg(list(a = param[1], b = param[2]), data$ssb)
        sum(dnorm(log(data$rec), pred, param[3], log = TRUE))
    }
}

nllik <- function(param, ...) -1 * llik(param, ...)

SmoothDamper1 <- function(x,Roof,Floor) {
  deltax <- 0.01
  if(Roof == Floor) return(x) 
  lb  <-  1.0 - deltax/2.0
  ub <- 1.0 + deltax/2.0
  if(x <= lb* Roof & x >= ub*Floor) return(x)
  if(x >= ub*Roof) return(Roof)
  if(x <= lb*Floor) return(Floor)
  if(x <= ub*Roof && x >= lb*Roof) {
    y <-  (x - ub*Roof);
    return(Roof - 0.5/deltax/Roof*y*y)
  }
  if(x >= lb*Floor && x <= ub*Floor) {
    y <-  (x - lb*Floor);
    return(Floor +0.5/deltax/Floor*y*y)
  }
}

SmoothDamper <- function(x,Roof,Floor) {
  x1 <- x*0
  for(i in 1:length(x)) x1[i] <- SmoothDamper1(x[i],Roof,Floor)
  return(x1)
}

Segreg <- function(ab,ssb){
  pr <- ab$a*ssb
  pr <- SmoothDamper(pr,ab$b*ab$a,0)
  return(log(pr))
}

rby <-
  fishvice::read_separ("/net/hafkaldi/export/u2/reikn/hoski/Mackerel/HCRSimulations/2periods/",
                       "HockeyAcf00",
                       run = "HockeyEstAcf",
                       mName = "sep")$rby %>% 
  select(year, rec = r, ssb) %>% 
  filter(year < 2016) %>% 
  as_tibble()


ndat <- nrow(rby)
initial <- function(df) {
    c(a = log(stats::median(df$rec) / min(df$ssb)), 
      b = log(min(df$ssb)),
      cv = 0)
}
fit <- nlminb(initial(rby), nllik, data = rby, logpar = TRUE)
exp(fit$par)

rby <-
  rby %>% 
  mutate(rec = rec/1e3,
         ssb = ssb/1e3)
fit <- nlminb(initial(rby), nllik, data = rby, logpar = TRUE)
exp(fit$par)
exp(initial(rby))
# However, easy to get wrong minimum
initial <- function(df) {
    c(a = log(stats::median(df$rec/df$ssb, na.rm = TRUE)), 
      b = log(stats::median(df$ssb)),
      cv = 0)
}
fit <- nlminb(initial(rby), nllik, data = rby, logpar = TRUE)