Setting up an SIR disease model

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)

Setting up and plotting an optimiation routine, fitting the SIR model to the data

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