Commit 04a9d98b authored by Mathilde Grimée's avatar Mathilde Grimée
Browse files

initial commit

parents
File added
This diff is collapsed.
This diff is collapsed.
####################
## Plot adjacency ##
## matrices ##
####################
#In this script, we plot the adjacency matrices for the first of each month,
#at baseline, after the IOM adjustment, and after the FB adjustment
load("../Data_preparation/adjacencies_BASE.rds")
load("../Data_preparation/adjacencies_FB.rds")
load("../Data_preparation/adjacencies_IOM.rds")
library(ggplot2)
library(gridExtra)
plotmatrix <- function(d2, title = NULL, leg.pos = "none"){
d2.df <- data.frame(x=rep(dimnames(d2)[[1]],each=ncol(d2)),y=rep(dimnames(d2)[[2]],times=nrow(d2)),z=as.vector(d2))
plot <- ggplot(data=d2.df,aes(x=x,y=y,fill=z))+
geom_tile() +
ggtitle(title) +
scale_fill_gradient(low = "#000000",
high = "#ffffff",
guide = "colourbar",
limits = c(0, 1))+
labs(y = " ", x= " ")+
theme(legend.position = leg.pos, axis.text.x = element_text(angle=-45))+
coord_fixed(ratio = 1)
return(plot)
}
p <- plotmatrix(adjacencies_BASE, "Adjacencies at baseline")
pdf("../Figures/Adjacencies_BASE.pdf")
p
dev.off()
p2 <- plotmatrix(adjacencies_BASE, "At baseline")
#p_fb_3 <- plotmatrix(adjacencies_FB[,,"2020-03-01"])
p_fb_3_2 <- plotmatrix(adjacencies_FB[,,"2020-03-01"], "Facebook adjustment")
p_fb_4<- plotmatrix(adjacencies_FB[,,"2020-04-01"])
p_fb_5 <- plotmatrix(adjacencies_FB[,,"2020-05-01"])
p_fb_6 <- plotmatrix(adjacencies_FB[,,"2020-06-01"])
p_fb_7 <- plotmatrix(adjacencies_FB[,,"2020-07-01"])
p_fb_8 <- plotmatrix(adjacencies_FB[,,"2020-08-01"])
#p_iom_3 <- plotmatrix(adjacencies_IOM[,,"2020-03-01"])
p_iom_3_2 <- plotmatrix(adjacencies_IOM[,,"2020-03-01"], "IOM adjustment")
p_iom_4 <- plotmatrix(adjacencies_IOM[,,"2020-04-01"])
p_iom_5 <- plotmatrix(adjacencies_IOM[,,"2020-05-01"])
p_iom_6 <- plotmatrix(adjacencies_IOM[,,"2020-06-01"])
p_iom_7 <- plotmatrix(adjacencies_IOM[,,"2020-07-01"])
p_iom_8 <- plotmatrix(adjacencies_IOM[,,"2020-08-01"])
pdf("../Figures/Adjacencies.pdf", width = 8, height = 18)
grid.arrange(p2+labs(y = "March"), p_fb_3_2, p_iom_3_2, p+labs(y = "April"), p_fb_4, p_iom_4, p+labs(y = "May"), p_fb_5, p_iom_5,
p+labs(y = "June"), p_fb_6, p_iom_6, p+labs(y = "July"), p_fb_7, p_iom_7, p+labs(y = "August"), p_fb_8, p_iom_8,
ncol = 3)
dev.off()
#################
## Model plots ##
#################
#In this script, we create plots of the fitted values
#for the three best performing models.
#We also create a plot of the lag distribution for the best performing model.
#dependencies
load("models.Rdata")
load("inputs.Rdata")
load("overview.Rdata")
#axis def
myaxis <- list(epochsAsDate = TRUE,
# how often to put ticks
xaxis.tickFreq = list("%m" = atChange, "%m" = atChange),
# labels
xaxis.labelFormat = "%b")#"%d-%b")
#3 best models
best_models <- list()
for(i in 1:3){
print(i)
best_i <- rownames(overview)[which(overview$BIC == sort(overview$BIC)[i])]
print(best_i)
best_models <- append(best_models, list(best_i))
}
best_inputs <- list()
for(modname in best_models){
print(modname)
best_input <- inputs[[paste("input", substr(modname, 4, nchar(modname)), sep = "")]]
best_inputs <- append(best_inputs, list(best_input))
}
#plot
pdf("../Figures/Plot_best1.pdf")
plotHHH4lag_fitted(profile_par_lag(cases_sts, best_inputs[[1]]), xaxis = myaxis,
units = NULL,
pch = 20,
hide0s = TRUE,
legend.args = list(x = "topleft",
legend = c("Epidemic", "Endemic")),
col = c("#eead3d", "white", "#4a89be"),
names = NULL, ylab = "COVID-19 cases")
dev.off()
pdf("../Figures/Plot_best2.pdf")
plot2 <- plotHHH4lag_fitted(profile_par_lag(cases_sts, best_inputs[[2]]), xaxis = myaxis,
units = NULL,
pch = 20,
hide0s = TRUE,
legend.args = list(x = "topleft",
legend = c("Epidemic", "Endemic")),
col = c("#eead3d", "white", "#4a89be"),
names = NULL, ylab = "COVID-19 cases")
dev.off()
pdf("../Figures/Plot_best3.pdf")
plot3 <- plotHHH4lag_fitted(profile_par_lag(cases_sts, best_inputs[[3]]), xaxis = myaxis,
units = NULL,
pch = 20,
hide0s = TRUE,
legend.args = list(x = "topleft",
legend = c("Epidemic", "Endemic")),
col = c("#eead3d", "white", "#4a89be"),
names = NULL, ylab = "COVID-19 cases")
dev.off()
#lag_distribution
best <- get_vals(best_inputs[[1]])
pdf("../Figures/lag_distr.pdf")
plot(best$lag_distr, main="Lag distribution",
xlab = "lag",
ylab = "weight")
lines(best$lag_distr)
dev.off()
########################################
## Plot the predicted number of cases ##
## under base vs real case counts ##
########################################
#In this script, we create plots comparing the predicted daily cases under the baseline
#scenario versus the observed values.
#We also create plots for the predicted values under the three scenarios (baseline, CF A and CF B)
load("../Data_preparation/parameters.Rdata")
load("mom_pred.Rdata")
load("../Data_preparation/cases.rds")
load("overview.Rdata")
load("inputs.Rdata")
load("models.Rdata")
library(gridExtra)
colors <- c("#f5006e", "#797676")
#baseline scenario
cases.df <- as.data.frame(cases)
df <- cases.df[which(rownames(cases.df) %in% as.character(seq(as.Date("2020-03-03"), as.Date("2020-08-04"), "day"))),]
df <- rbind(df, mom_base$mu_matrix)
df$date <- rep(as.character(seq(as.Date("2020-03-03"), as.Date("2020-08-04"), "day")), 2)
df$date <- as.Date(df$date)
df$`Cases under` <- c(rep("reality", nrow(df)/2), rep("baseline scenario", nrow(df)/2))
p1 <- ggplot(data = df, aes(x=date, y=CH01, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("CH01")+
labs(y = "cases")
p2 <- ggplot(data = df, aes(x=date, y=CH02, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("CH02")+
labs(y = "cases")+
theme(legend.position = "none")
p3 <- ggplot(data = df, aes(x=date, y=CH03, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("CH03")+
labs(y = "cases")+
theme(legend.position = "none")
p4 <- ggplot(data = df, aes(x=date, y=CH04, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("CH04")+
labs(y = "cases")+
theme(legend.position = "none")
p5 <- ggplot(data = df, aes(x=date, y=CH05, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("CH05")+
labs(y = "cases")+
theme(legend.position = "none")
p6 <- ggplot(data = df, aes(x=date, y=CH06, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("CH06")+
labs(y = "cases")+
theme(legend.position = "none")
p7<-ggplot(data = df, aes(x=date, y=CH07, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("CH07")+
labs(y = "cases")+
theme(legend.position = "none")
p8 <- ggplot(data = df, aes(x=date, y=ITC1, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("ITC1")+
labs(y = "cases")+
theme(legend.position = "none")
p9<- ggplot(data = df, aes(x=date, y=ITC2, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("ITC2")+
labs(y = "cases")+
theme(legend.position = "none")
p10 <- ggplot(data = df, aes(x=date, y=ITC4, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("ITC4")+
labs(y = "cases")+
theme(legend.position = "none")
p11 <- ggplot(data = df, aes(x=date, y=ITH1, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors) +
ggtitle("ITH1")+
labs(y = "cases")+
theme(legend.position = "none")
g_legend<-function(a.gplot){
tmp <- ggplot_gtable(ggplot_build(a.gplot))
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
mylegend<-g_legend(p1)
pdf("../Figures/pred_vs_real.pdf")
grid.arrange(p1+theme(legend.position="none"), p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, ncol = 3,
top = "Comparing real case counts with baseline prediction",
mylegend)
dev.off()
#scenarios A and B and baseline
df_cf <- rbind(as.data.frame(mom_base$mu_matrix),
as.data.frame(mom_B$mu_matrix),
as.data.frame(mom_A$mu_matrix))
df_cf$`Cases under` <- c(rep("Baseline scenario", nrow(df_cf)/3),
rep("Scenario B", nrow(df_cf)/3),
rep("Scenario A", nrow(df_cf)/3))
df_cf$date <- rep(as.character(seq(as.Date("2020-03-03"), as.Date("2020-08-04"), "day")), 3)
df_cf$date <- as.Date(df_cf$date)
colors_cf <- c("#f5006e", "#000080", "#ffc125")
p1_cf <- ggplot(data = df_cf, aes(x=date, y=CH01, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("CH01")+
labs(y = "cases")
p2_cf <- ggplot(data = df_cf, aes(x=date, y=CH02, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("CH02")+
labs(y = "cases")+
theme(legend.position = "none")
p3_cf <- ggplot(data = df_cf, aes(x=date, y=CH03, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("CH03")+
labs(y = "cases")+
theme(legend.position = "none")
p4_cf <- ggplot(data = df_cf, aes(x=date, y=CH04, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("CH04")+
labs(y = "cases")+
theme(legend.position = "none")
p5_cf <- ggplot(data = df_cf, aes(x=date, y=CH05, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("CH05")+
labs(y = "cases")+
theme(legend.position = "none")
p6_cf <- ggplot(data = df_cf, aes(x=date, y=CH06, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("CH06")+
labs(y = "cases")+
theme(legend.position = "none")
p7_cf <- ggplot(data = df_cf, aes(x=date, y=CH07, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("CH07")+
labs(y = "cases")+
theme(legend.position = "none")
p8_cf <- ggplot(data = df_cf, aes(x=date, y=ITC1, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("ITC1")+
labs(y = "cases")+
theme(legend.position = "none")
p9_cf <- ggplot(data = df_cf, aes(x=date, y=ITC2, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("ITC2")+
labs(y = "cases")+
theme(legend.position = "none")
p10_cf <- ggplot(data = df_cf, aes(x=date, y=ITC4, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("ITC4")+
labs(y = "cases")+
theme(legend.position = "none")
p11_cf <- ggplot(data = df_cf, aes(x=date, y=ITH1, group = `Cases under`))+
geom_line(aes(color=`Cases under`)) +
scale_color_manual(values=colors_cf) +
ggtitle("ITH1")+
labs(y = "cases")+
theme(legend.position = "none")
mylegend_cf<-g_legend(p1_cf)
pdf("../Figures/Scen_A_B_base.pdf")
grid.arrange(p1_cf+theme(legend.position="none"), p2_cf, p3_cf, p4_cf,
p5_cf, p6_cf, p7_cf, p8_cf, p9_cf, p10_cf, p11_cf, ncol = 3,
top = "Predictions under scenarios A, B and baseline",
mylegend_cf)
dev.off()
######################
## Plot diff ##
## and RI over time ##
######################
#In this script we plot the difference in daily/cumulative cases predicted under the counterfactual
#scenarios versus under the baseline scenario. We also plot the daily relative increase
#in cases between the counterfactual scenarios and the baseline scenario.
#The releative increase is defined as
#(counterfactual scenario-baseline scenario)/baseline scenario.
load("results.rds")
load("mom_pred.Rdata")
load("../Data_preparation/parameters.Rdata")
library(gridExtra)
library(ggplot2)
#CF scenario B
#sum up over Switzerland
mom_B$total_CH <-
mom_B$mu_matrix[,1] + mom_B$mu_matrix[,2] +
mom_B$mu_matrix[,3] + mom_B$mu_matrix[,4] +
mom_B$mu_matrix[,5] + mom_B$mu_matrix[,6] +
mom_B$mu_matrix[,7]
mom_B$total_CH_cum <- NA
for(i in 1:length(mom_B$total_CH)){
mom_B$total_CH_cum[i] <- sum(mom_B$total_CH[1:i])
}
#sum up over Italy
mom_B$total_IT <-
mom_B$mu_matrix[,8] + mom_B$mu_matrix[,9] +
mom_B$mu_matrix[,10] + mom_B$mu_matrix[,11]
mom_B$total_IT_cum <- NA
for(i in 1:length(mom_B$total_IT)){
mom_B$total_IT_cum[i] <- sum(mom_B$total_IT[1:i])
}
#sum up all
mom_B$total <- mom_B$total_CH + mom_B$total_IT
mom_B$total_cum <- NA
for(i in 1:length(mom_B$total)){
mom_B$total_cum[i] <- sum(mom_B$total[1:i])
}
#CF scenario A
#sum up over Switzerland
mom_A$total_CH <-
mom_A$mu_matrix[,1] + mom_A$mu_matrix[,2] +
mom_A$mu_matrix[,3] + mom_A$mu_matrix[,4] +
mom_A$mu_matrix[,5] + mom_A$mu_matrix[,6] +
mom_A$mu_matrix[,7]
mom_A$total_CH_cum <- NA
for(i in 1:length(mom_A$total_CH)){
mom_A$total_CH_cum[i] <- sum(mom_A$total_CH[1:i])
}
#sum up over Italy
mom_A$total_IT <-
mom_A$mu_matrix[,8] + mom_A$mu_matrix[,9] +
mom_A$mu_matrix[,10] + mom_A$mu_matrix[,11]
mom_A$total_IT_cum <- NA
for(i in 1:length(mom_A$total_IT)){
mom_A$total_IT_cum[i] <- sum(mom_A$total_IT[1:i])
}
#sum up all
mom_A$total <- mom_A$total_CH + mom_A$total_IT
mom_A$total_cum <- NA
for(i in 1:length(mom_A$total)){
mom_A$total_cum[i] <- sum(mom_A$total[1:i])
}
#base scenario
#sum up over Switzerland
mom_base$total_CH <-
mom_base$mu_matrix[,1] + mom_base$mu_matrix[,2] +
mom_base$mu_matrix[,3] + mom_base$mu_matrix[,4] +
mom_base$mu_matrix[,5] + mom_base$mu_matrix[,6] +
mom_base$mu_matrix[,7]
mom_base$total_CH_cum <- NA
for(i in 1:length(mom_base$total_CH)){
mom_base$total_CH_cum[i] <- sum(mom_base$total_CH[1:i])
}
#sum up over Italy
mom_base$total_IT <-
mom_base$mu_matrix[,8] + mom_base$mu_matrix[,9] +
mom_base$mu_matrix[,10] + mom_base$mu_matrix[,11]
mom_base$total_IT_cum <- NA
for(i in 1:length(mom_base$total_IT)){
mom_base$total_IT_cum[i] <- sum(mom_base$total_IT[1:i])
}
#sum up all
mom_base$total <- mom_base$total_CH + mom_base$total_IT
mom_base$total_cum <- NA
for(i in 1:length(mom_base$total)){
mom_base$total_cum[i] <- sum(mom_base$total[1:i])
}
colors <- c("#87a7ca", "#83b695")
#CREATE DF FOR GGPLOT: new cases
#TOTAL
df_total <- data.frame(date = sts_dates[9:length(sts_dates)],
difference = mom_B$total - mom_base$total,
RI = (mom_B$total-mom_base$total)/mom_base$total,
Scenario = "B")
df_total <- rbind(df_total, data.frame(date = sts_dates[9:length(sts_dates)],
difference= mom_A$total - mom_base$total,
RI = (mom_A$total-mom_base$total)/mom_base$total,
Scenario = "A"))
p1 <- ggplot(data = df_total, aes(x=date, y=difference, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
p2 <- ggplot(data = df_total, aes(x=date, y=RI, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
pdf("Plot_RI_new_cases_TOT.pdf", width = 12, height = 8)
grid.arrange(p1 , p2, ncol = 1,
top = "Change in daily new cases\nfor counterfactual scenarios compared to baseline\n(TOTAL)")
dev.off()
#CH
df_CH <- data.frame(date = sts_dates[9:length(sts_dates)],
difference = mom_B$total_CH - mom_base$total_CH,
RI = (mom_B$total_CH-mom_base$total_CH)/mom_base$total_CH,
Scenario = "B")
df_CH <- rbind(df_CH, data.frame(date = sts_dates[9:length(sts_dates)],
difference= mom_A$total_CH - mom_base$total_CH,
RI = (mom_A$total_CH-mom_base$total_CH)/mom_base$total_CH,
Scenario = "A"))
p3 <- ggplot(data = df_CH, aes(x=date, y=difference, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
p4 <- ggplot(data = df_CH, aes(x=date, y=RI, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
pdf("Plot_RI_new_cases_CH.pdf", width = 12, height = 8)
grid.arrange(p3, p4, ncol = 1,
top = "Change in daily new cases\nfor counterfactual scenarios compared to baseline\n(Switzerland)")
dev.off()
#IT
df_IT <- data.frame(date = sts_dates[9:length(sts_dates)],
difference = mom_B$total_IT - mom_base$total_IT,
RI = (mom_B$total_IT - mom_base$total_IT)/mom_base$total_IT,
Scenario = "B")
df_IT <- rbind(df_IT, data.frame(date = sts_dates[9:length(sts_dates)],
difference= mom_A$total_IT - mom_base$total_IT,
RI = (mom_A$total_IT - mom_base$total_IT)/mom_base$total_IT,
Scenario = "A"))
p5 <- ggplot(data = df_IT, aes(x=date, y=difference, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
p6 <- ggplot(data = df_IT, aes(x=date, y=RI, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
pdf("Plot_RI_new_cases_IT.pdf", width = 12, height = 8)
grid.arrange(p5, p6, ncol = 1,
top = "Change in daily new cases\nfor counterfactual scenarios compared to baseline\n(Italy)")
dev.off()
#CREATE DF FOR GGPLOT: cumulative cases
#TOTAL
df_total_cum <- data.frame(date = sts_dates[9:length(sts_dates)],
difference = mom_B$total_cum - mom_base$total_cum,
RI = (mom_B$total_cum-mom_base$total_cum)/mom_base$total_cum,
Scenario = "B")
df_total_cum <- rbind(df_total_cum, data.frame(date = sts_dates[9:length(sts_dates)],
difference= mom_A$total_cum - mom_base$total_cum,
RI = (mom_A$total_cum-mom_base$total_cum)/mom_base$total_cum,
Scenario = "A"))
p1_cum <- ggplot(data = df_total_cum, aes(x=date, y=difference, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
p2_cum <- ggplot(data = df_total_cum, aes(x=date, y=RI, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
pdf("Plot_RI_cum_cases_TOT.pdf", width = 12, height = 8)
grid.arrange(p1_cum, p2_cum, ncol = 1,
top = "Change in cumulative cases\nfor counterfactual scenarios compared to baseline\n(TOTAL)")
dev.off()
#CH
df_CH_cum <- data.frame(date = sts_dates[9:length(sts_dates)],
difference = mom_B$total_CH_cum - mom_base$total_CH_cum,
RI = ( mom_B$total_CH_cum-mom_base$total_CH_cum)/mom_base$total_CH_cum,
Scenario = "B")
df_CH_cum <- rbind(df_CH_cum, data.frame(date = sts_dates[9:length(sts_dates)],
difference= mom_A$total_CH_cum - mom_base$total_CH_cum,
RI = (mom_A$total_CH_cum-mom_base$total_CH_cum)/mom_base$total_CH_cum,
Scenario = "A"))
p3_cum <- ggplot(data = df_CH_cum, aes(x=date, y=difference, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
p4_cum <- ggplot(data = df_CH_cum, aes(x=date, y=RI, group = Scenario))+
geom_line(aes(color=Scenario)) +
scale_color_manual(values=colors)
pdf("Plot_RI_cum_cases_CH.pdf", width = 12, height = 8)
grid.arrange(p3_cum, p4_cum, ncol = 1,
top = "Change in cumulative cases\nfor counterfactual scenarios compared to baseline\n(Switzerland)")