AdmArwiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== ADMAR-ADMB Framework ======
===== 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:\\
^N<sub>i,t</sub>      ^= ^(N<sub>i,t-1</sub>-C<sub>i,t-1</sub>)e<sup>-M</sup>  ^[1] ^
^f(N<sub>i,0</sub>,q<sub>k</sub>,s,M) ^= ^ sum<sub>i,t,k</sub>{-(2s<sup>2</sup>)<sup>-1</sup>[ln(I<sub>i,t,k</sub>)-ln(q<sub>k</sub>N<sub>i,t</sub>)]<sup>2</sup>-ln(s)} ^[2] ^

^Parameter       ^ Description        ^ Dimension         ^
^ N<sub>i,t</sub>    | Number of indivduals in cohort i, at age t  | m-by-(maxage+1)                 |
^ M    | Natural mortality          | Scalar               |
^ C<sub>i,t</sub>    | Catch of indivduals in cohort i, at age t           | m-by-(maxage+1)           |
^ I<sub>i,t,k</sub>    | 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                 |
^ q<sub>k</sub>    | 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 N<sub>i,0</sub> as well as q<sub>k</sub>, s and M.

==== Assumptions ====

The statistical assumption we make is that the survey data is log-normally distributed, that is:
ln(I<sub>i,t,k</sub>) %%~%% N(ln(q<sub>k</sub>N<sub>i,t</sub>),s<sup>2</sup>)

==== Implementation and Coding ====

We implement this as an ADMB template file:

<file C++ 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);
</file>

==== 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. ({{:admb_real:cod_1985_3cohorts_dat.txt|.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.

| {{:admb_real:simple_model85_87.png?200|}} | <file txt 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
</file> |

===== 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.

^N<sub>i,t</sub>      ^= ^(N<sub>i,t-1</sub>-C<sub>i,t-1</sub>)e<sup>-M</sup>  for t<7 ^[1] ^
^N<sub>i,t</sub>      ^= ^(N<sub>i,t-1</sub>-C<sub>i,t-1</sub>)e<sup>-(M+Mvec)</sup> for t >= 7 ^::: ^
^f(N<sub>i,0</sub>,q<sub>k</sub>,s,M,Mvec) ^= ^ sum<sub>i,t,k</sub> {-(2s<sup>2</sup>)<sup>-1</sup>[ln(I<sub>i,t,k</sub>)-ln(q<sub>k</sub>N<sub>i,t</sub>)]<sup>2</sup>-ln(s)} -0.5 Mvec<sup>2</sup>^[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(N<sub>i,0</sub>,q<sub>k</sub>,s,M) | Joint function [2] with Mvec as a constant, and N<sub>i,0</sub>,q<sub>k</sub>,s,M as variables |
| h(Mvec) | Joint function [2] with N<sub>i,0</sub>,q<sub>k</sub>,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(N<sub>i,0</sub>,q<sub>k</sub>,s,M) + 0.5 log( det( H(__Mvec__) ) ) ^ [3] ^

In other words, for each function evaluation with the variables N<sub>i,0</sub>,q<sub>k</sub>,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(I<sub>i,t,k</sub>) %%~%% N(ln(q<sub>k</sub>N<sub>i,t</sub>),s<sup>2</sup>)

==== Implementation ====

The .tpl file then reads:

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

<code txt 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);
</code>
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:
<code txt cod.pin>
# N0
1 1 1
# q
0.40 
# log(s)
0.1
# M
0.2
# Mvec
0
</code>

==== 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 ({{:admb_real:cod_1985_3cohorts_dat.txt|.dat-file}}), we get the trajectories, 95% confidence intervals and output data below:

| {{:admb_real:twom_model85_87.png?200|}} | <file txt 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
</file> |

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 ({{:admb_real:cod_1985_1cohorts_dat.txt|.dat-file}}), the trajectory is significantly different from the corresponding trajectory when estimating three cohorts, as can be seen in the figure below:

| {{:admb_real:twom_model85.png?200|}} |<file txt 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
 </file>| 

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

2 cohorts: ({{:admb_real:cod_1985_2cohorts_dat.txt|.dat-file}})

| {{:admb_real:twom_model85_86.png?200|}}| <file txt 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
 </file>| 

4 cohorts: ({{:admb_real:cod_1985_4cohorts_dat.txt|.dat-file}})

| {{:admb_real:twom_model85_88.png?200|}}| <file txt 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
 </file>| 

5 cohorts: ({{:admb_real:cod_1985_5cohorts_dat.txt|.dat-file}})
| {{:admb_real:twom_model85_89.png?200|}}| <file txt 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
</file>|

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 {{:admb_real:sigmoid_model_rename_to_cod_tpl.txt|this .tpl file}}. Note that the parameter structure is now different, so a .pin file for three cohorts should look like:
<file txt cod.pin>
# N0
1.000000  1.000000  1.000000  

# q  alpha beta
0.60 1 1

# log(s)
0.1

# M
0.2
</file> 

We run the model on the same dataset as before, and get:
| {{:admb_real:sigmoid_model85_87.png?200|}} | <file txt 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

</file> |
| {{:admb_real:sigmoid_q85_87.png?200|}}| ::: |

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:

|{{:admb_real:sigmoid_model_q05_85_87.png?200|}} |<file txt 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
</file> |
|{{:admb_real:sigmoid_q_q05_85_87.png?200|}} | ::: |

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 {{:admb_real:cod_two_q_rename_to_cod_tpl.txt|this .tpl file.}} A .pin file for three cohorts may read:
<file txt cod.pin>
# N0
3.07 2.95 2.9
# q02
0.20 
# q3plus
0.2
# log(s)
0.1
# M
0.1
</file>

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

|{{:admb_real:two_q_model_85_87.png?200|}}| <file txt 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
</file> |

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 {{:admb_real:cod_two_q_two_m_rename_to_cod_tpl.txt|here.}} A .pin file should look like
<file txt cod.pin>
# N0
3.07 2.95 2.9
# q02
0.20 
# q3plus
0.2
# log(s)
0.1
# M
0.1
# Mvec
0
</file>

We run it on the standard dataset and get:

| {{:admb_real:two_q_two_m_model_85_87.png?200|}} | <file txt 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
</file>|

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