This document is for the course Fisheries Technology and Evaluation of Resources (Univ. Algarve).
Load needed library (for plotting)
library(ggplot2)
dat <- data.frame(
length = c(19.5,20.5,21.5,22.5,23.5,24.5,
25.5,26.5,27.5,28.5,29.5,30.5),
codend = c(1,1,2,10,13,10,14,7,1,1,2,1),
cover = c(4,25,69,151,190,192,146,102,62,24,4,8)
)
# Total catch at length
dat$total <- dat$codend + dat$cover
# Observed selectivity r(l) = codend / (codend + cover)
dat$sel_obs <- dat$codend / dat$total
# look/check the dataset
dat
## length codend cover total sel_obs
## 1 19.5 1 4 5 0.20000000
## 2 20.5 1 25 26 0.03846154
## 3 21.5 2 69 71 0.02816901
## 4 22.5 10 151 161 0.06211180
## 5 23.5 13 190 203 0.06403941
## 6 24.5 10 192 202 0.04950495
## 7 25.5 14 146 160 0.08750000
## 8 26.5 7 102 109 0.06422018
## 9 27.5 1 62 63 0.01587302
## 10 28.5 1 24 25 0.04000000
## 11 29.5 2 4 6 0.33333333
## 12 30.5 1 8 9 0.11111111
ggplot(dat, aes(x = length, y = sel_obs)) +
geom_point(size = 2) +
labs(
x = "Length (cm)",
y = "Observed selectivity (proportion in codend)",
title = "Codend–cover experiment: observed selectivity"
) +
theme_minimal(base_size = 13)
# Input some plausible biological parameters
L50_concept <- 30 # length at 50% retention
SR_concept <- 6 # selection range = L75 - L25
# For a logistic curve, SR = 2 * log(3) / b -> b = 2 * log(3) / SR
b_concept <- 2 * log(3) / SR_concept
a_concept <- -b_concept * L50_concept
x_concept <- seq(L50_concept - 4*SR_concept,
L50_concept + 4*SR_concept,
length.out = 200)
p_concept <- 1 / (1 + exp(-(a_concept + b_concept * x_concept)))
concept_df <- data.frame(length = x_concept, p = p_concept)
# plot
ggplot(concept_df, aes(x = length, y = p)) +
geom_line(linewidth = 1.2) +
geom_hline(yintercept = 0.5, linetype = "dashed") +
geom_vline(xintercept = L50_concept, linetype = "dotted") +
coord_cartesian(ylim = c(0, 1)) +
labs(
x = "Length (cm)",
y = "Retention probability",
title = "Typical logistic selectivity curve (general example)"
) +
theme_minimal(base_size = 13)
# Model:
# codend_l ~ Binomial(total_l, p_l)
# logit(p_l) = a + b * length_l
fit_glm <- glm(
cbind(codend, cover) ~ length,
data = dat,
family = binomial(link = "logit")
)
#summary of the model
summary(fit_glm)
##
## Call:
## glm(formula = cbind(codend, cover) ~ length, family = binomial(link = "logit"),
## data = dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.79258 1.60242 -2.367 0.0179 *
## length 0.04297 0.06509 0.660 0.5091
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13.019 on 11 degrees of freedom
## Residual deviance: 12.587 on 10 degrees of freedom
## AIC: 51.656
##
## Number of Fisher Scoring iterations: 4
# Extract coefficients
a_glm <- coef(fit_glm)[1]; a_glm
## (Intercept)
## -3.792585
b_glm <- coef(fit_glm)[2]; b_glm
## length
## 0.04297322
# For logistic:
# logit(p) = log(p / (1-p)) = a + b * L
# For p = 0.5 -> logit(0.5) = 0 -> 0 = a + b * L50 -> L50 = -a/b
L50_glm <- -a_glm / b_glm
# For p = 0.25 and 0.75:
logit25 <- log(0.25 / 0.75)
logit75 <- log(0.75 / 0.25)
L25_glm <- (logit25 - a_glm) / b_glm
L75_glm <- (logit75 - a_glm) / b_glm
SR_glm <- L75_glm - L25_glm
# see values
cat("\nLogistic GLM (binomial, maximum likelihood):\n")
##
## Logistic GLM (binomial, maximum likelihood):
cat(" L25 =", round(L25_glm, 2), "cm\n")
## L25 = 62.69 cm
cat(" L50 =", round(L50_glm, 2), "cm\n")
## L50 = 88.25 cm
cat(" L75 =", round(L75_glm, 2), "cm\n")
## L75 = 113.82 cm
cat(" SR =", round(SR_glm, 2), "cm\n\n")
## SR = 51.13 cm
# We will use ggplot, that creates nice looking plots
# Prediction grid around L50 to see the most critical section of the S-shape
x_grid <- seq(min(dat$length), max(dat$length), length.out = 200)
p_grid <- 1 / (1 + exp(-(a_glm + b_glm * x_grid)))
# create a dataframe with the predicted data
pred_df <- data.frame(length = x_grid, sel = p_grid)
ggplot() +
# observed points
geom_point(data = dat,
aes(x = length, y = sel_obs),
size = 2) +
# fitted GLM logistic curve
geom_line(data = pred_df,
aes(x = length, y = sel),
linewidth = 1.1) +
# reference lines at 0.5 and L50
geom_hline(yintercept = 0.5, linetype = "dashed") +
geom_vline(xintercept = L50_glm, linetype = "dotted") +
annotate("text", x = L50_glm, y = 0.05,
label = paste0("L50 = ", round(L50_glm, 1), " cm"),
angle = 90, vjust = -0.5, hjust = 0) +
coord_cartesian(ylim = c(0, 1)) +
labs(
x = "Length (cm)",
y = "Retention probability",
title = "Codend selectivity – logistic GLM fit"
) +
theme_minimal(base_size = 13)
# new experimental data - smaller mesh size
dat2 <- data.frame(
length = c(19.5,20.5,21.5,22.5,23.5,24.5,
25.5,26.5,27.5,28.5,29.5,30.5),
# codend catches
codend = c( 4, 3, 9, 14, 22, 30,
44, 55, 63, 58, 41, 30),
# cover catches
cover = c(76, 72, 66, 64, 58, 50,
46, 35, 27, 22, 17, 10)
)
# see data
dat2
## length codend cover
## 1 19.5 4 76
## 2 20.5 3 72
## 3 21.5 9 66
## 4 22.5 14 64
## 5 23.5 22 58
## 6 24.5 30 50
## 7 25.5 44 46
## 8 26.5 55 35
## 9 27.5 63 27
## 10 28.5 58 22
## 11 29.5 41 17
## 12 30.5 30 10