[R-sig-dyn-mod] Failed attempt to use modFit- part 2

[Thomas Petzoldt]

Hello,

thanks for your detailed feedback and analysis. I rewrote your example a 
little bit using with() which may cost a little bit CPU time but 
increases readability.

The reason for your difficulties is the high collinearity. It can be 
seen when plotting simulation and data (see last two lines of example 
below).

The essential problem lies in the poor representation of the saturation 
K by the given data, because it is questionable how much the population 
will grow after the last point, so it is a model-specification problem 
(inadequate model given the data) and not just a technical issue.


Thomas P.





#####################################################################

library("simecol")
library("FME")

censo <- read.table("censos.csv",header=TRUE)

modelo_generalizado <- new("odeModel",
   main = function(time,init,parms, ...) {
       with(as.list(c(init, parms)), {
         if (pob < K)
             dpob <- pob ^ alfa * r *(1-(pob/K)^beta)^gamma
         else
             dpob <- 0

         list(dpob)
       })
   },
   init = c(pob=1732411),
   parms = c(r=1.488546e08, K=64672787,
     alfa=1.943699, beta=1,gamma=4.223137),
   times = seq(from=1873, to=2011, by=1),
   solver = "lsoda"
)

#función de costo
fCosto <- function(p) {
     parms(modelo_generalizado) <- p
     out <- out(sim(modelo_generalizado))
     return(modCost(out, censo, weight="std"))
}

#ajuste del modelo
resultado <- modFit(fCosto, parms(modelo_generalizado) )
parms(modelo_generalizado) <- resultado$par
parms(modelo_generalizado)
cat("SC : ",ssqOdeModel(NULL, modelo_generalizado, censo$time, 
censo[2]), "\n")

# plot the model
modelo_generalizado <- sim(modelo_generalizado)
plot(modelo_generalizado, obs=censo)