In datalimited/datalimited: Stock Assessment Methods for Data-limited Fisheries

Test the R package code against the original code in Dropbox.

Load the package functions:

library(dplyr)
devtools::load_all("..")

source("FunctionsCOM.R")
ct <- blue_gren$ct

set.seed(1)
system.time({
orig_comsir_priors <- Priors(Nsim = 1e4, 
  MyPar = list(k = 800, r = 0.6, x = 0.5, a = 0.8), Catch = ct, 
  LogisticModel = TRUE, Obs = FALSE, start.r = c(0.2,1.0), NormK = FALSE, Normr = FALSE, Norma = FALSE, Normx = FALSE, logK = TRUE)
})

set.seed(1)
system.time({
pkg2_comsir_priors <- comsir_priors(ct = ct,
   k = 800, r = 0.6, x = 0.5, a = 0.8, start_r = c(0.2, 1),
   mink = max(ct), maxk = max(ct) * 100, nsim = 1e4, logistic_model = TRUE, obs = FALSE, norm_k = FALSE, norm_r = FALSE, 
  logk = TRUE, norm_x = FALSE, norm_a = FALSE)
})
p_pkg <- pkg2_comsir_priors
p_orig <- orig_comsir_priors

p_pkg <- p_pkg[p_pkg$like != 0, ]
p_orig <- p_orig[p_orig$Like != 0, ]
stopifnot(identical(nrow(p_pkg), nrow(p_orig)))
stopifnot(identical(p_pkg$biomass, p_orig$B))
stopifnot(identical(p_pkg$prop, p_orig$Prop))
stopifnot(identical(p_pkg$n1, p_orig$N1))
stopifnot(identical(as.numeric(as.matrix(p_pkg)), as.numeric(as.matrix(p_orig))))

Compare:

# par(mfrow = c(10, 2), cex = 0.5)
# for(i in 1:ncol(pkg_comsir_priors)) {
  # hist(pkg2_comsir_priors[,i], main = paste("pkg", names(pkg_comsir_priors)[i]))
  # hist(orig_comsir_priors[,i], main = paste("orig", names(pkg_comsir_priors)[i]))
# }

Looks good now.

Check the "estimation" phase:

set.seed(1)
system.time({est_orig <- Estimate(p_orig, Catch = ct, CV = 0.4, 
  LogisticModel = TRUE, NormalL = TRUE)})
set.seed(1)
system.time({est_pkg <- comsir_est(n1 = p_pkg$n1, k = p_pkg$k, r = p_pkg$r, a = p_pkg$a, 
  x = p_pkg$x, h = p_pkg$h, z = p_pkg$z, like = p_pkg$like, ct = ct, 
  logistic_model = TRUE, normal_like = TRUE, cv = 0.4)})

identical(est_orig$B, est_pkg$B) %>% stopifnot
identical(est_orig$K, est_pkg$k) %>% stopifnot
identical(est_orig$Prop, est_pkg$prop) %>% stopifnot
identical(est_orig$Like, est_pkg$like) %>% stopifnot

Check the resampling phase:

set.seed(1)
r_pkg <- with(est_pkg, comsir_resample(k = k, r = r, a = a, x = x, h = h, 
  like = like, yr = blue_gren$yr, n_posterior = 10L, ct = blue_gren$ct, logistic_model = TRUE))
set.seed(1)
r_orig <- MYsumpars(ParVals = est_orig, Npost = 10L, Catch = blue_gren$ct, LogisticModel = TRUE)

r_orig$posterior %>% head
r_pkg$posterior %>%  head
identical(r_pkg$like, r_orig$like) %>% stopifnot
identical(r_pkg$Biomass, r_orig$B) %>% stopifnot

Check effort dynamics code:

e_pkg <- effortdyn(h = 0, k = 1000, r = 0.5, x = 0.4, a = 0.9, yr = 1:3, 
  ct = c(500, 200, 305), logistic_model = TRUE)
e_orig <- effortdyn_orig(TrueP = c(h = 0, K = 1000, r = 0.5, x = 0.4, a = 0.9), 
  LogisticModel = TRUE, Catch = c(500, 200, 305), What = 1:5)
identical(as.numeric(e_pkg[,-6]), as.numeric(e_orig)) %>% stopifnot()

Checking the whole function:

set.seed(123)
system.time({x <- comsir(ct = blue_gren$ct, yr = blue_gren$yr, k = 800, r = 0.6, 
  nsim = 1e5, a = 0.8, x = 0.5, n_posterior = 5e2, start_r = c(0.2, 1.0), 
  logistic_model = TRUE, obs = FALSE, norm_k = FALSE, norm_r = FALSE, 
  logk = TRUE, norm_x = FALSE, norm_a = FALSE, normal_like = FALSE)})

# x <- qq
bbmsy <- reshape2::dcast(x$quantities, sample_id ~ yr,
    value.var = "bbmsy")[,-1] # convert long to wide format
bbmsy_out <- summarize_bbmsy(bbmsy, log = TRUE)
bbmsy_out$year <- seq_len(nrow(bbmsy_out))
library("ggplot2")
ggplot(bbmsy_out, aes(year, bbmsy_q50)) + geom_line() + geom_ribbon(aes(ymin = bbmsy_q25, ymax = bbmsy_q75), alpha = 0.2) + geom_hline(yintercept = 1, lty = 2)

Compare the original code:

set.seed(123)
system.time({out <- DoProject(TruePar = as.data.frame(t(c(N1 = 800, K = 800, 
  r = 0.6, z = 1.0, a = 0.8, x = 0.5, h = 0))), MyData = data.frame(catch = ct), 
  Nsim = 1e5, Npost = 5e2, LogisticModel = TRUE, Obs = FALSE, EstLogisticM = TRUE, start.r = c(0.2,1.0), NormK = FALSE, Normr = FALSE, Norma = FALSE, Normx = FALSE, logK = TRUE)})

bbmsy_orig_out <- summarize_bbmsy(t(out$BoverBmsy), log = TRUE)
bbmsy_orig_out$year <- seq_len(nrow(bbmsy_orig_out))
ggplot(bbmsy_orig_out, aes(year, bbmsy_q50)) + geom_line() + geom_ribbon(aes(ymin = bbmsy_q25, ymax = bbmsy_q75), alpha = 0.2) + geom_hline(yintercept = 1, lty = 2)

plot(bbmsy_out$bbmsy_q50, bbmsy_orig_out$bbmsy_q50);abline(a = 0, b = 1)

Check another stock where the likelihood is getting too close to zero:

set.seed(1)
system.time({orig <- DoProject(TruePar = as.data.frame(t(c(N1 = 800, K = 800, 
  r = 0.6, z = 1.0, a = 0.8, x = 0.5, h = 0))), MyData = data.frame(catch = gom$c_touse), 
  Nsim = 1e4, Npost = 5e2, NormalL=F, LogisticModel = TRUE, Obs = FALSE, EstLogisticM = TRUE, start.r = c(0.2,1.0), NormK = FALSE, Normr = FALSE, Norma = FALSE, Normx = FALSE, logK = TRUE)})

Plot the two:

names(out)
names(x)

median(x$posterior$r)
median(out$posterior$r)

median(x$quantities$bbmsy)
median(out$BoverBmsy)

po <- readRDS("priors-orig.rds")
pn <- readRDS("prior.rds")
median(pn$biomass)
median(po$B)
median(pn$k)
median(po$K)
median(pn$prop)
median(po$Prop)

eo <- readRDS("est-orig.rds")
en <- readRDS("est.rds")

median(eo$B)
median(en$B)

median(eo$Prop)
median(en$prop)


# MYsumpars<- function(ParVals,Npost=1000,Catch,Plot=F, LogisticModel=FALSE)



fake_est_orig <- list(N1 = en$n1, K = en$k, r = en$r, z = en$z, a = en$a, x = en$x, h = rep(0, length(length(en$x))), B = en$B, Prop = en$prop, LogLike = en$loglike, Like = en$like)
qq <- MYsumpars(ParVals = fake_est_orig, Npost = 5e3, Catch = blue_gren$ct, LogisticModel = TRUE)

function(k, r, a, x, h, like, yr, ct, n_posterior = 1000L,
  logistic_model = TRUE) {

qq_p <- with(en, comsir_resample(k = k, r = r, a = a, x = en$x, h = h, like = en$like, yr = blue_gren$yr, ct = blue_gren$ct, n_posterior = 5e3))