[R] Maximum likelihood estimation of parameters make no biological sense

Luis Ridao Cruz luisr at hav.fo
Thu Sep 24 12:51:22 CEST 2009


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.




More information about the R-help mailing list