HEADS - HIDA - COMPSTAT: Assignment 3

Instructions

Consider the SIR model discussed in class. Using and adapting the code provided, deliver a report on different scenarios simulated using parameters that could be estimated from:

1. Portuguese data (all)

2. Other countries’ data (all)

3. Selected periods of data (e.g. using landmark windows or sliding windows in data used to fit beta and gamma)

Submitted a PDF file with the report, explaining and interpreting the simulated scenario(s).

Preparing the data

For simplicity the code used for fetching the cases_dataset is not printed, but is available in the embedded code.

First we will filter the data so we start modelling from the first recovered case.

N <- 10260893 # Portugal population, estimative for 2020

start <- min(which(cases_dataset$recuperados > 0))
end <- nrow(cases_dataset)

data <- cases_dataset[start:end, ]
data <- data %>%
  mutate(
    time = seq_len(end - start + 1) - 1, S = as.numeric(N - (casos_confirmados + recuperados + n_obitos)),
    I = as.numeric(casos_confirmados), R = as.numeric(recuperados + n_obitos)
  ) %>%
  select(time, S, I, R)

Initial SIR Model

Next, we compute the initial beta and gamma for the SIR Model.

S0 <- data$S[1]
I0 <- data$I[1]
R0 <- data$R[1]

S1 <- data$S[2]
I1 <- data$I[2]
R1 <- data$R[2]

S0_ <- S1 - S0
# S0_ = -beta*I0*S0/N
beta0 <- -(S0_ * N) / (I0 * S0)


I0_ <- I1 - I0
# I0_ = beta*I0*S0/N-gamma*I0
# gamma0 = (beta0*(I0*S0)/N -I0_)/I0

### OR

R0_ <- R1 - R0
# R0_  = gamma*I0
gamma0 <- R0_ / I0

beta0

[1] 0.4556289

gamma0

[1] 0.00591716

This SIR function is a set of three differential equations that relates the S, I and R variables in respect with time:

\(\frac{dS}{dt} = - \frac{\beta I S}{N}\)

\(\frac{dI}{dt} = \frac{\beta I S}{N}- \gamma I\)

\(\frac{dR}{dt} = \gamma I\)

and for that we will use the ode() solver from package deSolve.

####### SIR Model #######
SIR <- function(time, state, pars) { # returns rate of change
  with(as.list(c(state, pars)), {
    dS <- -beta * S * I / N
    dI <- beta * S * I / N - gamma * I
    dR <- gamma * I
    return(list(c(dS, dI, dR)))
  })
}

##### Initial Parameters ########
state <- c(S = S0, I = I0, R = R0)
pars <- list(beta = beta0, gamma = gamma0)
time <- seq(0, 200, by = 1)

###### differential equations ########
out <- ode(y = state, times = time, func = SIR, parms = pars)

###### Initial Plot #############
oldpar <- par(las = 1, tcl = 0.5, mgp = c(3, 0.8, 0), bty = "l")
matplot(out[, 1], out[, -1],
  type = "l", lty = 1:3, lwd = 2,
  col = 1:3, xlab = "time, days", ylab = ""
)
legend("topright", c("S", "I", "R"),
  lty = 1:3, lwd = 2, col = 1:4
)

minor.tick(2, 2, 0.7)
minor.tick(10, 10, 0.5)

par(oldpar)

Incremental SIR Model

Now we create a function that will fit the model to the available data using modFit() from package FME.

sir_step <- function(data, beta, gamma) {

  state <- c(S = data$S[1], I = data$I[1], R = data$R[1])

  time <- seq(0, 200, by = 1) ## output times
    
  Objective <- function(x, parset = names(x)) {
    pars[parset] <- x
    out <- ode(y = state, times = time, func = SIR, parms = pars, method = method)
    ## Model cost
    return(modCost(out, data))
  }
  
  Fit <- modFit(Objective, c(beta = beta0, gamma = gamma0))
  
  Fit$par
  R_0 <- Fit$par[1] / Fit$par[2]
  R_0
  
  out2 <- ode(y = state, times = time, func = SIR, parms = Fit$par, method = method)

  return(list(beta = Fit$par[1], gamma = Fit$par[2], re = R_0, forecast = out2))
}

Here we simulate as if we get new data, and change de model accordingly.

res <- list(beta = beta0, gamma = gamma0)

invisible(saveGIF( # save all graphics in an animated gif
  for (i in seq(2, nrow(data), by = 1)) {
    be <- res$beta
    ga <- res$gamma
    
    slice <- 1:i
    
    data_slice <- data[slice,]

    res <- sir_step(data_slice, 0.5, 0.5)
    
    if(i == 2) {
      out <- res$forecast[1:150, ]
    } else {
      out[min(slice):150, ] <- res$forecast[min(slice):150, ]
    }
    
    oldpar <- par(las = 1, tcl = 0.5, mgp = c(3, 0.8, 0), bty = "l")
    matplot(out[, 1], out[, -1],
      type = "l", lty = 1:3, lwd = 2,
      main = "Animated SIR",
      col = 1:3, xlab = paste0("Time\n", 
                               sprintf("R %6.2f, beta %.3f, gamma %.3f", res$re, res$beta, res$gamma)), ylab = ""
    )
    legend("topright", c("S", "I", "R"),
      lty = 1:3, lwd = 2, col = 1:4
    )
    minor.tick(2, 2, 0.7)
    minor.tick(10, 10, 0.5)
    
    points(data_slice$R, pch = 16, col = "green")
    points(data_slice$I, pch = 16, col = "red")
    points(data_slice$S, pch = 16, col = "black")
    
    par(oldpar)
  },
  "ani_sir.gif",
  interval = 0.5, autobrowse = FALSE, ani.res = 96, ani.height = 500, ani.width = 800
))

For better view, here is a zoomed plot

res <- list(beta = beta0, gamma = gamma0)

invisible(saveGIF( # save all graphics in an animated gif
  for (i in seq(2, nrow(data), by = 1)) {
    be <- res$beta
    ga <- res$gamma
    
    slice <- 1:i
    
    data_slice <- data[slice,]

    res <- sir_step(data_slice, 0.5, 0.5)
    
    if(i == 2) {
      out <- res$forecast[1:150, ]
    } else {
      out[min(slice):150, ] <- res$forecast[min(slice):150, ]
    }
    
    oldpar <- par(las = 1, tcl = 0.5, mgp = c(3, 0.8, 0), bty = "l")
    matplot(out[, 1], out[, -1],
      type = "l", lty = 1:3, lwd = 2,
      main = "Animated SIR Zoom",
      col = 1:3, xlab = paste0("Time\n", 
                               sprintf("R %6.2f, beta %.3f, gamma %.3f", res$re, res$beta, res$gamma)), ylab = "",
      ylim = c(0, max(data$I)*1.1), xlim = c(0, nrow(data)*1.1)
    )
    legend("topright", c("S", "I", "R"),
      lty = 1:3, lwd = 2, col = 1:4
    )
    minor.tick(2, 2, 0.7)
    minor.tick(10, 10, 0.5)
    
    points(data_slice$R, pch = 16, col = "green")
    points(data_slice$I, pch = 16, col = "red")
    points(data_slice$S, pch = 16, col = "black")
    
    par(oldpar)
  },
  "ani_sir_zoom.gif",
  interval = 0.5, autobrowse = FALSE, ani.res = 96, ani.height = 500, ani.width = 800
))