CHAPTER 7 - ESTIMATION OF PARAMETERS

In the previous chapters, several models used in stock assessment were analysed, the respective parameters having been defined. In the corresponding exercises, it was not necessary to estimate the values of the parameters because they were given. In this chapter, several methods of estimating parameters will be analysed. In order to estimate the parameters, it is necessary to know the sampling theory and statistical inference.

This manual will use one of the general methods most commonly used in the estimation of parameters - the least squares method. In many cases this method uses iterative processes, which require the adoption of initial values. Therefore, particular methods will also be presented, which obtain estimates close to the real values of the parameters. In many situations, these initial estimates also have a practical interest. These methods will be illustrated with the estimation of the growth parameters and the S-R stock-recruitment relation.

The least squares method is presented under the forms of Simple linear Regression, multiple linear model and non linear models (method of Gauss-Newton).

Subjects like residual analysis, sampling distribution of the estimators (asymptotic or empiric Bookstrap and jacknife), confidence limits and intervals, etc., are important. However, these matters would need a more extensive course.

7.1 SIMPLE LINEAR REGRESSION - LEAST SQUARES METHOD

Model

Consider the following variables and parameters:

The response variable is linear with the parameters

Y = A+BX

Objective

The objective of the method is to estimate the parameters of the model, based on the observed pairs of values and applying a certain criterium function (the observed pairs of values are constituted by selected values of the auxiliary variable and by the corresponding observed values of the response variable), that is:

Object function (or criterium function)

Estimation method

In the least squares method the estimators are the values of A and B which minimize the object function. Thus, one has to calculate the derivatives ∂Φ/∂A e ∂Φ/∂B, equate them to zero and solve the system of equations in A and B.

The solution of the system can be presented as:

Notice that the observed values y, for the same set of selected values of X, depend on the collected sample. For this reason, the problem of the simple linear regression is usually presented in the form:

y = A + BX + ε

where ε is a random variable with expected value equal to zero and variance equal to σ².

So, the expected value of y will be Y or A+BX and the variance of y will be equal to the variance of ε.

The terms deviation and residual will be used in the following ways:

Deviation is the difference between y_observed and y_mean () i.e., deviation = (y-)

while

Residual is the difference between y_observed and Y_estimated (), i.e., residual =.

To analyse the adjustment of the model to the observed data, it is necessary to consider the following characteristics:

Sum of squares of the residuals:

This quantity indicates the residual variation of the observed values in relation to the estimated values of the response variable of the model, which can be considered as the variation of the observed values that is not explained by the model.

Sum of squares of the deviations of the estimated values of the response variable of the model:

This quantity indicates the variation of the estimated values of the response variable of the model in relation to its mean, that is the variation of the response variable explained by the model.

Total sum of squares of the deviations of the observed values equal to:

This quantity indicates the total variation of the observed values in relation to the mean

It is easy to verify the following relation:

SQ_total = SQ_model + SQ_residual

or

or

1 = r² + (1 - r²)

where

r² (coefficient of determination) is the percentage of the total variation that is explained by the model and

1-r² is the percentage of the total variation that is not explained by the model.

7.2 MULTIPLE LINEAR REGRESSION - LEAST SQUARES METHOD

Model

Consider the following variables and parameters:

The response variable is linear with the parameters

Y = B₁X₁+B₂X₂+... + B_kX_k = Σ B_jX_j

Objective

The objective of the method is to estimate the parameters of the model, based on the observed n sets of values and by applying a certain criterium function (the observed sets of values are constituted by selected values of the auxiliary variable and by the corresponding observed values of the response variable), that is:

Observed values x_1,i x_2,i,..., x_j,i,.., x_k,i and y_i for each set i, where i=1,2,...,i,...n

Values to be estimated B₁,B₂,...,B_j,...,B_k et (Y₁,Y₂,..., Y_i,..., Y_n)

The estimated values can be represented by:

(ou b₁,b₂,...,b_j,...,b_k) et

Object function (or criterium function)

Estimation method

In the least squares method the estimators are the values of B_j which minimize the object function.

As with the simple linear model, the procedure of minimization requires equating the partial derivatives of Φ to zero in order to each parameter, B_j, where j=1, 2,..., k. The system is preferably solved using matrix calculus.

Matrix version

Matrix X_(n,k) = Matrix of the n observed values of each of the k auxiliary variables
Vector y_(n,1) = Vector of the n observed values of the response variable
Vector Y_(n,1) = Vector of the values of the response variable given by the model (unknown)
Vector B_(k,1) = Vector of the parameters
Vector or b_(k,1) = Vector of the estimators of the parameters

Model

Y_(n,1) = X_(n,k). B_(k,1) ou Y=X.B+ε

Object function

Φ_(1,1) = (y-Y)^T.(y-Y) ou Φ_(1,1) = (y-X.B)^T.(y-X.B)

To calculate the least squares estimators it will suffice to put the derivative dΦ/dB of Φ in order to vector B, equal to zero. dΦ/dB is a vector with components ∂Φ/∂B₁, ∂Φ/∂B₂,..., ∂Φ/∂B_k. Thus:

dΦ/dB_(k,1) = -2.X^T.(y-X.B) = 0

or X^Ty - (X^T.X). B = 0

and b = = (X^T.X)^-1. X^Ty

The results can be written as:

b_(k,1) = (X^T.X)^-1.X^Ty

= X.b or = X (X^T.X)^-1.X^T y

residuals_(n,1) = (y-)

Comments

In statistical analysis it is convenient to write the estimators and the sums of the squares using idempotent matrices. Then the idempotent matrices L, (I - L) and (I - M) with L_(n,n) = X (X^T. X)^-1. X^T, I = unity matrix and M_(n,n) = mean_(n,1) matrix = 1/n [1] where [1] is a matrix with all its elements equal to one, are used.

It is also important to consider the sampling distributions of the estimators assuming that the variables ε_i are independent and have a normal distribution.

A summary of the main properties of the expected value and variance of the estimators is presented:

6 - Sum of squares

6.1 - Residual Sum of squares = SQ residual_(1.1) = (y-)^T(y-) = y^T (I-L)y

This quantity indicates the residual variation of the observed values in relation to the estimated values of the model, that is, the variation not explained by the model.

6.2 - Sum of squares of the deviation of the model = SQ model_(1.1) = (-)^T(-) = y^T (L-M)y

This quantity indicates the variation of the estimated response values of the model in relation to the mean, that is, the variation explained by the model.

6.3 - Total Sum of the squares of the deviations = SQ total_(1.1) = (y-)^T(y-) = y^T (I-M) y

This quantity indicates the total variation of the observed values in relation to the mean.

It is easy to verify the following relation:

SQ_total = SQ_model + SQ_residual or

or 1 = R² + (1 - R²)

where:

R² is the percentage of the total variation that is explained by the model. In matrix terms it will be:

R² = [y^T(L - M)y].[ (y^T(I - M)y]^-1

1-R² is the percentage of the total variation that is not explained by the model.

The ranks of the matrices (I-L), (I-M) and (L-M) respectively equal to (n-k), (n-1) and (k-1), are the degrees of freedom associated with the respective sums of squares.

7.3 NON-LINEAR MODEL - METHOD OF GAUSS-NEWTON - LEAST SQUARES METHOD

Model

Consider the following variables and parameters:

The response variable is non-linear with the parameters

Y = f(X;B) where B is a vector with the components B₁,B₂,...,B_j,...,B_k

Objective

The objective of the method is to estimate the parameters of the model, based on the n observed pairs of values and by applying a certain criterium function (the observed sets of values are constituted by selected values of the auxiliary variable and by the corresponding observed values of the response variable), that is:

Observed values x_i and y_i for each pair i, where i=1,2,...,i,...n

Values to be estimated B₁,B₂,...,B_j,..,B_k and (Y₁,Y₂,...,Y_i,...,Y_n) form the n pairs of observed values.

(Estimates = or b₁,b₂,...,b_j,...,b_k and )

Object function or criterium function

Estimation criterium

The estimators will be the values of B_j for which the object function is minimum.

(This criterium is called the least squares method).

Matrix version

It is convenient to present the problem using matrices.

So:

Vector X_(n,1) = Vector of the observed values of the auxiliary variable
Vector y_(n,1) = Vector of the observed values of the response variable
Vector Y_(n,1) = Vector of the values of the response variable given by the model
Vector B_(k,1) = Vector of the parameters
Vector b_(k,1) = Vector of the estimators of the parameters

Model

Y_(n,1) = f(X; B)

Object function

Φ_(1,1) = (y-Y)^T.(y-Y)

In the case of the non linear model, it is not easy to solve the system of equations resulting from equating the derivative of the function Φ in order to the vector B, to zero. Estimation by the least squares method can, based on the Taylor series expansion of function Y, use iterative methods.

Revision of the Taylor series expansion of a function

Here is an example of the expansion of a function in the Taylor series in the case of a function with one variable.

The approximation of Taylor means to expand a function Y = f(x) around a selected point, x_0, in a power series of x:

Y = f(x) = f(x₀) +(x-x₀).f’(x₀)/1! + (x-x₀)²f’’(x₀)/2! +... + (x- x₀)ⁱ f⁽ⁱ⁾(x₀)/i!+...

where

f⁽ⁱ⁾(x₀) = i^th derivatives of f(x) in order to x, at the point x₀.

The expansion can be approximated to the desired power of x. When the expansion is approximated to the power 1 it is called a linear approximation, that is,

Y ≅ f(x₀) + (x-x₀).f’(x₀)

The Taylor expansion can be applied to functions with more than one variable. For example, for a function Y = f(x₁,x₂) of two variables, the linear expansion would be:

which may be written, in matrix notation, as

Y = Y₍₀₎+A₍₀₎.(x-x₍₀₎)

where Y₍₀₎ is the value of the function at the point x₍₀₎,with components x₁₍₀₎ and x₂₍₀₎,and A₍₀₎ is the matrix of derivatives whose elements are equal to the partial derivatives of f(x₁,x₂) in order to x₁,x₂ at the point (x₁₍₀₎, x₂₍₀₎).

To estimate the parameters, the Taylor series expansion of function Y is made in order to the parameters B and not to the vector X.

For example, the linear expansion of Y = f(x,B) in B₁, B₂,..., B_k, would be:

Y = f(x;B) = f(x; B₍₀₎) + (B₁-B₁₍₀₎) f /B₁ (x;B₍₀₎) +..... +
(B₂-B₂₍₀₎)f /B₂ (x;B₍₀₎) +...... +..........+ (B_k-B_k(0)) f /B_k (x;B₍₀₎)

or, in matrix notation, it would be:

Y_(n,1) = Y_{(0) (n,1)} + A_{(0) (n,k).} ΔB_{(0) (k,1)}

where

A = matrix of order (n,k) of the partial derivatives of the matrix f(x;B) in order to the vector B at the point B₍₀₎ and

ΔB₍₀₎ = vector (B - B₍₀₎).

Then, the object function will be:

Φ = (y-Y)^T.(y-Y) = (y-Y₍₀₎ - A₍₀₎. ΔB₍₀₎)^T(y-Y₍₀₎ - A_(0). ΔB₍₀₎)

To obtain the minimum of this function it is more convenient to differentiate Φ in order to the vector ΔB than in relation to vector B and put it equal to zero. Thus:

or

Therefore:

If ΔB₍₀₎ is "equal to zero" then the estimate of B is equal to B₍₀₎.

(In practice, when we say "equal to zero" in this process, we really mean smaller than the approximation vector one has to define beforehand).

If ΔB₍₀₎ is not "equal to zero" then the vector B₍₀₎ will be replaced by:

B₍₁₎ =B₍₀₎ + ΔB₍₀₎

And the process will be repeated, that is, there will be another iteration with B₍₀₎ replaced by B₍₁₎ (and A₍₀₎ replaced by A₍₁₎). The iterative process will go on until the convergence at the desired level of approximation is reached.

Comments

1. It is not guaranteed that the process always converges. Sometimes it does not, some other times it is too slow (even for computers!) and some other times it converges to another limit!!

2. The above described method is the Gauss-Newton method which is the basis of many other methods. Some of those methods introduce modifications in order to obtain a faster convergence like the Marquardt method (1963), which is frequently used in fisheries research. Other methods use the second order Taylor expansion (Newton-Raphson method), looking for a better approximation. Some others, combine the two modifications.

3. These methods need the calculation of the derivatives of the functions. Some computer programs require the introduction of the mathematical expressions of the derivatives, while others use sub-routines with numerical approximations of the derivatives.

4. In fisheries research, there are methods to calculate the initial values of the parameters, for example in growth, mortality, selectivity or maturity analyses.

5. It is important to point out that the convergence of the iterative methods is faster and more likely to approach the true limit when the initial value of the vector B₍₀₎ is close to the real value.

7.4 ESTIMATION OF GROWTH PARAMETERS

The least squares method (non-linear regression) allows the estimation of the parameters K, L_∞ and t_o of the individual growth equations.

The starting values of K, L_∞ and t₀ for the iterative process of estimation can be obtained by simple linear regression using the following methods:

Ford-Walford (1933-1946) and Gulland and Holt (1959) Methods

The Ford-Walford and Gulland and Holt expressions, which were presented in Section 3.4, are already in their linear form, allowing the estimation of K and L_∞ with methods of simple linear regression on observed L_i and T_i. The Gulland and Holt expression allows the estimation of K and L_∞ even when the intervals of time T_i are not constant. In this case, it is convenient to re-write the expression as:

Stamatopoulos and Caddy Method (1989)

These authors also present a method to estimate K, L_∞ and t_o (or L_o) using the simple linear regression. In this case the von Bertalanffy equation should be expressed as a linear relation of L_t against e^-Kt.

Consider n pairs of values t_i, L_i where t_i is the age and L_i the length of the individual i where i=1,2,...., n.

The von Bertalanffy equation, in its general form is (as previously seen):

L_∞ - L_t = (L_∞- L_a). e^-K(t-ta)

It can be written as:

L_t = L_∞ - (L_∞- L_a). e^+Kta. e^-Kt

The equation has the simple linear form, y = a + bx, where:

If one takes L_a = 0, then t_a=t_o, but, if one considers t_a = 0, then L_a = L_o.

The parameters to estimate from a and b will be L_∞, t_o or L_o.

The authors propose adopting an initial value K₍₀₎, of K, and estimating a₍₀₎, b₍₀₎ and r²₍₀₎ by simple linear regression between y (= L_t) and x(=e^k₍₀₎). The procedure may be repeated for several values of K, that is, K₍₁₎ K₍₂₎,.... One can then adopt the regression that results in the larger value of r², to which K_max, a_max and b_max correspond. From the values of a_max, b_max and K_max one can obtain the values of the remaining parameters.

One practical process towards finding K_max can be:

(i). To select two extreme values of K which include the required value, for example K= 0 and K=2 (for practical difficulties, use K = 0.00001 instead of K = 0).

(ii). Calculate the 10 regressions for equally-spaced values of K between those two values in regular intervals.

(iii). The corresponding 10 values of r² will allow one to select two new values of K which determine another interval, smaller than the one in (i), containing another maximum value of r².

(iv). The steps (ii) and (iii) can be repeated until an interval of values of K with the desired approximation is obtained. Generally, the steps do not need many repetitions.

7.5 ESTIMATION OF M - NATURAL MORTALITY COEFFICIENT

Several methods were proposed to estimate M, and they are based on the association of M with other biological parameters of the resource. These methods can produce approximate results.

7.5.1 RELATION OF M WITH THE LONGEVITY,

Longevity: Maximum mean age t_λ of the individuals in a non-exploited population.

Duration of the exploitable life: t_λ - t_r = λ (Figure 7.1)

Figure 7.1 Duration of the exploitable life

Tanaka (1960) proposes "NATURAL" Survival Curves (Figure 7.2) to obtain the values of M from longevity.

A cohort practically vanishes when only a fraction, p, of the recruited individuals survives. In that case, N_λ = R · e^-M·λ, and it can be written:

Different values of the survival fraction produce different survival curves of M in function of λ.

Figure 7.2 Survival curves by Tanaka

Any value of p can be chosen, for instance, p = 5%, (i.e. one in each twenty recruits survives until the age t_λ) as variable value of the survival curves.

7.5.2 RELATION BETWEEN M AND GROWTH

Beverton and Holt Method (1959)

Gulland (1969) mentions that Beverton and Holt verified that species with a larger mortality rate M also presented larger values of K. Looking for a simple relation between these two parameters, they concluded approximately that:

Pauly Method (1980)

Based on the following considerations:

1. Resources with a high mortality rate cannot have a very big maximum size;
2. In warmer waters, the metabolism is accelerated, so the individuals can grow up to a larger size and reach the maximum size faster than in colder waters.

Based on data of 175 species, Pauly adjusted multiple linear regressions of transformed values of M against the corresponding transformed values of K, L_∞ and temperature, T, and selected one that was considered to have a better adjustment, that is, the following empirical relation:

with the parameters expressed in the following units:

M = year^-1
L_∞ = cm of total length
K = year^-1
T° = surface temperature of the waters in °C

Pauly highlights the application of this expression to small pelagic fishes and crustaceans. The Pauly relation uses decimal logarithms to present the first coefficient different from the value -0.0152 which was given in the previous expression, written with natural logarithms.

7.5.3 RELATION BETWEEN M AND REPRODUCTION

Rikhter and Efanov Method (1976)

These authors analysed the dependency between M and the age of first (or 50 percent) maturity. They used data from short, mean and long life species, and suggested the following relation of M with the, t_mat, age of 1^st maturity:

Gundersson Method (1980)

Based on the assumption that the natural mortality rate should be related to the investment of the fish in reproduction, beyond the influence of other factors, Gundersson established several relations between M and those factors.

He proposed, however, the following simple empirical relation, using the Gonadosomatic Index (GSI) (estimated for mature females in the spawning period) in order to calculate M:

M = 4.64xGSI - 0.37

7.5.4 KNOWING THE STOCK AGE STRUCTURE, AT BEGINNING AND END OF YEAR, AND CATCHES IN NUMBER, BY AGE, DURING THE YEAR

The natural mortality coefficients M_i, at age i can be calculated from the catch, C_i, in numbers, and the survival numbers, N_i and N_i+1 at the beginning and end of a year, by following the steps:

calculate
calculate
calculate

The several values of M obtained in each age could be combined to calculate a constant value, M, for all ages.

Paloheimo Method (1961)

Let us consider the supposition that F_i is proportional to f_i for several years i, that is

for T_i = 1 year, F_i = q · f_i,

Z_i = q · f_i + M

So, the linear regression between Z_i and f_i has a slope b = q and an intercept a = M.

7.6 ESTIMATION OF Z - TOTAL MORTALITY COEFFICIENT

There are several methods of estimating the total mortality coefficient, Z, assumed to be constant during a certain interval of ages or years.

It is convenient to group the methods, according to the basic data, into those using ages or those using lengths.

7.6.1 METHODS USING AGE DATA

The different methods are based on the general expression of the number of survivors of a cohort, at the instant t, submitted to the total mortality, Z, during an interval of time, that is:

Z is supposed to be constant in the interval of time (t_a,t_b).

Taking logarithms and re-arranging the terms, the expression will be:

lnN_t = Cte - Z.t

where Cte is a constant (= ln N_a+Zt_a).

This expression shows that the logarithm of the number of survivors is linear with the age, being the slope equal to -Z.

Any constant expression which does not affect the determination of Z will be referred to as Cte.

1. If Z can be considered constant inside the interval (t_a,t_b) and, having available abundance data, N_i, or indices of abundance in number, U_i in several ages, i, then, the application of the simple linear regression allows one to estimate the total mortality coefficient Z.

In fact

so . Constant

and, as

then, by substitution:

(T_i = const = 1 year)

and also

The simple linear regression between and t_i allows the estimation of Z (notice that the constant, Cte is different from the previous one. In this case only the slope matters to estimate Z).

2. If ages are not at constant intervals, the expression could be approximated and expressed in terms of the t_centrali. For T_i variable, it will be:

3. When using indices U_i, the situation is similar because U_i = q. N_i, with q constant, and then, also:

The simple linear regression between and t_i allows one to estimate Z.

4. If the intervals are not constant, the expression should be modified to:

Simple linear regression can be applied to obtain Z, from catches, C_i, and ages, t_i, supposing that F_i is constant.

and so, when T_i is constant. So:

lnC_i = Cte - Z. t_i

5. If the intervals are not constant, the expression should be modified to:

lnC_i/T_i ≈ Cte - Z. t_centrali

6. Let V_i be the cumulative catch from t_i until the end of the life, then:

V_i = Σ C_k = Σ F_k N_kcum,

Where the sum goes from the last age until age i,

As F_k and Z_k are supposed to be constant ΣN_kcum = N_i/Z and so:

Therefore:

ln V_i = Cte - Z. t_i

7. Following Beverton and Holt (1956), Z can be expressed as:

Then, it is possible to estimate Z from the mean age t

This expression was derived, considering the interval (t_a, t_b) as (t_a, ∞).

7.6.2 METHODS USING LENGTH DATA

When one has available data by length classes instead of by age, the methods previously referred to can still be applied. For that purpose, it is convenient to define the relative age.

Using the von Bertalanffy equation one can obtain the age t in function of the length, as:

(the expression is written in the general form in relation to t_a and not to t₀)

t = t_a - (1/K).ln[(L_∞- L_t)/(L_∞- L_a)]

or

(This equation is referred to by some authors as the inverse von Bertalanffy equation).

The difference t-ta is called relative age, t^*,.

So: t^* =-(1/K).ln[(L_∞- L_t)/(L_∞- L_a)] or t^* =-(1/K)ln[1-(L_t-L_a)/ (L_∞- L_a)]

For t_a = t_o, L_a = 0 and:

t^* is called a relative age because the absolute ages, t, are related to a constant age, t_a.

In this way, the duration of the interval T_i can either be calculated by the difference of the absolute ages or by the difference of the relative ages at the extremes of the interval:

T_i = t_i+1 -t_i = t^*_{i +1} - t^*_i

Also:

t*_centrali = t_centrali + Cte

So, the previous expressions still hold when the absolute ages are replaced by the relative ages:

ln N_i = Cte - Z. t^*_centrali
ln U_i = Cte - Z. t^*_centrali
ln V_i = Cte - Z. t^*_i
ln C_i/T_i = Cte - Z. t^*_centrali

Finally, the expression would also be:

Beverton and Holt (1957) proved that:

must be calculated as the mean of the lengths weighted with abundances (or their indices) or with the catches in numbers.

Comments

1. The application of any of these methods must be preceeded by the graphical representation of the corresponding data, in order to verify if the assumptions of the methods are acceptable or not and also to determine the adequate interval, (t_a, t_b).

2. These formulas are proved with the indications that were presented, but it is a good exercise to develop the demonstrations as they clarify the methods.

3. It is useful to estimate a constant Z, even when it is not acceptable, because it gives a general orientation about the size of the values one can expect.

4. The methods are sometimes referred to by the names of the authors. For example, the expression ln V_i = Cte - Z.t^*_i is called the Jones and van Zalinge method (1981).

5. The mean age as well as the mean length in the catch can be calculated from the following expressions:

with C_i = catch in number in the age class i

where C_i = catch in number in the length class i

with C_i = catch in number in the age class.

The relative age should be t^* = - (1/K).ln[(L_∞- L_t)/(L_∞- L_a)]

Summary of the Methods to Estimate the Total Mortality Coefficient, Z

7.7 ESTIMATION OF THE PARAMETERS OF THE STOCK-RECRUITMENT (S-R) RELATION

The least squares method (non-linear model) can be used to estimate the parameters, α and k, of any of the S-R models.

The initial values of the Beverton and Holt model (1957) can be obtained by re-writing the equation as:

and estimating the simple linear regression between y (= S/R) and x (=S) which will give the estimations of 1/α and 1/(αk). From these values, it will then be possible to estimate the parameters α and k. These values can be considered as the initial values in the application of the non-linear model.

In the Ricker model (1954) the parameters can be obtained by re-writing the equation as:

and applying the simple linear regression between y (= ln R/S) and x (=S) to estimate ln α and (-1/k). From these values, it will be possible to estimate the parameters (α and k) of the model, which can be considered as the initial values in the application of the non-linear model.

It is useful to represent the graph of y against x in order to verify if the marked points are adjustable to a straight line before applying the linear regression in any of these models.

In the models with the flexible parameter, c, like for example, the Deriso model (1980), the equation can be re-written as:

For a given value of c the linear regression between y (= (R/S)^c) and x (=S) allows the estimation of the parameters α and k.

One can try several values of c to verify which one will have a better adjustment with the line y against x; for example, values of c between -1 and 1.

The values thus obtained for α, k and c, can be considered as initial values in the application of the iterative method, to estimate the parameters α, k and c of the non-linear Deriso model.

7.8 ESTIMATION OF THE MATRIX [F] AND OF THE MATRIX [N] -COHORT ANALYSIS - AC AND LCA

7.8.1 COHORT ANALYSIS BY AGE- (AC)

The cohort analysis is a method to estimate the fishing mortality coefficients, F_i, and the number of survivors, N_i, at the beginning of each age, from the annual structures of the stock catches, in number, over a period of years.

More specifically, consider a stock where the following is known:

Data

age, i, where i = 1,2,...,k
year, j, where j = 1,2,...,n
Matrix of catches [C] with
C_i,j = Annual catch, in number, of the individuals with the age i and during the year j
Matrix of natural mortality [M] with
M_i,j = natural mortality coefficient, at the age i and in the year j.
Vector [T] where
T_i = Size of the age interval i (in general, T_i=T=1 year)

Objective

To estimate

matrix [F]

and

matrix [N].

In the resolution of this problem, it is convenient to consider these estimations separately; one interval of age i (part 1); all the ages during the life of a cohort (part 2); and finally, all the ages and years (part 3).

PART 1 (INTERVAL T_I)

Consider that the following characteristics of a cohort, in an interval T_i are known:

C_i = Catch in number
M_i = Natural mortality coefficient
T_i = Size of the interval

Adopting a value of F_i, it is then possible to estimate the number of survivors at the beginning, N_i, and at the end, N_i+1, of the interval.

In fact, from the expression:

one can calculate N_i which is the only unknown variable in the expression.

To calculate N_i+1 one can use the expression where the values N_i, F_i and M_i were previously obtained.

PART 2 (DURING THE LIFE)

Suppose now that the catches C_i of each age i, of a cohort during its life, the values of M_i and the sizes of the interval T_i are known.

Adopting a certain value, F_final, for the Fishing Mortality Coefficient in the last class of ages, it is possible, as mentioned in part 1, to estimate all the parameters (related to numbers) in that last age group. In this way, one will know the number of survivors at the beginning and end of the last age.

The number at the beginning of that last class of ages, is also the number N_last at the end of the previous class, that is, N_final is the initial number of survivors of the class before last.

Using the C_i expression, resulting from the combination of the two expressions above:

one can estimate F_i in the previous class, which is the only unknown variable in the expression. The estimation may require iterative methods or trial and error methods.

Finally, to estimate the number N_i of survivors at the beginning of the class i, the following expression can be used:

Repeating this process for all previous classes, one will successively obtain the parameters in all ages, until the first age.

In the case of a completely caught cohort, the number at the end of the last class is zero and the catch C has to be expressed as:

Pope Method

Pope (1972) presented a simple method to estimate the number of survivors at the beginning of each age of the cohort life, starting from the last age.

It is enough to apply successively in a backward way, the expression:

N_i ≈ (N_i+1 e ^MT/2 + C_i).e ^MT/2

Pope indicates that the approximation is good when MT ≤ 0.6

Pope’s expression is obtained, supposing that the catch is made exactly at the central point of the interval T_i (Figure 7.3).

Figure 7.3 Number of survivors during the interval T_i = t_i+1 - t_i with the catch extracted at the central point of the interval

Proceeding from the end to the beginning one calculates successively:

N" = N_i+1e^+MTi/2
N’ = N" + C_i
N_i = N’.e^+MT/2

substituting N’ by N"+C_i, the expression will be:

N_i = (N" + C_i).e ^MT/2

Finally, substituting N" by N_i+1.e ^+MTi/2, it will be:

Part 3 (period of years)

Let us suppose now that the Catch matrix [C], the natural mortality [M] matrix and the vector size of the intervals [T], are known for a period of years.

Let us also assume that the values of F in the last age of all the years represented in the matrices and the values of F of all the ages of the last year were adopted. These values will be designated by F_terminal (Figure 7.4)

Figure 7.4 Matrix of catch, [C], with F_terminal in the last line and in the last column of the matrix C. The shadowed zones exemplify the catches of a cohort

Notice that in this matrix the elements of the diagonal correspond to values of the same cohort, because one element of a certain age and a certain year will be followed, in the diagonal, by the element that is a year older.

From parts 1 and 2 it will then be possible to estimate successively Fs and Ns for all the cohorts present in the catch matrix.

Comments

1. The values of M_i,j are considered constant and equal to M, when there is no information to adopt other values.

2. When data is referred to ages, the values T_i will be equal to 1 year.

3. The last age group of each year is, sometimes grouped ages(+). The corresponding catches are composed of individuals caught during those years, with several ages. So, the cumulative values do not belong to the same cohorts, but are survivors of several previous cohorts with different recruitments and submitted to different fishing patterns. It would not be appropriate to use the catch of a group (+) and to apply cohort analysis. Despite this fact, the group (+) is important in order to calculate the annual totals of the catches in weight, Y, of total biomasses, B, and the spawning stock biomass. So, it is usual to start with the cohort analysis on the age immediately before the group (+) and use the group (+) only to calculate the annuals Y, B and (SP). The value of F in that group (+) in each year, can be estimated as being the same fishing mortality coefficient as the previous age or, in some cases, as being a reasonable value in relation to the values of F_i in the year that is being considered.

4. A difficulty in the technical application appears when the number of ages is small or when the years are few. In fact, in those cases, the cohorts have few age classes represented in the Matrix [C] and the estimations will be very dependent on the adopted values of F_terminals.

5. The cohort analysis (CA) has also been designated as: VPA (Virtual Population Analysis), Derzhavin method, Murphy method, Gulland method, Pope method, Sequential Analysis, etc. Sometimes, CA is referred to when the Pope formula and the VPA are used in other cases. Megrey (1989) presents a very complete revision about the cohort analyses.

6. It is also possible to estimate the remaining parameters in an age i, related to numbers, that is, N_cumi, N_i, D_i, Z_i and E_i. When the information on initial individual or mean weights matrices [w] or [w] are available, one can also calculate the matrices of annual catch in weight [Y], of biomasses at the beginning of the years, [B], and of mean biomasses during the years [B]. If one has information on maturity ogives in each year, for example at the beginning of the year, spawning biomasses [SP] can also be calculated. Usually, only the total catches Y, the stock biomasses (total and spawning) at the beginning and the mean biomasses of the stock (total and spawning) in each year are estimated.

7. The elements on the first line of the matrix [N] can be considered estimates of the recruitment to the fishery in each year.

8. The fact that the F_terminals are adopted and that these values have influence on the resulting matrix [F] and matrix [N], forces the selection of values of F_terminals to be near the real ones. The agreement between the estimations of the parameters mentioned in the points 6. and 7. and other independent data or indices (for example, estimations by acoustic methods of recruitment or biomasses, estimations of abundance indices or cpue´s, of fishing efforts, etc) must be analysed.

9. The hypothesis that the exploitation pattern is constant from year to year, means that the fishing level and the exploitation pattern can be separated, or F_sepi = F_j x s_i. This hypothesis can be tested based on the matrix [ F ] obtained from the cohort analysis.

It is usual to call this separation VPA-Separable (SVPA).

Then, if F_ij = F_j.s_i one can prove that .

If the estimated values of F_ij are the same as the previous Fsep_ij = F_j.s_i then the hypothesis is verified. This comparison can be carried out in two different ways, the simplest is to calculate the quotients Fsep_ij /F_ij. If the hypothesis is true this quotient is equal to one. If the hypothesis is not verified it is always possible to consider other hypotheses with the annual vector [s] constant in some years only, mainly the last years.

10. It is usual to consider an interval of ages, where it can be assumed that the individuals caught are "completely recruited". In that case, the interval of ages corresponds to exploitation pattern constant (for the remaining ages, not completely recruited, the exploitation pattern should be smaller). For that interval of ages, the means of the values of F_i,j in each year are then calculated. Those means, F_j, are considered as fishing levels in the respective years. The exploitation pattern in each cell, would then be the ratio F_i,j / F_j. The s_i, for the period of years considered, can be taken as the mean of the relative pattern of exploitation calculated before. Alternatively, they can also be taken as referring to s_i of an age chosen for reference.

7.8.2 LENGTH COHORT ANALYSIS - (LCA)

The technique of the cohort ansalysis, applied to the structure of the catches of a cohort during its life, can be made with non constant intervals of time, T_i,. This means that the length classes structure of the catches of a cohort during its life, can also be analysed.

The methods of analysis of the cohort in those cases is called the LCA (Length Cohort Analysis). The same techniques; Pope method, iterative method, etc, of the CA for the ages, can be applied to the LCA analysis (the intervals T_i´s can be calculated from the relative ages).

One way to apply the LCA to the length annual catch compositions, will be: to group the catches of length classes belonging to the same age interval in each year. The technique CA can then be applied directly to the resulting age composition of the catches by age of the matrix [C]. This technique is known as "slicing" the length compositions. To "slice", one usually inverts the von Bertalanffy length growth equation and estimates the age t_i for each length L_i (sometimes using the relative ages t^*_i) (Figure 7.5). It is possible that when grouping the length classes of the respective age interval, there are length classes composed by elements that belong to two consecutive age groups. In these cases, it will be necessary to "break" the catch of these extreme classes into two parts and distribute them to each of those ages. In the example of Figure 7.5, the catches of the length class (24-26] belong to age 0 and to age 1. So, it is necessary to distribute that catch to the two ages. One simple method is to attribute to age 0 the fraction (1.00 - 0.98)/(1.06 - 0.98) = 0.25 of the annual catch of that length class and to age 1 the fraction (1.06 - 1.00)/(1.06 - 0.98) = 0.75. The method may not be the most appropriate one, because it is based on the assumption that, in the length classes, the distribution of the individuals by length is uniform. So, it is necessary to use the smallest possible interval of length classes, when applying this distribution technique.

Another way to do the length cohort analysis is to use the catches in the length classes of the same age group. It is possible to follow the cohorts in the matrix [C], through the length classes belonging to a same age, in a certain year, with the length classes of the next age, in the following year, etc. In this way, the different cohorts existing in the matrix will be separated and the evolution of each one of them will be by length classes, not by age (see Figure 7.5).

Figure 7.5 Example of a matrix [C] with the catches of the cohort shadowed, written in bold, recruited at year 2000, "sliced" by length classes,

The LCA R. Jones method (1961), of analysing a length composition during the life of a cohort can then be applied. The different values of T_i are calculated as T_i = t_i+1*-t_i*, where t_i* e t_i+1* are the relative ages corresponding to the extremes of the length interval i. The vector [N] can also be obtained as the number of initial survivors in each length class of the cohort, and in each age class.

Comments on cohort analyses

1. Certain models, called integrated models, combine all the available information (catches, data collected on research cruises, effort and cpue data, etc) with the matrix [C], and integrate in a unique model, in order to optimize the previously defined criterium function. A model integrating CA and the hypothesis of constant exploitation pattern was developed and called SVPA, separable VPA, because the Fishing level and Exploitation pattern are "separable".

2. Fry (1949) considered the cumulative catches of a cohort by age during its life, from the end to the beginning, as an image of the number of survivors at the beginning of each age (which the author designated as "virtual population"):

In the fishery that Fry studied, M was practically equal to zero.

If M is different from zero it can also be said that the number N_i of survivors at the beginning of the interval i will be

where D_k represents the number of total deaths at the interval k.

Adopting the initial values, E_k(0), for the exploitation rates, E, in all the classes, one can calculate the total deaths:

D_k(0) = C_k/E_k(0).

N _i(0) can be calculated as the cumulative total deaths from the last class up to the i^th class, that is:

Then the expression will be:

Z_i(1).T_i = ln(N_i+1(0)/N_i(0))

and:

F_i(1).T_i = E_i(0).Z_i(1).T_i

Comparing E_i(1) with E_i(0), the new values of E will be:

E_i(1) = F_i(1).T_i/ F_i(1).T_i + M_i.T_i

One can then estimate the values of E with the desired approximation by an iterative method, repeating the five calculations (of D_i, N_i, Z_iT_i, F_iT_i and E_i,) using E_i(1) instead of E_i(0).

In the last class, the number, N_last, can be taken as equal to the number of deaths, D_last, and in this case, N_last will be calculated as:

N_last = D_last = C_last / E_last

3. Finally, the results of CA and of LCA give a perspective view of the stock in the previous years. That information is useful for the short and long-term projections. Usually, data concerning the catches is not available for the year in which the assessment is done and so it is necessary to project the catches and the biomasses to the beginning of the present year before calculating the short-term projection.

4. When the relative ages are calculated, it is usual to adopt zero as the age t_a corresponding to the value of L_a, taken as the lower limit of the first length class represented in the catches.