Go back to Flu Mortality

data2 = c(0.8760308, 14.6437410, 28.4196368, 47.2210290, 73.0648590, 136.9674535, 166.9442884, 163.0097641, 125.1769971, 70.4576305, 44.8616645, 71.3973124, 34.0708812, 76.8866802, 36.8469589, 36.9518742, 22.1994882, 27.5857962, 4.1047848)

t2 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18)

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)))
})

}

incidence = function(t, par) { state1 = c(S=par[3]-par[4], I = par[4]) pars1 = c(beta = par[1], alpha = par[2], N = par[3]) # Junling: use ode45 method (Runge-Kutta 4/5) out1 = ode(y = state1, times = t, func = SIR, parms = pars1, method = 'ode45')

# Junling: use diff function to simplify code	

return (-diff(out1[,“S”])) }

logl = function(par, t, y) { print(par)

Beta1 = par[1] Alpha1 = par[2] Ne = par[3] I1 = par[4] Beta2 = par[5] Alpha2 = par[6] I2 = par[7]

par1 = c(Beta1, Alpha1, Ne, I1) par2 = c(Beta2, Alpha2, Ne, I2)

lamb1 = incidence(t, par1) lamb2 = incidence(t, par2)

loglik = sum(y*log(lamb1 + lamb2) - lamb1 - lamb2)

return(loglik) }

fit_t = c(0:19)

p = optim(c(0.3, 0.2, 1377, 1, 0.3, 0.2, 1), logl, t = fit_t, y = data2, method = “BFGS”, control = list(fnscale = -1, maxit = 300), hessian = TRUE)

if (p$convergence==0) { # print fitted value print(p$par) # print the covariance matrix of parameters print(solve(-p$hessian)) # print the maximum log-likelihood print(p$value)

para = p$par

Beta1 = para[1] Alpha1 = para[2] Ne = para[3] I1 = para[4] Beta2 = para[5] Alpha2 = para[6] I2 = para[7]

par1 = c(Beta1, Alpha1, Ne, I1) par2 = c(Beta2, Alpha2, Ne, I2)

lamb1 = incidence(fit_t, par1) lamb2 = incidence(fit_t, par2)

plot(t2, data2, col='red', ylab='incidence', main='comparison') lines(t2, data2, col='red') points(t2, lamb1, col = “purple”) lines(t2, lamb1, col = “purple”) points(t2, lamb2, col = “green”) lines(t2, lamb2, col = “green”) points(t2, lamb1 + lamb2, col = “blue”) lines(t2, lamb1 + lamb2, col = “blue”)

} else print(“Maximization did not converge\n”)