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:

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 N_i,0 as well as q_k, s and M.

The statistical assumption we make is that the survey data is log-normally distributed, that is: ln(I_i,t,k) ~ N(ln(q_kN_i,t),s²)

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;
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(1000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2000000);
  gradient_structure::set_MAX_NVAR_OFFSET(40205);

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

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

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:

Then the Laplace approximation (our objective function) is:

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

As for the simple model, the statistical assumption we make is that the survey data is log-normally distributed, that is: ln(I_i,t,k) ~ N(ln(q_kN_i,t),s²)

The .tpl file then reads:

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;
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(1000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2000000);
  gradient_structure::set_MAX_NVAR_OFFSET(40205);

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

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.

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

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.

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.

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.