[R-sig-dyn-mod] Fitting SIR model with mle2 using a binomial distribution

[Thomas Petzoldt]

Simeon,

it seems that mle2 assigns unrealistic parameters during the fit. If we 
introduce a cat() function in sir.nll:


sir.nll <- function(b, g) {
   cat("b=", b, ", g=", g, "\n")  # <---------
   parms <- c(exp(b), exp(g))
   out <- as.data.frame(ode(y = c(S0, I0, R0),
     times = times, closed.sir.model, parms = parms, hmax = 1/120))
   nll <- -sum(dbinom(x = round(dat$prev*dat$N), size = dat$N,
     prob = out[,3]/500, log = TRUE))
}

times <- 1:14
S0 <- 499; I0 <- 1; R0 <- 0
params <- list(b=-6, g=-0.7)
fit <- mle2(sir.nll, start = params)

## output:
b= -6 , g= -0.7
b= -5.999 , g= -0.7
b= -6.001 , g= -0.7
b= -6 , g= -0.699
b= -6 , g= -0.701
b= 4757.967 , g= -1823.75
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
       such that in the machine, T + H = T on the next step
      (H = step size). Solver will continue anyway.
In above message, R =
[1] 1 0
DINTDY-  T (=R1) illegal


... we see that the differential equation solver breaks because mle2 
assigned completely unrealistic values to b and g.

I would therefore recommend to constrain the parameters within realistic 
ranges (i.e. an optimizer supporting 'lower' and 'upper'). The last 
example on the ?mle2 page shows one possibility how to do it, and 
another possibility would be to use transformation formulae like this:

http://cran.r-project.org/web/packages/FME/FME.pdf#26



Hope it helps

Thomas P.


-- 
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

E-Mail: thomas.petzoldt at tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt