Go back to Flu Mortality

data = c(3, 6, 25, 73, 222, 294, 258, 237, 191, 125, 69, 27, 11, 4)

t = c(1:14)

model_SIR = function(t, theta) {

SIR = function(Time, State, Pars) {
with(as.list(c(State, Pars)), {

	DecreaseS = -1 * beta * S * I / N
	DecreaseI = -1 * alpha * I

	dS = DecreaseS
	dI = -1 * DecreaseS + DecreaseI

	return(list(c(dS, dI)))
	})
}	

with(as.list(theta), {
	
	state = c(S = N - I, I = I)
	pars = c(beta = beta, alpha = alpha, N = N)
	
	out = ode(y = state, times = t, func = SIR, parms = pars, 		method = 'ode45')

	return (-diff(out[,"S"]))
})

}

loglik_poisson = function(data, mean) {

sum(dpois(x=data, lambda=mean, log=TRUE))

}

exp_fit_given_logL = function(times, data, theta0, model, logL) {

# construct a general negative log-likelyhood function

L <- function() 
{
	# reconstruct the vector of named pairs of 
	# parameters from the list of arguments
	l=length(theta0)
	theta=c()
	vars=names(theta0)
	for (i in 1:l) {
		item=c(x=get(vars[i]))
		names(item) <- vars[i]
		theta <-c(theta, item)
	}

	# compute the mean from the model

	mean = model(t, theta)

	# compute the negative log-likelihood

	L = -logL(data, mean)

	# if not a number, return a very small likelihood

	if (is.na(L)) return(1e10)
	return(L)
}

# replace the input arguments of L by the list of parameters

formals(L)←as.list(theta0)

# mle2 # sometimes confint may find a better solution. # In this case we need to refit the model with the better # solution as a starting paoint, because the parameter # names in the returned better fit can be wrong

fit = mle2(L, method = “BFGS”, start = as.list(theta0), control = list(maxit = 1e6*length(theta0)))

# fit=mle2(L, method='BFGS', start=as.list(theta0), control=list(maxit=1e5*length(theta0), reltol=1e-10)) # if there is only one parameter to be fitted # the returned confidence interval is a vector # change it to a matrix

return(fit) }

p = exp_fit_given_logL(c(1:15), data, c(beta = 0.3, alpha = 0.2, N = 1800, I = 1), model_SIR, loglik_poisson)

lamb1 = model_SIR(c(1:15), c(beta = 1.6632951, alpha = 0.7186123, N = 1797.4594924, I = 1.9795067))

plot(t, lamb1, type = “b”, col='blue', ylab='incidence', main='comparison') points(t, data, type = “b”,col = “red”)