# go back to Flu Mortality
#Data data = c(3, 6, 25, 73, 222, 294, 258, 237, 191, 125, 69, 27, 11, 4)
#Time t = c(1:14)
plot(t, data, col = “red”) lines(t, data, col = “red”)
# simple_sir_model_example_data.pdf
# Junling: Here I used standard incidence (beta SI/N), so that beta # is easier to interpret (per individual contact rate). # R0=beta/alpha 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))) })
}
#Junling: define the incidence curve 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"]))
}
#Here I used beta and alpha as the parameters to be maximized #Junling: Here I added the population N as a parameter to be passed # to the SIR model logl = function(par, t, y) {
lamb1 = incidence(t, par) # Junling: do not remove NA, because that implies # lamb1<0, which is caused by wrong parameter value loglik = sum(y*log(lamb1) - lamb1) # Junling: Simply return log-likelyhood, # and use optim() to maximize # optim can handle NA values return(loglik)
}
# Junling: because we use standard incidence now, to have an epidemic, # we need beta>alpha so that R0>1. # We use the initial values beta = 0.3 and alpha = 0.2, for example # Junling: here I added fnscale=-1, which actually make optim to maximize # note we need one more data points than data to take differences fit_t = 1:15 p = optim(c(0.3, 0.2, 1800, 1), logl, t = fit_t, y = data, 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) #plot lamb1 = incidence(t=fit_t, par=p$par) plot(t, lamb1, col='blue', ylab='incidence', main='comparison') lines(t, lamb1, col='blue') plot(t, data, col = "red") lines(t, data, col = "red")
} else
print("Maximization did not converge\n")