Recall that the rows of matrix mu contain the values for \(S\), \(I\), and \(R\), respectively. The columns correspond to the values at all the periods.
S <- mu[1, ]I <- mu[2, ]R <- mu[3, ]N <- S + I + RD <-100- Ndf_lockdown <-tibble(S = S, I = I, R = R, t = time, type ="no-lockdown")df_lockdown
Now, let us create a function that will run the SIR model with a lockdown. The code from the SIR model without lockdown is slightly modified, to introduce the lockdown and its severity. The following function allows us to easily make a simulation depending on the start and end of the lockdown, as well as its severity:
#' Simulate the SIR model with a lockdown#' #' @param lock_start start of the lockdown (number of days from the outbreak)#' @param lock_end end of the lockdown (number of days from the outbreak)#' @param lock_severity severity of the lockdown (in [0,1], 0 for no lockdown)simulate_sir_lockdown <-function(lock_start, lock_end, lock_severity, type) {# Matrix with 3 rows, each corresponding to S, I and R, resp. mu <-matrix(rep(0, 3* (Nt +1)), nrow =3) S <- I <- R <- N <- F <- Re <-rep(0, Nt)# Transition matrix A_t <-vector(mode="list", length = Nt)# First column: initial conditions mu[, 1] <-matrix(c(S0, I0, R0), nrow =3) lockdown <-rep(0, Nt)for (n in1:(Nt -1)) { S[n] <- mu[1, n] I[n] <- mu[2, n] R[n] <- mu[3, n]if (time[n] >= lock_start & time[n] <= lock_end) {# Lockout lockdown[n] <- lock_severity } Re[n] <- reprod_number * (1- lockdown[n]) * S[n] A_t[[n]] <-matrix(c(-beta * (1-lockdown[n]) * I[n], beta * (1- lockdown[n]) * I[n], 0,0, -gamma, gamma,0, 0, 0 ),byrow =TRUE,ncol =3)if (implicit ==0) { mu[, n+1] <- dt * (t(A_t[[n]]) %*% mu[, n]) + mu[, n] } else { mu[n +1, ] <- (I(3) - dt * (t(A_t[[n]])) %/% mu[, n]) } } mu <- mu[, 1:Nt] S_tight <- mu[1, ] I_tight <- mu[2, ] R_tight <- mu[3, ] N_tight <- S_tight + I_tight + R_tight D_tight <-100- N_tighttibble(S = S_tight, I = I_tight, R = R_tight, t = time)}
The datasets for each scenario need to be reshaped to be used in ggplot2().
#' Reshape the data from two scenarios#' #' @param df_scenario_1 data for first scenario#' @param df_scenario_2 data for second scenario#' @param name_scenario_1 name of the scenario in `df_scenario_1`#' @param label_scenario_1 desired label for the name of the first scenario#' @param name_scenario_2 name of the scenario in `df_scenario_2`#' @param label_scenario_2 desired label for the name of the second scenarioreshape_data_graph <-function(df_scenario_1, df_scenario_2, name_scenario_1, label_scenario_1, name_scenario_2, label_scenario_2) { df_scenario_1 |>bind_rows(df_scenario_2) |>pivot_longer(cols =c("S", "I", "R")) |>filter(t <= Nt -1) |>mutate(name =factor(name, levels =c("S", "I", "R")),type =factor( type, levels =c(name_scenario_1, name_scenario_2), labels =c(label_scenario_1, label_scenario_2)))}
We define a theme function to make the graphs pretty (at least for us!).
And we define a function to plot the evolution of the two scenarios with respect to time.
#' Plots the evolution of two scenarios#' #' @param df_plot data with the two scenarios (obtained from #' `reshape_data_graph()`)#' @param lock_start start period of the lockdown#' @param lock_end end period of the lockdown#' @param lock_severity severity of the lockdownplot_scenario <-function(df_plot, lock_start, lock_end, lock_severity) {ggplot(data = df_plot) +geom_rect(data =tibble(x1 = lock_start, x2 = lock_end,y1 =0, y2 = lock_severity),mapping =aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2),fill ="grey", alpha = .8 ) +geom_line(aes(x = t, y = value, linetype = type, colour = name)) +scale_linetype_discrete(NULL) +scale_colour_manual(NULL, values =c("S"="#1b9e77", "I"="#d95f02", "R"="#7570b3")) +labs(x ="Days from the outbreak", y ="Proportion of cases") +theme_paper() +coord_cartesian(xlim =c(0, 150))}
The phase diagrams can also be plotted. Once again, the data need to be reshaped for use in ggplot2.
#' @param df_scenario_1 data for first scenario#' @param df_scenario_2 data for second scenario#' @param name_scenario_1 name of the scenario in `df_scenario_1`#' @param label_scenario_1 desired label for the name of the first scenario#' @param name_scenario_2 name of the scenario in `df_scenario_2`#' @param label_scenario_2 desired label for the name of the second scenarioreshape_data_phase <-function(df_scenario_1, df_scenario_2, name_scenario_1, label_scenario_1, name_scenario_2, label_scenario_2) { df_scenario_1 |>filter(I >0.001) |>group_by(S) |>slice(1) |>ungroup() |>bind_rows( df_scenario_2 |>filter(I >0.001) |>group_by(S) |>slice(1) |>ungroup() ) |>mutate(type =factor( type, levels =c(name_scenario_1, name_scenario_2), labels =c(label_scenario_1, label_scenario_2)) )}
If we want to display some arrows to highlight the direction on the diagrams, we can create a function that will do so for a single scenario.
#' Gives the coordinates to draw a small arrow on the phase diagram of one #' scenario#' #' @param df_plot table with the data of the diagram of one scenario ready for #' use in ggplot2#' @param val_S value for Screate_df_arrow <-function(df_plot, val_S) { df_arrows <- df_plot |>arrange(desc(S)) |>group_by(type) |>filter(S <!!val_S) |>slice(1:2) |>select(S, I, type) |>ungroup() |>mutate(toto =rep(c("beg", "end"), 2)) df_arrows |>select(-I) |>pivot_wider(names_from = toto, values_from = S, names_prefix ="S_") |>left_join( df_arrows |>select(-S) |>pivot_wider(names_from = toto, values_from = I, names_prefix ="I_"),by ="type" )}
Then, let us create a function that will plot the phase diagrams for two scenarios.
