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

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)

Hi,

here another ODE model that fits your data rather well. The approach 
uses deSolve without simecol. It applies a two-equation growth model
that splits the population in a non-reproducing and a an actively 
reproducing fraction.

Note also that the parameters are again highly correlated and that 
versions with fixed parameter "a" would also result in a plausible 
curves, but with very different r and K.


Thomas P.


##===================================================================

library(FME)
library(deSolve)

## model function
model <- function(t, y, p) {
   with(as.list(c(y, p)), {
     dy1 <- -a * y1               # "inactive" part of population
     dy2 <-  a * y1 + r * y2 * (1 - y2 / K)  # growing population
     list(c(dy1, dy2), pob = y1 + y2)        # total pop = y1 + y2
   })
}

## cost function
fCost <- function(p) {
   parms[names(p)] <- p    # "names" allow to fit parameter subsets
   out <- ode(init, times = censo$time, model, parms = parms)
   cost <- modCost(out, censo)
   cat(parms, " : ", cost$model, "\n")
   cost
}

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

init  <- c(y1 = 1732411, y2 = 0)
parms <- c(r = 0.05, K = 3.0e7, a = 0.001)
times <- 1873:2100

out <- ode(init, times, model, parms)

plot(censo, xlim = c(1850, 2100), ylim = c(0, 5e7))
lines(out[, 1], out[, 4]) ## initial parameters

##--------------------------------------------------------------------


## fit model (all parameters)
result <- modFit(fCost, p = parms[c("r", "K", "a")])

## alternative: leave "a" untouched
#result <- modFit(fCost, p = parms[c("r", "K")])

## retrieve and display parameters
(parms[names(coef(result))] <- coef(result))
deviance(result)

## simulate fitted model and plot results
out <- ode(init, times = 1873:2100, model, parms)
lines(out[, 1], out[, 4], col = "red")
legend("topright", lty = 1, col = c("black", "red"),
   legend = c("initial parms", "fitted"))


-------------- next part --------------
A non-text attachment was scrubbed...
Name: censos.csv
Type: application/vnd.ms-excel
Size: 259 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-dynamic-models/attachments/20121229/4b50b2f8/attachment.xlb>