library(package = "deSolve")
SIR <- function(t, y, parms) {
with(as.list(c(y, parms)), {
dS <- -beta*S*I
dI <- beta*S*I - gamma*I
dR <- gamma*I
return(list(c(dS, dI, dR)))
})
}
p_vec <- c(beta = 0.5, gamma = 0.1)
t_vec <- seq(from = 0, to = 15, by = 0.1)
y0 <- c(S = 0.99, I = 0.01, R = 0)
sir_out <- ode(y = y0, times = t_vec, parms = p_vec, func = SIR)
plot(las = 1, x = sir_out[,1], y = sir_out[,2], xlab = "Time", ylab = "Desity", type = "l", ylim = c(0, 1))
lines(x = sir_out[,1], y = sir_out[,3], col = "red")
lines(x = sir_out[,1], y = sir_out[,4], col = "blue")
legend(x = "topright", bty = "n", lty = 1, col = c("black", "red", "blue"), legend = c("S", "I", "R"))
## Entering and vizualizing real data
flu <- c(3, 8, 28, 76, 222, 293, 257, 237, 192, 126, 70, 28, 12, 5)
day <- 0:(length(flu)-1)
plot(x = day, y = flu, pch = 19, xlab = "Time (days)", ylab = "Density", las = 1)
lines(x = day, y = flu)
flu_data <- data.frame(day, flu)
parm0 <- c(beta = 0.01, gamma = 0.01)
y0 <- c(S = 762, I = 1, R = 0)
sse_sir <- function(parms, init, model, data) {
cases <- data[,2]
out <- ode(y = init, times = data[,1], func = model, parms = parms, hmax = 0.01)
sse <- sum((out[,3] - cases)^2)
}
fit_out <- optim(par = parm0, fn = sse_sir, model = SIR, init = y0, data = flu_data)
p_fit <- c(beta = fit_out$par[["beta"]], gamma = fit_out$par[["gamma"]])
fit_ode <- ode(y = y0, parms = p_fit, func = SIR, times = t_vec)
plot(las = 1, x = fit_ode[,1], y = fit_ode[,2], xlab = "Time", ylab = "Desity", type = "l", ylim = c(0, 763))
lines(x = fit_ode[,1], y = fit_ode[,3], col = "red")
lines(x = fit_ode[,1], y = fit_ode[,4], col = "blue")
points(x = day, y = flu, col = "purple", pch = 19)
lines(x = day, y = flu, col = "purple")
legend(x = "right", bty = "n", lty = 1, col = c("black", "red", "blue", "purple"), legend = c("S", "I", "R", "data"), pch = c(NA, NA, NA, 19))