## Simple Model

The ADMB framework is first tested on cod (Gadus morhua) data. In the first step, we start with a simple model [1] and an objective function [2], which is the log-likelihood function:

Ni,t = (Ni,t-1-Ci,t-1)e-M [1]
f(Ni,0,qk,s,M) = sumi,t,k{-(2s2)-1[ln(Ii,t,k)-ln(qkNi,t)]2-ln(s)} [2]
Parameter Description Dimension
Ni,t Number of indivduals in cohort i, at age t m-by-(maxage+1)
M Natural mortality Scalar
Ci,t Catch of indivduals in cohort i, at age t m-by-(maxage+1)
Ii,t,k Survey index for cohort i, at age t, by survey k At most m-by-(maxage+1)-by-num_surveys
s Standard deviation for I Scalar
qk Catchability of survey k 1-by-num_surveys

The number of cohorts and the number of surveys is determined by the input data. The free variables of the likelihood function are the m initial cohort sizes Ni,0 as well as qk, s and M.

### Assumptions

The statistical assumption we make is that the survey data is log-normally distributed, that is: ln(Ii,t,k) ~ N(ln(qkNi,t),s2)

### Implementation and Coding

We implement this as an ADMB template file:

cod.tpl
```// File: herring.tpl.
// Nedstrippet versjon til bruk ifm torsk på HI, 2010.
//
// Opprinnelig:
// Date: 29.09.02
// Author: Hans Julius Skaug (skaug@imr.no)
//
// Description: ADMB implementation of Sigurd's herring assessment.
//

DATA_SECTION

// Fixed parameters in the model
init_int 	m			// Number of cohorts

!! cout << "Number of cohorts: " << m << endl;

// Starting and ending year of cohort observation. Cohort is born in the
// starting year.
init_vector n_start_year(1,m)
init_vector n_end_year(1,m)

// Keep track of largest index of each cohort (it is not
// a given that cohorts are followed for the same amount of years)
// note that this index starts at zero
vector n_largest_year_index(1,m)
!! n_largest_year_index = n_end_year - n_start_year;
!! cout << "n_largest_year_index " << n_largest_year_index << endl;

// MATRIX WITH CATCH DATA, one cohort per row, one year per column
// NOTE: The size of C is such that catch in the last year
// (which won't be used for any computation) MUST be present!!
// Setting it to zero is fine
init_matrix C(1,m,0,n_largest_year_index)

// Survey data
init_int 	n_surv			// Number surveys

!! cout << "Number of surveys: " << n_surv << endl;

init_int 	n_I			// Number survey indices
init_matrix I(1,n_I,1,4)		// "cohort","age","survey","I"
// end of survey data

// Number of calendar years our model spans. Computed from start and end years in data.
int n_span
!! n_span = max(n_end_year) - min(n_start_year) + 1;

!! cout << "n_span = max(n_end_year) - min(n_start_year) + 1 = " << n_span << endl;
!! cout << "Ignoring survey data for which there is no catch data!" << endl;

PARAMETER_SECTION

// Recruitement
init_bounded_vector   N0(1,m,0.01,10000);

// Population trajectory, one cohort per row, one year per column.
// (Entry i,j is year j in the life span of cohort i, which is calendar year
// n_start_year(i)+j.) The width corresponds to the largest
// life span we follow, so not all columns will be used unless all spans are the same.
matrix 		N(1,m,0,n_span-1)

// Catchability (survey specific)
init_bounded_vector 	q(1,n_surv,.05,1.0)
// Single q, from meeting with Sam 19Apr2010
// init_bounded_number q0(.05,1.0)
init_bounded_number 	logs(-5,3)		// Log standard deviation of survey errors
number 	s				// Standard deviation of survey errors
init_bounded_number   M(0,.5)                 // Mortality

number tmp

objective_function_value l

PRELIMINARY_CALCS_SECTION

GLOBALS_SECTION

// Helper function to keep population nonnegative
#include <fvar.hpp>
dvariable posfun(_CONST dvariable&x,const double eps)
{
if (x>=eps)
{
return x;
}
else
{
dvariable y=1.0-x/eps;
return eps*(1./(1.+y+y*y));
}
}

// Removed by LF 12032010
// #include <df1b2fun.h>

//#include <df1b2loc.h>

PROCEDURE_SECTION

int i,j;

s= exp(logs);

// initialize likelihood value
l = 0;

// Initialize N
for(i=1;i<=(int) m;i++)
N(i,0) = N0(i);

// Compute population trajectory, note use of helper function posfun
for(i=1;i<=m;i++) {
for(j=1;j<=(int) n_end_year(i)-n_start_year(i);j++) {

N(i,j) = posfun(N(i,j-1)-C(i,j-1),.01)*exp(-M);

}
}

// Likelihood contributions from surveys
for(i=1;i<=(int) n_I;i++)
{
// NEW --- LF, 19 Sept 2009
// Only include survey index if it comes from a year where N actually has a meaningful value
// That is, that it has been recursively computed
if (  I(i,2) <= n_largest_year_index(I(i,1))  ) {
tmp = (log(I(i,4))  - log(q(I(i,3)) *N(I(i,1),I(i,2)))) /s;
// Single q, from meeting with Sam 19Apr2010
// tmp = (log(I(i,4))  - log(q0 *N(I(i,1),I(i,2)))) /s;

l += -.5*tmp*tmp - logs;
}
}

// Reverse objective function sign since ADMB does minimization
l *= -1;

cout << setprecision(4);

REPORT_SECTION
report << setprecision(4);
//cout << "l="<< l << " N0=" << N0 << " q="<< q << " M=" << M  << " s=" << s << endl;

report << N << endl;
TOP_OF_MAIN_SECTION
arrmblsize = 50000000L;

### Preliminary Results

To test the initial model on real data, we use table 3.7 from the 2009 ICES AFWG report, with the catch for 0-2 year olds set to zero. The survey data is taken from table A.13 in the same document. (.dat-file) Population trajectories with 95% confidence intervals for the initial population size, for three cohorts starting in 1985, are given below. On the right are the contents of the file ending with .cor output by ADMB, which contains the estimates, standard deviations, and correlation matrix.

 cod.cor ``` The logarithm of the determinant of the hessian = 43.7225 index name value std dev 1 2 3 4 5 6 1 N0 3.4712e-01 2.8425e-01 1.0000 2 N0 1.5502e-01 9.2802e-02 0.8905 1.0000 3 N0 1.3625e-01 9.4196e-02 0.8854 0.8779 1.0000 4 q 1.1603e-01 1.0183e-01 -0.9477 -0.9396 -0.9343 1.0000 5 logs -9.5234e-02 1.3868e-01 0.0000 0.0000 0.0000 -0.0000 1.0000 6 M 3.3178e-07 2.3460e-04 0.0047 0.0053 0.0049 -0.0034 0.0004 1.0000```

## Model with two mortality parameters

We now extend the model by introducing a latent variable, Mvec, which adjusts the mortality for age classes 7 and up.

Ni,t = (Ni,t-1-Ci,t-1)e-M for t<7 [1]
Ni,t = (Ni,t-1-Ci,t-1)e-(M+Mvec) for t >= 7
f(Ni,0,qk,s,M,Mvec) = sumi,t,k {-(2s2)-1[ln(Ii,t,k)-ln(qkNi,t)]2-ln(s)} -0.5 Mvec2[2]

NOTE: It is not recommended to optimize over both the regular and the latent variable(s). Instead we remove the latter by integration, specifically by a Laplace approximation to the integral. Define:

 g(Ni,0,qk,s,M) Joint function [2] with Mvec as a constant, and Ni,0,qk,s,M as variables h(Mvec) Joint function [2] with Ni,0,qk,s,M as constants, and Mvec as the variable Mvec The maximum of h(Mvec) (for a given set of constants). H(Mvec) The Hessian of h(Mvec), evaluated at Mvec.

Then the Laplace approximation (our objective function) is:

g(Ni,0,qk,s,M) + 0.5 log( det( H(Mvec) ) ) [3]

In other words, for each function evaluation with the variables Ni,0,qk,s,M, as the input, an inner optimization problem must be solved to determine Mvec. This is done automatically by ADMB.

### Assumptions

As for the simple model, the statistical assumption we make is that the survey data is log-normally distributed, that is: ln(Ii,t,k) ~ N(ln(qkNi,t),s2)

### Implementation

Note that the N0 variables are scaled with a factor 1000 in this version.

cod.tpl
```DATA_SECTION
!! cout << "***************************************************************" << endl;
!! cout << "* VERSION WITH SCALING, INPUT N0 WILL BE MULTIPLIED WITH 1000 *" << endl;
!! cout << "***************************************************************" << endl;
!! cout << endl;

// Fixed parameters in the model
init_int      m                       // Number of cohorts

!! cout << "Number of cohorts: " << m << endl;

// Starting and ending year of cohort observation. Cohort is born in the
// starting year.
init_vector n_start_year(1,m)
init_vector n_end_year(1,m)

// Keep track of largest index of each cohort (it is not
// a given that cohorts are followed for the same amount of years)
// note that this index starts at zero
vector n_largest_year_index(1,m)
!! n_largest_year_index = n_end_year - n_start_year;
!! cout << "n_largest_year_index " << n_largest_year_index << endl;

// MATRIX WITH CATCH DATA, one cohort per row, one year per column
// NOTE: The size of C is such that catch in the last year
// (which won't be used for any computation) MUST be present!!
// Setting it to zero is fine
init_matrix C(1,m,0,n_largest_year_index)

// Survey data
init_int      n_surv                  // Number surveys

!! cout << "Number of surveys: " << n_surv << endl;

init_int      n_I                     // Number survey indices
init_matrix I(1,n_I,1,4)              // "cohort","age","survey","I"
// end of survey data

// Number of calendar years our model spans. Computed from start and end years in data.
int n_span
!! n_span = max(n_end_year) - min(n_start_year) + 1;

!! cout << "n_span = max(n_end_year) - min(n_start_year) + 1 = " << n_span << endl;
!! cout << "Ignoring survey data for which there is no catch data!" << endl;

PARAMETER_SECTION

// Recruitement
init_bounded_vector   N0(1,m,0.01,10000);

// Population trajectory, one cohort per row, one year per column.
// (Entry i,j is year j in the life span of cohort i, which is calendar year
// n_start_year(i)+j.) The width corresponds to the largest
// life span we follow, so not all columns will be used unless all spans are the same.
matrix                N(1,m,0,n_span-1)

// Catchability (survey specific)
init_bounded_vector   q(1,n_surv,.05,1.0)
// Single a, from meeting with Sam 19Apr2010
//init_bounded_number q0(.05,1.0)
init_bounded_number   logs(-5,3)              // Log standard deviation of survey errors
number        s                               // Standard deviation of survey errors
init_bounded_number   M(0,.5)                 // Mortality

random_effects_vector Mvec(1,1)

number tmp

objective_function_value l

PRELIMINARY_CALCS_SECTION

GLOBALS_SECTION

// Helper function to keep population nonnegative
#include <fvar.hpp>
dvariable posfun(_CONST dvariable&x,const double eps)
{
if (x>=eps)
{
return x;
}
else
{
dvariable y=1.0-x/eps;
return eps*(1./(1.+y+y*y));
}
}

// Removed by LF 12032010
// #include <df1b2fun.h>

//#include <df1b2loc.h>

PROCEDURE_SECTION

int i,j;

s= exp(logs);

// initialize likelihood value
l = 0;

// Penalize random effects
l -= 0.5 * norm2(Mvec);

// Initialize N
for(i=1;i<=(int) m;i++)
N(i,0) = 1000*N0(i);

// Compute population trajectory, note use of helper function posfun
for(i=1;i<=m;i++) {
for(j=1;j<=(int) n_end_year(i)-n_start_year(i);j++) {

// NEW - mortality has two components, base mortality for age 0-6, and
// a random effects correction term for ages 7 and up
if (j>6) {
N(i,j) = posfun(N(i,j-1)-C(i,j-1),.01)*exp(-M-Mvec(1));
} else {
N(i,j) = posfun(N(i,j-1)-C(i,j-1),.01)*exp(-M);
}
}
}

// Likelihood contributions from surveys
for(i=1;i<=(int) n_I;i++)
{
// NEW --- LF, 19 Sept 2009
// Only include survey index if it comes from a year where N actually has a meaningful value
// That is, that it has been recursively computed
if (  I(i,2) <= n_largest_year_index(I(i,1))  ) {
tmp = (log(I(i,4))  - log(q(I(i,3)) *N(I(i,1),I(i,2)))) /s;
// Single q, from meeting with Sam 19apr2010
//  tmp = (log(I(i,4))  - log(q0 *N(I(i,1),I(i,2)))) /s;
l += -.5*tmp*tmp - logs;
}
}

// Reverse objective function sign since ADMB does minimization
l *= -1;

cout << setprecision(4);

REPORT_SECTION
report << setprecision(4);
//cout << "l="<< l << " N0=" << N0 << " q="<< q << " M=" << M  << " s=" << s << endl;

report << N << endl;
TOP_OF_MAIN_SECTION
arrmblsize = 50000000L;

The corresponding .pin file must have an initial value for the new parameter and may (if there are three cohorts and one survey) look like:

cod.pin
```# N0
1 1 1
# q
0.40
# log(s)
0.1
# M
0.2
# Mvec
0```

### Preliminary Results

NOTE: For all the experiments in this section, the estimate of q is 0.05, its lower bound.

If we run the same experiment as for the simple model (.dat-file), we get the trajectories, 95% confidence intervals and output data below:

 cod.cor ``` The logarithm of the determinant of the hessian = 73.3174 index name value std dev 1 2 3 4 5 6 1 N0 9.4596e-01 2.8898e-01 1.0000 2 N0 3.6858e-01 1.0512e-01 0.1834 1.0000 3 N0 3.1155e-01 9.1028e-02 0.1012 0.1101 1.0000 4 q 5.0000e-02 8.1311e-09 -0.0000 -0.0000 -0.0000 1.0000 5 logs -1.2616e-01 1.4133e-01 -0.0065 -0.0053 -0.0039 0.0000 1.0000 6 M 1.2339e-08 3.8799e-05 0.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 7 Mvec 3.1205e-01 1.7955e-01 -0.4104 -0.4444 -0.2465 0.0000 0.0152 0.0000 1.0000```

The red curve in particular illustrates the effect of two different mortalities, since the slope changes significantly after seven years.

We now perform the same experiment, but instead of following three cohorts we follow 1, 2, 4 and 5 cohorts, to see how the output behaves.

For the case of only one cohort (.dat-file), the trajectory is significantly different from the corresponding trajectory when estimating three cohorts, as can be seen in the figure below:

 cod.cor ``` The logarithm of the determinant of the hessian = 39.9846 index name value std dev 1 2 3 4 1 N0 3.1453e+00 1.6824e+00 1.0000 2 q 5.0000e-02 9.1511e-08 -0.0000 1.0000 3 logs -5.1246e-01 2.4933e-01 -0.0129 0.0000 1.0000 4 M 2.9277e-01 1.2740e-01 0.9053 -0.0000 -0.0182 1.0000 5 Mvec -1.2934e-01 3.1406e-01 0.5776 0.0000 -0.0224 0.8012 1.0000  ```

For 2, 4, and 5 cohorts, the trajectories do not change much, and are consistent with the case of three cohorts.

2 cohorts: (.dat-file)

 cod.cor ``` The logarithm of the determinant of the hessian = 62.9246 index name value std dev 1 2 3 4 5 1 N0 9.4914e-01 2.8341e-01 1.0000 2 N0 3.6989e-01 1.0341e-01 0.2107 1.0000 3 q 5.0000e-02 4.7730e-08 -0.0000 -0.0000 1.0000 4 logs -1.6446e-01 1.7132e-01 -0.0091 -0.0075 0.0000 1.0000 5 M 1.2348e-08 9.2555e-05 0.0000 0.0000 -0.0000 0.0000 1.0000 6 Mvec 3.1706e-01 1.8830e-01 -0.4403 -0.4756 0.0000 0.0196 0.0000 1.0000  ```

4 cohorts: (.dat-file)

 cod.cor ``` The logarithm of the determinant of the hessian = 75.5664 index name value std dev 1 2 3 4 5 6 7 1 N0 9.4633e-01 2.8648e-01 1.0000 2 N0 3.6870e-01 1.0416e-01 0.1784 1.0000 3 N0 3.1162e-01 9.0387e-02 0.0982 0.1069 1.0000 4 N0 8.7346e-01 2.8254e-01 0.0325 0.0353 0.0195 1.0000 5 q 5.0000e-02 1.4173e-08 -0.0000 -0.0000 -0.0000 -0.0000 1.0000 6 logs -1.3245e-01 1.2494e-01 -0.0055 -0.0045 -0.0033 -0.0013 0.0000 1.0000 7 M 1.2341e-08 2.9666e-05 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 8 Mvec 3.1254e-01 1.7543e-01 -0.4046 -0.4384 -0.2426 -0.0803 0.0000 0.0132 0.0000 1.0000  ```

5 cohorts: (.dat-file)

 cod.cor ``` The logarithm of the determinant of the hessian = 77.7625 index name value std dev 1 2 3 4 5 6 7 8 1 N0 9.4714e-01 2.6779e-01 1.0000 2 N0 3.6894e-01 9.7359e-02 0.1788 1.0000 3 N0 3.1177e-01 8.4444e-02 0.0985 0.1072 1.0000 4 N0 8.7364e-01 2.6385e-01 0.0326 0.0354 0.0195 1.0000 5 N0 3.1106e+00 1.0306e+00 0.0000 0.0000 0.0000 0.0000 1.0000 6 q 5.0000e-02 7.1536e-09 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 1.0000 7 logs -2.0112e-01 1.1466e-01 -0.0048 -0.0039 -0.0029 -0.0011 -0.0000 0.0000 1.0000 8 M 1.2339e-08 2.7192e-05 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 9 Mvec 3.1373e-01 1.6409e-01 -0.4052 -0.4391 -0.2430 -0.0805 -0.0000 0.0000 0.0113 0.0000 1.0000```

Note that the scale of the y-axis is different for the case of five cohorts.

## Model with sigmoid function for catchability

Since the estimates for the catchability q are very low, we try to model it with a sigmoid function. That is,

Catchability = q0 exp(a x + b)/(1+exp(a x + b))

where the x represents the age of the cohort.

NOTE this means that instead of one q per survey, there are now three parameters, q, a, and b (called q0, q_alpha and q_beta in the code), regardless of the number of surveys.

The model is implemented in this .tpl file. Note that the parameter structure is now different, so a .pin file for three cohorts should look like:

cod.pin
```# N0
1.000000  1.000000  1.000000

# q  alpha beta
0.60 1 1

# log(s)
0.1

# M
0.2```

We run the model on the same dataset as before, and get:

 cod.cor ``` The logarithm of the determinant of the hessian = 25.2341 index name value std dev 1 2 3 4 5 6 7 8 1 N0 3.0186e+01 1.1253e+01 1.0000 2 N0 1.1064e+01 3.9272e+00 0.5094 1.0000 3 N0 9.4634e+00 3.7190e+00 0.5202 0.5138 1.0000 4 q0 5.0006e-02 5.8935e-03 -0.2905 -0.2920 -0.2677 1.0000 5 q_alpha 9.3271e-01 2.0696e-01 -0.3959 -0.3923 -0.4035 -0.0115 1.0000 6 q_beta -4.5843e+00 5.2739e-01 -0.1483 -0.1427 -0.1564 -0.0168 -0.6298 1.0000 7 logs -2.3144e-01 1.3868e-01 -0.0001 -0.0001 -0.0001 0.0003 -0.0000 -0.0000 1.0000 8 M 5.0000e-01 8.6228e-04 0.0185 0.0192 0.0171 -0.0000 0.0012 -0.0115 -0.0004 1.0000```

As one can see, the effective catchability is very low, due to the fact that q0 is at its lower bound of 0.05. Therefore we set the lower bound to be 0.5, (which requires editing of the .tpl file and recompilation), and then get:

 cod.cor ``` The logarithm of the determinant of the hessian = 39.2739 index name value std dev 1 2 3 4 5 6 7 8 1 N0 5.4677e+00 1.6838e+00 1.0000 2 N0 2.5934e+00 4.7455e-01 0.3641 1.0000 3 N0 2.0265e+00 5.4008e-01 0.4213 0.3511 1.0000 4 q0 5.0000e-01 2.0633e-03 -0.0064 -0.0065 -0.0054 1.0000 5 q_alpha 8.1589e-01 1.4679e-01 -0.4167 -0.3754 -0.4119 -0.0023 1.0000 6 q_beta -5.1365e+00 4.4860e-01 -0.0600 -0.0124 -0.0439 -0.0045 -0.7160 1.0000 7 logs -2.2158e-01 1.3868e-01 -0.0000 -0.0000 -0.0000 0.0004 -0.0000 0.0000 1.0000 8 M 5.0000e-01 1.0834e-03 0.0268 0.0424 0.0288 -0.0000 0.0050 -0.0165 -0.0004 1.0000```

As one can see, q is still at its lower bound, which seems to be the most important parameter for this model on this dataset.

## Model with two separate catchabilities

This model has two q variables for each survey, one which is applied to age groups 0-2, and one which is applied to age group 3 and up. It is implemented in this .tpl file. A .pin file for three cohorts may read:

cod.pin
```# N0
3.07 2.95 2.9
# q02
0.20
# q3plus
0.2
# log(s)
0.1
# M
0.1```

We run it on the usual data set, and get:

 cod.cor ``` The logarithm of the determinant of the hessian = 50.7902 index name value std dev 1 2 3 4 5 6 7 1 N0 2.0581e-01 8.4094e-02 1.0000 2 N0 1.1073e-01 2.3668e-02 0.7116 1.0000 3 N0 9.2034e-02 2.6383e-02 0.7050 0.6679 1.0000 4 q02 1.1075e-01 4.9119e-02 -0.5731 -0.5228 -0.5390 1.0000 5 q3plus 2.5265e-01 1.3221e-01 -0.8607 -0.8178 -0.8072 0.5695 1.0000 6 logs -1.5297e-01 1.3868e-01 0.0000 0.0000 0.0000 -0.0000 -0.0000 1.0000 7 M 2.0589e-06 1.4559e-03 0.0336 0.0514 0.0403 -0.0223 -0.0140 0.0004 1.0000```

As one can see the catchability for the classes age 3 and up is larger than for 0-2 year olds, which is consistent with our expectations.

## Model with two separate catchabilities and two mortalities

This model is a combination of the models with two Ms and two qs. The .tpl file can be found here. A .pin file should look like

cod.pin
```# N0
3.07 2.95 2.9
# q02
0.20
# q3plus
0.2
# log(s)
0.1
# M
0.1
# Mvec
0```

We run it on the standard dataset and get:

 cod.cor ``` The logarithm of the determinant of the hessian = 82.4993 index name value std dev 1 2 3 4 5 6 7 1 N0 5.9029e-01 2.4390e-01 1.0000 2 N0 2.4613e-01 8.6909e-02 0.6200 1.0000 3 N0 2.0132e-01 7.4850e-02 0.5779 0.5697 1.0000 4 q02 5.0000e-02 7.1973e-09 -0.0000 -0.0000 -0.0000 1.0000 5 q3plus 9.9498e-02 4.1441e-02 -0.7242 -0.6946 -0.7115 0.0000 1.0000 6 logs -1.7842e-01 1.4154e-01 0.0122 0.0161 0.0102 0.0000 -0.0092 1.0000 7 M 1.2339e-08 5.4534e-05 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 1.0000 8 Mvec 3.7216e-01 1.7757e-01 -0.2653 -0.3269 -0.1560 0.0000 -0.0631 -0.0195 0.0000 1.0000```

Note that the catchability for 0-2 year olds, q02, is at its lower bound of 0.05, whereas the catchability q3plus is not.