Go back to Flu Mortality

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), {
	
state1 = c(S = N - I1, I = I1)	
pars1 = c(beta = beta1, alpha = alpha1, N = N)

state2 = c(S = N - I2, I = I2)	
pars2 = c(beta = beta2, alpha = alpha2, N = N)

state3 = c(S = N - I3, I = I3)	
pars3 = c(beta = beta3, alpha = alpha3, N = N)
	
out1 = ode(y = state1, times = t, func = SIR, parms = pars1, 	method = 'ode45')

out2 = ode(y = state2, times = t, func = SIR, parms = pars2, 	method = 'ode45')

out3 = ode(y = state3, times = t, func = SIR, parms = pars3, 	method = 'ode45')

	return(cbind(-diff(out1[,"S"]), -diff(out2[,"S"]), -		diff(out3[,"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)
	}
	
	l = model(times, theta)
	
	#"base" is the vector of baseline values
	mean = l[,1] + l[,2] + l[,3] + base
	
	# 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) }

p1 = exp_fit_given_logL(c(0:19), data2, c(beta1 = 0.3, alpha1 = 0.2, N = 1377, I1 = 1, beta2 = 0.3, alpha2 = 0.2, I2 = 1, beta3 = 0.3, alpha3 = 0.2, I3 = 1), model_SIR, loglik_poisson)