R help - Sep 2009 - Maximum likelihood estimation of parameters make no biological sense

R-help,

I'm trying to estimate some parameters using the Maximum Likehood method.
The model describes fish growth using a sigmoidal-type of curve:

fn_w <- function(params) {
        Winf   <-  params[1]
        k       <-  params[2]
        t0      <-  params[3]
        b       <-  params[4]
       sigma <-  params[5]
        what  <-  Winf * (1-exp(- k *(tt - t0)))^b
        logL   <-  -sum(dnorm(log(wobs),log(what),sqrt(sigma),TRUE))
        return(logL)
                    }

tt      <- 4:14
wobs <- c(1.545, 1.920, 2.321 ,2.591, 3.676, 4.425 ,5.028, 5.877, 6.990,
6.800 ,6.900)

An then the optimization method:

OPT <-optim(c(8, .1, 0, 3, 1), fn_w, method="L-BFGS-B"
,lower=c(0.0, 0.001, 0.001,0.001, 0.01), upper = rep(Inf, 5), hessian=TRUE,
control=list(trace=1))

which gives:

$par     Winf               k                     t0                     b      
sigma
[1] 24.27206813  0.04679844  0.00100000  1.61760492  0.01000000

$value
[1] -11.69524

$counts
function gradient 
     143      143 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

$hessian
              [,1]          [,2]         [,3]          [,4]          [,5]
[1,]  1.867150e+00  1.262763e+03    -7.857719 -5.153276e+01 -1.492850e-05
[2,]  1.262763e+03  8.608461e+05 -5512.469266 -3.562137e+04  9.693180e-05
[3,] -7.857719e+00 -5.512469e+03    41.670222  2.473167e+02 -5.356813e+01
[4,] -5.153276e+01 -3.562137e+04   247.316675  1.535086e+03 -1.464370e-03
[5,] -1.492850e-05  9.693180e-05   -53.568127 -1.464370e-03  1.730462e+04

after iteration number 80.
>From the biological point of view Winf =24(hipothesized asimptotical maximum
weight)
makes no sense while the b parameter is no nearly close to b=3 leading to a
non-sigmoidal
curve.

However using the least-squares method provide with more sensible parameter
estimates

$par     Winf              k                  t0                b
[1] 10.3827256  0.0344187  3.1751376  2.9657368

$value
[1] 2.164403

$counts
function gradient 
      48       48 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Is there anything wrong with my MLE function and parameters?
I have tried with distinct initial parameters also.

Can anyone give me a clue on this?

Thanks in advance.

Hi,

Your results are do to using an unstable parameterization
of the Von Bertalanffy growth curve, combined with the unreliable
optimization methods supplied with R.  I coded up your model in
AD Model Builder which supplies exact derivatives through
AD.

I used your starting values and ran the model with no optimization steps
just to se that we had the same value for the -log-likelihood
Results are
# Number of parameters = 5  Objective function value = -11.6954  Maximum
gradien
t component = 0.00000
# winf:
24.2720681300
# k:
0.0467984400000
# t0:
0.00100000000000
# vhat:
0.0100000000000
# b:
1.61760492000

However the R routine is stuck. When I let the ADMB code run it produced

# Number of parameters = 5  Objective function value = -13.8515  Maximum
gradient component = 9.41643e-05
# winf:
15.7188821203
# k:
0.118198731245
# t0:
-32.9089295327
# vhat:
0.00471832483493
# b:
184.999879271

Note that b--> infinity. I have it bounded at 185.
t0--> -infinity  so that the model is only using a small part of the
growth curve which happens to fit the data better.

The estimated correlation matrix for the parameter estimates tells the story

index name  value    std dev       1       2       3       4       5
 1   winf 1.5719e+01 5.1252e+00   1.0000
 2   k    1.1820e-01 2.7849e-02  -0.9832  1.0000
 3   t0  -3.2909e+01 7.6867e+00  -0.9748  0.9990  1.0000
 4   vhat 4.7183e-03 2.0119e-03   0.0000  0.0000  0.0000  1.0000
 5   b    1.8500e+02 1.6374e+00  -0.0002  0.0003 -0.0094  0.0000  1.0000

You can see that several of the parameters are highly confounded.

Also the eigenvalues of the Hessian are

  0.01691149331  0.02045399106 963.2994413  2255.900979  4225373.963

So you have a condition number of about 10^8. Very difficult to work
with such a function with only approximate derivatives.

I think the moral of the story is that you should use a more stable
parameterization or an industrial strength estimation system or maybe
both.

  Cheers,

   Dave

   Cheers,

    Dave







-- 
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com