What order are shark fisheries exploited in? Brainstorm for Leo's dissertation.

sharkexploitation.R

library(tidyverse)
set.seed(123)

ex_countries <- tibble(
  country = c("Canada", "USA", "Mexico", "Brazil", "Chile"),
  ruleoflaw = c(0.8, 0.7, 0.41, 0.5, 0.66),
  ruleoflawZ = (ruleoflaw - mean(ruleoflaw)) / sd(ruleoflaw)
) 
sharks <- tibble(
  commonname = c("Tiger", "Bull", "Blacktip", "Nurse"),
  # Irschick and Hammerschlag (2015) doi:10.1111/bij.12404
  dorsalarea = c(287, 344, 150, 242),
  dorsalareaZ = (dorsalarea - mean(dorsalarea)) / sd(dorsalarea)
)

# Define coefficients
beta_country <- rnorm(5, 0, 0.25)
names(beta_country) <- ex_countries$country
beta_ruleoflaw <- -0.25
beta_commonname <- rnorm(4, 0, 0.5)
names(beta_commonname) <- sharks$commonname
beta_dorsalareaZ <- 0.5

# Simulate data
ex_grid <- expand_grid(ex_countries, sharks) %>% 
  mutate(eta = beta_country[country] + 
           beta_ruleoflaw * ruleoflawZ +
           beta_commonname[commonname] +
           beta_dorsalareaZ * dorsalareaZ,
         yearstarted = NA)

softmax <- function(x) {
  exp(x) / sum(exp(x))
}
for (t in seq(nrow(ex_grid) - 1)) {
  remaining <- which(is.na(ex_grid$yearstarted))
  remaining_etas <- ex_grid$eta[remaining]
  fishery_probs <- softmax(remaining_etas)
  new_fishery <- sample(remaining, 1, prob = fishery_probs)
  ex_grid$yearstarted[new_fishery] <- t
}
ex_grid$yearstarted[is.na(ex_grid$yearstarted)] <- nrow(ex_grid)

# What's the likelihood function?
L <- function(theta, data) {
  data$eta <- with(data,
                   theta$beta_country[country] + 
                     theta$beta_ruleoflaw * ruleoflawZ +
                     theta$beta_commonname[commonname] +
                     theta$beta_dorsalareaZ * dorsalareaZ)
  negloglik <- 0
  for (t in seq(nrow(data))) {
    remaining <- which(data$yearstarted >= t)
    remaining_etas <- data$eta[remaining]
    remaining_yearstarted <- data$yearstarted[remaining]
    fishery_probs <- softmax(remaining_etas)
    negloglik <- negloglik + -log(fishery_probs[remaining_yearstarted == t])
  }
  names(negloglik) <- NULL
  negloglik
}
L(theta = list(beta_country = beta_country,
               beta_ruleoflaw = beta_ruleoflaw,
               beta_commonname = beta_commonname,
               beta_dorsalareaZ = beta_dorsalareaZ),
  data = ex_grid)

# Max likelihood estimate?
shark_mle <- optim(
  rep(0, 5 + 1 + 4 + 1),
  function(pars) {
    theta <- list(beta_country = c(Canada = pars[1],
                                   USA = pars[2],
                                   Mexico = pars[3],
                                   Brazil = pars[4],
                                   Chile = pars[5]),
                  beta_ruleoflaw = pars[6],
                  beta_commonname = c(Tiger = pars[7],
                                      Bull = pars[8],
                                      Blacktip = pars[9],
                                      Nurse = pars[10]),
                  beta_dorsalareaZ = pars[11])
    L(theta, ex_grid)
  } 
)
tibble(par = c("beta_country_Canada",
               "beta_country_USA",
               "beta_country_Mexico",
               "beta_country_Brazil",
               "beta_country_Chile",
               "beta_ruleoflaw",
               "beta_commonname_Tiger",
               "beta_commonname_Bull",
               "beta_commonname_Blacktip",
               "beta_commonname_Nurse",
               "beta_dorsalareaZ"),
       sim_value = c(beta_country, beta_ruleoflaw, beta_commonname, beta_dorsalareaZ),
       mle_value = shark_mle$par
)

# What order did the dominos fall?
ggplot(ex_grid, aes(yearstarted, eta)) +
  geom_line() +
  geom_text(aes(label = paste(country, commonname))) +
  theme_minimal()
cum_softmax <- function(x) {
  result <- numeric(length(x))
  for (i in seq_along(x)) {
    result[i] <- softmax(x[i:length(x)])[1]
  }
  result
}
ex_grid %>% 
  arrange(yearstarted) %>% 
  mutate(prob = cum_softmax(eta)) %>% 
  ggplot(aes(yearstarted, prob)) +
  geom_line() +
  geom_text(aes(label = paste(country, commonname))) +
  theme_minimal()