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

Thomas Petzoldt Thomas.Petzoldt at TU-Dresden.de
Sat Dec 29 18:28:32 CET 2012


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>


More information about the R-sig-dynamic-models mailing list