Skip to content

Commit 8ae05ce

Browse files
s-mishras-mishra
authored
Merge pull request #107 from ImperialCollegeLondon/Brazil
Brazil
2 parents 193d1ea + 8bdea04 commit 8ae05ce

16 files changed

+317291
-0
lines changed

‎.gitignore‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,14 @@ Italy/results/*.png
5858
Italy/figures/*.png
5959
Italy/figures/*.pdf
6060
Italy/results/*.RDS
61+
Brazil/code/stan-models/*.rds
62+
Brazil/results/*.Rdata
63+
Brazil/results/*.pdf
64+
Brazil/results/*.png
65+
Brazil/figures/*.png
66+
Brazil/figures/*.pdf
67+
Brazil/figures/*.txt
68+
Brazil/results/*.RDS
6169

6270

6371

‎Brazil/code/plot-3-panel.r‎

Lines changed: 233 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,233 @@
1+
make_data_plot <- function(filename){
2+
load(paste0('Brazil/results/',filename,'-stanfit.Rdata'))
3+
interventions <- read.csv2(paste0("Brazil/data/brazil-interventions.csv"), sep=";")
4+
interventions[,2] <- dmy(as.character(interventions[,2]))
5+
interventions[,3] <- dmy(as.character(interventions[,3]))
6+
interventions[,4] <- dmy(as.character(interventions[,4]))
7+
interventions[,5] <- dmy(as.character(interventions[,5]))
8+
colnames(interventions) = c("region","Emergency","Retail and Service","Transport","School Closing")
9+
table_paper = NULL
10+
for(i in 1:length(countries)){
11+
print(i)
12+
N <- length(dates[[i]])
13+
country <- countries[[i]]
14+
15+
predicted_cases <- colMeans(prediction[,1:N,i])
16+
predicted_cases_li <- colQuantiles(prediction[,1:N,i], probs=.025)
17+
predicted_cases_ui <- colQuantiles(prediction[,1:N,i], probs=.975)
18+
predicted_cases_li2 <- colQuantiles(prediction[,1:N,i], probs=.25)
19+
predicted_cases_ui2 <- colQuantiles(prediction[,1:N,i], probs=.75)
20+
21+
estimated_deaths <- colMeans(estimated.deaths[,1:N,i])
22+
estimated_deaths_li <- colQuantiles(estimated.deaths[,1:N,i], probs=.025)
23+
estimated_deaths_ui <- colQuantiles(estimated.deaths[,1:N,i], probs=.975)
24+
estimated_deaths_li2 <- colQuantiles(estimated.deaths[,1:N,i], probs=.25)
25+
estimated_deaths_ui2 <- colQuantiles(estimated.deaths[,1:N,i], probs=.75)
26+
27+
out <- rstan::extract(fit)
28+
rt <- colMeans(out$Rt_adj[,1:N,i])
29+
rt_li <- colQuantiles(out$Rt_adj[,1:N,i],probs=.025)
30+
rt_ui <- colQuantiles(out$Rt_adj[,1:N,i],probs=.975)
31+
rt_li2 <- colQuantiles(out$Rt_adj[,1:N,i],probs=.25)
32+
rt_ui2 <- colQuantiles(out$Rt_adj[,1:N,i],probs=.75)
33+
34+
data_country <- data.frame("time" = as_date(as.character(dates[[i]])),
35+
"country" = rep(country, length(dates[[i]])),
36+
"reported_cases" = reported_cases[[i]],
37+
"reported_cases_c" = cumsum(reported_cases[[i]]),
38+
"predicted_cases_c" = cumsum(predicted_cases),
39+
"predicted_min_c" = cumsum(predicted_cases_li),
40+
"predicted_max_c" = cumsum(predicted_cases_ui),
41+
"predicted_cases" = predicted_cases,
42+
"predicted_min" = predicted_cases_li,
43+
"predicted_max" = predicted_cases_ui,
44+
"predicted_min2" = predicted_cases_li2,
45+
"predicted_max2" = predicted_cases_ui2,
46+
"deaths" = deaths_by_country[[i]],
47+
"deaths_c" = cumsum(deaths_by_country[[i]]),
48+
"estimated_deaths_c" = cumsum(estimated_deaths),
49+
"death_min_c" = cumsum(estimated_deaths_li),
50+
"death_max_c"= cumsum(estimated_deaths_ui),
51+
"estimated_deaths" = estimated_deaths,
52+
"death_min" = estimated_deaths_li,
53+
"death_max"= estimated_deaths_ui,
54+
"death_min2" = estimated_deaths_li2,
55+
"death_max2"= estimated_deaths_ui2,
56+
"rt" = rt,
57+
"rt_min" = rt_li,
58+
"rt_max" = rt_ui,
59+
"rt_min2" = rt_li2,
60+
"rt_max2" = rt_ui2)
61+
62+
aux = data_country[,c("reported_cases","predicted_cases","predicted_min2","predicted_max2","deaths","death_min2","death_max2")]
63+
aux2 = data.frame(as.character(country),
64+
tail(apply(aux, 2, cumsum),1),
65+
(df_pop[which(df_pop$region==country),2])
66+
)
67+
table_paper = rbind.data.frame(table_paper,aux2)
68+
data_cases_95 <- data.frame(data_country$time, data_country$predicted_min,
69+
data_country$predicted_max)
70+
names(data_cases_95) <- c("time", "cases_min", "cases_max")
71+
data_cases_95$key <- rep("nintyfive", length(data_cases_95$time))
72+
data_cases_50 <- data.frame(data_country$time, data_country$predicted_min2,
73+
data_country$predicted_max2)
74+
names(data_cases_50) <- c("time", "cases_min", "cases_max")
75+
data_cases_50$key <- rep("fifty", length(data_cases_50$time))
76+
data_cases <- rbind(data_cases_95, data_cases_50)
77+
levels(data_cases$key) <- c("ninetyfive", "fifty")
78+
make_plots(data_country, data_cases, country,filename, interventions)
79+
}
80+
colnames(table_paper) = c("State","reported_cases","predicted_cases","predicted_min2","predicted_max2",
81+
"deaths","death_min2","death_max2","pop")
82+
make_table(table_paper)
83+
}
84+
make_plots <- function(data_country, data_cases, country, filename, interventions){
85+
p1 <- ggplot(data_country) +
86+
geom_bar(data = data_country, aes(x = time, y = reported_cases),
87+
fill = "coral4", stat='identity', alpha=0.5) +
88+
geom_ribbon(data = data_cases,
89+
aes(x = time, ymin = cases_min, ymax = cases_max, fill = key)) +
90+
xlab("") +
91+
ylab("Daily number of infections\n") +
92+
scale_x_date(date_breaks = "2 weeks", labels = date_format("%e %b")) +
93+
scale_y_continuous(expand = c(0, 0), labels = comma) +
94+
scale_fill_manual(name = "", labels = c("50%", "95%"),
95+
values = c(alpha("deepskyblue4", 0.55),
96+
alpha("deepskyblue4", 0.45))) +
97+
theme_pubr() +
98+
theme(axis.text.x = element_text(angle = 45, hjust = 1),
99+
legend.position = "None") + ggtitle(df_region_codes[which(df_region_codes[,1]==country),2]) +
100+
guides(fill=guide_legend(ncol=1))
101+
102+
data_deaths_95 <- data.frame(data_country$time, data_country$death_min,
103+
data_country$death_max)
104+
names(data_deaths_95) <- c("time", "death_min", "death_max")
105+
data_deaths_95$key <- rep("nintyfive", length(data_deaths_95$time))
106+
data_deaths_50 <- data.frame(data_country$time, data_country$death_min2,
107+
data_country$death_max2)
108+
names(data_deaths_50) <- c("time", "death_min", "death_max")
109+
data_deaths_50$key <- rep("fifty", length(data_deaths_50$time))
110+
data_deaths <- rbind(data_deaths_95, data_deaths_50)
111+
levels(data_deaths$key) <- c("ninetyfive", "fifty")+ coord_fixed(ratio = 10)
112+
113+
p2 <- ggplot(data_country, aes(x = time)) +
114+
geom_bar(data = data_country, aes(y = deaths, fill = "reported"),
115+
fill = "coral4", stat='identity', alpha=0.5) +
116+
geom_ribbon(
117+
data = data_deaths,
118+
aes(ymin = death_min, ymax = death_max, fill = key)) +
119+
scale_x_date(date_breaks = "2 weeks", labels = date_format("%e %b")) +
120+
scale_y_continuous(expand = c(0, 0), labels = comma) +
121+
scale_fill_manual(name = "", labels = c("50%", "95%"),
122+
values = c(alpha("deepskyblue4", 0.55),
123+
alpha("deepskyblue4", 0.45))) +
124+
ylab("Daily number of deaths\n") +
125+
xlab("") +
126+
theme_pubr() +
127+
theme(axis.text.x = element_text(angle = 45, hjust = 1),
128+
legend.position = "None") +
129+
guides(fill=guide_legend(ncol=1))
130+
131+
# Plotting interventions
132+
data_rt_95 <- data.frame(data_country$time,
133+
data_country$rt_min, data_country$rt_max)
134+
names(data_rt_95) <- c("time", "rt_min", "rt_max")
135+
data_rt_95$key <- rep("nintyfive", length(data_rt_95$time))
136+
data_rt_50 <- data.frame(data_country$time, data_country$rt_min2,
137+
data_country$rt_max2)
138+
names(data_rt_50) <- c("time", "rt_min", "rt_max")
139+
data_rt_50$key <- rep("fifty", length(data_rt_50$time))
140+
data_rt <- rbind(data_rt_95, data_rt_50)
141+
levels(data_rt$key) <- c("ninetyfive", "fifth")
142+
143+
# interventions
144+
# # delete these 2 lines
145+
covariates_country <- interventions[which(interventions$region == country),-1]
146+
covariates_country_long <- gather(covariates_country, key = "key",
147+
value = "value")
148+
covariates_country_long$x <- rep(NULL, length(covariates_country_long$key))
149+
un_dates <- unique(covariates_country_long$value)
150+
151+
for (k in 1:length(un_dates)){
152+
idxs <- which(covariates_country_long$value == un_dates[k])
153+
max_val <- round(max(data_country$rt_max)) + 0.3
154+
for (j in idxs){
155+
covariates_country_long$x[j] <- max_val
156+
max_val <- max_val - 0.3
157+
}
158+
}
159+
160+
covariates_country_long$value <- as_date(covariates_country_long$value)
161+
covariates_country_long$country <- rep(country,
162+
length(covariates_country_long$value))
163+
164+
plot_labels <- c("Emergency","Retail and Service","School Closing","Transport")
165+
166+
p3 <- ggplot(data_country) +
167+
geom_ribbon(data = data_rt, aes(x = time, ymin = rt_min, ymax = rt_max,
168+
group = key,
169+
fill = key)) +
170+
geom_hline(yintercept = 1, color = 'black', size = 0.1) +
171+
geom_segment(data = covariates_country_long,
172+
aes(x = value, y = 0, xend = value, yend = max(x)),
173+
linetype = "dashed", colour = "grey", alpha = 0.75) +
174+
geom_point(data = covariates_country_long, aes(x = value,
175+
y = x,
176+
group = key,
177+
shape = key,
178+
col = key), size = 2) +
179+
xlab("") +
180+
ylab(expression(R[t])) +
181+
scale_fill_manual(name = "", labels = c("50%", "95%"),
182+
values = c(alpha("seagreen", 0.75), alpha("seagreen", 0.5))) +
183+
scale_shape_manual(name = "Interventions", labels = plot_labels,
184+
values = c(21, 22, 23, 24, 25, 12)) +
185+
scale_colour_discrete(name = "Interventions", labels = plot_labels) +
186+
scale_x_date(date_breaks = "weeks", labels = date_format("%e %b"),
187+
limits = c(data_country$time[1],
188+
data_country$time[length(data_country$time)])) +
189+
scale_y_continuous(expand = expansion(mult=c(0,0.1))) +
190+
theme_pubr() +
191+
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
192+
theme(legend.position="right")
193+
194+
195+
ptmp <- plot_grid(p1, p2, p3, ncol = 3, rel_widths = c(0.75, 0.75, 1))
196+
#print(ptmp)
197+
#save_plot(filename = paste0("Brazil/figures/", country, "_three_pannel_", filename2, ".png"), p, base_width = 14)
198+
#ggsave(ptmp, file=paste0("Brazil/figures/", country, "_three_pannel_", JOBID,'-',filename2, ".png"), width = 14)
199+
ggsave(ptmp,
200+
file=paste0("Brazil/figures/",country,"-three-pannel-", filename, ".png"), width = 14, height = 5)
201+
202+
}
203+
204+
make_table <- function(table_paper){
205+
weight_fatality<-read.csv(paste0("Brazil/data/IFRS-all.csv"))[c("X","State","IFR_Peru_poorer")]
206+
cfr.by.country<-weight_fatality
207+
colnames(cfr.by.country)<-c(" ","region","weighted_fatality")
208+
ifrs=cfr.by.country[c("region","weighted_fatality")]
209+
colnames(ifrs) = c("State","IFR")
210+
table=merge(table_paper,ifrs,x.by="State")
211+
table_f = table[c("State","IFR","pop","deaths","predicted_cases","predicted_min2","predicted_max2")]
212+
table_f$death_per_capita = table_f$deaths/as.double(table$pop)
213+
table_f$attackrate = 100*table_f$predicted_cases/table_f$pop
214+
table_f$attackrate_min2 = 100*table_f$predicted_min2/table_f$pop
215+
table_f$attackrate_max2 = 100*table_f$predicted_max2/table_f$pop
216+
table_f$predicted_cases = signif(table_f$predicted_cases, 3)
217+
table_f$predicted_min2 = signif(table_f$predicted_min2, 3)
218+
table_f$predicted_max2 = signif(table_f$predicted_max2, 3)
219+
table_f$IFR = signif(100*table_f$IFR, 2)
220+
table_f$death_per_capita = round(signif(1000000*table_f$death_per_capita, 3),3)
221+
table_f$attackrate = signif(table_f$attackrate, 3)
222+
table_f$attackrate_min2 = signif(table_f$attackrate_min2, 3)
223+
table_f$attackrate_max2 = signif(table_f$attackrate_max2, 3)
224+
table_f = table_f[order(table_f$deaths,decreasing = T),]
225+
table_paper_final=table_f[c("State","IFR","pop","deaths","death_per_capita",
226+
"predicted_cases","predicted_min2","predicted_max2",
227+
"attackrate","attackrate_min2","attackrate_max2")]
228+
colnames(table_paper_final)=c("State","IFR","Pop.","Deaths","Deaths ppm",
229+
"cases","case5%","case95%",
230+
"AR","AR5%","AR95%")
231+
print(table_paper_final)
232+
write.csv(table_paper_final, file=paste0("Brazil/figures/TABLE_FINAL_counterfactual_",filename, ".txt"),row.names=FALSE)
233+
}

0 commit comments

Comments
 (0)