Hi together,
 few examples in the literature* showed how
to fit a closed deterministic SIR model to data
on an flu outbreak in a British boarding school.

In both cited examples the maximum likelihood fit 
was done using a Poisson distribution for the 
negative log-likelihood.
This seems to be an appropriate assumption, 
however if the sample size of each observation
differs (and therewith the confidence interval)
a binomial distribution would be better?!

I tried to slightly transform the flue example, adding
a sample size and gave prevalence level instead
of positive counts.

I assume that the prevalence data are subsamples
from a population (N=500) and therefore used the
total population size for the SIR simulation.

BUT, I cannot get it to run (see example below)!

Any suggestions or conceptual criticism is
very much appreciated.

Simeon

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

# start example

library(deSolve)
library(bbmle)

dat <- data.frame(time=1:14, N=round(runif(14,350,500)),
	          prev=c(0.006,0.015,0.05,0.15,0.45,0.6,0.5,0.45,0.4,0.25,0.15,0.05,0.025,0.01))

plot(dat$t, dat$prev*100, type="o", col="red", xlab="time", ylab="prevalence")
	axis(3,at=1:14,labels=paste("N=",round(dat$N,0),sep=""))

closed.sir.model <- function(t, x, params) {
	S <- x[1]
	I <- x[2]
	R <- x[3]
	
	beta <- params[1]
	gamma <- params[2]
	
	dS <- -beta*S*I
	dI <- beta*S*I - gamma*I
	dR <- gamma*I
	
	list(c(dS,dI,dR))
}

sir.nll <- function(b, g){
	times <- seq(0,14)
	parms <- c(exp(b), exp(g))
	
	S0 <- 499
	I0 <- 1
	R0 <- 0
	
	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))
}

params <- list(b=-6, g=-0.7)
   fit <- mle2(sir.nll, start=params0); fit

# end example
	[[alternative HTML version deleted]]

