Stochastic life table analysis applied to shark population dynamics
Matheus de Barros
2022-07-22
Introduction
These are sets of functions that, using basic life-history information, allows the user to stochastically execute demographic simulations for shark populations using single-sex age-structured life tables (Rockwood, 2015) and project the effect of increases in fishing pressure in crucial demographic parameters.
Rockwood, L. L. (2015). Introduction to population ecology. John Wiley & Sons.
The base function “life.table” uses information on growth and reproductive parameters to execute basic life table calculations such as survival and age-specific fecundity/reproductive output to estimate basic demographic parameters. It also includes options to incorporate age-independent and age-specific fishing mortality in the survival calculations to estimate its effects on population growth. The parameter “MHA” (which stands for “Maximum Harvestable Age”) allows the user to specify a Maximum Age Limit for harvest, so the age classes < MHA will have its survival calculated by M + F, while age classes > MHA by just M.
life.table <- function (K = NULL,
Linf = NULL,
A50 = NULL,
fecund = NULL,
ob.max.age = NULL,
L0 = NULL,
MHA = ob.max.age,
FM = 0,
plot = FALSE) {
options(warn = 0)
if(MHA == 0) {
stop("ALC must be either set to 2 or higher, remembering that if ALC = 2, only age classes < 2 are harvested. If not set (= NULL), the function applies fishing mortality to all age classes")
} else if (is.null(K)) {
stop("Growth rate (K) not given")
} else if (is.null(Linf)) {
stop("Asymptotic length (Linf) not given")
}
K.dist <- K
Linf.dist <- Linf
t0 <- exp(-0.3922 - 0.2752*log(Linf) - 1.038*log(K.dist)) # Linf must be in cm
tmax <- ceiling(ob.max.age)
Ages <- 0:tmax #vector for the age range based on specified maximum age
# Mortality
M.frisk1 <- exp(0.42*log(K.dist) - 0.83)
M.frisk2 <- 1/(0.44*A50 + 1.87)
M.jensen <- 1.65/(A50)
M.HH <- -log(0.0178)/tmax
M.Dureuil <- exp(1.583 - 1.087*log(tmax))
M.random <- cbind(M.frisk1, M.frisk2, M.HH, M.Dureuil, M.jensen)
Mortality <- M.random[,sample(x=ncol(M.random),
size = 1)]
Survival <- NULL
for (i in 3:tmax) {
if (Ages[i] <= MHA) {
Survival[1] <- 1
Survival[2] <- exp(-(Mortality))-(1-exp(-(FM)))
Survival[i] <- exp(-(Mortality+FM))*Survival[i-1]
Survival[tmax + 1] <- exp(-(Mortality+FM))*Survival[tmax]
} else {
Survival[i] <- exp(-(Mortality))*Survival[i-1]
Survival[tmax + 1] <- exp(-(Mortality))*Survival[tmax]
}
}
lx <- NULL # probability of surviving from a given age class to another
for (i in 0:tmax) {
lx[1] <- (Survival[2])
lx[i] <- (Survival[i] + Survival[i +1])/2
lx[tmax + 1] <- 0
}
mx <- NULL
for (i in 0:tmax) {
ifelse(i >= A50 + 1,
mx[i] <- fecund*Survival[i],
mx[i] <- 0)
mx[tmax + 1] <- fecund*Survival[tmax]
}
life.tab <- cbind(Ages, Survival, lx, mx)
life.tab <- as.data.frame(life.tab)
life.tab$lxmx <- life.tab$lx*life.tab$mx
if (plot == TRUE) {
par(mfrow = c(1,1))
plot(life.tab$Survival ~ life.tab$Ages, pch = 20, xlab = "Age class", ylab = "Survival")
lines(life.tab$Survival ~ life.tab$Ages)
}
GRR <- sum(life.tab$mx)
R0 <- sum(life.tab$lxmx)
G.time <- (sum(life.tab$Ages*life.tab$lxmx)/R0)
r.growth <- (log(R0)/G.time)
lambda <- exp(r.growth)
output <- list(life.tab, GRR, R0, G.time, r.growth, lambda, FM, Mortality, tmax)
names(output) <- c('Life.table',
'GRR',
'R0',
'G',
'r',
'lambda',
'FM',
'M',
"max.age")
return(output)
}“life.table” outputs:
# FM = 0, age-independent fishing mortality
life.table(K = 0.08,
Linf = 260,
A50 = 11,
fecund = 26,
FM = 0,
ob.max.age = 20,
L0 = 60,
plot = TRUE,
MHA = )
## $Life.table
## Ages Survival lx mx lxmx
## 1 0 1.00000000 0.86154215 0.000000 0.00000000
## 2 1 0.86154215 0.80189851 0.000000 0.00000000
## 3 2 0.74225488 0.69086937 0.000000 0.00000000
## 4 3 0.63948386 0.59521308 0.000000 0.00000000
## 5 4 0.55094230 0.51280116 0.000000 0.00000000
## 6 5 0.47466002 0.44179982 0.000000 0.00000000
## 7 6 0.40893961 0.38062916 0.000000 0.00000000
## 8 7 0.35231871 0.32792807 0.000000 0.00000000
## 9 8 0.30353742 0.28252385 0.000000 0.00000000
## 10 9 0.26151028 0.24340621 0.000000 0.00000000
## 11 10 0.22530213 0.20970471 0.000000 0.00000000
## 12 11 0.19410728 0.18066945 5.046789 0.91180064
## 13 12 0.16723161 0.15565434 4.348022 0.67678847
## 14 13 0.14407708 0.13410278 3.746004 0.50234955
## 15 14 0.12412848 0.11553520 3.227340 0.37287140
## 16 15 0.10694191 0.09953844 2.780490 0.27676562
## 17 16 0.09213497 0.08575656 2.395509 0.20543063
## 18 17 0.07937816 0.07388289 2.063832 0.15248189
## 19 18 0.06838763 0.06365323 1.778078 0.11318042
## 20 19 0.05891882 0.05483994 1.531889 0.08400872
## 21 20 0.05076105 0.00000000 1.531889 0.00000000
##
## $GRR
## [1] 28.44984
##
## $R0
## [1] 3.295677
##
## $G
## [1] 13.21913
##
## $r
## [1] 0.09021861
##
## $lambda
## [1] 1.094414
##
## $FM
## [1] 0
##
## $M
## M.frisk2
## 0.1490313
##
## $max.age
## [1] 20# FM = 0.05, age-independet fishing mortality
life.table(K = 0.08,
Linf = 260,
A50 = 11,
fecund = 26,
FM = 0.05,
ob.max.age = 20,
L0 = 60,
plot = TRUE)
## $Life.table
## Ages Survival lx mx lxmx
## 1 0 1.000000000 0.780162851 0.0000000 0.000000000
## 2 1 0.780162851 0.697662918 0.0000000 0.000000000
## 3 2 0.615162985 0.550111303 0.0000000 0.000000000
## 4 3 0.485059622 0.433765990 0.0000000 0.000000000
## 5 4 0.382472358 0.342027028 0.0000000 0.000000000
## 6 5 0.301581698 0.269690318 0.0000000 0.000000000
## 7 6 0.237798938 0.212652398 0.0000000 0.000000000
## 8 7 0.187505857 0.167677663 0.0000000 0.000000000
## 9 8 0.147849469 0.132214821 0.0000000 0.000000000
## 10 9 0.116580173 0.104252162 0.0000000 0.000000000
## 11 10 0.091924151 0.082203441 0.0000000 0.000000000
## 12 11 0.072482732 0.064817896 1.8845510 0.122152632
## 13 12 0.057153059 0.051109291 1.4859795 0.075947361
## 14 13 0.045065523 0.040299976 1.1717036 0.047219626
## 15 14 0.035534429 0.031776767 0.9238952 0.029358401
## 16 15 0.028019106 0.025056167 0.7284967 0.018253336
## 17 16 0.022093229 0.019756935 0.5744240 0.011348857
## 18 17 0.017420641 0.015578459 0.4529367 0.007056055
## 19 18 0.013736277 0.012283706 0.3571432 0.004387042
## 20 19 0.010831135 0.009685775 0.2816095 0.002727606
## 21 20 0.008540414 0.000000000 0.2816095 0.000000000
##
## $GRR
## [1] 8.142349
##
## $R0
## [1] 0.3184509
##
## $G
## [1] 12.51699
##
## $r
## [1] -0.09141872
##
## $lambda
## [1] 0.9126355
##
## $FM
## [1] 0.05
##
## $M
## M.Dureuil
## 0.1876154
##
## $max.age
## [1] 20# FM = 0.05, but only applied to age classes < 3
life.table(K = 0.08,
Linf = 260,
A50 = 11,
fecund = 26,
FM = 0.05,
ob.max.age = 20,
L0 = 60,
plot = TRUE,
MHA = 3)
## $Life.table
## Ages Survival lx mx lxmx
## 1 0 1.00000000 0.81193740 0.000000 0.00000000
## 2 1 0.81193740 0.73834776 0.000000 0.00000000
## 3 2 0.66475812 0.60450802 0.000000 0.00000000
## 4 3 0.54425792 0.50635252 0.000000 0.00000000
## 5 4 0.46844713 0.43582166 0.000000 0.00000000
## 6 5 0.40319618 0.37511517 0.000000 0.00000000
## 7 6 0.34703417 0.32286462 0.000000 0.00000000
## 8 7 0.29869508 0.27789216 0.000000 0.00000000
## 9 8 0.25708924 0.23918400 0.000000 0.00000000
## 10 9 0.22127876 0.20586757 0.000000 0.00000000
## 11 10 0.19045639 0.17719186 0.000000 0.00000000
## 12 11 0.16392733 0.15251045 4.262111 0.65001641
## 13 12 0.14109356 0.13126696 3.668433 0.48154400
## 14 13 0.12144036 0.11298252 3.157449 0.35673657
## 15 14 0.10452468 0.09724496 2.717642 0.26427695
## 16 15 0.08996523 0.08369951 2.339096 0.19578118
## 17 16 0.07743379 0.07204084 2.013279 0.14503827
## 18 17 0.06664788 0.06200612 1.732845 0.10744699
## 19 18 0.05736436 0.05336916 1.491473 0.07959869
## 20 19 0.04937396 0.04593526 1.283723 0.05896816
## 21 20 0.04249656 0.00000000 1.283723 0.00000000
##
## $GRR
## [1] 23.94977
##
## $R0
## [1] 2.339407
##
## $G
## [1] 13.20987
##
## $r
## [1] 0.06433808
##
## $lambda
## [1] 1.066453
##
## $FM
## [1] 0.05
##
## $M
## M.jensen
## 0.15
##
## $max.age
## [1] 20Stochastically estimating demographic parameters
The function “stochastic.LT” is used to stochastically estimate demographic parameters using means and standard errors in life-history traits:
stochastic.LT <- function(K = NULL,
K.se = NULL,
Linf = NULL,
Linf.se = NULL,
L0 = NULL,
L0.se = NULL,
Fec = NULL,
Fec.se = NULL,
A50 = NULL,
A50.se = NULL,
FM = NULL,
FM.se = NULL,
MHA = ob.max.age,
ob.max.age = NULL,
N = 1000, plot = FALSE) {
prior <- data.frame(rnorm(N, K, K.se),
rnorm(N, Linf, Linf.se),
rnorm(N, L0, L0.se),
rnorm(N, Fec, Fec.se),
rnorm(N, A50, A50.se),
rnorm(N, FM, FM.se),
rnorm(N, MHA, 0),
rnorm(N, ob.max.age, 0))
names(prior) <- c('K', 'Linf', 'L0', 'Fec', 'A50', 'FM', 'MHA', 'Agemax')
# create empty vectors to store the parameter distributions
GRR.out <- as.vector(rep(NA, length(prior$K)))
R0.out <- as.vector(rep(NA, length(prior$K)))
G.out <- as.vector(rep(NA, length(prior$K)))
r.out <- as.vector(rep(NA, length(prior$K)))
lambda.out <- as.vector(rep(NA, length(prior$K)))
FM.out <- as.vector(rep(NA, length(prior$K)))
M.out <- as.vector(rep(NA, length(prior$K)))
Survival.df <- NULL # empty object to store survival curves for each run
N <- length(prior$K) # number of runs is the number of parameter replicates in the priors
agemax <- ifelse(ob.max.age > (1/K)*log10((Linf - L0)/((1 - 0.99)*Linf)),
ob.max.age,
(1/K)*log10((Linf - L0)/((1 - 0.99)*Linf)))
agemax <- ceiling(agemax)
for (i in 1:N) {
LT.out <- life.table(K = prior$K[i],
Linf = prior$Linf[i],
L0 = prior$L0[i],
A50 = prior$A50[i],
fecund = prior$Fec[i],
FM = prior$FM[i],
MHA = prior$MHA[i],
ob.max.age = agemax,
plot = FALSE)
GRR.out[i] <- LT.out$GRR
R0.out[i] <- LT.out$R0
G.out[i] <- LT.out$G
r.out[i] <- LT.out$r
lambda.out[i] <- LT.out$lambda
FM.out[i] <- LT.out$FM
M.out[i] <- LT.out$M
Survival.df <- cbind(Survival.df,
LT.out$Life.table$Survival,
deparse.level = 0)
}
par.dist <- cbind(GRR.out,
R0.out,
G.out,
r.out,
lambda.out,
FM.out,
M.out)
par.dist <- as.data.frame(par.dist)
Mortalities <- as.data.frame(M.out)
Survival.mean <- apply(Survival.df, 1, mean)
Survival.up <- apply(Survival.df, 1, quantile, probs = c(.95))
Survival.low <- apply(Survival.df, 1, quantile, probs = c(.05))
if (plot == TRUE) {
par(mfrow = c(2,2))
hist(par.dist$lambda.out,
main = "Finite population growth rate",
xlab = "λ", breaks = 60, col = 0, cex.main = 0.8)
mtext(paste("Median λ =", round(median(lambda.out),3)),
side = 3, cex = 0.7)
hist(par.dist$G.out,
main = "Generation time",
xlab = "G", breaks = 60, col = 0, cex.main = 0.8)
mtext(paste("Median G =", round(median(G.out),3)),
side = 3, cex = 0.7)
hist(par.dist$R0.out,
main = "Net reproductive rate", breaks = 80,
col = 0, cex.main = 0.8, xlab = "R0")
mtext(paste("Median R0 =", round(median(R0.out),3)),
side = 3, cex = 0.7)
plot(Survival.mean ~ LT.out$Life.table$Ages,
xlab = "Age class", ylab = "Survival",
main = "Survival curve", cex.main = 0.8, type = "b")
points(Survival.mean ~ LT.out$Life.table$Ages, pch = 20)
lines(Survival.up ~ LT.out$Life.table$Ages)
lines(Survival.low ~ LT.out$Life.table$Ages)
}
lambda.median <- median(par.dist$lambda.out, na.rm = TRUE)
G.median <- median(par.dist$G.out, na.rm = TRUE)
R0.median <- median(par.dist$R0.out, na.rm = TRUE)
output <- list(lambda.median, G.median, R0.median, par.dist, Survival.df, agemax)
names(output) <- c("Median_Lambda",
"Median_G",
"Median_R0",
"Par_dists",
"Survival_table",
"agemax")
return(output)
}“stochastic.LT” output:
stoch.LT.out <- stochastic.LT(K = 0.08,
K.se = 0.007,
Linf = 260,
Linf.se = 20,
A50 = 11,
A50.se = 1,
Fec = 26,
Fec.se = 5,
FM = 0,
FM.se = 0,
ob.max.age = 20,
L0 = 60,
L0.se = 0,
plot = TRUE,
N = 10000)
# age-independent F
stoch.LT.out <- stochastic.LT(K = 0.08,
K.se = 0.007,
Linf = 260,
Linf.se = 20,
A50 = 11,
A50.se = 1,
Fec = 26,
Fec.se = 5,
FM = 0.05,
FM.se = 0,
ob.max.age = 20,
L0 = 60,
L0.se = 0,
plot = TRUE,
N = 10000)
stoch.LT.out <- stochastic.LT(K = 0.08,
K.se = 0.007,
Linf = 260,
Linf.se = 20,
A50 = 11,
A50.se = 1,
Fec = 26,
Fec.se = 5,
FM = 0.05,
FM.se = 0,
MHA = 3,
ob.max.age = 20,
L0 = 60,
L0.se = 0,
plot = TRUE,
N = 10000)