Merge pull request #101 from ImperialCollegeLondon/Italy

Italy

.gitignore

@@ -51,6 +51,13 @@ figures/*.pdf
 Rplots.pdf
 web/*
 results/*.csv
+Italy/code/stan-models/*.rds
+Italy/results/*.Rdata
+Italy/results/*.pdf
+Italy/results/*.png
+Italy/figures/*.png
+Italy/figures/*.pdf
+Italy/results/*.RDS
 
 
 

Italy/code/plotting/format-data-plotting.r

@@ -0,0 +1,85 @@
+library(ggplot2)
+library(tidyr)
+library(dplyr)
+library(rstan)
+library(data.table)
+library(lubridate)
+library(gdata)
+library(EnvStats)
+library(matrixStats)
+library(scales)
+library(gridExtra)
+library(bayesplot)
+library(cowplot)
+
+
+#---------------------------------------------------------------------------
+format_data <- function(i, dates, countries, estimated_cases_raw, estimated_deaths_raw, 
+ reported_cases, reported_deaths, out, forecast=0, SIM = FALSE){
+
+ N <- length(dates[[i]])
+ if(forecast > 0) {
+ dates[[i]] = c(dates[[i]], max(dates[[i]]) + 1:forecast)
+ N = N + forecast
+ reported_cases[[i]] = c(reported_cases[[i]],rep(NA,forecast))
+ reported_deaths[[i]] = c(reported_deaths[[i]],rep(NA,forecast))
+ }
+
+ country <- countries[[i]]
+
+ estimated_cases <- colMeans(estimated_cases_raw[,1:N,i])
+ estimated_cases_li <- colQuantiles(estimated_cases_raw[,1:N,i], probs=.025)
+ estimated_cases_ui <- colQuantiles(estimated_cases_raw[,1:N,i], probs=.975)
+ estimated_cases_li2 <- colQuantiles(estimated_cases_raw[,1:N,i], probs=.25)
+ estimated_cases_ui2 <- colQuantiles(estimated_cases_raw[,1:N,i], probs=.75)
+
+ estimated_deaths <- colMeans(estimated_deaths_raw[,1:N,i])
+ estimated_deaths_li <- colQuantiles(estimated_deaths_raw[,1:N,i], probs=.025)
+ estimated_deaths_ui <- colQuantiles(estimated_deaths_raw[,1:N,i], probs=.975)
+ estimated_deaths_li2 <- colQuantiles(estimated_deaths_raw[,1:N,i], probs=.25)
+ estimated_deaths_ui2 <- colQuantiles(estimated_deaths_raw[,1:N,i], probs=.75)
+
+ rt <- colMeans(out$Rt_adj[,1:N,i])
+ rt_li <- colQuantiles(out$Rt_adj[,1:N,i],probs=.025)
+ rt_ui <- colQuantiles(out$Rt_adj[,1:N,i],probs=.975)
+ rt_li2 <- colQuantiles(out$Rt_adj[,1:N,i],probs=.25)
+ rt_ui2 <- colQuantiles(out$Rt_adj[,1:N,i],probs=.75)
+
+ if (SIM == FALSE){
+ mu <- mean(out$mu[,i])
+ mu_li <- quantile(out$mu[,i], probs=.025)
+ mu_ui <- quantile(out$mu[,i], probs=.975)
+ }
+
+
+
+ data_state_plotting <- data.frame("date" = dates[[i]],
+ "country" = rep(country, length(dates[[i]])),
+ "reported_cases" = reported_cases[[i]], 
+ "predicted_cases" = estimated_cases,
+ "cases_min" = estimated_cases_li,
+ "cases_max" = estimated_cases_ui,
+ "cases_min2" = estimated_cases_li2,
+ "cases_max2" = estimated_cases_ui2,
+ "reported_deaths" = reported_deaths[[i]],
+ "estimated_deaths" = estimated_deaths,
+ "deaths_min" = estimated_deaths_li,
+ "deaths_max"= estimated_deaths_ui,
+ "deaths_min2" = estimated_deaths_li2,
+ "deaths_max2"= estimated_deaths_ui2,
+ "rt" = rt,
+ "rt_min" = rt_li,
+ "rt_max" = rt_ui,
+ "rt_min2" = rt_li2,
+ "rt_max2" = rt_ui2)
+
+ if (SIM == FALSE){
+ data_state_plotting$mu_rep = rep(mu, length(dates[[i]]))
+ data_state_plotting$mu_li = rep(mu_li, length(dates[[i]]))
+ data_state_plotting$mu_ui = rep(mu_ui, length(dates[[i]]))
+ }
+
+ return(data_state_plotting)
+
+}
+

Italy/code/plotting/make-plots.r

@@ -0,0 +1,86 @@
+library(lubridate)
+library(ggplot2)
+source("Italy/code/utils/read-data-subnational.r")
+source("Italy/code/plotting/format-data-plotting.r")
+source("Italy/code/plotting/make-three-panel-plots.r")
+source("Italy/code/plotting/make-rt-plot.r")
+
+make_plots_all <- function(filename, SIM=FALSE, label = "", last_date_data, ext = ".png"){
+ print("In subnational code")
+ load(filename)
+ countries <- states
+
+ out <- rstan::extract(fit)
+
+ rt_data_long <- NULL
+ rt_data_wide <- NULL
+
+ interventions <- read_interventions()
+
+ covariates <- interventions
+ covariates$Country <- as.factor(covariates$Country)
+ for (i in 1:length(countries)){
+ print(countries[i])
+
+ data_country_plot <- format_data(i = i, dates = dates, countries = countries, 
+ estimated_cases_raw = estimated_cases_raw, 
+ estimated_deaths_raw = estimated_deaths_raw, 
+ reported_cases = reported_cases,
+ reported_deaths = reported_deaths,
+ out = out, SIM = SIM)
+ # Cuts data on last_data_date
+ data_country_plot <- data_country_plot[which(data_country_plot$date <= last_date_data),]
+ # Read in covariates
+
+ covariates_long <- gather(covariates[which(covariates$Country == countries[i]), 
+ 2:ncol(covariates)], 
+ key = "key", value = "value")
+ covariates_long$x <- rep(NA, length(covariates_long$key))
+ un_dates <- unique(covariates_long$value)
+
+ for (k in 1:length(un_dates)){
+ idxs <- which(covariates_long$value == un_dates[k])
+ max_val <- ceiling(max(data_country_plot$rt_max) +0.3)
+ for (k in idxs){
+ covariates_long$x[k] <- max_val
+ max_val <- max_val - 0.3
+ }
+ }
+
+ print(sprintf("Last line of data: %s", data_country_plot$date[length(data_country_plot$rt)]))
+
+
+ len <- length(data_country_plot$rt)
+ rt_data_state_long <- data.frame("state" = c(as.character(countries[i]), as.character(countries[i])),
+ "x" = c("start", "end"),
+ "rt" = c(data_country_plot$rt[1], 
+ mean(data_country_plot$rt[(len-6):len])),
+ "rt_min" = c(data_country_plot$rt_min[1],
+ mean(data_country_plot$rt_min[(len-6):len])),
+ "rt_max" = c(data_country_plot$rt_max[1],
+ mean(data_country_plot$rt_max[(len-6):len])))
+ rt_data_long <- rbind(rt_data_long, rt_data_state_long)
+
+ # Make the three panel plot
+ make_three_panel_plots(data_country_plot, jobid = JOBID, country = countries[i], 
+ covariates_long = covariates_long, label = label)
+ }
+
+ print("Making rt plot")
+ rt_data_long$x <- factor(rt_data_long$x, levels = c("start", "end"))
+
+ region_to_macro=rbind(
+ data.frame(country= c("Aosta","Liguria","Lombardy","Piedmont"), macro="NorthWest"),
+ data.frame(country=c("Emilia-Romagna","Friuli-Venezia_Giulia","Trento","Bolzano","Veneto"),macro="NorthEast"),
+ data.frame(country=c("Lazio","Marche","Tuscany","Umbria"),macro="Centre"),
+ data.frame(country=c("Abruzzo","Apulia","Basilicata","Calabria","Campania","Molise"),macro="South"),
+ data.frame(country=c("Sardinia","Sicily"),macro="Islands")
+ )
+ names(region_to_macro) <- c("state", "macro")
+ rt_data_long = rt_data_long %>%inner_join(region_to_macro,)
+
+
+ make_rt_point_plot(rt_data_long, JOBID = JOBID, label = label)
+
+
+}

Italy/code/plotting/make-rt-plot.r

@@ -0,0 +1,55 @@
+# Make arrow RT plot
+library(ggplot2)
+library(ggrepel)
+
+make_rt_point_plot <- function(rt_data_long, JOBID, label= ""){
+ # Only choose end
+ rt_data_end <- rt_data_long[which(rt_data_long$x == "end"),]
+ rt_data_end <- rt_data_end[order(-rt_data_end$rt),]
+ rt_data_end$state <- factor(rt_data_end$state, levels = rt_data_end$state)
+
+ rt_data_long$state[which(rt_data_long$state=="Friuli-Venezia_Giulia")]<-"Friuli-Venezia Giulia"
+ rt_data_end$state[which(rt_data_end$state=="Friuli-Venezia_Giulia")]<-"Friuli-Venezia Giulia"
+
+ p1 <- ggplot(rt_data_end) + 
+ geom_point(aes(x = state, y = rt, col=macro), stat="identity") + 
+ geom_errorbar(aes(x = state, ymin = rt_min, ymax = rt_max, col=macro), width=0) + 
+ geom_hline(aes(yintercept=1)) + 
+ xlab("Region") + ylab(expression(R[t])) + 
+ scale_y_continuous(expand = c(0, 0)) + 
+ scale_colour_discrete(name = "") + 
+ coord_flip() + 
+ theme(axis.line = element_line(colour = "black"), 
+ panel.grid.major = element_blank(), 
+ panel.grid.minor = element_blank(),
+ panel.background = element_blank(),
+ legend.key = element_blank(),legend.position = "bottom")+guides(colour = guide_legend(nrow = 2))
+ #ggtitle("Final Rt") 
+ p1
+ ggsave(paste0("Italy/figures/rt_point_final", "_", label, "_", JOBID, ".pdf"), 
+ p1, width = 5, height=10)
+
+ # Only choose start
+ rt_data_start <- rt_data_long[which(rt_data_long$x == "start"),]
+ rt_data_start <- rt_data_start[order(-rt_data_start$rt),]
+ rt_data_start$state <- factor(rt_data_start$state, levels = rt_data_start$state)
+
+ rt_data_start$state[which(rt_data_start$state=="Friuli-Venezia_Giulia")]<-"Friuli-Venezia Giulia"
+
+ p2 <- ggplot(rt_data_start) + 
+ geom_point(aes(x = state, y = rt, col = macro), stat="identity") + 
+ geom_errorbar(aes(x = state, ymin = rt_min, ymax = rt_max, col=macro), width = 0) + 
+ geom_hline(aes(yintercept=1)) + 
+ xlab("Region") + ylab(expression(R[t])) + 
+ scale_y_continuous(expand = c(0, 0)) + 
+ scale_colour_discrete(name = "") +
+ coord_flip() + 
+ theme(axis.line = element_line(colour = "black"), 
+ panel.grid.major = element_blank(), 
+ panel.grid.minor = element_blank(),
+ panel.background = element_blank(),
+ legend.key = element_blank(),legend.position = "bottom") +guides(colour = guide_legend(nrow = 2))
+ #ggtitle("Inital Rt")
+ ggsave(paste0("Italy/figures/rt_point_start_", label, "_", JOBID, ".png"), p2, width = 5, height=10)
+
+}

Italy/code/plotting/make-scenario-plots-top7.r

@@ -0,0 +1,110 @@
+library(dplyr)
+library(lubridate)
+library(grid)
+library(gtable)
+source("Italy/code/plotting/format-data-plotting.r")
+
+make_scenario_comparison_plots_mobility <- function(JOBID, StanModel, len_forecast, last_date_data, 
+ baseline=FALSE, mobility_increase = 20,top=7){
+ print(paste0("Making scenario comparision plots for ", mobility_increase , "%"))
+ load(paste0('Italy/results/sim-constant-mob-', StanModel, '-', len_forecast, '-0-', JOBID, '-stanfit.Rdata'))
+
+ countries <- states
+
+ out <- rstan::extract(fit)
+
+ mob_data <- NULL
+ for (i in 1:length(countries)){
+ data_state_plot <- format_data(i = i, dates = dates, countries = countries, 
+ estimated_cases_raw = estimated_cases_raw, 
+ estimated_deaths_raw = estimated_deaths_raw, 
+ reported_cases = reported_cases,
+ reported_deaths = reported_deaths,
+ out = out, forecast = 0, SIM = TRUE)
+ # Cuts data on last_data_date
+ data_state_plot <- data_state_plot[which(data_state_plot$date <= last_date_data),]
+
+ subset_data <- select(data_state_plot, country, date, reported_deaths, estimated_deaths, 
+ deaths_min, deaths_max)
+ subset_data$key <- rep("Constant mobility", length(subset_data$country))
+ mob_data <- rbind(mob_data, subset_data)
+ } 
+
+ if (baseline == TRUE){
+ load(paste0('Italy/results/sim-increase-mob-baseline-', StanModel, '-', len_forecast, '-', mobility_increase, '-', JOBID, 
+ '-stanfit.Rdata'))
+ out <- rstan::extract(fit)
+ } else {
+ load(paste0('Italy/results/sim-increase-mob-current-', StanModel, '-', len_forecast, '-', mobility_increase, '-', 
+ JOBID, '-stanfit.Rdata'))
+ out <- rstan::extract(fit)
+ }
+ for (i in 1:length(countries)){
+ data_state_plot <- format_data(i = i, dates = dates, countries = countries, 
+ estimated_cases_raw = estimated_cases_raw, 
+ estimated_deaths_raw = estimated_deaths_raw, 
+ reported_cases = reported_cases,
+ reported_deaths = reported_deaths,
+ out = out, forecast = 0, SIM = TRUE)
+ # Cuts data on last_data_date
+ data_state_plot <- data_state_plot[which(data_state_plot$date <= last_date_data),]
+ subset_data <- select(data_state_plot, country, date, reported_deaths, estimated_deaths, 
+ deaths_min, deaths_max)
+ subset_data$key <- rep("Increased mobility", length(subset_data$country))
+ mob_data <- rbind(mob_data, subset_data)
+ } 
+
+ data_half <- mob_data[which(mob_data$key == "Increased mobility"),]
+ mob_data$key <- factor(mob_data$key)
+ data_half$key <- factor(data_half$key)
+
+ #nametrans <- read.csv("Subnational_Analysis/Italy/province_name_translation.csv")
+
+ # To do top 7:
+ if(top==7){
+ mob_data <- mob_data %>% filter(country %in% c("Lombardy","Marche","Veneto","Tuscany","Piedmont","Emilia-Romagna","Liguria")) %>%
+ droplevels()
+
+ data_half <- data_half %>% filter(country %in% c("Lombardy","Marche","Veneto","Tuscany","Piedmont","Emilia-Romagna","Liguria")) %>%
+ droplevels()
+ }
+ if(top==8){
+ # # To do all others"
+ mob_data <- mob_data %>% filter((country %in% c("Abruzzo","Basilicata","Calabria","Campania","Friuli-Venezia_Giulia","Lazio","Molise"))) %>%
+ droplevels()
+ data_half <- data_half %>% filter((country %in% c("Abruzzo","Basilicata","Calabria","Campania","Friuli-Venezia_Giulia","Lazio","Molise"))) %>%
+ droplevels()
+ }
+ if(top==9){
+ # # To do all others"
+ mob_data <- mob_data %>% filter((country %in% c("Bolzano","Trento","Apulia","Sardinia","Sicily","Umbria","Aosta"))) %>%
+ droplevels()
+ data_half <- data_half %>% filter((country %in% c("Bolzano","Trento","Apulia","Sardinia","Sicily","Umbria","Aosta"))) %>%
+ droplevels()
+ }
+
+ last_date_data<-mob_data$date[nrow(mob_data)]
+
+ #mob_data$label <- mob_data$key %>% str_replace_all(" ", "_") %>% recode( Constant_Mobility= "Mobility held constant", Increased_Mobility = "Increased mobility: ",mobility_increase,"% return to pre-lockdown level")
+
+ levels(mob_data$key)=c("Mobility held constant",paste0("Increased mobility: ",mobility_increase,"% return to pre-lockdown level"))
+
+ p <- ggplot(mob_data) +
+ geom_bar(data = mob_data, aes(x = date, y = reported_deaths), stat='identity') +
+ geom_ribbon(aes(x = date, ymin = deaths_min, ymax = deaths_max, group = key, fill = key), alpha = 0.5) +
+ #geom_line(aes(date,deaths_max),color="black",size=0.2)+
+ #geom_line(aes(date,deaths_min),color="black",size=0.2)+
+ #geom_line(aes(date,estimated_deaths),group = key,size=0.5)+
+ geom_line(aes(date,estimated_deaths, group = key, color = key),size = 1) +scale_colour_manual(values= c("skyblue","red"))+
+ #geom_ribbon(aes(x = date, ymin = deaths_min, ymax = deaths_max, fill = "ICL"), alpha = 0.5) +
+ scale_fill_manual(name = "", labels = c("Mobility held constant", paste0("Increased mobility: ",mobility_increase,"% return to pre-lockdown level")), values = c("skyblue","red")) + 
+ scale_x_date(date_breaks = "2 weeks", labels = date_format("%e %b"), limits = c(as.Date("2020-03-02"), last_date_data)) + 
+ #facet_wrap(~country, scales = "free",nrow=7) + 
+ facet_grid(country ~key, scales = "free_y")+
+ xlab("") + ylab("Daily number of deaths") +
+ theme_minimal() + 
+ theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 26), axis.title = element_text( size = 26 ),axis.text = element_text( size = 26),
+ legend.position = "none",strip.text = element_text(size = 26),legend.text=element_text(size=26))
+ ggsave(paste0("Italy/figures/scenarios_increase_baseline-", len_forecast, '-', mobility_increase, '-', JOBID, "top_",top,".png"), p, height = 30, width = 20)
+
+}