[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