Title: | Fit Sex-Specific Life History Models with Missing Classifications |
---|---|
Description: | Fits sex-specific life-history models for fish and other taxa where some of the individuals have unknown sex. |
Authors: | Coilin Minto [aut, cre], John Hinde [aut], Rui Coelho [ctb, dtc] |
Maintainer: | Coilin Minto <[email protected]> |
License: | GPL (>= 3) |
Version: | 0.1.0.9000 |
Built: | 2025-03-08 02:59:22 UTC |
Source: | https://github.com/mintoc/lhmixr |
Growth data for deepwater smooth lantern shark Etmopterus pusillus from Coelho and Erzini (2007). Data are cross-sectional with one observation (row) per individual.
Epusillus
Epusillus
A data frame with five variables:
species
Full species name
sex
Sex of the animal: Female (F
) or Male (M
)
age
Age in years
length
Total length in centimetres
maturity
Maturation status: immature
or mature
Coelho, R. and Erzini, K. (2007). Population parameters of the smooth lantern shark, Etmopterus pusillus, in southern Portugal (NE Atlantic). Fisheries Research, 86, 42–57.
Growth data for deepwater velvet belly lantern shark Etmopterus spinax from Coelho and Erzini (2008). Data are cross-sectional with one observation (row) per individual.
Espinax
Espinax
A data frame with five variables:
species
Full species name
sex
Sex of the animal: Female (F
) or Male (M
)
age
Age in years
length
Total length in centimetres
maturity
Maturation status: immature
or mature
Coelho, R. and Erzini, K. (2008). Life history of a wide-ranging deepwater lantern shark in the north-east Atlantic, Etmopterus spinax (Chondrichthyes: Etmopteridae), with implications for conservation. Journal of Fish Biology, 73, 1419–1443.
get_growth_post_prob
returns the probability of
the observation(s) arising from the female component given a set of growth parameters
and an assumed distribution (normal or lognormal). The component probability is
given by Bayes' theorem. Used internally.
get_growth_post_prob(mixprop, muF, muM, sigmaF, sigmaM, data, distribution)
get_growth_post_prob(mixprop, muF, muM, sigmaF, sigmaM, data, distribution)
mixprop |
Numeric scalar of mixing proportion (overall sex ratio) |
muF |
Numeric vector with predicted female lengths |
muM |
Numeric vector with predicted male lengths |
sigmaF |
Numeric scalar for female residual standard deviation |
sigmaM |
Numeric scalar for male residual standard deviation |
data |
A data.frame with column "length". Note predicted means "muF" and "muM" must come from corresponding ages. |
distribution |
Character with options: "normal" or "lognormal". |
Numeric vector of the posterior probability of being female.
Minto, C., Hinde, J. and Coelho, R. (2017). Including unsexed individuals in sex-specific growth models. Canadian Journal of Fisheries and Aquatic Sciences. DOI: 10.1139/cjfas-2016-0450.
get_growth_post_prob(mixprop = 0.5, muF = 4, muM = 6, sigmaF = 1, sigmaM = 1, data = data.frame(length = 4.5), distribution = "normal")
get_growth_post_prob(mixprop = 0.5, muF = 4, muM = 6, sigmaF = 1, sigmaM = 1, data = data.frame(length = 4.5), distribution = "normal")
sim_vb_data
simulates sex-specific growth data according to
the von Bertalanffy growth model and a logistic model governing maturity.
sim_vb_data(nfemale, nmale, mean_ageF, mean_ageM, growth_parF, growth_parM, mat_parF, mat_parM, distribution)
sim_vb_data(nfemale, nmale, mean_ageF, mean_ageM, growth_parF, growth_parM, mat_parF, mat_parM, distribution)
nfemale |
Numeric scalar for number of female observations. |
nmale |
Numeric scalar for number of male observations. |
mean_ageF |
Numeric scalar for female mean age - used to generate ages from |
mean_ageM |
Numeric scalar for male mean age - used to generate ages from |
growth_parF |
Named ("linf", "k", "t0", "sigma") numeric vector with female growth parameters |
growth_parM |
Named ("linf", "k", "t0", "sigma") numeric vector with male growth parameters |
mat_parF |
Named ("A50", "MR") numeric vector with female maturation parameters A50 is the age at 50% maturity, MR is age range between 25% and 75% mature. |
mat_parM |
Named ("A50", "MR") numeric vector with male maturation parameters. |
distribution |
Character with options: "normal" or "lognormal" for simulated length-at-age distributon. |
data.frame with columns "age", "length", "true.sex", "obs.sex" (observed sex assuming immature animals are unclassified), "maturity" (binary: 1 if mature; 0 if immature).
sim.dat <- sim_vb_data(nfemale = 30, nmale = 30, mean_ageF = 3, mean_ageM = 3, growth_parF = c(linf = 30, k = 0.2, t0 = -1, sigma = 0.1), growth_parM = c(linf = 25, k = 0.2, t0 = -1, sigma = 0.1), mat_parF = c(A50 = 3, MR = 1), mat_parM = c(A50 = 2, MR = 1), distribution = "lognormal") plot(jitter(sim.dat$age), sim.dat$length, xlim=c(0, max(sim.dat$age)), ylim = c(0, max(sim.dat$length)), col = c("red", "blue", "grey")[match(sim.dat$obs.sex, c("female", "male", "unclassified"))], pch = 19, xlab = "age", ylab = "Length")
sim.dat <- sim_vb_data(nfemale = 30, nmale = 30, mean_ageF = 3, mean_ageM = 3, growth_parF = c(linf = 30, k = 0.2, t0 = -1, sigma = 0.1), growth_parM = c(linf = 25, k = 0.2, t0 = -1, sigma = 0.1), mat_parF = c(A50 = 3, MR = 1), mat_parM = c(A50 = 2, MR = 1), distribution = "lognormal") plot(jitter(sim.dat$age), sim.dat$length, xlim=c(0, max(sim.dat$age)), ylim = c(0, max(sim.dat$length)), col = c("red", "blue", "grey")[match(sim.dat$obs.sex, c("female", "male", "unclassified"))], pch = 19, xlab = "age", ylab = "Length")
vb_bind_gr
returns the parameter gradients of negative log-likelihood for the von Bertalanffy model. Equality constraints across sexes can be implemented for any combination of parameters using the binding
argument.
vb_bind_gr(theta, binding, data, distribution)
vb_bind_gr(theta, binding, data, distribution)
theta |
A parameter vector of the same length as the maximum of |
binding |
A (4x2) parameter index matrix with rows named (in order): "lnlinf", "lnk", "lnnt0", "lnsigma" and the left column for the female parameter index and right column for male parameter index. Used to impose arbitrary equality constraints across the sexes (see Examples). |
data |
A data.frame with columns: "age", "length" and "weights". "weights" are set to 1 or 0 for known females or males, respectively; proportions otherwise. |
distribution |
Character with options: "normal" or "lognormal". |
Vector of parameter gradients:
## Unconstrained model binding <- matrix(c(1:8), ncol = 2, byrow = TRUE) rownames(binding) <- c("lnlinf", "lnk", "lnnt0", "lnsigma") colnames(binding) <- c("female", "male") ## starting values start.par <- c(rep(log(25), 2), rep(log(0.2), 2), rep(log(3), 2), rep(log(1), 2)) vb_bind_gr(theta = start.par, binding = binding, data = data.frame(age = rep(1, 2), length = rep(10, 2), weights = c(1, 0)), distribution = "lognormal")
## Unconstrained model binding <- matrix(c(1:8), ncol = 2, byrow = TRUE) rownames(binding) <- c("lnlinf", "lnk", "lnnt0", "lnsigma") colnames(binding) <- c("female", "male") ## starting values start.par <- c(rep(log(25), 2), rep(log(0.2), 2), rep(log(3), 2), rep(log(1), 2)) vb_bind_gr(theta = start.par, binding = binding, data = data.frame(age = rep(1, 2), length = rep(10, 2), weights = c(1, 0)), distribution = "lognormal")
vb_bind_nll
returns the negative log-likelihood for the von Bertalanffy model. Equality constraints across sexes can be implemented for any combination of parameters using the binding
argument.
vb_bind_nll(theta, binding, data, distribution)
vb_bind_nll(theta, binding, data, distribution)
theta |
A parameter vector of the same length as the maximum of |
binding |
A (4x2) parameter index matrix with rows named (in order): "lnlinf", "lnk", "lnnt0", "lnsigma" and the left column for the female parameter index and right column for male parameter index. Used to impose arbitrary equality constraints across the sexes (see Examples). |
data |
data.frame with columns: "age", "length" and "weights". "weights" are set to 1 or 0 for known females or males, respectively; proportions otherwise. |
distribution |
Character with options: "normal" or "lognormal" |
Complete data negative log-likelihood:
## Unconstrained model binding <- matrix(c(1:8), ncol = 2, byrow = TRUE) rownames(binding) <- c("lnlinf", "lnk", "lnnt0", "lnsigma") colnames(binding) <- c("female", "male") ## starting values start.par <- c(rep(log(25), 2), rep(log(0.2), 2), rep(log(3), 2), rep(log(1), 2)) vb_bind_nll(theta = start.par, binding = binding, data = data.frame(age = rep(1, 2), length = rep(10, 2), weights = c(1, 0)), distribution = "normal")
## Unconstrained model binding <- matrix(c(1:8), ncol = 2, byrow = TRUE) rownames(binding) <- c("lnlinf", "lnk", "lnnt0", "lnsigma") colnames(binding) <- c("female", "male") ## starting values start.par <- c(rep(log(25), 2), rep(log(0.2), 2), rep(log(3), 2), rep(log(1), 2)) vb_bind_nll(theta = start.par, binding = binding, data = data.frame(age = rep(1, 2), length = rep(10, 2), weights = c(1, 0)), distribution = "normal")
vb_growth_mix
fits sex-specific growth models where some of the animals are of unknown sex. Optimisation is via the Expectation-Maximisation algorithm. Equality constraints across sexes can be implemented for any combination of parameters using the binding
argument.
vb_growth_mix(start.list, data, binding, maxiter.em = 1000, reltol = 1e-08, plot.fit = FALSE, verbose = TRUE, optim.method = "BFGS", estimate.mixprop = TRUE, distribution)
vb_growth_mix(start.list, data, binding, maxiter.em = 1000, reltol = 1e-08, plot.fit = FALSE, verbose = TRUE, optim.method = "BFGS", estimate.mixprop = TRUE, distribution)
start.list |
A list with a list called par containing starting values for: "mixprop", "growth.par" (see Examples). |
data |
A data.frame with columns: "age", "length" and "obs.sex". "obs.sex" must have values "female", "male", "unclassified". |
binding |
A (4x2) parameter index matrix with rows named (in order): "lnlinf", "lnk", "lnnt0", "lnsigma" and the left column for the female parameter index and right column for male parameter index. Used to impose arbitrary equality constraints across the sexes (see Examples). |
maxiter.em |
Integer for maximum number of EM iterations (1e3 default). |
reltol |
Relative tolerance for EM observed data log likelihood convergence (1e-8 default). |
plot.fit |
Logical, if TRUE fit plotted per iteration. Red and blue circles are used for known females and males, respectively. Unclassified animals are plotted as triangle with the colour indicating the expected probability of being female or male (FALSE default). |
verbose |
Logical, if TRUE iteration and observed data log-likelihood printed. |
optim.method |
Character, complete data optimisation method to use in |
estimate.mixprop |
Logical, if TRUE the mixing proportion is estimated, otherwise fixed at the starting value. |
distribution |
Character with options: "normal" or "lognormal". |
List containing the components:
logLik.vec |
Observed data log-likelihood at each iteration. |
logLik |
Observed data log-likelihood on the last EM iteration. |
complete_data |
Data frame of the data (re-ordered) with component probabilities (tau). |
coefficients |
Parameter estimates (on the real line) and associated standard errors on the real line. |
vcov |
Estimated variance covariance matrix of the parameters estimated on the real line. Can be used to obtain parameter standard errors on the natural scale. |
convergence |
Binary with a "0" denoting convergence of the EM algorithm. |
Minto, C., Hinde, J. and Coelho, R. (2017). Including unsexed individuals in sex-specific growth models. Canadian Journal of Fisheries and Aquatic Sciences. DOI: 10.1139/cjfas-2016-0450.
set.seed(1010) sim.dat <- sim_vb_data(nfemale = 50, nmale = 50, mean_ageF = 4, mean_ageM = 4, growth_parF = c(linf = 30, k = 0.5, t0 = -1, sigma = 0.1), growth_parM = c(linf = 25, k = 0.5, t0 = -1, sigma = 0.1), mat_parF = c(A50 = 5, MR = 2), mat_parM = c(A50 = 3, MR = 2), distribution = "lognormal") ## Model fit with contrained Brody's growth coefficient ## Set up the constraint binding <- matrix(c(1:2, rep(3, 2), 4:7), ncol = 2, byrow = TRUE) rownames(binding) <- c("lnlinf", "lnk", "lnnt0", "lnsigma") colnames(binding) <- c("female", "male") ## note: lnnt0 is the natural logarithm of the negative of t0 (t0 < 0) ## starting values start.par <- c(c(log(30), log(25)), rep(log(0.3), 1), rep(log(1), 2), rep(log(.1), 2)) start.list <- list(par = list(mixprop = 0.5, growth.par = start.par)) vb.bind.fit <- vb_growth_mix(data = sim.dat, start.list = start.list, binding = binding, distribution = "lognormal", reltol = 1e-6)
set.seed(1010) sim.dat <- sim_vb_data(nfemale = 50, nmale = 50, mean_ageF = 4, mean_ageM = 4, growth_parF = c(linf = 30, k = 0.5, t0 = -1, sigma = 0.1), growth_parM = c(linf = 25, k = 0.5, t0 = -1, sigma = 0.1), mat_parF = c(A50 = 5, MR = 2), mat_parM = c(A50 = 3, MR = 2), distribution = "lognormal") ## Model fit with contrained Brody's growth coefficient ## Set up the constraint binding <- matrix(c(1:2, rep(3, 2), 4:7), ncol = 2, byrow = TRUE) rownames(binding) <- c("lnlinf", "lnk", "lnnt0", "lnsigma") colnames(binding) <- c("female", "male") ## note: lnnt0 is the natural logarithm of the negative of t0 (t0 < 0) ## starting values start.par <- c(c(log(30), log(25)), rep(log(0.3), 1), rep(log(1), 2), rep(log(.1), 2)) start.list <- list(par = list(mixprop = 0.5, growth.par = start.par)) vb.bind.fit <- vb_growth_mix(data = sim.dat, start.list = start.list, binding = binding, distribution = "lognormal", reltol = 1e-6)
vb_lengths
returns the predicted length-at-age for given named set of parameters for the von Bertalanffy growth function:
vb_lengths(theta, age)
vb_lengths(theta, age)
theta |
A numeric vector with named values "linf", "k", "t0". |
age |
A numeric vector of ages. |
Predicted length-at-age.
vb_lengths(theta = c("linf" = 30,"k" = 0.2,"t0" = -1), age = 0:10)
vb_lengths(theta = c("linf" = 30,"k" = 0.2,"t0" = -1), age = 0:10)