Package 'ef'

Modelling Framework for the Estimation of Salmonid Abundance

Description

A set of functions to estimate capture probabilities and densities from multipass removal (electrofishing) data. The main function is efp, which fits a model for capture probability, and predictions of capture probability are available from predict.efp.

Author(s)

Maintainer: Colin Millar [email protected]

Authors:

• Rob Fryer

• Iain Malcolm

• Ross Glover

Other contributors:

• Marine Scotland Science [copyright holder]

References

Millar, C. P., Fryer, R. J., Millidine, K. J. and Malcolm, I. A. (2016) Modelling capture probability of Atlantic salmon (Salmo salar) from a diverse national electrofishing dataset: Implications for the estimation of abundance. Fisheries Research, 177, 1-12. doi:10.1016/j.fishres.2016.01.001

Glover, R. S., Fryer, R. J., Soulsby, C., Bacon, P. J. and Malcolm, I. A. (2018) Incorporating estimates of capture probability and river network covariance in novel habitat - abundance models: Assessing the effects of conservation stocking on catchment-scale production of juvenile Atlantic salmon (Salmo salar) from a long-term electrofishing dataset. Ecological Indicators, 93, 302-315. doi:10.1016/j.ecolind.2018.05.013

See Also

Useful links:

• https://faskally.github.io/ef/

• https://github.com/Faskally/ef

• Report bugs at https://github.com/Faskally/ef/issues

---

Conversion of an efp fit to a gam fit for plotting

Description

Converts an efp fit to a gam object, so that the plotting and prediction tools in the mgcv package can be used. The coefficients and their variances are those estimated by efp, but the smoothing parameters and effective degrees of freedom in the gam object are not.

Usage

as.gam(object)

Arguments

object

an efp fitted model

Details

The conversion fits a model using mgcv, which can take a very long time if the model has many random effects. Predictions do not need the conversion, see predict.efp.

Value

an object of class gam, see gamObject

Examples

data(ef_data)
ef_data$sampleID <- factor(paste(ef_data$date, ef_data$siteID,
                                 ef_data$lifestage, ef_data$species))
m1 <- efp(count ~ lifestage + s(doy, k = 3),
          data = ef_data, pass = pass, id = sampleID)

m1_gam <- as.gam(m1)
plot(m1_gam)

---

Counts of salmon fry and parr from electrofishing.

Description

A dataset containing counts from multipass electrofishing samples over three sites and 5 years.

Usage

ef_data

Format

A data frame with 90 rows and 10 variables:

siteID

Unique identifier for each electrofishing site

year

year of data collection

date

date of data collections

doy

day of the year of data collection, a continuous covariate which can be used to fit smooth terms

EFpasscount

the total number of electrofishing passes

pass

the electrofishing pass on which the fish were caught

species

fish species, Salmon

lifestage

fish lifestage, Fry or Parr

count

the number of fish caught per pass for each site visit and species, etc.

area

the area of river fished

Source

https://www.gov.scot/publications/scotland-river-temperature-monitoring-network-srtmn/

---

Estimate capture probabilities from electrofishing data

Description

This function uses the marginal likelihood of capture probabilities to estimate model parameters

Usage

efp(
  formula,
  data = NULL,
  pass = pass,
  id = id,
  offset = NULL,
  verbose = FALSE,
  hessian = TRUE,
  fit = TRUE,
  control = list(),
  sparse = TRUE
)

Arguments

formula	a formula object. Random effects are specified using s(f, bs = "re") terms, where f is a factor, and smooths using other s() terms, see Details.
data	a data.frame containing all relevant info
pass	a vector of integers giving the pass number of the observation
id	a vector of integers identifying an observation (a set of electrofishing passes)
offset	an optional offset for the linear predictor of capture probability
verbose	if TRUE optimiser messages are printed to the screen
hessian	if TRUE the hessian is computed and the covariance matrix of the parameters is returned via Vb
fit	if TRUE model is fitted if FALSE the data that would be passed to the optimiser is returned
control	a list of control commands to be passed to optim
sparse	if TRUE (the default) penalised smooths are set up so that sparse bases and penalties stay sparse, see Details. If FALSE the constrained basis and penalties from mgcv are used, which are dense and so slower to fit.

Details

Random effects on the logit capture probability are included by adding s(f, bs = "re") terms to the formula, where f is a factor in data. For example, s(visitID, bs = "re") adds between visit variation and s(sampleID, bs = "re") adds between sample variation. Each term has its own standard deviation, and the random effects are integrated out of the likelihood using the Laplace approximation, so that logLik and AIC are based on the marginal likelihood. Fitted values are conditional on the estimated random effects.

Smooth terms can also be included in the formula. These are penalised, and the amount of smoothing is estimated: a smooth is fitted as a random effect, in which the part of the smooth that is not penalised (for example a linear trend) is a fixed effect, and the standard deviation of the random part controls how smooth the term is. A standard deviation close to zero means that only the unpenalised part of the smooth is needed. As for other random effects, logLik and AIC are based on the marginal likelihood.

The recommended smooths are first and second order random walks, using the Markov random field basis of mgcv with a penalty from rw_penalty: s(f, bs = "mrf", xt = list(penalty = rw_penalty(levels(f)))), where f is a factor whose levels are the steps of the walk. The basis and the penalty are sparse, and are kept sparse when fitting, so these models are quick to fit even when there are many levels. The standard deviation is that of the steps of the random walk, and the starting values of the walk (the first level for a first order walk, the first two for a second order walk) are fixed effects: this means that the intercept is the logit capture probability at the first level of the walk, and that a second order walk has an additional fixed effect for its linear trend, which is on an arbitrary scale.

Any other penalised smooth from mgcv can be used, such as s(x) or te(x, z), but their penalties are dense so they will take longer to fit, particularly if the basis dimension k is large. To fit an unpenalised regression spline, whose flexibility is set entirely by k, use fx = TRUE, for example s(x, k = 3, fx = TRUE).

The fixed effects are estimated by maximum likelihood, so the likelihood depends slightly on how the unpenalised part of a smooth is separated from the penalised part. With sparse = TRUE one coefficient of the smooth is fixed for each function that is not penalised, as described above for random walks. With sparse = FALSE, and for smooths with by variables or models without an intercept, the identifiability constraints of mgcv are applied to the smooth first, which makes the basis and the penalty dense. The two approaches give very similar, but not identical, fits.

Value

If fit = TRUE, an object of class efp, a list which inherits from glm, containing among other things the fixed effect estimates (coefficients), the random effect estimates (random), the standard deviations of the random effects and smooths (sigma), the fitted capture probabilities (fitted), the log likelihood (llik) and the variance matrix of the fixed effects (Vb). See efp-methods for the methods available. If fit = FALSE, a list containing the data that would be passed to the optimiser.

References

Millar, C. P., Fryer, R. J., Millidine, K. J. and Malcolm, I. A. (2016) Modelling capture probability of Atlantic salmon (Salmo salar) from a diverse national electrofishing dataset: Implications for the estimation of abundance. Fisheries Research, 177, 1-12. doi:10.1016/j.fishres.2016.01.001

Examples

# create two electrofishing site visits with 3 and 4 passes and 2 lifestages
ef_data <- data.frame(n      = c(100, 53, 24, 50, 26, 12,
                                 100, 53, 24, 50, 26, 12),
                      pass   = c(  1,  2,  3,  1,  2,  3,
                                   1,  2,  3,  1,  2,  3),
                      stage  = c(  1,  1,  1,  2,  2,  2,
                                   1,  1,  1,  2,  2,  2),
                      sample = c(  1,  1,  1,  2,  2,  2,
                                   3,  3,  3,  4,  4,  4))

ef_data2 <- data.frame(n      = c(100, 53, 24, 50, 26, 12,
                                  100, 53, 24, 12, 50, 26, 12, 6),
                       pass   = c(  1,  2,  3,  1,  2,  3,
                                    1,  2,  3,  4,  1,  2,  3, 4),
                       stage  = c(  1,  1,  1,  2,  2,  2,
                                    1,  1,  1,  1,  2,  2,  2, 2),
                       sample = c(  1,  1,  1,  2,  2,  2,
                                    3,  3,  3,  3,  4,  4,  4, 4))

ef_data3 <- data.frame(n      = c(100, 53, 24, 50, 26, 12, 40,
                                  100, 53, 24, 12, 50, 26, 12, 6, 40),
                       pass   = c(  1,  2,  3,  1,  2,  3, 1,
                                    1,  2,  3,  4,  1,  2,  3, 4, 1),
                       stage  = c(  1,  1,  1,  2,  2,  2, 1,
                                    1,  1,  1,  1,  2,  2,  2, 2, 2),
                       sample = c(  1,  1,  1,  2,  2,  2, 5,
                                    3,  3,  3,  3,  4,  4,  4, 4, 6))

# Fit a simple model
m2 <- efp(n ~ 1 + factor(stage), data = ef_data, pass = pass, id = sample)
cbind(ef_data, fit = fitted(m2))
m3 <- efp(n ~ 1 + factor(stage), data = ef_data2, pass = pass, id = sample)
cbind(ef_data2, fit = fitted(m3))
m4 <- efp(n ~ 1 + factor(stage), data = ef_data3, pass = pass, id = sample)
cbind(ef_data3, fit = fitted(m4))

# create two electrofishing site visits with 3 and 4 passes and 2 lifestages
ef_data <- data.frame(n      = c(200, 53, 24, 100, 26, 12,
                                 200, 53, 24, 100, 26, 12),
                      pass   = c(  1,  2,  3,  1,  2,  3,
                                   1,  2,  3,  1,  2,  3),
                      stage  = c(  1,  1,  1,  2,  2,  2,
                                   1,  1,  1,  2,  2,  2),
                      sample = c(  1,  1,  1,  2,  2,  2,
                                   3,  3,  3,  4,  4,  4))
# Fit a simple model
m2 <- efp(n ~ 1 + factor(stage) + factor(replace(pass, pass> 2, 2)),
          data = ef_data, pass = pass, id = sample)
out <- cbind(ef_data, p = fitted(m2, type = "p"))
out

# a smooth effect of a continuous covariate, day of the year, using a
# second order random walk
data(ef_data)
ef_data$sampleID <- factor(paste(ef_data$date, ef_data$siteID,
                                 ef_data$lifestage, ef_data$species))
ef_data$pass12 <- factor(replace(ef_data$pass, ef_data$pass > 2, 2))
days <- min(ef_data$doy):max(ef_data$doy)
ef_data$fdoy <- factor(ef_data$doy, levels = days)
rw2 <- rw_penalty(days, order = 2)
m5 <- efp(count ~ siteID + pass12 + lifestage +
            s(fdoy, bs = "mrf", xt = list(penalty = rw2)),
          data = ef_data, pass = pass, id = sampleID)
m5

# predicted capture probability over the smooth, with standard errors
newdata <- data.frame(siteID = "Site1", pass12 = "1", lifestage = "Fry",
                      fdoy = factor(days, levels = days))
pred <- predict(m5, newdata = newdata, type = "response", se.fit = TRUE)
cbind(newdata, fit = pred$fit, se = pred$se.fit)

# an unpenalised regression spline
m6 <- efp(count ~ siteID + pass12 + lifestage + s(doy, k = 3, fx = TRUE),
          data = ef_data, pass = pass, id = sampleID)
m6

---

Methods for efp fits

Description

Methods to summarise and extract information from a model fitted using efp.

Usage

## S3 method for class 'efp'
fitted(object, type = "p", ...)

## S3 method for class 'efp'
vcov(object, ...)

## S3 method for class 'efp'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'efp'
logLik(object, ...)

## S3 method for class 'efp'
formula(x, ...)

## S3 method for class 'efp'
summary(object, ...)

## S3 method for class 'summary.efp'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'efp'
residuals(object, ...)

## S3 method for class 'efp'
anova(object, ...)

## S3 method for class 'efp'
plot(x, ...)

## S3 method for class 'efp'
model.matrix(object, ...)

## S3 method for class 'efp'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

object, x	an efp fitted model
type	the type of fitted value: "p" for the capture probability in each pass, "pi" for the probability that a fish is first caught in each pass, "both" for a list containing both and "lpmatrix" for the design matrix.
...	further arguments, ignored unless stated
parm	the parameters to compute confidence intervals for, either a vector of numbers or names. The default is all fixed effects.
level	the confidence level required
digits	the number of significant digits to print

Details

The variance matrix returned by vcov is for the fixed effects only, the joint variance matrix of the fixed and random effects is available as object$Vp. Confidence intervals are Wald intervals based on this variance matrix.

The log likelihood is the marginal log likelihood when the model contains random effects or penalised smooths, and the degrees of freedom are the number of fixed effects plus the number of standard deviations.

Residuals, analysis of deviance tables and diagnostic plots are not available for efp fits, and the corresponding methods return an error.

Value

fitted returns a vector, list or matrix depending on type; model.matrix a matrix; vcov a matrix; confint a matrix with columns giving the lower and upper confidence limits; logLik an object of class logLik; formula the model formula; summary an object of class summary.efp, a list containing the coefficient table (coefficients) and the table of random effect standard deviations (sigma). The print methods return their argument invisibly.

Examples

data(ef_data)
ef_data$sampleID <- factor(paste(ef_data$date, ef_data$siteID,
                                 ef_data$lifestage, ef_data$species))
m1 <- efp(count ~ lifestage, data = ef_data, pass = pass, id = sampleID)

summary(m1)
vcov(m1)
confint(m1)
logLik(m1)
head(fitted(m1))

---

Inverse logistic transformation

Description

This function transforms values on the logistic scale to values on the probability scale.

Usage

invlogit(x)

Arguments

x	a numeric vector

Value

vector of values between 0 and 1

Examples

invlogit(c(-2, 0, 2))

---

Logistic transformation

Description

This function transforms values on the probability scale to values on the logistic scale.

Usage

logit(p)

Arguments

p	a numeric vector with values between 0 and 1

Value

vector of values between -Inf and Inf

Examples

logit(c(0.1, 0.5, 0.9))

---

Predict capture probabilities from an efp fit

Description

Predictions of capture probability from a fitted efp model, for the data used to fit the model or for new data.

Usage

## S3 method for class 'efp'
predict(
  object,
  newdata = NULL,
  type = c("link", "response", "lpmatrix"),
  se.fit = FALSE,
  exclude = NULL,
  ...
)

Arguments

object	an efp fitted model
newdata	a data.frame containing the covariates in the model formula. If NULL, predictions are for the data used to fit the model.
type	the type of prediction: "link" for the logit capture probability, "response" for the capture probability and "lpmatrix" for the design matrix which gives the logit capture probability when multiplied by the model parameters.
se.fit	if TRUE standard errors of the predictions are also returned
exclude	a character vector of terms to leave out of the prediction, named as in the model output, for example "s(visitID)". This is the way to predict for levels of a random effect that were not in the data used to fit the model. Variables that are only used by excluded smooth or random effect terms need not be supplied in newdata.
...	further arguments, ignored

Details

Predictions are conditional on the estimated random effects, unless these are removed using exclude. For a random walk smooth, predictions can be made for any level in the penalty, including levels with no data, which are interpolated. Standard errors on the response scale are calculated using the delta method.

The offset supplied to efp is included in the predictions only when newdata is NULL.

Value

a vector of predictions, or if se.fit = TRUE a list with elements fit and se.fit. If type = "lpmatrix" a matrix.

Examples

data(ef_data)
ef_data$sampleID <- factor(paste(ef_data$date, ef_data$siteID,
                                 ef_data$lifestage, ef_data$species))
m1 <- efp(count ~ lifestage, data = ef_data, pass = pass, id = sampleID)

newdata <- data.frame(lifestage = c("Fry", "Parr"))
predict(m1, newdata = newdata, type = "response", se.fit = TRUE)

---

Random walk penalties for smooth terms

Description

Creates a first or second order random walk penalty matrix for the levels of an ordered covariate, for use with the Markov random field (mrf) basis of mgcv. This is the recommended way to include smooth terms in efp: the basis and the penalty are both sparse, so models are quick to fit even when there are many levels.

Usage

rw_penalty(levels, order = 1)

Arguments

levels	the levels of the covariate, in order, and assumed to be equally spaced. A level is needed for every step of the random walk, including any for which there are no data.
order	the order of the random walk, 1 or 2

Details

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. In both cases the standard deviation estimated by efp is the standard deviation of the steps of the random walk, on the logit scale.

To use the penalty, the covariate must be a factor with the same levels as the penalty, see the examples.

Value

a square matrix with one row and column for each level, with the levels as dimnames. The matrix has an attribute, null.space, which is a basis for the functions which are not penalised.

Examples

rw_penalty(1:5)
rw_penalty(1:5, order = 2)

# a first order random walk over day of the year
data(ef_data)
ef_data$sampleID <- factor(paste(ef_data$date, ef_data$siteID,
                                 ef_data$lifestage, ef_data$species))
days <- min(ef_data$doy):max(ef_data$doy)
ef_data$fdoy <- factor(ef_data$doy, levels = days)
rw1 <- rw_penalty(days)

m1 <- efp(count ~ lifestage + s(fdoy, bs = "mrf", xt = list(penalty = rw1)),
          data = ef_data, pass = pass, id = sampleID)
m1

---

Utility function to convert parameters to probabilities

Description

The matrix G should be of dimension n x p, and the parameter vector should be length p

Usage

transpar(par, G)

Arguments

par	fitted model parameters
G	The design matrix for a model

Value

a vector of probabilities of length n

Examples

G <- cbind(1, c(0, 0, 1, 1))
transpar(c(0.5, -0.2), G)

---