Go back to Flu Mortality
library(bbmle)
library(deSolve)
data1 = c(106, 121, 136, 156, 183, 248, 279, 276, 239, 185, 160, 187, 150, 193, 153, 153, 138, 143, 119)
t1 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18)
base1 = c(105.1240, 106.3563, 107.5804, 108.7790, 109.9351, 111.0325, 112.0557, 112.9902, 113.8230, 114.5424, 115.1383, 115.6027, 115.9291, 116.1133, 116.1530, 116.0481, 115.8005, 115.4142, 114.8952)
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)
out1 = ode(y = state1, times = t, func = SIR, parms = pars1, method = 'ode45')
out2 = ode(y = state2, times = t, func = SIR, parms = pars2, method = 'ode45')
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))
}
loglik_poisson = function(theta, data, mean)
{
sum(dpois(x=data, lambda=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)
}
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 theta0 contains no r,
if (is.na(theta0['r']))
theta0 = c(theta0, r = 10000)
while (TRUE)
{
no_r = is.na(theta0['r']);
if (!no_r)
{
p = SIR_fit_given_logL(times, data, base, theta0, model, loglik_negbin)
print(p)
}
# if the fitted shape parameter r is too large (> 1000)
# then it is approximately poisson
# in this case we use poisson
r1 = coef(p)[“r”]
if (no_r || r1"r" > 1000) {
# we need to remove r from theta
if (!no_r) {
old_theta0 = coef(p)
theta0 = c()
for (i in 1:length(old_theta0))
if (names(old_theta0[i]) != 'r')
theta0 = c(theta0, old_theta0[i])
}
p = SIR_fit_given_logL(times, data, base, theta0, model, loglik_poisson)
print(p)
}
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))
print(theta0)
if (is.na(theta0['r']))
theta0 = c(theta0, r = 10000)
}
return(list(fit=p, conf = conf, r=coef(p)[“r”]))
}
c = SIR_fit(c(0:19), data1, base1, c(beta1 = 0.3, alpha1 = 0.2, N = 1377, I1 = 1, beta2 = 0.3, alpha2 = 0.2, I2 = 1), model_SIR)
print©
aicc = AICc(c"fit", nobs = length(data1))
print(aicc)