# 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")