Skip to content

Commit

Permalink
Update Figure.R
Browse files Browse the repository at this point in the history
  • Loading branch information
@wxwx1993
wxwx1993 committed Sep 19, 2020
1 parent 307cb2b commit e7e5ee8
Showing 1 changed file with 102 additions and 99 deletions.
201 changes: 102 additions & 99 deletions Figure.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,123 +11,127 @@ library(gridExtra)
library("lme4")

# Figure 1 the US PM2.5 and COVID-19 death maps
us<-map_data('state')
us <- map_data('state')

#Plot prevalence of COVID-19
states <- st_read("./cb_2017_us_county_500k/cb_2017_us_county_500k.shp")

aggregate_pm_census_cdc_test_beds$fips = str_pad(aggregate_pm_census_cdc_test_beds$fips, 5, pad = "0")
aggregate_pm_census_cdc_test_beds$fips <- str_pad(aggregate_pm_census_cdc_test_beds$fips, 5, pad = "0")
covid_us <- mutate(aggregate_pm_census_cdc_test_beds,
STATEFP = str_sub(fips, 1, 2),
COUNTYFP = str_sub(fips, 3, 5))
str(covid_us)
str(states)
states$STATEFP=as.character(states$STATEFP)
states$COUNTYFP=as.character(states$COUNTYFP)
states$STATEFP <- as.character(states$STATEFP)
states$COUNTYFP <- as.character(states$COUNTYFP)
statesCOVID <- left_join(states, covid_us, by = c("STATEFP", "COUNTYFP"))
statesCOVID$mortality = statesCOVID$Deaths/statesCOVID$population*10^6
statesCOVID$logmortality = log10(statesCOVID$Deaths/statesCOVID$population*10^6)
statesCOVID$logmortality[statesCOVID$logmortality < 0] = -1
statesCOVID$logmortality[is.na(statesCOVID$logmortality)] = -1

g1<-ggplot(statesCOVID)+
xlim(-125,-65)+ylim(25,50)+
statesCOVID$mortality <- statesCOVID$Deaths/statesCOVID$population * 10^6
statesCOVID$logmortality <- log10(statesCOVID$Deaths / statesCOVID$population * 10^6)
statesCOVID$logmortality[statesCOVID$logmortality < 0] <- (-1)
statesCOVID$logmortality[is.na(statesCOVID$logmortality)] <- (-1)

g1 <- ggplot(statesCOVID) +
xlim(-125, -65) +
ylim(25, 50) +
# geom_sf(aes(fill = PD_p),color=NA,size=0.025)+
geom_sf(aes(fill = logmortality),color='grey',size=0.005)+
geom_sf(aes(fill = logmortality), color = 'grey', size = 0.005) +
# scale_fill_viridis_c(option="magma",begin=0.4)+
scale_fill_gradient2(expression(paste("# COVID-19 deaths per 1 million")),low = "#1e90ff", mid="#ffffba", high = "#8b0000",midpoint = 1,
breaks = c(-1,0,1,2,3),
labels = c("0","1","10","100","1000+"),limits = c(-1,3.4) , na.value = "white") +
scale_fill_gradient2(expression(paste("# COVID-19 deaths per 1 million")),
low = "#1e90ff",
mid = "#ffffba",
high = "#8b0000",
midpoint = 1,
breaks = c(-1, 0, 1, 2, 3),
labels = c("0", "1", "10", "100", "1000+"),
limits = c(-1, 3.4),
na.value = "white") +
# labs(title = expression(paste("Cumulative Deaths Related to COVID-19 until March 30, 2020"))) +
theme_minimal() +
theme(plot.title = element_text(size = 24*2,hjust = 0.5),
theme(plot.title = element_text(size = 24 * 2, hjust = 0.5),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
line = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.direction = "horizontal",
legend.text = element_text(angle = 60, size = 20*2),
legend.text = element_text(angle = 60, size = 20 * 2),
legend.text.align = 0.75,
legend.title = element_text(size = 18*2),
legend.key.width = unit(150*2, "points"),
legend.title = element_text(size = 18 * 2),
legend.key.width = unit(150 * 2, "points"),
panel.grid.major = element_line(colour = "transparent"))

png("county_covid.jpeg", height = 1024*0.6*2, width = 1024*2)
png("county_covid.jpeg", height = 1024 * 0.6 * 2, width = 1024 * 2)
g1
dev.off()

county_pm_aggregated = county_pm %>%
county_pm_aggregated <- county_pm %>%
group_by(fips) %>%
summarise(mean_pm25 = mean(pm25))

county_pm_aggregated$fips = str_pad(county_pm_aggregated$fips, 5, pad = "0")
county_pm_aggregated$fips <- str_pad(county_pm_aggregated$fips, 5, pad = "0")
county_pm_aggregated <- mutate(county_pm_aggregated,
STATEFP = str_sub(fips, 1, 2),
COUNTYFP = str_sub(fips, 3, 5))
str(county_pm_aggregated)
str(states)
states$STATEFP=as.character(states$STATEFP)
states$COUNTYFP=as.character(states$COUNTYFP)
statesPM <- left_join(states, county_pm_aggregated, by = c("STATEFP", "COUNTYFP"))

g2<-ggplot(statesPM)+
xlim(-125,-65)+ylim(25,50)+
geom_sf(aes(fill = mean_pm25),color='grey',size=0.005)+
states$STATEFP <- as.character(states$STATEFP)
states$COUNTYFP <- as.character(states$COUNTYFP)
statesPM <- left_join(states,
county_pm_aggregated,
by = c("STATEFP", "COUNTYFP"))

g2 <- ggplot(statesPM) +
xlim(-125, -65) +
ylim(25, 50) +
geom_sf(aes(fill = mean_pm25), color='grey', size = 0.005) +
# scale_fill_viridis_c(option="magma",begin=0.4)+
# scale_fill_viridis_c(option="magma",begin=0,direction = -1, breaks = c(0,4,8,12,16,20))+
scale_fill_gradient2(expression(paste("PM"[2.5])),low = "#1e90ff", mid="#ffffba", high = "#8b0000",midpoint = 9,
breaks = c(0,3,6,9,12),
labels = c("0","3","6","9","12+"),limits = c(0,15) , na.value = "white") +
scale_fill_gradient2(expression(paste("PM"[2.5])),
low = "#1e90ff",
mid = "#ffffba",
high = "#8b0000",
midpoint = 9,
breaks = c(0, 3, 6, 9, 12),
labels = c("0", "3", "6", "9", "12+"),
limits = c(0, 15),
na.value = "white") +
# labs(title = expression(paste("Annual Average of PM"[2.5]," per ",mu,g/m^3," in 2000-2016"))) +
theme_minimal() +
theme(plot.title = element_text(size = 24*2,hjust = 0.5),
theme(plot.title = element_text(size = 24 * 2,hjust = 0.5),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
line = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.direction = "horizontal",
legend.text = element_text(angle = 60, size = 20*2),
legend.text = element_text(angle = 60, size = 20 * 2),
legend.text.align = 0.75,
legend.title = element_text(size = 24*2),
legend.key.width = unit(150*2, "points"),
legend.title = element_text(size = 24 * 2),
legend.key.width = unit(150 * 2, "points"),
panel.grid.major = element_line(colour = "transparent"))

png("county_pm.jpeg", height = 1024*0.6*2, width = 1024*2)
png("county_pm.jpeg", height = 1024 * 0.6 * 2, width = 1024 * 2)
g2
dev.off()

# Figure 2 main analyses
data <- data.frame(method = c("Main Analysis", "Omit # Hospital Beds", "Omit # Tested","Omit BRFSS",
"Omit Weather", "Exclude NY State", "Exclude Confirmed < 10"),
RR = exp(c(0.14, 0.06, 0.14, 0.12, 0.10, 0.12, 0.11)
),
lower_CI = exp(c(0.05,-0.01,0.05,0.03,0.03,0.03,0.01)
),
upper_CI = exp(c(0.22,0.13,0.22,0.20,0.18,0.20,0.20)
),
RR = exp(c(0.14, 0.06, 0.14, 0.12, 0.10, 0.12, 0.11)),
lower_CI = exp(c(0.05, -0.01, 0.05, 0.03, 0.03, 0.03, 0.01)),
upper_CI = exp(c(0.22, 0.13, 0.22, 0.20, 0.18, 0.20, 0.20)),
Methods = c(1:7))

bracketsGrob <- function(...){
l <- list(...)
e <- new.env()
e$l <- l
grid:::recordGrob( {
do.call(grid.brackets, l)
}, e)
}

pdf("RR1.pdf",width = 28, height = 11)
pdf("RR1.pdf", width = 28, height = 11)
#jpeg("RR1.jpeg", height = 1024*0.85, width = 1024*2)
p1 <- ggplot(data[1:7,], aes(x=Methods, y=RR),size= 5) +
p1 <- ggplot(data[1:7, ], aes(x = Methods, y = RR), size = 5) +
ylab("Mortality Rate Ratios") +
geom_point(aes(size = 1)) +
geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=.2,size=1) +
geom_hline(yintercept = 1.0, linetype="dashed", size=0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1-12/100, label = data$method, size = 5*2) +
geom_errorbar(aes(ymin = lower_CI, ymax = upper_CI), width = 0.2 , size = 1) +
geom_hline(yintercept = 1.0, linetype = "dashed", size = 0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1 - 12 / 100, label = data$method, size = 5 * 2) +
#annotate(geom = "text", x = c(3), y = 1+21/100, label = c("Entire Medicare Enrollees (2000-2016)"), size = 8) +
#annotate(geom = "text", x = c(c(1.5+5 * (0),4+5 * (0))), y = 1-1.75/100, label = c(rep(c("Regression"),1),rep(c("Causal"),1)), size = 8) +
coord_cartesian(ylim = c(0.90, 1.4), xlim = c(0.5, 7.5), expand = FALSE, clip = "off") +
Expand All @@ -136,21 +140,21 @@ p1 <- ggplot(data[1:7,], aes(x=Methods, y=RR),size= 5) +
theme(plot.margin = unit(c(3, 2, 4, 1), "lines"),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(size=16*2),
axis.title.y = element_text(size = 20*2),
axis.text.y = element_text(size=16 * 2),
axis.title.y = element_text(size = 20 * 2),
legend.position = "none")
p1
dev.off()


# Figure S1
data <- data.frame(method = c(rep(c("Main Analysis", "Omit # Hospital Beds", "Omit # Tested","Omit BRFSS","Omit Weather", "Exclude NY State", "Exclude Confirmed < 10"),4)),
data <- data.frame(method = c(rep(c("Main Analysis", "Omit # Hospital Beds", "Omit # Tested","Omit BRFSS","Omit Weather", "Exclude NY State", "Exclude Confirmed < 10"), 4)),
RR = exp(c(0.14, 0.06, 0.14, 0.12, 0.10, 0.12, 0.11,
0.12, 0.03, 0.12, 0.10, 0.09, 0.12, 0.09,
0.12, 0.05, 0.12, 0.11, 0.08, 0.09, 0.08,
0.15, 0.09, 0.15, 0.15, 0.10, 0.12, 0.13)
),
lower_CI = exp(c(0.05, -0.01, 0.05, 0.03, 0.03, 0.03,0.01,
lower_CI = exp(c(0.05, -0.01, 0.05, 0.03, 0.03, 0.03, 0.01,
0.02, -0.05, 0.02, 0.01, 0.00, 0.02, -0.03,
0.03, -0.02, 0.03, 0.02, 0.00, -0.01, -0.02,
0.03, -0.02, 0.03, 0.03, -0.01, -0.01, -0.01)
Expand All @@ -160,89 +164,88 @@ data <- data.frame(method = c(rep(c("Main Analysis", "Omit # Hospital Beds", "Om
0.28, 0.13, 0.21, 0.19, 0.16, 0.18, 0.19,
0.27, 0.19, 0.27, 0.26, 0.20, 0.24, 0.27)
),
Methods = c(1:7,1:7,1:7,1:7),
Exposure = c(rep("17-year mean (2000-2016) by (1)",7),rep("the most recent (2016) by (1)",7),
rep("17-year mean (2000-2016) by (2)",7),rep("the most recent (2016) by (2)",7)))
Methods = c(1:7, 1:7, 1:7, 1:7),
Exposure = c(rep("17-year mean (2000-2016) by (1)", 7), rep("the most recent (2016) by (1)", 7),
rep("17-year mean (2000-2016) by (2)", 7), rep("the most recent (2016) by (2)", 7)))

pdf("RR_sup.pdf",width = 40, height = 44)
pdf("RR_sup.pdf", width = 40, height = 44)

p1 <- ggplot(data[1:7,], aes(x=Methods, y=RR),size= 5) +
p1 <- ggplot(data[1:7, ], aes(x = Methods, y = RR), size = 5) +
ylab("Mortality Rate Ratios") +
geom_point(aes(size = 1)) +
geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=.2,size=1) +
geom_hline(yintercept = 1.0, linetype="dashed", size=0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1-12/100, label = data$method, size = 6.5*2) +
annotate(geom = "text", x = c(4), y = 1+42/100, label = c(expression(paste("P-1"))), size = 9*2) +
geom_errorbar(aes(ymin = lower_CI, ymax = upper_CI), width = 0.2, size = 1) +
geom_hline(yintercept = 1.0, linetype = "dashed", size = 0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1 - 12 / 100, label = data$method, size = 6.5 * 2) +
annotate(geom = "text", x = c(4), y = 1 + 42 / 100, label = c(expression(paste("P-1"))), size = 9 * 2) +
#annotate(geom = "text", x = c(c(1.5+5 * (0),4+5 * (0))), y = 1-1.75/100, label = c(rep(c("Regression"),1),rep(c("Causal"),1)), size = 8) +
coord_cartesian(ylim = c(0.90, 1.4), xlim = c(0.5, 7.5), expand = FALSE, clip = "off") +
#geom_segment(aes(x = 2.5, y = 1-10/100, xend = 2.5, yend = 1),colour="black") +
theme_bw() +
theme(plot.margin = unit(c(4, 2, 4, 2), "lines"),
plot.title = element_text(hjust=0.5),
plot.title = element_text(hjust = 0.5),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(size=16*2),
axis.title.y = element_text(size = 18*2),
axis.text.y = element_text(size = 16 * 2),
axis.title.y = element_text(size = 18 * 2),
legend.position = "none")


p2 <- ggplot(data[8:14,], aes(x=Methods, y=RR),size= 5) +
p2 <- ggplot(data[8:14, ], aes(x = Methods, y = RR), size = 5) +
ylab("Mortality Rate Ratios") +
geom_point(aes(size = 1)) +
geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=.2,size=1) +
geom_hline(yintercept = 1.0, linetype="dashed", size=0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1-12/100, label = data$method, size = 6.5*2) +
annotate(geom = "text", x = c(4), y = 1+42/100, label = c(expression(paste("P-2"))), size = 9*2) +
geom_errorbar(aes(ymin = lower_CI, ymax = upper_CI), width = 0.2, size = 1) +
geom_hline(yintercept = 1.0, linetype = "dashed", size = 0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1 - 12 / 100, label = data$method, size = 6.5 * 2) +
annotate(geom = "text", x = c(4), y = 1 + 42 / 100, label = c(expression(paste("P-2"))), size = 9 * 2) +
#annotate(geom = "text", x = c(c(1.5+5 * (0),4+5 * (0))), y = 1-1.75/100, label = c(rep(c("Regression"),1),rep(c("Causal"),1)), size = 8) +
coord_cartesian(ylim = c(0.90, 1.4), xlim = c(0.5, 7.5), expand = FALSE, clip = "off") +
#geom_segment(aes(x = 2.5, y = 1-10/100, xend = 2.5, yend = 1),colour="black") +
theme_bw() +
theme(plot.margin = unit(c(4, 2, 4, 2), "lines"),
plot.title = element_text(hjust=0.5),
plot.title = element_text(hjust = 0.5),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(size=16*2),
axis.title.y = element_text(size = 18*2),
axis.text.y = element_text(size = 16 * 2),
axis.title.y = element_text(size = 18 * 2),
legend.position = "none")

p3 <- ggplot(data[15:21,], aes(x=Methods, y=RR),size= 5) +
p3 <- ggplot(data[15:21, ], aes(x = Methods, y = RR), size = 5) +
ylab("Mortality Rate Ratios") +
geom_point(aes(size = 1)) +
geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=.2,size=1) +
geom_hline(yintercept = 1.0, linetype="dashed", size=0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1-12/100, label = data$method, size = 6.5*2) +
annotate(geom = "text", x = c(4), y = 1+42/100, label = c(expression(paste("P-3"))), size = 9*2) +
geom_errorbar(aes(ymin = lower_CI, ymax = upper_CI), width = 0.2, size = 1) +
geom_hline(yintercept = 1.0, linetype = "dashed", size = 0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1-12 / 100, label = data$method, size = 6.5 * 2) +
annotate(geom = "text", x = c(4), y = 1 + 42 / 100, label = c(expression(paste("P-3"))), size = 9 * 2) +
#annotate(geom = "text", x = c(c(1.5+5 * (0),4+5 * (0))), y = 1-1.75/100, label = c(rep(c("Regression"),1),rep(c("Causal"),1)), size = 8) +
coord_cartesian(ylim = c(0.90, 1.4), xlim = c(0.5, 7.5), expand = FALSE, clip = "off") +
#geom_segment(aes(x = 2.5, y = 1-10/100, xend = 2.5, yend = 1),colour="black") +
theme_bw() +
theme(plot.margin = unit(c(4, 2, 4, 2), "lines"),
plot.title = element_text(hjust=0.5),
plot.title = element_text(hjust = 0.5),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(size=16*2),
axis.title.y = element_text(size = 18*2),
axis.text.y = element_text(size = 16 * 2),
axis.title.y = element_text(size = 18 * 2),
legend.position = "none")

p4 <- ggplot(data[22:28,], aes(x=Methods, y=RR),size= 5) +
p4 <- ggplot(data[22:28, ], aes(x = Methods, y = RR), size = 5) +
ylab("Mortality Rate Ratios") +
geom_point(aes(size = 1)) +
geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=.2,size=1) +
geom_hline(yintercept = 1.0, linetype="dashed", size=0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1-12/100, label = data$method, size = 6.5*2) +
annotate(geom = "text", x = c(4), y = 1+42/100, label = c(expression(paste("P-4"))), size = 9*2) +
geom_errorbar(aes(ymin = lower_CI, ymax = upper_CI), width = 0.2, size = 1) +
geom_hline(yintercept = 1.0, linetype = "dashed", size = 0.8) +
annotate(geom = "text", x = seq_len(nrow(data)), y = 1 - 12 / 100, label = data$method, size = 6.5 * 2) +
annotate(geom = "text", x = c(4), y = 1 + 42 / 100, label = c(expression(paste("P-4"))), size = 9 * 2) +
#annotate(geom = "text", x = c(c(1.5+5 * (0),4+5 * (0))), y = 1-1.75/100, label = c(rep(c("Regression"),1),rep(c("Causal"),1)), size = 8) +
coord_cartesian(ylim = c(0.90, 1.4), xlim = c(0.5, 7.5), expand = FALSE, clip = "off") +
#geom_segment(aes(x = 2.5, y = 1-10/100, xend = 2.5, yend = 1),colour="black") +
theme_bw() +
theme(plot.margin = unit(c(4, 2, 4, 2), "lines"),
plot.title = element_text(hjust=0.5),
plot.title = element_text(hjust = 0.5),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(size=16*2),
axis.title.y = element_text(size = 18*2),
axis.text.y = element_text(size = 16 * 2),
axis.title.y = element_text(size = 18 * 2),
legend.position = "none")

grid.arrange(p1, p2,p3,p4, nrow = 4)

grid.arrange(p1, p2, p3, p4, nrow = 4)
dev.off()

0 comments on commit e7e5ee8

Please sign in to comment.