9 October 2026
the model
For a sample s (a unique combination of site visit, species and lifestage) fished with n_s passes, let y_{si} be the number of fish caught on pass i and p_{si} the capture probability on that pass. The probability that a fish is first caught on pass i is
\pi_{si} = p_{si} \prod_{j < i} (1 - p_{sj})
Conditional on the total catch, T_s = \sum_i y_{si}, the counts are multinomial with probabilities \pi_{si} / P_s, where P_s = \sum_i \pi_{si} is the cumulative capture probability. This gives the log likelihood
\ell = \sum_s \left( \sum_i y_{si} \log \pi_{si} - T_s \log P_s \right)
which does not depend on the unknown number of fish at the site. Capture probability is modelled on the logit scale
\text{logit}(p_{si}) = \mathbf{x}_{si}^\top \boldsymbol{\beta} + \mathbf{z}_{si}^\top \mathbf{b}, \qquad b_k \sim N(0, \sigma_{g(k)}^2)
where \boldsymbol{\beta} are fixed effects and \mathbf{b} are optional random effects, for example for visit and sample, with one standard deviation \sigma_g for each random effect term. The random effects are integrated out of the likelihood using the Laplace approximation, and the fixed effects and standard deviations are estimated by maximum likelihood using RTMB.
Fish density is then estimated by dividing the total catch by the area fished and the cumulative capture probability: T_s / (A_s P_s).
installation
To install the ef package from r-universe run
install.packages("ef", repos = c('https://faskally.r-universe.dev', 'https://cloud.r-project.org'))or to install the latest development version from GitHub run
remotes::install_github("faskally/ef")Create Required Fields for Analysis
Create a factor for pass with separate levels for the first and subsequent passes. NOTE: Assuming the second and third pass have the same capture probability avoids problems with model identifiability (Millar et al., 2016)
Create a sampleID for grouping EF counts together. A sample is a unique combination of site visit x species x lifestage and is required by the EF function
ef_data$sampleID <-
as.numeric(
factor(
paste(ef_data$date, ef_data$siteID, ef_data$lifestage, ef_data$species)
)
)Create a visitID for unique site visits. This is used to include between visit variation as a random effect (see below).
Create a factor for year
ef_data$fyear <- factor(ef_data$year)Fit a Model
Fit a simple model. This contains a factor for pass (2 levels), lifestage (2 levels), site (3 levels) and year (5 levels). Continuous variables can be added as linear or smooth terms, see below and Millar et al., 2016 for further details.
m1 <-
efp(
count ~ siteID + fyear + pass12 + lifestage,
data = ef_data, pass = pass, id = sampleID
)##
## Call:
## efp(formula = count ~ siteID + fyear + pass12 + lifestage, data = ef_data,
## pass = pass, id = sampleID)
##
## Formula:
## count ~ siteID + fyear + pass12 + lifestage
##
## No random effects
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z^2|)
## (Intercept) 0.652301 0.1723 3.78656 1.527e-04
## siteIDSite2 0.074540 0.1354 0.55059 5.819e-01
## siteIDSite3 -0.023719 0.1473 -0.16098 8.721e-01
## fyear2012 -0.209949 0.1736 -1.20936 2.265e-01
## fyear2013 0.046976 0.1855 0.25317 8.001e-01
## fyear2014 0.358431 0.1834 1.95389 5.071e-02
## fyear2015 -0.004013 0.2114 -0.01899 9.849e-01
## pass122 -0.239160 0.1311 -1.82436 6.810e-02
## lifestageParr 0.590641 0.1210 4.88208 1.050e-06
##
## Degrees of Freedom: 90 Total (i.e. Null); 81 Residual
## AIC: 3188Smooth Terms
Smooth terms are penalised, and the amount of smoothing is estimated. A smooth is fitted as a random effect: the part of the smooth that is not penalised is a fixed effect, and the standard deviation of the random part controls how smooth the term is.
The recommended smooths are first and second order random walks. These use the Markov random field (mrf) basis from mgcv with a penalty created by rw_penalty(). Both the basis and the penalty are sparse, so these models are quick to fit even when the walk has many steps.
- a first order random walk (
order = 1) penalises the differences between neighbouring levels and shrinks towards a constant - a second order random walk (
order = 2) penalises the second differences and shrinks towards a linear trend
The covariate must be a factor whose levels are the steps of the walk, in order. Here the year factor is replaced by a first order random walk over years.
rw1 <- rw_penalty(levels(ef_data$fyear), order = 1)
m1_smooth <-
efp(
count ~ siteID + s(fyear, bs = "mrf", xt = list(penalty = rw1)) +
pass12 + lifestage,
data = ef_data, pass = pass, id = sampleID
)##
## Call:
## efp(formula = count ~ siteID + s(fyear, bs = "mrf", xt = list(penalty = rw1)) +
## pass12 + lifestage, data = ef_data, pass = pass, id = sampleID)
##
## Formula:
## count ~ siteID + s(fyear, bs = "mrf", xt = list(penalty = rw1)) +
## pass12 + lifestage
##
## Random effects:
## Estimate Std. Error
## sigma s(fyear) 0.192 0.1239
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z^2|)
## (Intercept) 0.59939 0.1565 3.8304 1.279e-04
## siteIDSite2 0.08302 0.1345 0.6171 5.371e-01
## siteIDSite3 -0.02424 0.1466 -0.1653 8.687e-01
## pass122 -0.24364 0.1309 -1.8614 6.269e-02
## lifestageParr 0.59448 0.1160 5.1239 2.993e-07
##
## Degrees of Freedom: 90 Total (i.e. Null); 84 Residual
## AIC: 3193The standard deviation, sigma s(fyear), is that of the steps of the random walk on the logit scale. The starting value of the walk is a fixed effect, so the intercept is for the first year.
The smooth effect of year on capture probability, here for the first pass for fry at Site1, with approximate 95% confidence intervals
newdata <-
data.frame(
siteID = "Site1", pass12 = "1", lifestage = "Fry",
fyear = factor(levels(ef_data$fyear), levels = levels(ef_data$fyear))
)
pred <- predict(m1_smooth, newdata = newdata, se.fit = TRUE)
year <- as.numeric(levels(ef_data$fyear))
plot(year, invlogit(pred$fit), type = "b", pch = 16, ylim = c(0, 1),
xlab = "Year", ylab = "Capture probability")
lines(year, invlogit(pred$fit - 2 * pred$se.fit), lty = 2)
lines(year, invlogit(pred$fit + 2 * pred$se.fit), lty = 2)
The models with a year factor and a random walk over years can be compared using AIC
AIC(m1, m1_smooth)A random walk can have steps with no data, for example a walk over every day of the year when sites were only fished on some of them. All the steps are included in the penalty, and the factor is given the same levels. Here a second order random walk is used for day of the year (doy).
days <- min(ef_data$doy):max(ef_data$doy)
ef_data$fdoy <- factor(ef_data$doy, levels = days)
rw2 <- rw_penalty(days, order = 2)
m1_doy <-
efp(
count ~ siteID + s(fdoy, bs = "mrf", xt = list(penalty = rw2)) +
pass12 + lifestage,
data = ef_data, pass = pass, id = sampleID
)##
## Call:
## efp(formula = count ~ siteID + s(fdoy, bs = "mrf", xt = list(penalty = rw2)) +
## pass12 + lifestage, data = ef_data, pass = pass, id = sampleID)
##
## Formula:
## count ~ siteID + s(fdoy, bs = "mrf", xt = list(penalty = rw2)) +
## pass12 + lifestage
##
## Random effects:
## Estimate Std. Error
## sigma s(fdoy) 3.937e-05 0.001866
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z^2|)
## (Intercept) 0.678300 0.1178 5.75842 8.491e-09
## siteIDSite2 0.001903 0.1448 0.01314 9.895e-01
## siteIDSite3 -0.158737 0.1653 -0.96058 3.368e-01
## pass122 -0.254290 0.1308 -1.94349 5.196e-02
## lifestageParr 0.607910 0.1135 5.35677 8.472e-08
## s(fdoy).null1 -0.550210 0.3837 -1.43391 1.516e-01
##
## Degrees of Freedom: 90 Total (i.e. Null); 83 Residual
## AIC: 3196In this case the standard deviation is very close to zero, so the smooth has been shrunk to its unpenalised part, a linear trend in day of the year (the fixed effect s(fdoy).null1).
newdata <-
data.frame(
siteID = "Site1", pass12 = "1", lifestage = "Fry",
fdoy = factor(days, levels = days)
)
pred <- predict(m1_doy, newdata = newdata, se.fit = TRUE)
plot(days, invlogit(pred$fit), type = "l", ylim = c(0, 1),
xlab = "Day of year", ylab = "Capture probability")
lines(days, invlogit(pred$fit - 2 * pred$se.fit), lty = 2)
lines(days, invlogit(pred$fit + 2 * pred$se.fit), lty = 2)
rug(ef_data$doy)
Any other penalised smooth from mgcv can be used, for example s(doy) or te(doy, area), but their penalties are not sparse so they take longer to fit. An unpenalised regression spline, whose flexibility is set by its basis dimension k, can be fitted using fx = TRUE, for example s(doy, k = 3, fx = TRUE).
Get Predicted Values for P
Predict P from the model to the full dataset. Note, if all the data is multi-pass then can just use fitted values. However, if you want P for both single and multi-pass data then predict to the new dataframe
ef_data$prob <- predict(m1, newdata = ef_data, type = "response")Standard errors of the predictions are available using se.fit = TRUE. To use the plotting and prediction tools in mgcv the model can also be converted to a gam object using as.gam(m1).
Obtain Cumulative Estimates of P (for estimation of densities later)
Aggregate counts by sampleID (i.e. counts of fry or parr per site visit)
counts <- aggregate(count ~ sampleID, ef_data, sum)Reshape to get capture probability by species x lifestage x site visit x pass. This facilitates the summary of P.
probs <- cast(ef_data, formula = sampleID ~ pass, value = "prob")Get cumulative capture probability by sampleID (see Glover et al., 2018)
Remove individual pass capture probabilities
probs <- select(probs, sampleID, prob)Remove pass and pass-specific data from dataframe
densdata <- unique(select(ef_data, -count, -prob, -pass, -pass12))Join data and cumulative counts by sampleID
densdata <- left_join(densdata, counts, by = "sampleID")Join data and cumulative capture probability by sampleID
densdata <- left_join(densdata, probs, by = "sampleID")Estimate Density
Density (N/m-2) is estimated as \frac{counts}{area*P}
densdata $ Density_Estimate <- densdata $ count / (densdata $ area *
densdata $ prob)If you want to model counts directly (Poisson model), you need to have cumulative P and area in the offset. For a Poisson model the offset is \log({area*cumulativeP}).
densdata $ offset <- log(densdata $ area * densdata $ prob)Plot to check
Comparing Models
Electrofishing count data are typically overdispersed: capture probability varies between visits to a site, and between samples within a visit, by more than the covariates can explain. This variation is included in the model using random effects, which are specified by adding s(f, bs = "re") terms to the formula, where f is a factor.
The sampleID needs to be a factor to be used as a random effect
ef_data$fsampleID <- factor(ef_data$sampleID)Fit a model with visit and sample random effects
mfull <- efp(count ~ siteID + year + lifestage + pass12 +
s(visitID, bs = "re") + s(fsampleID, bs = "re"),
data = ef_data, pass = pass, id = sampleID)
mfull##
## Call:
## efp(formula = count ~ siteID + year + lifestage + pass12 + s(visitID,
## bs = "re") + s(fsampleID, bs = "re"), data = ef_data, pass = pass,
## id = sampleID)
##
## Formula:
## count ~ siteID + year + lifestage + pass12 + s(visitID, bs = "re") +
## s(fsampleID, bs = "re")
##
## Random effects:
## Estimate Std. Error
## sigma s(visitID) 0.01161 0.06545
## sigma s(fsampleID) 0.07751 0.16450
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z^2|)
## (Intercept) 0.648112 0.1778 3.64471 2.677e-04
## siteIDSite2 0.079936 0.1421 0.56244 5.738e-01
## siteIDSite3 -0.016579 0.1556 -0.10658 9.151e-01
## year2012 -0.208267 0.1798 -1.15807 2.468e-01
## year2013 0.042128 0.1931 0.21814 8.273e-01
## year2014 0.361853 0.1899 1.90504 5.678e-02
## year2015 -0.003164 0.2172 -0.01457 9.884e-01
## lifestageParr 0.588162 0.1253 4.69505 2.665e-06
## pass122 -0.236059 0.1318 -1.79112 7.327e-02
##
## Degrees of Freedom: 90 Total (i.e. Null); 79 Residual
## AIC: 3192The estimated standard deviations of the visit and sample random effects (on the logit scale) are
mfull$sigmaBecause the random effects are part of the likelihood, models can be compared directly using AIC
- Full model
AIC(mfull)- Model without lifestage
mlife <- efp(count ~ siteID + year + pass12 +
s(visitID, bs = "re") + s(fsampleID, bs = "re"),
data = ef_data, pass = pass, id = sampleID)
AIC(mlife)contributing
After making changes to the code, run:
devtools::document()to update help files and the NAMESPACE
then
devtools::check()to check for errors, and fix warnings and notes
Finally, the readme can be rebuilt using
rmarkdown::render("README.Rmd")then changes can be commited and pushed.