In cttedwards/lhm: Life-history module

knitr::opts_chunk$set(fig.path = 'fig/bdm-examples-', fig.width = 6, tidy = TRUE, tidy.opts = list(blank = TRUE, width.cutoff = 95), message = FALSE, warning = FALSE, collapse = TRUE, comment = "#>")

Development of a prior for $r$

If life history data are available we can populate an object of the lhm class, which is obtained from the lhm (life-history module) package, and is used to store and manipulate life-history data. Monte-Carlo samples are generated, and application of the lhm::rCalc function to this class of object calcuates values of the intrinsic growth rate $r$ for each iteration, producing an object of the prior class. The prior class contains a numeric vector with an additional slot for holding parameters of the associated distribution. Currently only the log-normal distribution is supported, with parameters stored in [e-mail].

library(lhm)

# initialise lhm data object for calculation of r with uncertainty
rdat <- lhm(ainf = 100, iter = 200)

# then life-history vectors can be assigned to each iteration
# with or without uncertainty
nmort(rdat)    <- list(mu = 0.18)
maturity(rdat) <- c(0.0,0.01,0.02,0.06,0.14,0.28,0.50,0.72,0.86,0.94,0.98,0.99,1.00)
size(rdat)     <- list(mu = list(Linf = 106.5, k = 0.229, t0 = 0.01))
mass(rdat)     <- list(mu = list(a = 1.88e-9, b = 3.305))
sr(rdat)       <- list(type = 'BH', mu = 0.90, cv = 0.10)

# calculate r prior and fit log-normal distribution
r <- rCalc(rdat)
plot(r)

Refine prior

It is possible to refine the prior distribution using an allometric scaling constant

print_rmax <- function(data) {

    cat('rmax:', round(quantile(data$rmax, c(0.5, 0.025, 0.975)), 3), '\n')
    cat('rmax (rT_exact):', with(subset(data, rT_exact), round(quantile(rmax, c(0.5, 0.025, 0.975)), 3)), '\n')
    cat('rmax (rT_adjusted):', with(subset(data, rT_adjusted), round(quantile(rmax, c(0.5, 0.025, 0.975)), 3)), '\n')
}

refine_rmax <- function(data, a = 1, sigma2 = 0.045, delta = 0.1) {

    # check headers
    stopifnot(all(c("rmax", "tmax") %in% names(data)))

    # estimated allometric scaling constant
    data$rT <- with(data, tmax * rmax)

    # trim using 'exact' method
    data$rT_exact <- with(data, abs(rT - a) < delta)

    # trim using 'adjusted' method
    asim <- rnorm(nrow(data), a, sqrt(sigma2))
    data$rT_adjusted <- with(data, abs(rT - asim) < delta)

    # save scaling constants
    data$a    <- a
    data$asim <- asim

    # return
    return(data)
}

dfr <- data.frame(rmax = r@.Data, tmax = [e-mail])
dfr <- refine_rmax(dfr)

print_rmax(dfr)

Alternative estimators

# natural mortality
M <- apply(nmort(rdat), 2, unique)

# spawning biomass per recruit
SBPR  <- apply(rdat@lhdat[['survivorship']] * rdat@lhdat[['mass']] * rdat@lhdat[['maturity']], 2, sum)
#cat('SBPR:',SBPR,'\n')

# recruits per unit of spawning biomass
alpha <- (4 * rdat@lhdat[['h']])/(SBPR * (1 - rdat@lhdat[['h']]))
#cat('alpha:',alpha,'\n')

# recruits per recruit
alpha_hat <- alpha * SBPR

# recruits per recruit per year
alpha_tilde <- alpha_hat * (1 - exp(-M))

# age at knife edge maturity
amat <- (1:rdat@ainf)[apply(maturity(rdat), 2, function(x) max(which(x <= 0.5)))]

dfr <- data.frame(l = exp(-M), 
                  h = as.numeric(rdat@lhdat[['h']]),
                  SBPR = as.numeric(SBPR), 
                  alpha = as.numeric(alpha), 
                  alpha_hat = as.numeric(alpha_hat), 
                  alpha_tilde = as.numeric(alpha_tilde), 
                  amat = as.numeric(amat))

# objective function
obj <- function(r, par) {

  amat        <- par[1]
  l           <- par[2]
  alpha_tilde <- par[3]

  return(exp(r)^amat - l * exp(r)^(amat - 1) - alpha_tilde * l)
}

dfr$r <- apply(dfr, 1, function(x) uniroot(obj, interval = c(0,10), par = c(x['amat'], x['l'], x['alpha_tilde']))$root)

plot(r@.Data, dfr$r)
abline(0,1)

cttedwards/lhm documentation built on May 6, 2022, 6:30 p.m.