From 0a560430fef15e2f69fedaacb43d3685ba4722af Mon Sep 17 00:00:00 2001 From: "Quang C. Huynh" Date: Mon, 5 Jun 2023 16:06:40 -0700 Subject: [PATCH 1/4] French labels to plot_maturity_months --- R/mat-freq-month.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/mat-freq-month.R b/R/mat-freq-month.R index 3fad398..33aac74 100644 --- a/R/mat-freq-month.R +++ b/R/mat-freq-month.R @@ -156,7 +156,7 @@ plot_maturity_months <- function(dat, theme_pbs() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + theme(panel.spacing = unit(-0.1, "lines")) + - labs(title = title, colour = "Sex", fill = "Sex") + labs(title = title, colour = en2fr("Sex"), fill = en2fr("Sex")) if (nrow(dat) >= 10) { g <- g + From 1204859a74d3e8da9c4693328990e2bcacc1378f Mon Sep 17 00:00:00 2001 From: "Quang C. Huynh" Date: Wed, 3 Aug 2022 15:13:21 -0700 Subject: [PATCH 2/4] Schnute growth fn --- NAMESPACE | 2 + R/growth.R | 333 +++++++++++++++++++++++++++--- inst/stan/schnute-nopriors.stan | 75 +++++++ inst/stan/schnute.stan | 75 +++++++ inst/tmb/schnute.R | 76 +++++++ inst/tmb/schnute.cpp | 58 ++++++ man/{fit_vb.Rd => fit_schnute.Rd} | 67 +++++- man/plot_growth.Rd | 10 +- 8 files changed, 657 insertions(+), 39 deletions(-) create mode 100644 inst/stan/schnute-nopriors.stan create mode 100644 inst/stan/schnute.stan create mode 100644 inst/tmb/schnute.R create mode 100644 inst/tmb/schnute.cpp rename man/{fit_vb.Rd => fit_schnute.Rd} (62%) diff --git a/NAMESPACE b/NAMESPACE index 9dfc452..6996e60 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -6,6 +6,7 @@ export(bind_samples) export(fit_cpue_index_glmmtmb) export(fit_length_weight) export(fit_mat_ogive) +export(fit_schnute) export(fit_survey_sets) export(fit_vb) export(plot_age_precision) @@ -23,6 +24,7 @@ export(plot_predictor_bubbles) export(plot_qres_histogram) export(plot_qres_qq) export(plot_sample_avail) +export(plot_schnute) export(plot_survey_index) export(plot_survey_sets) export(plot_vb) diff --git a/R/growth.R b/R/growth.R index ec355b3..c8b968f 100644 --- a/R/growth.R +++ b/R/growth.R @@ -1,6 +1,9 @@ -#' Fit a von Bertalanffy growth model +#' Fit a growth model for length-at-age #' -#' For use with data for a single species. +#' For use with data for a single species, fit either a Schnute growth curve or +#' a von Bertalanffy growth curve. `fit_vb(method = "tmb")` will use `fit_schnute()` +#' with `p = 1`. Otherwise, separate Stan files are used. The Richards function +#' is also nested in the Schnute function with `0 < p < 1`. #' #' @param dat Input data frame. Should be from [gfdata::get_survey_samples()] or #' [gfdata::get_commercial_samples()]. @@ -8,6 +11,14 @@ #' @param method `"mpd"` for the mode of the posterior distribution (with #' [rstan::optimizing()]) or `"mcmc"` for full MCMC sampling with Stan (with #' [rstan::sampling()]). `"tmb"` for a TMB model. +#' @param growth Whether to estimate the inflection parameter p (`"schnute"`) or +#' fix the parameter p = 1 (`"vb"`) to reduce to the von Bertalanffy equation. +#' @param a1 The first age parameter in the Schnute function, the corresponding +#' length `"L1"` will be estimated. You may want to specify the smallest age +#' class in the data. +#' @param a2 The second age parameter in the Schnute function, the corresponding +#' length `"L2"` will be estimated. You may want to specify the largest age +#' class in the data. Set `"a2" = 999` to effectively estimate `linf`. #' @param downsample If not `Inf` this represents a number of fish specimens to #' sample prior to model fitting. Can be useful for large data sets that you #' want to fit with MCMC for testing purposes. @@ -34,8 +45,25 @@ #' @param ... Any other arguments to pass on to [rstan::sampling()] or #' [rstan::optimizing()]. #' @importFrom stats median quantile rlnorm runif median +#' @references +#' Schnute, J. 1981. A Versatile Growth Model with Statistically Stable Parameters. +#' Canadian Journal of Fisheries and Aquatic Sciences, 38(9), 1128–1140. +#' #' #' @export +#' @details +#' Linf and t0 are derived from the Schnute equations as: +#' +#' \deqn{ +#' linf = \left(\dfrac{\exp(ka_2) L_2^p - \exp(ka_1)L_1^p}{\exp(ka_2) - \exp(ka_1)}\right)^{1/p} +#' } +#' +#' \deqn{ +#' t_0 = a_1 + a_2 - \dfrac{1}{k}\log\left(\dfrac{\exp(ka_2) L_2^p - \exp(ka_1)L_1^p}{L_2^p - L_1^p}\right) +#' } +#' +#' `t_0` is undefined when `p < 0`. +#' #' @examples #' \dontrun{ #' # with `rstan::optimizing()` for the mode of the posterior density: @@ -45,6 +73,11 @@ #' model_f$model #' model_f$predictions #' +#' # Schnute functions +#' model_fs <- fit_schnute(pop_samples, sex = "female") +#' model_ms <- fit_schnute(pop_samples, sex = "male") +#' plot_schnute(model_fs, model_ms) +#' #' # You can also fit both sexes combined if you want. #' # Just note that you need to specify the colours explicitly in the plot. #' model_all <- fit_vb(pop_samples, sex = "all") @@ -65,7 +98,196 @@ #' obj <- fit_vb(pop_samples[1:2,]) #' plot_vb(obj, obj) #' } +fit_schnute <- function(dat, + sex = c("female", "male", "all"), + method = c("tmb", "mpd", "mcmc"), + growth = c("schnute", "vb"), + a1 = 0, + a2 = 999, + downsample = Inf, + chains = 4L, + iter = 1000L, + cores = parallel::detectCores(), + allow_slow_mcmc = FALSE, + est_method = median, + min_samples = 50L, + too_high_quantile = 1.0, + uniform_priors = FALSE, + ageing_method_codes = NULL, + usability_codes = c(0, 1, 2, 6), + check_convergence_tmb = TRUE, + tmb_inits = list(k = 0.5, L1 = 10, L2 = 40, log_sigma = log(0.1), p = 1), + ...) { + + method <- match.arg(method) + growth <- match.arg(growth) + if (growth == "vb" && method != "tmb") { + stop("Use fit_vb(method = \"", method, "\")") + } + + if ("species_common_name" %in% names(dat)) { + if (length(unique(dat$species_common_name)) > 1L) { + stop("Multiple species detected via the `species_common_name` column. ", + "fit_schnute() is for use with a single species. Filter the data yourself ", + "first.", + call. = FALSE + ) + } + } + + if (!is.null(usability_codes)) { + dat <- filter(dat, .data$usability_code %in% usability_codes) + } + + if (!is.null(ageing_method_codes)) { + dat <- filter(dat, .data$ageing_method %in% ageing_method_codes) + } + dat <- dat[!duplicated(dat$specimen_id), , drop = FALSE] + dat <- filter(dat, !is.na(.data$sex), !is.na(.data$length), !is.na(.data$age)) + ql <- quantile(dat$length, probs = too_high_quantile) + dat <- filter(dat, length <= ql) + + sex <- match.arg(sex) + dat <- switch(sex, + "female" = filter(dat, sex == 2L), + "male" = filter(dat, sex == 1L), + "all" = dat + ) + + if (nrow(dat) < min_samples) { + return(list( + predictions = tibble( + ages = NA, length = NA + ), + pars = list(k = NA, L1 = NA, L2 = NA, log_sigma = NA, p = NA, linf = NA, t0 = NA), + data = dat, + model = NA + )) + } + + if (nrow(dat) > downsample) { + dat <- dat[sample(seq_len(nrow(dat)), downsample), , drop = FALSE] + } + + if (method %in% c('mpd', 'mcmc')) { + rstan::rstan_options(auto_write = TRUE) + if (uniform_priors) { + model_file <- system.file("stan", "schnute-nopriors.stan", package = "gfplot") + .f <- "schnute-nopriors.stan" + } else { + model_file <- system.file("stan", "schnute.stan", package = "gfplot") + .f <- "schnute.stan" + } + + if (!file.exists(.f)) { + file.copy(model_file, to = .f) + } + + schnute_mod_gfplot <- rstan::stan_model(.f) + assign("schnute_mod_gfplot", schnute_mod_gfplot, envir = globalenv()) + } + + if (method == "mpd") { + mpd_init <- list() + mpd_init$k <- 0.2 # wild guess + mpd_init$L1 <- 10 # wild guess + mpd_init$L2 <- 40 # wild guess + mpd_init$sigma <- 0.1 + mpd_init$p <- 1 + + m <- suppressMessages(rstan::optimizing(schnute_mod_gfplot, + data = list( + N = nrow(dat), age = dat$age, length = dat$length, + len_upper_sd = quantile(dat$length, 0.99)[[1]] * 2, + a1 = a1, + a2 = a2 + ), + init = mpd_init, hessian = TRUE, ... + )) + + if (m$return_code != 0L && check_convergence_tmb) { + stop("Schnute growth model did not converge!") + } + + pars <- as.list(m$par) + } + if (method == "mcmc") { + if (nrow(dat) > 50000 && !allow_slow_mcmc) { + stop("There are > 50,000 aged fish. ", + "MCMC sampling may take a long time. Set `allow_slow_mcmc = TRUE` ", + "if you still want to sample from the posterior.", + call. = FALSE + ) + } + m <- rstan::sampling(schnute_mod_gfplot, + data = list( + N = nrow(dat), age = dat$age, length = dat$length, + len_upper_sd = quantile(dat$length, 0.99)[[1]] * 2, + a1 = a1, + a2 = a2 + ), + chains = chains, iter = iter, cores = cores, ... + ) + + pars <- lapply(rstan::extract(m), est_method) + } + + if (method == "tmb") { + dlls <- getLoadedDLLs() + if (!any(vapply(dlls, function(x) x[["name"]] == "schnute", FUN.VALUE = TRUE))) { + lib.f <- system.file("tmb", "schnute.cpp", package = "gfplot") + .f <- "schnute_gfplot.cpp" + if (!file.exists(.f)) { + file.copy(lib.f, to = .f) + } + TMB::compile(.f) + dyn.load(TMB::dynlib("schnute_gfplot")) + } + data <- list(len = dat$length, age = dat$age, a1 = a1, a2 = a2) + map <- list() + if (growth == "vb") map$p <- factor(NA) + + obj <- TMB::MakeADFun(data, tmb_inits, DLL = "schnute_gfplot", map = map, silent = TRUE) + opt <- stats::nlminb(obj$par, obj$fn, obj$gr) + if (opt$convergence != 0L && check_convergence_tmb) + stop("Schnute growth model did not converge!") + pars <- as.list(opt$par) + pars$linf <- obj$report(obj$env$last.par.best)$linf + pars$t0 <- obj$report(obj$env$last.par.best)$t0 + if (growth == "vb") pars$p <- tmb_inits$p + m <- obj + } + + if (is.na(pars[["t0"]])) warning("Derived t0 is NA.") + + ages <- seq(min(dat$age), max(dat$age), length.out = 200L) + pred <- schnute(ages, par = pars, a1 = a1, a2 = a2) + + list( + predictions = tibble(age = ages, length = pred), + pars = pars, data = as_tibble(dat), model = m + ) +} + +schnute <- function(ages, par, a1, a2) { + L1 <- par[["L1"]] + L2 <- par[["L2"]] + k <- par[["k"]] + p <- par[["p"]] + if(is.null(p)) p <- 1 + + tmp1 <- L1^p + tmp2 <- (L2^p - tmp1)/(1 - exp(-k * (a2 - a1))) + + tmp <- tmp1 + tmp2 * (1 - exp(-k * (ages - a1))) + invp <- 1/p + + tmp^invp +} + +#' @rdname fit_schnute +#' @export fit_vb <- function(dat, sex = c("female", "male", "all"), method = c("tmb", "mpd", "mcmc"), @@ -83,6 +305,41 @@ fit_vb <- function(dat, check_convergence_tmb = TRUE, tmb_inits = list(k = 0.5, linf = 40, log_sigma = log(0.1), t0 = -1), ...) { + + method <- match.arg(method) + + if (method == "tmb") { + + m <- fit_schnute( + dat = dat, + sex = sex, + method = method, + growth = "vb", + a1 = 0, + a2 = 999, + downsample = downsample, + chains = chains, + iter = iter, + cores = cores, + allow_slow_mcmc = allow_slow_mcmc, + est_method = est_method, + min_samples = min_samples, + too_high_quantile = too_high_quantile, + uniform_priors = uniform_priors, + ageing_method_codes = ageing_method_codes, + usability_codes = usability_codes, + check_convergence_tmb = check_convergence_tmb, + tmb_inits = list( + k = tmb_inits$k, L1 = 0.1 * tmb_inits$linf, + L2 = tmb_inits$linf, log_sigma = tmb_inits$log_sigma, + p = 1 + ), + ... + ) + + return(m) + } + if ("species_common_name" %in% names(dat)) { if (length(unique(dat$species_common_name)) > 1L) { stop("Multiple species detected via the `species_common_name` column. ", @@ -192,28 +449,6 @@ fit_vb <- function(dat, pars <- lapply(rstan::extract(m), est_method) } - if (method == "tmb") { - dlls <- getLoadedDLLs() - if (!any(vapply(dlls, function(x) x[["name"]] == "vb", FUN.VALUE = TRUE))) { - lib.f <- system.file("tmb", "vb.cpp", package = "gfplot") - .f <- "vb_gfplot.cpp" - if (!file.exists(.f)) { - file.copy(lib.f, to = .f) - } - TMB::compile(.f) - dyn.load(TMB::dynlib("vb_gfplot")) - } - data <- list(len = dat$length, age = dat$age) - - obj <- TMB::MakeADFun(data, tmb_inits, DLL = "vb_gfplot", silent = TRUE) - opt <- stats::nlminb(obj$par, obj$fn, obj$gr) - if (opt$convergence != 0L && check_convergence_tmb) - stop("VB growth model did not converge!") - pars <- as.list(opt$par) - m <- obj - } - - vb <- function(ages, linf, k, t0) linf * (1 - exp(-k * (ages - t0))) ages <- seq(min(dat$age), max(dat$age), length.out = 200L) pred <- vb(ages, linf = pars$linf, k = pars$k, t0 = pars$t0) @@ -223,6 +458,10 @@ fit_vb <- function(dat, ) } +vb <- function(ages, linf, k, t0) linf * (1 - exp(-k * (ages - t0))) + + + #' Fit a length-weight model #' #' For use with data for a single species. @@ -386,9 +625,10 @@ fit_length_weight <- function(dat, #' plot_vb(object_female = model_f, object_male = model_m) #' } + plot_growth <- function(object_female, object_male, object_all, - type = c("vb", "length-weight"), + type = c("schnute", "length-weight"), downsample = 2000L, pt_alpha = 0.2, xlab = "Age (years)", @@ -401,8 +641,11 @@ plot_growth <- function(object_female, object_male, col = c("Female" = "black", "Male" = "grey40"), french = FALSE, jitter = FALSE) { - xvar <- if (type[[1]] == "vb") "age" else "length" - yvar <- if (type[[1]] == "vb") "length" else "weight" + + if (type[1] == "vb") type <- "schnute" + type <- match.arg(type) + xvar <- if (type == "schnute") "age" else "length" + yvar <- if (type == "schnute") "length" else "weight" if (missing(object_all)) { line_dat <- bind_rows( @@ -437,8 +680,8 @@ plot_growth <- function(object_female, object_male, } no_pts <- if (nrow(pt_dat) == 0L) TRUE else FALSE - xdat <- if (type[[1]] == "vb") pt_dat$age else pt_dat$length - ydat <- if (type[[1]] == "vb") pt_dat$length else pt_dat$weight + xdat <- if (type == "schnute") pt_dat$age else pt_dat$length + ydat <- if (type == "schnute") pt_dat$length else pt_dat$weight if (nrow(pt_dat) > downsample) { if (!is.null(seed)) set.seed(seed) @@ -482,7 +725,7 @@ plot_growth <- function(object_female, object_male, ), size = 1.0) } - ann_func <- if (type[[1]] == "vb") ann_vb else ann_lw + ann_func <- if (type == "schnute") ann_schnute else ann_lw if (missing(object_all)) { if (!is.na(object_female$pars[[1]])) { @@ -522,6 +765,13 @@ plot_vb <- function(..., type = "vb") { ggplot2::ggtitle("Growth") } +#' @export +#' @rdname plot_growth +plot_schnute <- function(..., type = "schnute") { + plot_growth(..., type = type) + + ggplot2::ggtitle("Growth") +} + #' @export #' @rdname plot_growth plot_length_weight <- function(..., type = "length-weight", xlab = "Length (cm)", @@ -537,16 +787,31 @@ plot_length_weight <- function(..., type = "length-weight", xlab = "Length (cm)" } # annotation helpers: -ann_vb <- function(gg, pars, title, col, x, y, gap, french = FALSE) { - gg + ggplot2::annotate("text", +ann_schnute <- function(gg, pars, title, col, x, y, gap, french = FALSE) { + gout <- gg + ggplot2::annotate("text", label = title, x = x, y = y, hjust = 0, col = col, size = 3 ) + ann("k", pars[["k"]], dec = 2, x, y - gap, french = french) + - ann("linf", pars[["linf"]], dec = 1, x, y - gap * 2, french = french) + - ann("t0", pars[["t0"]], dec = 2, x, y - gap * 3, french = french) + ann("linf", pars[["linf"]], dec = 1, x, y - gap * 2, french = french) + + gap_counter <- 3 + if (!is.null(pars[["t0"]]) && !is.na(pars[["t0"]])) { + gout <- gout + + ann("t0", pars[["t0"]], dec = 2, x, y - gap * gap_counter, french = french) + gap_counter <- gap_counter + 1 + } + + if (!is.null(pars[["p"]]) && !is.na(pars[["p"]]) && pars[["p"]] != 1) { + gout <- gout + + ann("p", pars[["p"]], dec = 2, x, y - gap * gap_counter, french = french) + } + + gout } +ann_vb <- ann_schnute + ann_lw <- function(gg, pars, title, col, x, y, gap, french = FALSE) { gg + ggplot2::annotate("text", label = title, diff --git a/inst/stan/schnute-nopriors.stan b/inst/stan/schnute-nopriors.stan new file mode 100644 index 0000000..4a4c475 --- /dev/null +++ b/inst/stan/schnute-nopriors.stan @@ -0,0 +1,75 @@ +functions { + +real get_linf(real k, data real a1, data real a2, real L1, real L2, real p) { + real linf_num; + real linf_den; + real linf; + + linf_num = exp(k * a2) * pow(L2, p); + linf_num -= exp(k * a1) * pow(L1, p); + + linf_den = exp(k * a2) - exp(k * a1); + linf = pow(linf_num/linf_den, 1/p); + return linf; +} + +real get_t0(real k, data real a1, data real a2, real L1, real L2, real p) { + real linf_num; + real t0_den; + real t0; + + linf_num = exp(k * a2) * pow(L2, p); + linf_num -= exp(k * a1) * pow(L1, p); + + t0_den = pow(L2, p) - pow(L1, p); + t0 = -log(linf_num/t0_den)/k; + t0 += a1 + a2; + return t0; +} + +vector get_eta(data int N, data vector age, real k, data real a1, data real a2, real L1, real L2, real p) { + + vector[N] eta; + vector[N] tmp; + + real tmp1 = pow(L1, p); + real tmp2 = pow(L2, p) - tmp1; + tmp2 /= 1 - exp(-k * (a2 - a1)); + + tmp = 1 - exp(-k * (age - a1)); + tmp *= tmp2; + tmp += tmp1; + + for (i in 1:N) eta[i] = pow(tmp[i], 1/p); + + return eta; +} + +} +data { + int N; + vector[N] length; + vector[N] age; + real len_upper_sd; + real a1; + real a2; +} +parameters { + real k; + real L1; + real L2; + real sigma; + real p; +} +model { + //k ~ normal(0, 2); + //L1 ~ normal(0, len_upper_sd); + //L2 ~ normal(0, len_upper_sd); + //sigma ~ student_t(3, 0, 2); + //p ~ normal(0, 5); + length ~ lognormal(log(get_eta(N, age, k, a1, a2, L1, L2, p)), sigma); +} +generated quantities { + real linf = get_linf(k, a1, a2, L1, L2, p); + real t0 = get_t0(k, a1, a2, L1, L2, p); +} diff --git a/inst/stan/schnute.stan b/inst/stan/schnute.stan new file mode 100644 index 0000000..9084799 --- /dev/null +++ b/inst/stan/schnute.stan @@ -0,0 +1,75 @@ +functions { + +real get_linf(real k, data real a1, data real a2, real L1, real L2, real p) { + real linf_num; + real linf_den; + real linf; + + linf_num = exp(k * a2) * pow(L2, p); + linf_num -= exp(k * a1) * pow(L1, p); + + linf_den = exp(k * a2) - exp(k * a1); + linf = pow(linf_num/linf_den, 1/p); + return linf; +} + +real get_t0(real k, data real a1, data real a2, real L1, real L2, real p) { + real linf_num; + real t0_den; + real t0; + + linf_num = exp(k * a2) * pow(L2, p); + linf_num -= exp(k * a1) * pow(L1, p); + + t0_den = pow(L2, p) - pow(L1, p); + t0 = -log(linf_num/t0_den)/k; + t0 += a1 + a2; + return t0; +} + +vector get_eta(data int N, data vector age, real k, data real a1, data real a2, real L1, real L2, real p) { + + vector[N] eta; + vector[N] tmp; + + real tmp1 = pow(L1, p); + real tmp2 = pow(L2, p) - tmp1; + tmp2 /= 1 - exp(-k * (a2 - a1)); + + tmp = 1 - exp(-k * (age - a1)); + tmp *= tmp2; + tmp += tmp1; + + for (i in 1:N) eta[i] = pow(tmp[i], 1/p); + + return eta; +} + +} +data { + int N; + vector[N] length; + vector[N] age; + real len_upper_sd; + real a1; + real a2; +} +parameters { + real k; + real L1; + real L2; + real sigma; + real p; +} +model { + k ~ normal(0, 2); + L1 ~ normal(0, len_upper_sd); + L2 ~ normal(0, len_upper_sd); + sigma ~ student_t(3, 0, 2); + p ~ normal(0, 5); + length ~ lognormal(log(get_eta(N, age, k, a1, a2, L1, L2, p)), sigma); +} +generated quantities { + real linf = get_linf(k, a1, a2, L1, L2, p); + real t0 = get_t0(k, a1, a2, L1, L2, p); +} diff --git a/inst/tmb/schnute.R b/inst/tmb/schnute.R new file mode 100644 index 0000000..4a88f07 --- /dev/null +++ b/inst/tmb/schnute.R @@ -0,0 +1,76 @@ +library(dplyr) +library(TMB) + +## Compile and load the model +compile("inst/tmb/schnute.cpp") +dyn.load(dynlib("inst/tmb/schnute")) + +## Data and parameters +dat <- pop_samples +dat <- readRDS("../gfsynopsis/report/data-cache/shortspine-thornyhead.rds") +dat <- dat$survey_samples +usability_codes = c(0, 1, 2, 6) +dat <- filter(dat, usability_code %in% usability_codes) +dat <- dat[!duplicated(dat$specimen_id), , drop = FALSE] +dat <- filter(dat, !is.na(sex), !is.na(length), !is.na(age)) +dat <- switch("female", + "female" = filter(dat, sex == 2L), + "male" = filter(dat, sex == 1L), + "all" = dat +) +data <- list(len=dat$length, age=dat$age, a1 = 0, a2 = 999) +parameters <- list(k = 0.1 * rlnorm(1, sdlog = 0.1), + L1 = 5 * rlnorm(1, sdlog = 0.1), + L2 = 40 * rlnorm(1, sdlog = 0.1), + log_sigma = rnorm(1, mean = 0.1, sd = 0.01), + p = 1 * rlnorm(1, sdlog = 0.1)) + +plot(dat$age, dat$length) + +# dat <- filter(dat, age < 90) +## Make a function object +obj <- MakeADFun(data, parameters, DLL = "schnute", silent = TRUE) + +## Call function minimizer +opt <- nlminb(obj$par, obj$fn, obj$gr) +opt$par +## Get parameter uncertainties and convergence diagnostics +sdr <- sdreport(obj) +sdr + + +library(gfplot) +model_f <- fit_schnute(pop_samples, sex = "male") +model_f$pars +model_f <- fit_schnute(pop_samples, sex = "male", uniform_priors = TRUE) +model_f$pars + +library(tmbstan) +m <- tmbstan::tmbstan(obj = obj, iter = 500, chains = 4, cores = 4, + control = list(adapt_delta = 0.99)) +m + +y <- sapply(1:100, function(x) { + model_f <- fit_schnute(pop_samples, sex = "male", uniform_priors = TRUE) + # model_f$pars$t0 + model_f$model$value +}) +table(round(y, 1)) + +y <- sapply(1:100, function(x) { + model_f <- fit_schnute(pop_samples, sex = "male") + # model_f$pars$t0 + model_f$model$value +}) +table(round(y, 2)) + +y <- sapply(1:100, function(x) { + opt <- nlminb(obj$par, obj$fn, obj$gr) + opt$objective +}) +table(round(y, 2)) + +y <- fit_schnute(pop_samples, sex = "male", method = "tmb") +y$pars +y <- fit_schnute(pop_samples, sex = "female", method = "tmb") +y$pars diff --git a/inst/tmb/schnute.cpp b/inst/tmb/schnute.cpp new file mode 100644 index 0000000..441c47d --- /dev/null +++ b/inst/tmb/schnute.cpp @@ -0,0 +1,58 @@ +#include + +template +Type dlnorm(Type x, Type meanlog, Type sdlog, int give_log=0){ + Type logres; + logres = dnorm( log(x), meanlog, sdlog, true) - log(x); + if(give_log)return logres; else return exp(logres); +} + +template +Type objective_function::operator() () +{ + DATA_VECTOR(len); + DATA_VECTOR(age); + DATA_SCALAR(a1); + DATA_SCALAR(a2); + + PARAMETER(k); + PARAMETER(L1); + PARAMETER(L2); + PARAMETER(log_sigma); + PARAMETER(p); // vB: p = 1; Richards: 0 < p < 1 with p = ratio of length of inflection to Linf + + int n = len.size(); + vector eta(n); + Type nll = 0.0; + + Type tmp1 = pow(L1, p); + Type tmp2 = pow(L2, p) - tmp1; + tmp2 /= 1 - exp(-k * (a2 - a1)); + Type sigma = exp(log_sigma); + + for(int i = 0; i < n; i++){ + Type tmp = 1 - exp(-k * (age(i) - a1)); + tmp *= tmp2; + tmp += tmp1; + eta(i) = pow(tmp, 1/p); + nll -= dlnorm(len(i), log(eta(i)) - 0.5 * sigma * sigma, sigma, true); + } + + Type linf_num = exp(k * a2) * pow(L2, p); + linf_num -= exp(k * a1) * pow(L1, p); + + Type linf_den = exp(k * a2) - exp(k * a1); + Type linf = pow(linf_num/linf_den, 1/p); + + Type t0_den = pow(L2, p) - pow(L1, p); + Type t0 = -log(linf_num/t0_den)/k; + t0 += a1 + a2; + + ADREPORT(linf); + ADREPORT(t0); + REPORT(linf); + REPORT(t0); + REPORT(eta); + + return nll; +} diff --git a/man/fit_vb.Rd b/man/fit_schnute.Rd similarity index 62% rename from man/fit_vb.Rd rename to man/fit_schnute.Rd index 3be3b57..339ad47 100644 --- a/man/fit_vb.Rd +++ b/man/fit_schnute.Rd @@ -1,9 +1,33 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/growth.R -\name{fit_vb} +\name{fit_schnute} +\alias{fit_schnute} \alias{fit_vb} -\title{Fit a von Bertalanffy growth model} +\title{Fit a growth model for length-at-age} \usage{ +fit_schnute( + dat, + sex = c("female", "male", "all"), + method = c("tmb", "mpd", "mcmc"), + growth = c("schnute", "vb"), + a1 = 0, + a2 = 999, + downsample = Inf, + chains = 4L, + iter = 1000L, + cores = parallel::detectCores(), + allow_slow_mcmc = FALSE, + est_method = median, + min_samples = 50L, + too_high_quantile = 1, + uniform_priors = FALSE, + ageing_method_codes = NULL, + usability_codes = c(0, 1, 2, 6), + check_convergence_tmb = TRUE, + tmb_inits = list(k = 0.5, L1 = 10, L2 = 40, log_sigma = log(0.1), p = 1), + ... +) + fit_vb( dat, sex = c("female", "male", "all"), @@ -34,6 +58,17 @@ fit_vb( \code{\link[rstan:stanmodel-method-optimizing]{rstan::optimizing()}}) or \code{"mcmc"} for full MCMC sampling with Stan (with \code{\link[rstan:stanmodel-method-sampling]{rstan::sampling()}}). \code{"tmb"} for a TMB model.} +\item{growth}{Whether to estimate the inflection parameter p (\code{"schnute"}) or +fix the parameter p = 1 (\code{"vb"}) to reduce to the von Bertalanffy equation.} + +\item{a1}{The first age parameter in the Schnute function, the corresponding +length \code{"L1"} will be estimated. You may want to specify the smallest age +class in the data.} + +\item{a2}{The second age parameter in the Schnute function, the corresponding +length \code{"L2"} will be estimated. You may want to specify the largest age +class in the data. Set \code{"a2" = 999} to effectively estimate \code{linf}.} + \item{downsample}{If not \code{Inf} this represents a number of fish specimens to sample prior to model fitting. Can be useful for large data sets that you want to fit with MCMC for testing purposes.} @@ -74,7 +109,23 @@ Set to \code{NULL} to include all samples.} \code{\link[rstan:stanmodel-method-optimizing]{rstan::optimizing()}}.} } \description{ -For use with data for a single species. +For use with data for a single species, fit either a Schnute growth curve or +a von Bertalanffy growth curve. \code{fit_vb(method = "tmb")} will use \code{fit_schnute()} +with \code{p = 1}. Otherwise, separate Stan files are used. The Richards function +is also nested in the Schnute function with \verb{0 < p < 1}. +} +\details{ +Linf and t0 are derived from the Schnute equations as: + +\deqn{ +linf = \left(\dfrac{\exp(ka_2) L_2^p - \exp(ka_1)L_1^p}{\exp(ka_2) - \exp(ka_1)}\right)^{1/p} +} + +\deqn{ +t_0 = a_1 + a_2 - \dfrac{1}{k}\log\left(\dfrac{\exp(ka_2) L_2^p - \exp(ka_1)L_1^p}{L_2^p - L_1^p}\right) +} + +\code{t_0} is undefined when \code{p < 0}. } \examples{ \dontrun{ @@ -85,6 +136,11 @@ plot_vb(model_f, model_m) model_f$model model_f$predictions +# Schnute functions +model_fs <- fit_schnute(pop_samples, sex = "female") +model_ms <- fit_schnute(pop_samples, sex = "male") +plot_schnute(model_fs, model_ms) + # You can also fit both sexes combined if you want. # Just note that you need to specify the colours explicitly in the plot. model_all <- fit_vb(pop_samples, sex = "all") @@ -106,3 +162,8 @@ obj <- fit_vb(pop_samples[1:2,]) plot_vb(obj, obj) } } +\references{ +Schnute, J. 1981. A Versatile Growth Model with Statistically Stable Parameters. +Canadian Journal of Fisheries and Aquatic Sciences, 38(9), 1128–1140. +\url{https://doi.org/10.1139/f81-153} +} diff --git a/man/plot_growth.Rd b/man/plot_growth.Rd index 560d62f..60bbed2 100644 --- a/man/plot_growth.Rd +++ b/man/plot_growth.Rd @@ -3,6 +3,7 @@ \name{plot_growth} \alias{plot_growth} \alias{plot_vb} +\alias{plot_schnute} \alias{plot_length_weight} \title{Plot von Bertalanffy or length-weight fits} \usage{ @@ -10,7 +11,7 @@ plot_growth( object_female, object_male, object_all, - type = c("vb", "length-weight"), + type = c("schnute", "length-weight"), downsample = 2000L, pt_alpha = 0.2, xlab = "Age (years)", @@ -21,11 +22,14 @@ plot_growth( lab_x_gap = 0.3, lab_y_gap = 0.06, col = c(Female = "black", Male = "grey40"), - french = FALSE + french = FALSE, + jitter = FALSE ) plot_vb(..., type = "vb") +plot_schnute(..., type = "schnute") + plot_length_weight( ..., type = "length-weight", @@ -68,6 +72,8 @@ plot_length_weight( \item{french}{Logical.} +\item{jitter}{Logical, whether to jitter data points.} + \item{...}{Arguments to pass to \code{\link[=plot_growth]{plot_growth()}}.} } \description{ From 83a6a4cc8227ee184efcdb1b4e242189a10c9ff7 Mon Sep 17 00:00:00 2001 From: "Quang C. Huynh" Date: Tue, 6 Jun 2023 12:24:34 -0700 Subject: [PATCH 3/4] Update fit_mat_ogive.Rd with link argument --- man/plot_mat_ogive.Rd | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/man/plot_mat_ogive.Rd b/man/plot_mat_ogive.Rd index cfad751..941d181 100644 --- a/man/plot_mat_ogive.Rd +++ b/man/plot_mat_ogive.Rd @@ -14,7 +14,8 @@ fit_mat_ogive( months = seq(1, 12), custom_maturity_at = NULL, ageing_method_codes = NULL, - usability_codes = c(0, 1, 2, 6) + usability_codes = c(0, 1, 2, 6), + link = c("logit", "probit", "cloglog", "cauchit", "log") ) plot_mat_ogive( @@ -70,6 +71,8 @@ See \code{\link[gfdata:lookup]{gfdata::get_age_methods()}}.} All usability codes not in this vector will be omitted. Set to \code{NULL} to include all samples.} +\item{link}{The link function for the binomial GLM (by default, \code{"logit"}).} + \item{object}{Output from \code{\link[=fit_mat_ogive]{fit_mat_ogive()}}.} \item{col}{A named character vector declaring the colors for F and M.} From ea7bb29d25ed2a9d86c3d84a8e577c3ea53440c6 Mon Sep 17 00:00:00 2001 From: Quang Huynh Date: Thu, 7 Sep 2023 14:16:33 -0700 Subject: [PATCH 4/4] Set maturity parameter bound to 1000 for large fish --- R/maturity.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/maturity.R b/R/maturity.R index 0e1ad54..30b85ac 100644 --- a/R/maturity.R +++ b/R/maturity.R @@ -633,9 +633,9 @@ extract_maturity_perc_re <- function(betas, random_intercepts, model) { out } -binomial_perc <- function(x, a, b, perc = 0.5, linkinv, ...) { +binomial_perc <- function(a, b, perc = 0.5, linkinv, ...) { f <- function(x) linkinv(a + b * x) - perc - uniroot(f, interval = c(-100, 100))$root + uniroot(f, interval = c(-100, 1000))$root } mat_par_delta_method <- function(model, perc = 0.5) {