[R-sig-dyn-mod] parameter fitting using fitOdeModel versus FME

[Soetaert, Karline]

Dear Jose,

Sorry that it took so long to answer.

The problem is not in FME, but in the solver that you use: rk4. 

This is a fixed-time-step runge kutta, and its time-step is chosen depending on the time differences in "times".
So when you put "times" = seq(0,36,0.01) as in your definition of your model, then rk4 will use a time step = 0.01 and all is ok; this time step is small enough.


Now, before you fit the model using FME and call modFit, you set: 

times(iwojima) <- obs$time. 

obs has resolution of a day, so rk4 will take a time step of 1 day. This is clearly very imprecise.  So, FME will fit a very imprecise model to the data. To create the FME "output" of the best-fit model though, you use again a time step of 0.01 day, so that is why the fit looks so bad (although it has much improved).

So.
1. To use modFit, it is not necessary to change "times" so that it matches the observation times. It means that you do not need to alter times(iwojima). 

2. *never* use rk4 unless for demonstration. The default integration routine in deSolve is usually the best, so I would not change it. If you do want to use a runge-kutta, then use rk.

Hope this helps,


Karline



-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of jose romero
Sent: Tuesday, February 22, 2011 2:57 PM
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] parameter fitting using fitOdeModel versus FME

Hello list:

I am using a Lanchester ODE model of the battle of IwoJima as an excercise in parameter fitting.  In this case, the parameters to estimate are the american and japanese effective rates.  There is day-to-day data of remaining combatant american troops during the 36-day battle.  However, from the Japanese side there is no such data- all that is known is that at day 0 there were approximately 21500 troops and by day 36, for all practical purposes, there were 0 japanese troops still fighting.  What this means is that the corresponding observations in the japanese column are all "NA".  Thus the problem is to fit parameters when some observations are missing.

If i try to fit the parameters by fitOdeModel using the data frame with its NA's, i get an error, as NA's are not numerical values.  By using only the american observations and ommiting the two known japanese observation values, i get a reasonable estimate for the effectiveness rates that seems to be much better than that given in the Engel (1954) paper.  

On the other hand, when I use FME i get parameter estimates that do not lead to a good fit between the model and the historical data.  What am I doing wrong? Could it be that this case is not well suited for FME as the error terms do not meet the assumptions of being normally distributed or homoscedastic?

Bellow i have included the data (which must be put into a file called "iwo_jima.csv" in the working directory), the iwo jima model using fitOdeModel and last, the same model using FME.

Thanks and Regards,
jose romero

-------------------------------------------------------------------------------------------------------------------
DATA (must create file called "iwo_jima.csv" in the working directory
-------------------------------------------------------------------------------------------------------------------
"time";"ame";"jap"
0;0;21500
1;52839;"NA"
2;50945;"NA"
4;56031;"NA"
5;53749;"NA"
6;53155;"NA"
7;65250;"NA"
8;64378;"NA"
9;62874;"NA"
10;62339;"NA"
11;61405;"NA"
12;60667;"NA"
13;59549;"NA"
14;59345;"NA"
15;59081;"NA"
16;58779;"NA"
17;58196;"NA"
18;57259;"NA"
19;56641;"NA"
20;54792;"NA"
21;55308;"NA"
22;54796;"NA"
23;54038;"NA"
24;53938;"NA"
25;53347;"NA"
26;53072;"NA"
27;52804;"NA"
28;52735;"NA"
29;52608;"NA"
30;52507;"NA"
31;52462;"NA"
32;52304;"NA"
33;52155;"NA"
34;52155;"NA"
35;52155;"NA"
36;52140;0

---------------------------------------------------------------------------------------------------------------------
fitOdeModel VERSION
---------------------------------------------------------------------------------------------------------------------
library(simecol)

iwojima <- new("odeModel",
    main = function (time, init, parms, ...) {
        x <- init
        p <- parms
        f <- approxTime1(inputs, time, rule = 2)["refuerzos"]
        dame <- f - p["j"] * x["jap"]
        djap <- - p["a"] * x["ame"]
        list(c(dame, djap))
    },
    parms = c(j=0.0577, a=0.0106),
    times = c(from=0, to=36, by=0.01),
    init = c(ame=0, jap=21500),
    inputs = as.matrix(data.frame(
                dia = c(0, 0.999, 1, 1.999, 2, 2.999, 3, 4.999,
                        5, 5.999, 6, 36),
                refuerzos = c(54000, 54000, 0, 0, 6000, 6000, 0, 0,
                            13000, 13000, 0, 0))),
    solver = "rk4"
)
#iwojima
iwojima.ajustado <- iwojima

#simular modelo original y registrar resultados
sim1 <- out(sim(iwojima))

#se lee la data observacional para ajustar los parametros del modelo obs <- read.csv2("iwo_jima.csv")

#ahora ajustamos el modelo
resultado <- fitOdeModel(iwojima,
                         obstime=obs$time,
                         yobs=obs[2],
                         fn=ssqOdeModel,
                         method="PORT",
                         scale=c(1/0.1, 1/0.01),
                         lower=c(j=0,a=0)) resultado

#actualizamos los parametros y volvemos a simular
parms(iwojima.ajustado) <- resultado$par
sim2 <- out(sim(iwojima.ajustado))

#finalmente graficamos
matplot(sim1$time, sim1[2:3],
        main = "Modelo1 : Engel(1954)",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"),
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")
matplot(sim2$time, sim2[2:3],
        main = "Modelo2 : fitOdeModel",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"),
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")

---------------------------------------------------------------------------------------------------------------------
FME VERSION
---------------------------------------------------------------------------------------------------------------------
library(simecol)

iwojima <- new("odeModel",
    main = function (time, init, parms, ...) {
        x <- init
        p <- parms
        f <- approxTime1(inputs, time, rule = 2)["refuerzos"]
        dame <- f - p["j"] * x["jap"]
        djap <- - p["a"] * x["ame"]
        list(c(dame, djap))
    },
    parms = c(j=0.0577, a=0.0106),
    times = c(from=0, to=36, by=0.01),
    init = c(ame=0, jap=21500),
    inputs = as.matrix(data.frame(
                dia = c(0, 0.999, 1, 1.999, 2, 2.999, 3, 4.999,
                        5, 5.999, 6, 36),
                refuerzos = c(54000, 54000, 0, 0, 6000, 6000, 0, 0,
                            13000, 13000, 0, 0))),
    solver = "rk4"
)
#iwojima
iwojima.ajustado <- iwojima

#simular modelo original y registrar resultados
sim1 <- out(sim(iwojima))

#se lee la data observacional para ajustar los parametros del modelo obs <- read.csv2("iwo_jima.csv")

#ahora ajustamos el modelo
library(FME)
fCosto <- function(p) {
    parms(iwojima) <- p
    out <- out(sim(iwojima))
    cost <- modCost(out, obs, weight="std")
    cost
}
times(iwojima) <- obs$time
ajuste <- modFit(f = fCosto, p = parms(iwojima))
summary(ajuste)

#actualizamos los parametros y volvemos a simular
parms(iwojima.ajustado) <- coef(ajuste)
sim2 <- out(sim(iwojima.ajustado))

#finalmente graficamos
matplot(sim1$time, sim1[2:3],
        main = "Modelo1 : Engel(1954)",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"),
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")
matplot(sim2$time, sim2[2:3],
        main = "Modelo2 : FME",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"),
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")





	[[alternative HTML version deleted]]