Go back to Flu Mortality

library(bbmle)

library(deSolve)

dyn.load(“mymod.so”)

data10 = c(117, 127, 141, 137, 155, 129, 149, 143, 147, 154, 154, 175, 160, 225, 292, 255, 269, 267, 262, 271, 312, 328, 305, 291, 243, 217, 187, 173, 158, 147, 142, 157, 134, 135)

t10 = c(0:33)

base10 = c(110.8113, 111.5707, 112.4397, 113.4061, 114.4563, 115.5756, 116.7482, 117.9574, 119.1863, 120.4174, 121.6333, 122.8168, 123.9512, 125.0205, 126.0096, 126.9047, 127.6931, 128.3638, 128.9076, 129.3171, 129.5866, 129.7128, 129.6941, 129.5313, 129.2273, 128.7867, 128.2164, 127.5250, 126.7231, 125.8226, 124.8370, 123.7810, 122.6702, 121.5213)

model_SIR = function(t, theta) {

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 = alpha1, N = N)
	
out1 = ode(y = state1, times = t, func = "derivs", parms = pars1, jacfunc = "jac", dllname = "mymod", initfunc = "initmod", nout = 1, outnames = "S")

out2 = ode(y = state2, times = t, func = "derivs", parms = pars2, jacfunc = "jac", dllname = "mymod", initfunc = "initmod", nout = 1, outnames = "S")

	return(cbind(-diff(out1[,"S"]), -diff(out2[,"S"])))
})

}

loglik_negbin = function(theta, data, mean) {

# for negative bimonial, mean = r p /(1-p)
r = theta[["r"]];

# if r < 0, return a tiny likelihood
if (r < 0)
{
	return(-1e10)
}

sum(dnbinom(x=data, size=r, mu=mean, log=TRUE))

}

SIR_fit_given_logL = function(times, data, base, 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)
	}
	
	N = theta[["N"]]
      
      if(N < sum(data))
	{
		return(1e10)
	}
	
	for(i in 1:length(theta))
	{
		if(theta[[i]] < 0)
		{
			return(1e10)
		}
	}
	
	l = model(times, theta)
	
	#"base" is the vector of baseline values
	mean = l[,1] + l[,2] + base
	
	# compute the negative log-likelihood
	L = -logL(theta, 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 = “Nelder-Mead”, start = as.list(theta0), control = list(maxit = 1e6*length(theta0), ndeps = 1e-4, 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) }

SIR_fit = function(times, data, base, theta0, model) {

# we start with negative binomial

if (is.na(theta0['r']))
{
	theta0 = c(theta0, r = 2)
}
	
while (TRUE) 
{
	p = SIR_fit_given_logL(times, data, base, theta0, model, loglik_negbin)
	
conf = confint(p)

if(is.vector(conf))
{
	conf = matrix(conf, nrow = 1)
}

if(is.matrix(conf))
{
	break
}

theta0 = coef(conf)
names(theta0) = names(coef(p))

}

return(list(fit=p, conf = conf, r=coef(p)["r"]))

}

c = SIR_fit(c(0:34), data10, base10, c(beta1 = 1.1, alpha1 = 1, N = (sum(data10)*1.166), I1 = 1, beta2 = 1.1, I2 = 1), model_SIR)

print©

aicc = AICc(c"fit", nobs = length(data10))

print(aicc)