Dear R users,
I'im trying to find the parameters of a dynamic biomass model using maximum
likelihood estimation. I used two approaches, one by hand, with optim()
function and the other using mle2() function from package bbmle. My problem
is that the results change a lot depending on the initial values... I can't
see what I am doing wrong...
Thank you for any help!
# Data
x <- 1995:2010
B <- c(3500, 3200, 3000, 2800, 2600, 3000, 3200, 3800, 4200, 4300, 4400,
4400, 4500, 4600, 5000, 4300)
Ct <- c(912, 767, 642, 482, 353, 331, 332, 309, 366, 402, 392, 478, 408,
434, 407, 637)
a <- c(0.539, 0.603, -0.948, 0.166, 1.895, 0.786, 0.901, 0.844, 0.337,
0.429, 0.304, 0.230, 1.001, 0.750, 0.507, 1.502)
Ag <- 0.55
# Function with quantity to minimize
modl <- function(par) {
ro <- par[1]
ko <- par[2]
n <- length(B)
Be <- rep(NA, n)
Be[1] <- ko * Ag
for ( k in 2:n)
Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1]
err <- (log(B) - log(Be))^2
ee <- sqrt( sum(err)/(n-2) )
LL <- (1/(sqrt(2*pi)*ee)) * exp( -(err/(2*ee^2) ) )
-crossprod(LL)
}
# Using function optim()
par.optim <- optim(par = list(ro=0.4, ko=8000), modl, method =
"BFGS")
ro <- par.optim$par[1]
ko <- par.optim$par[2]
# estimated values of "B"
n <- length(B)
Be <- rep(NA, n)
Be[1] <- ko * Ag
for ( k in 2:n)
Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1]
# Plot, estimation of "B" seems reasonable....
plot(x, B, ylim=c(1000, 7000))
lines(x, Be, col="blue", lwd=2)
# ... but it is very sensible to initial values...
par.optim2 <- optim(par = list(ro=0.4, ko=10000), modl, method =
"BFGS")
ro2 <- par.optim2$par[1]
ko2 <- par.optim2$par[2]
Be2 <- rep(NA, n)
Be2[1] <- ko2 * Ag
for ( k in 2:n)
Be2[k] <- Be2[k-1] + ro2 * a[k-1] * Be2[k-1] * (1 - Be2[k-1]/ko2) -
Ct[k-1]
lines(x, Be2, col="blue", lwd=2, lty=3)
# Uing mle2 function
library(bbmle)
LL <- function(ro, ko, mu, sigma) {
n <- length(B)
Be <- rep(NA, n)
Be[1] <- ko * Ag
for ( k in 2:n)
Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1]
err <- log(B) - log(Be)
R <- (dnorm(err, mu, sigma, log=TRUE))
-sum(R)
}
Bc.mle <- mle2(LL, start = list(ro=0.4, ko=8000, mu=0, sigma=1))
summary(Bc.mle)
ro3 <- coef(Bc.mle)[1]
ko3 <- coef(Bc.mle)[2]
Be3 <- rep(NA, n)
Be3[1] <- ko3 * Ag
for ( k in 2:n)
Be3[k] <- Be3[k-1] + ro3 * a[k-1] * Be3[k-1] * (1 - Be3[k-1]/ko3) -
Ct[k-1]
lines(x, Be3, col="red", lwd=2)
--
H?ctor Villalobos <Hector.Villalobos.O at gmail.com>
Depto. de Pesquer?as y Biolog?a Marina
Centro Interdisciplinario de Ciencias Marinas-Instituto Polit?cnico Nacional
CICIMAR-I.P.N.
A.P. 592. Colonia Centro
La Paz, Baja California Sur, M?XICO. 23000
Tels.: (+52 612) 122 53 44; 123 46 58; 123 47 34 ext.: 81602
Fax: (+52 612) 122 53 22
[[alternative HTML version deleted]]
I made no attempt to examine your details for problems, but in general, My problem> is that the results change a lot depending on the initial values... I can't > see what I am doing wrong... > > This is a symptom of an overparameterized model: The parameter estimates > are unstable even though the predictions may not change much. In other > words, your model may be too complex for your data.Whether that is true here, you or others will have to determine. Try simplifying your model as a start. -- Bert> > > # Data > x <- 1995:2010 > B <- c(3500, 3200, 3000, 2800, 2600, 3000, 3200, 3800, 4200, 4300, 4400, > 4400, 4500, 4600, 5000, 4300) > Ct <- c(912, 767, 642, 482, 353, 331, 332, 309, 366, 402, 392, 478, 408, > 434, 407, 637) > a <- c(0.539, 0.603, -0.948, 0.166, 1.895, 0.786, 0.901, 0.844, 0.337, > 0.429, 0.304, 0.230, 1.001, 0.750, 0.507, 1.502) > Ag <- 0.55 > > # Function with quantity to minimize > modl <- function(par) { > ro <- par[1] > ko <- par[2] > n <- length(B) > Be <- rep(NA, n) > Be[1] <- ko * Ag > for ( k in 2:n) > Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1] > err <- (log(B) - log(Be))^2 > ee <- sqrt( sum(err)/(n-2) ) > LL <- (1/(sqrt(2*pi)*ee)) * exp( -(err/(2*ee^2) ) ) > -crossprod(LL) > } > > # Using function optim() > par.optim <- optim(par = list(ro=0.4, ko=8000), modl, method = "BFGS") > ro <- par.optim$par[1] > ko <- par.optim$par[2] > > # estimated values of "B" > n <- length(B) > Be <- rep(NA, n) > Be[1] <- ko * Ag > for ( k in 2:n) > Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1] > > # Plot, estimation of "B" seems reasonable.... > plot(x, B, ylim=c(1000, 7000)) > lines(x, Be, col="blue", lwd=2) > > > # ... but it is very sensible to initial values... > par.optim2 <- optim(par = list(ro=0.4, ko=10000), modl, method = "BFGS") > ro2 <- par.optim2$par[1] > ko2 <- par.optim2$par[2] > > Be2 <- rep(NA, n) > Be2[1] <- ko2 * Ag > for ( k in 2:n) > Be2[k] <- Be2[k-1] + ro2 * a[k-1] * Be2[k-1] * (1 - Be2[k-1]/ko2) - > Ct[k-1] > > lines(x, Be2, col="blue", lwd=2, lty=3) > > > > # Uing mle2 function > library(bbmle) > LL <- function(ro, ko, mu, sigma) { > n <- length(B) > Be <- rep(NA, n) > Be[1] <- ko * Ag > for ( k in 2:n) > Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1] > err <- log(B) - log(Be) > R <- (dnorm(err, mu, sigma, log=TRUE)) > -sum(R) > } > > Bc.mle <- mle2(LL, start = list(ro=0.4, ko=8000, mu=0, sigma=1)) > summary(Bc.mle) > > ro3 <- coef(Bc.mle)[1] > ko3 <- coef(Bc.mle)[2] > > Be3 <- rep(NA, n) > Be3[1] <- ko3 * Ag > for ( k in 2:n) > Be3[k] <- Be3[k-1] + ro3 * a[k-1] * Be3[k-1] * (1 - Be3[k-1]/ko3) - > Ct[k-1] > > lines(x, Be3, col="red", lwd=2) > > > -- > > H?ctor Villalobos <Hector.Villalobos.O at gmail.com <javascript:;>> > > Depto. de Pesquer?as y Biolog?a Marina > > Centro Interdisciplinario de Ciencias Marinas-Instituto Polit?cnico > Nacional > > CICIMAR-I.P.N. > > A.P. 592. Colonia Centro > > La Paz, Baja California Sur, M?XICO. 23000 > > Tels.: (+52 612) 122 53 44; 123 46 58; 123 47 34 ext.: 81602 > > Fax: (+52 612) 122 53 22 > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org <javascript:;> mailing list -- To UNSUBSCRIBE and > more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.-- Bert Gunter "Data is not information. Information is not knowledge. And knowledge is certainly not wisdom." -- Clifford Stoll [[alternative HTML version deleted]]
Your model is producing -Inf entries in the vector Be (in function modl and LL) at some stage during the optimization process. You should first do something about that before anything else. Berend> On 17 Oct 2015, at 03:01, Bert Gunter <bgunter.4567 at gmail.com> wrote: > > I made no attempt to examine your details for problems, but in general, > > My problem >> is that the results change a lot depending on the initial values... I can't >> see what I am doing wrong... >> >> This is a symptom of an overparameterized model: The parameter estimates >> are unstable even though the predictions may not change much. In other >> words, your model may be too complex for your data. > > > Whether that is true here, you or others will have to determine. Try > simplifying your model as a start. > > -- Bert > > > >> >> >> # Data >> x <- 1995:2010 >> B <- c(3500, 3200, 3000, 2800, 2600, 3000, 3200, 3800, 4200, 4300, 4400, >> 4400, 4500, 4600, 5000, 4300) >> Ct <- c(912, 767, 642, 482, 353, 331, 332, 309, 366, 402, 392, 478, 408, >> 434, 407, 637) >> a <- c(0.539, 0.603, -0.948, 0.166, 1.895, 0.786, 0.901, 0.844, 0.337, >> 0.429, 0.304, 0.230, 1.001, 0.750, 0.507, 1.502) >> Ag <- 0.55 >> >> # Function with quantity to minimize >> modl <- function(par) { >> ro <- par[1] >> ko <- par[2] >> n <- length(B) >> Be <- rep(NA, n) >> Be[1] <- ko * Ag >> for ( k in 2:n) >> Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1] >> err <- (log(B) - log(Be))^2 >> ee <- sqrt( sum(err)/(n-2) ) >> LL <- (1/(sqrt(2*pi)*ee)) * exp( -(err/(2*ee^2) ) ) >> -crossprod(LL) >> } >> >> # Using function optim() >> par.optim <- optim(par = list(ro=0.4, ko=8000), modl, method = "BFGS") >> ro <- par.optim$par[1] >> ko <- par.optim$par[2] >> >> # estimated values of "B" >> n <- length(B) >> Be <- rep(NA, n) >> Be[1] <- ko * Ag >> for ( k in 2:n) >> Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1] >> >> # Plot, estimation of "B" seems reasonable.... >> plot(x, B, ylim=c(1000, 7000)) >> lines(x, Be, col="blue", lwd=2) >> >> >> # ... but it is very sensible to initial values... >> par.optim2 <- optim(par = list(ro=0.4, ko=10000), modl, method = "BFGS") >> ro2 <- par.optim2$par[1] >> ko2 <- par.optim2$par[2] >> >> Be2 <- rep(NA, n) >> Be2[1] <- ko2 * Ag >> for ( k in 2:n) >> Be2[k] <- Be2[k-1] + ro2 * a[k-1] * Be2[k-1] * (1 - Be2[k-1]/ko2) - >> Ct[k-1] >> >> lines(x, Be2, col="blue", lwd=2, lty=3) >> >> >> >> # Uing mle2 function >> library(bbmle) >> LL <- function(ro, ko, mu, sigma) { >> n <- length(B) >> Be <- rep(NA, n) >> Be[1] <- ko * Ag >> for ( k in 2:n) >> Be[k] <- Be[k-1] + ro * a[k-1] * Be[k-1] * (1 - Be[k-1]/ko) - Ct[k-1] >> err <- log(B) - log(Be) >> R <- (dnorm(err, mu, sigma, log=TRUE)) >> -sum(R) >> } >> >> Bc.mle <- mle2(LL, start = list(ro=0.4, ko=8000, mu=0, sigma=1)) >> summary(Bc.mle) >> >> ro3 <- coef(Bc.mle)[1] >> ko3 <- coef(Bc.mle)[2] >> >> Be3 <- rep(NA, n) >> Be3[1] <- ko3 * Ag >> for ( k in 2:n) >> Be3[k] <- Be3[k-1] + ro3 * a[k-1] * Be3[k-1] * (1 - Be3[k-1]/ko3) - >> Ct[k-1] >> >> lines(x, Be3, col="red", lwd=2) >> >> >> -- >> >> H?ctor Villalobos <Hector.Villalobos.O at gmail.com <javascript:;>> >> >> Depto. de Pesquer?as y Biolog?a Marina >> >> Centro Interdisciplinario de Ciencias Marinas-Instituto Polit?cnico >> Nacional >> >> CICIMAR-I.P.N. >> >> A.P. 592. Colonia Centro >> >> La Paz, Baja California Sur, M?XICO. 23000 >> >> Tels.: (+52 612) 122 53 44; 123 46 58; 123 47 34 ext.: 81602 >> >> Fax: (+52 612) 122 53 22 >> >> [[alternative HTML version deleted]] >> >> ______________________________________________ >> R-help at r-project.org <javascript:;> mailing list -- To UNSUBSCRIBE and >> more, see >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. > > > > -- > Bert Gunter > > "Data is not information. Information is not knowledge. And knowledge is > certainly not wisdom." > -- Clifford Stoll > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.