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 SR stockrecruitment relation.
The least squares method is presented under the forms of Simple linear Regression, multiple linear model and non linear models (method of GaussNewton).
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.
Model
Consider the following variables and parameters:
Response or dependent variable 
= Y 
Auxiliary or independent variable 
= X 
Parameters 
= A, B 
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:
Observed values 
x_{i} and y_{i} for each pair i, where i=1,2,...,i,...n 
Values to be estimated 
A and B_{ }and (Y_{1},Y_{2},...,Y_{i},...,Y_{n})_{ }for the n observed pairs of values 
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:
b = Sxy/Sxx 
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 σ^{2}.
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^{2} + (1  r^{2})
where
r^{2} (coefficient of determination) is the percentage of the total variation that is explained by the model and
1r^{2} is the percentage of the total variation that is not explained by the model.
Model
Consider the following variables and parameters:
Response or dependent variable 
= Y 
Auxiliary or independent variables 
= X_{1}, X_{2},..., X_{j},..., X_{k} 
Parameters 
= B_{1}, B_{2},..., B_{j},..., B_{k} 
The response variable is linear with the parameters
Y = B_{1}X_{1}+B_{2}X_{2}+... + B_{k}X_{k} = Σ B_{j}X_{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_{1},B_{2},...,B_{j},...,B_{k }et (Y_{1},Y_{2},..., Y_{i},..., Y_{n})
The estimated values can be represented by:
(ou b_{1},b_{2},...,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)} = (yY)^{T}.(yY) ou Φ_{(1,1)} = (yX.B)^{T}.(yX.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_{1}, ∂Φ/∂B_{2},..., ∂Φ/∂B_{k}. Thus:
dΦ/dB_{(k,1)} = 2.X^{T}.(yX.B) = 0
or X^{T}y  (X^{T}.X). B = 0
and b = = (X^{T}.X)^{1.} X^{T}y
The results can be written as:
b_{(k,1)} = (X^{T}.X)^{1}.X^{T}y
= 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:

E[c_{1}+c_{2}.u] = c_{1}+c_{2}.E[u] 
V[c_{1}+c_{2}.u] = c_{2}.V[u].c_{2}^{T} 
1  
Random variable, ε 
ε_{n.} (independent) 

Expected value of ε 
E[ε] = 0. 

Variance of ε 
V[ε]_{(n.n) }= E[ε.ε^{T}]=I. σ^{2} 
2  
Observed response variable y 
y = Y+ε 

Expected value of y 
E[y] = Y = X.B. 

Variance of y 
V[y]_{(n.n) }= V[ε]_{(n.n) }= I. σ^{2} 
3  
Estimator of B 
= (X^{T}.X)^{1}.X^{T}.y 

Expected value of 
E[] = B 

Variance of 
V[]_{(k.k) }= (X^{T}.X)^{1}. σ^{2} 
4  
Estimator of Y of the model 
= X. = L.y 

Expected value of 
E[] = Y. 

Variance of 
V[] = L. σ^{2} 
5  
Residual e 
e = y= (IL).y 

Expected value of e 
E[e] = 0 

Variance of e 
V[e] = (IL). σ^{2} 
6  Sum of squares
6.1  Residual Sum of squares = SQ residual_{(1.1)} = (y)^{T}(y) = y^{T }(IL)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 }(LM)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 }(IM) 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^{2} + (1  R^{2})
where:
R^{2} is the percentage of the total variation that is explained by the model. In matrix terms it will be:
R^{2 }= [y^{T}(L  M)y].[ (y^{T}(I  M)y]^{1}
1R^{2} is the percentage of the total variation that is not explained by the model.
The ranks of the matrices (IL), (IM) and (LM) respectively equal to (nk), (n1) and (k1), are the degrees of freedom associated with the respective sums of squares.
Model
Consider the following variables and parameters:
Response or dependent variable 
= Y 
Auxiliary or independent variable 
= X 
Parameters 
= B_{1},B_{2},...,B_{j},...,B_{k} 
The response variable is nonlinear with the parameters
Y = f(X;B) where B is a vector with the components B_{1},B_{2},...,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_{1},B_{2},...,B_{j},..,B_{k }and (Y_{1},Y_{2},...,Y_{i},...,Y_{n})_{ }form the n pairs of observed values.
(Estimates = or b_{1},b_{2},...,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)} = (yY)^{T}.(yY)
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_{0}) +(xx_{0}).f’(x_{0})/1! + (xx_{0})^{2}f’’(x_{0})/2! +... + (x x_{0})^{i} f^{(i)}(x_{0})/i!+...
where
f^{(i)}(x_{0}) = i^{th} derivatives of f(x) in order to x, at the point x_{0}.
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_{0}) + (xx_{0}).f’(x_{0})
The Taylor expansion can be applied to functions with more than one variable. For example, for a function Y = f(x_{1},x_{2}) of two variables, the linear expansion would be:
which may be written, in matrix notation, as
Y = Y_{(0)}+A_{(0)}.(xx_{(0)})
where Y_{(0)} is the value of the function at the point x_{(0)},with components x_{1(0)} and x_{2(0)},and A_{(0)} is the matrix of derivatives whose elements are equal to the partial derivatives of f(x_{1},x_{2}) in order to x_{1},x_{2} at the point (x_{1(0)}, x_{2(0)}).
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_{1}, B_{2},..., B_{k}, would be:
Y = f(x;B) = f(x; B_{(0)}) + (B_{1}B_{1(0)}) f /B_{1 }(x;B_{(0)}) +..... +
(B_{2}B_{2(0)})f /B_{2}^{ }(x;B_{(0)}) +...... +..........+ (B_{k}B_{k(0)}) f /B_{k}^{ }(x;B_{(0)})
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_{(0)} and
ΔB_{(0)} = vector (B  B_{(0)}).
Then, the object function will be:
Φ = (yY)^{T}.(yY) = (yY_{(0) } A_{(0)}. ΔB_{(0)})^{T}(yY_{(0) } A_{(0).} ΔB_{(0)})
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_{(0) }is "equal to zero" then the estimate of B is equal to B_{(0)}.
(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_{(0) }is not "equal to zero" then the vector B_{(0)} will be replaced by:
B_{(1) }=B_{(0) }+ ΔB_{(0)}
And the process will be repeated, that is, there will be another iteration with B_{(0) }replaced by B_{(1) }(and A_{(0) }replaced by A_{(1)}). 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 GaussNewton 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 (NewtonRaphson 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 subroutines 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_{(0)} is close to the real value.
The least squares method (nonlinear 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_{0} for the iterative process of estimation can be obtained by simple linear regression using the following methods:
FordWalford (19331946) and Gulland and Holt (1959) Methods
The FordWalford 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 rewrite 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(tta)}
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:
y = L_{t } 
a = L_{∞} 
b =  (L_{∞} L_{a}). e^{+Kta} 
x = e^{Kt} 


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_{(0)}, of K, and estimating a_{(0)}, b_{(0)} and r^{2}_{(0)} by simple linear regression between y (= L_{t}) and x(=e^{k}_{(0)}). The procedure may be repeated for several values of K, that is, K_{(1)} K_{(2)},.... One can then adopt the regression that results in the larger value of r^{2}, 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 equallyspaced values of K between those two values in regular intervals.
(iii). The corresponding 10 values of r^{2} 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^{2}.
(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.
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 nonexploited 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:
and so 
M = (1/λ).ln p 
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:
for small pelagic fishes 

for demersal fishes 
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:

(Units) 
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
then: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.
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 rearranging 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:

N_{i} ≈ N_{i.} e^{ZTi/2} 
and, as 
N_{i }= N_{a.} e ^{Z.(tita)} 
it will be 
N_{i} ≈ Cte. e^{Ztcentrali} 
and finally: 
ln N_{i }≈ Cte  Z. t_{centrali} 
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:
V_{i }= FN/Z 
and 
lnV_{i }= C^{te }+ lnN_{i} 
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_{0})
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 tta 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
Assumption: Z is constant in the interval of ages, (t_{a}, t_{b}) 

T Constant 





ln C_{i} = Cte  Z · t_{i} 

ln V_{i} = Cte  Z · t_{i} 

T_{i} variable 





ln V_{i} = Cte  Z · t_{i} 

(t_{b} = ∞) (Beverton and Holt equation of Z) 

Supposition: Z is constant in the interval of lengths, (L_{a}, L_{b}) 

Relative age 



(Gulland and Holt equation) 

(Jones and van Zalinge equation) 

(Beverton and Holt equation of Z) 
The least squares method (nonlinear model) can be used to estimate the parameters, α and k, of any of the SR models.
The initial values of the Beverton and Holt model (1957) can be obtained by rewriting 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 nonlinear model.
In the Ricker model (1954) the parameters can be obtained by rewriting 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 nonlinear 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 rewritten 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 nonlinear Deriso model.
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+1}e^{+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

Years 


Ages 
2000 
2001 
2002 
2003 

1 
C 
C 
C 
C 
F_{terminal} 
2 
C 
C 
C 
C 
F_{terminal} 
3 
C 
C 
C 
C 
F_{terminal} 

F_{terminal} 
F_{terminal} 
F_{terminal} 
F_{terminal} 

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 VPASeparable (SVPA).
We have 

and 

and 
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 (2426] 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,
Group 


Years 

Age 
Relative age 
Length 
2000 
2001 
2002 
2003 
0 
1.03 
20 
41 
30 
17 
49 
1.54 
22 
400 
292 
166 
472 

1.98 
24 
952 
699 
400 
1127 

1 
2.06 
26 
1766 
1317 
757 
2108 
2.30 
28 
2222 
1702 
985 
2688 

2.74 
30 
2357 
1872 
1093 
2902 

2.88 
32 
2175 
1091 
1067 
2739 

2 
3.00 
34 
1817 
948 
1416 
1445 
3.42 
36 
1529 
812 
1270 
1250 

3.64 
38 
1251 
684 
980 
1053 

3.83 
40 
1003 
560 
702 
710 

3.96 
42 
787 
290 
310 
558 

3 
4.01 
44 
595 
226 
179 
834 
4.25 
46 
168 
70 
71 
112 







Cohort of the year 2000 
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_{i}T_{i}, F_{i}T_{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 longterm 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 shortterm 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.