Hello list:

A while back I asked a question about fitting population growth models to Venezuelan census data.  This post to the list is intended as a follow up to those questions, in case anyone is interested.  Once again, the census data and the R/simecol/FME code to fit to the generalized logistic growth model (as given in Petzoldt's "Konstruktion ökologischer Mo delle mit der Open Source Software R", page 55):

censos.csv-------------------------------------------------
"time"    "pob"
1873    1732411
1881    2005139
1891    2221572
1920    2479525
1926    2814131
1936    3364347
1941    3850771
1950    5034838
1961    7523999
1971    10721522
1981    14516735
1990    18105265
2001    23232553
2011    27150095
-----------------------------------------------------

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

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

modelo_generalizado <- new("odeModel",
        main=function(time,init,parms, ...) {
            x <- init
            p <- parms
            dpob <- x["pob"]^p["alfa"]*p["r"]*(1-(x["pob"]/p["K"])^p["beta"])^p["gamma"]
            list(dpob)},
        init=c(pob=1732411),
        parms=c(r=1.488546e-08,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")

---------------------------------------

At that time, I had asked the reasons as to why the attempts to use FME's modFit had failed, giving the following error message: "need at least 2 non NA values to interpolate with approx".   simecol's fitOdeModel (using the PORT algorithm) worked only with specific choice of parameters and the rk family of solvers.

In a response, Thomas Petzoldt suggested running a test for simultanaeous identifiability of parameters and fixing some while fitting the others.  This I did and indeed, for the generalized logistic growth models, the collinearity is much much greater than 20, even when fixing the K parameter (asymptotic population level) based upon knowledge (guesswork, rather) of the population at inflexion point.  Still, I could not used modFit.

After banging my head against this problem for sometime, and even running the ode Model with the Euler method without using simecol, I realized that the problem was in this line:

            dpob <- x["pob"]^p["alfa"]*p["r"]*(1-(x["pob"]/p["K"])^p["beta"])^p["gamma"]

For some choice of parameters (while searching for the optimum parameter set), the current population at time t (x["pob"]) would exceed the K parameter and thus (1-(x["pob"]/p["K"])^p["beta"]) could take on negative values, which, when raised to a certain exponent p["gamma"]<1, would mean trouble, thus accounting for those NA values.

Of course, in theory x["pob"] would gradually increase towards K and would never surpass that level but when using numerical differential equation solvers ("lsoda", "rk", and the like) there was no guarantee that x["pob"] would not shoot past K.  So I patched the following fix in the main equation simecol slot:

            if (x["pob"]<=p["K"])
                dpob <- x["pob"]^p["alfa"]*p["r"]*(1-(x["pob"]/p["K"])^p["beta"])^p["gamma"]
            else
                dpob <- 0

Although now modFit works, its still a patch and it seems "artificial" because it's not the way the model was intended to be.  So my question now is: Is there any way to "tell" the deSolver in simecol that x["pob"] cannot shoot past K in the dpob differential equation?

Thanks and Regards,

jose loreto romero palma

P.S. Season's greetings and happy holidays.

	[[alternative HTML version deleted]]

