rm(list = ls(all.names = TRUE)) library(Amelia) library(bestNormalize) library(openxlsx) library(lubridate) library(purrr) library(dplyr) library(nlme) library(rsq) library(lm.beta) library(tidyverse) library(caret) library(leaps) library(ggplot2) library(reghelper) library(ggpubr) library(car) library(missMDA) library(psych) library(knitr) setwd("C:/Users/user/Dropbox/BCG-COVID-19/Ver 3/TF_SG") CV19 <- read.xlsx("Table_S1-9.xlsx") DATA <- CV19[, which (colnames(CV19) %in% c( "Country_code", "BCG.TB", "cases_per_million_UP", "deaths_per_million_UP", #"Positive.cases.per.test", "temp", "BCG.type","Current.BCG.vaccination", "prevalence_overweight_2012", "Fraction_urban_population", "Population_ages_65_and_above_2018_percent", "TB.incidents.per.10.thosand.mean_2000_2018", "TB.Incidents.100"))] DATA_COMPLETE <- na.omit(DATA) dim(DATA_COMPLETE) DATA_COMPLETE$cpm_N <- predict(orderNorm(DATA_COMPLETE$cases_per_million_UP)) hist(DATA_COMPLETE$cpm_N) DATA_COMPLETE$cpm_Z <- scale(DATA_COMPLETE$cases_per_million_UP)[,1] hist(DATA_COMPLETE$cpm_Z) DATA_COMPLETE$dpm_N <- predict(orderNorm(DATA_COMPLETE$deaths_per_million_UP)) hist(DATA_COMPLETE$dpm_N) DATA_COMPLETE$dpm_Z <- scale(DATA_COMPLETE$deaths_per_million_UP)[,1] hist(DATA_COMPLETE$dpm_Z) DATA_COMPLETE$TB_N <- predict(orderNorm(DATA_COMPLETE$TB.incidents.per.10.thosand.mean_2000_2018)) hist(DATA_COMPLETE$TB_N) DATA_COMPLETE$TB_Z <- scale(DATA_COMPLETE$TB.incidents.per.10.thosand.mean_2000_2018)[,1] DATA_COMPLETE$Ob_N <- predict(orderNorm(DATA_COMPLETE$prevalence_overweight_2012)) hist(DATA_COMPLETE$Ob_N) DATA_COMPLETE$temp_N <- predict(orderNorm(DATA_COMPLETE$temp)) hist(DATA_COMPLETE$temp_N) DATA_COMPLETE$frac_N <- predict(orderNorm(DATA_COMPLETE$Fraction_urban_population)) hist(DATA_COMPLETE$frac_N) DATA_COMPLETE$age_N <- predict(orderNorm(DATA_COMPLETE$Population_ages_65_and_above_2018_percent)) hist(DATA_COMPLETE$age_N) #DATA_COMPLETE$positives_N <- predict(orderNorm(DATA_COMPLETE$Positive.cases.per.test)) #hist(DATA_COMPLETE$positives_N) DATA_COMPLETE$Yy <- ifelse(DATA_COMPLETE$BCG.TB == "Y_y", 1, 0) DATA_COMPLETE$Yn <- ifelse(DATA_COMPLETE$BCG.TB == "Y_n", 1, 0) DATA_COMPLETE$Nn <- ifelse(DATA_COMPLETE$BCG.TB == "N_n", 1, 0) ###Analysis### #cases cases_lm <- lm(cpm_N ~ Current.BCG.vaccination + TB_Z + Ob_N + temp_N + age_N + frac_N , data = DATA_COMPLETE) summary(cases_lm, corr=T) shapiro.test(resid(cases_lm)) #List of residuals plot(density(resid(cases_lm))) #A density plot qqnorm(rstandard(cases_lm)) qqline(rstandard(cases_lm)) rsq.partial(cases_lm, adj=TRUE) lm.beta(cases_lm) Anova(cases_lm) #selection of predictors models_cases <- regsubsets(cpm_N ~ Current.BCG.vaccination + TB_Z + Ob_N + temp_N + age_N + frac_N , data = DATA_COMPLETE, nvmax = 9) summary(models_cases) res.sum <- summary(models_cases) data.frame( Adj.R2 = which.max(res.sum$adjr2), CP = which.min(res.sum$cp), BIC = which.min(res.sum$bic) ) coef(models_cases, 1:5) vcov(models_cases, 5) cases_lm_4 <- lm(cpm_N ~ Current.BCG.vaccination + TB_Z + age_N + frac_N, data = DATA_COMPLETE) summary(cases_lm_4, corr=T) rsq.partial(cases_lm_4, adj=TRUE) lm.beta(cases_lm_4) shapiro.test(resid(cases_lm_4)) #List of residuals #deaths deaths_lm <- lm(dpm_N ~ Current.BCG.vaccination + TB_Z + Ob_N + temp_N + age_N + frac_N , data = DATA_COMPLETE) summary(deaths_lm, corr=T) rsq.partial(deaths_lm, adj=TRUE) lm.beta(deaths_lm) Anova(deaths_lm) #residual checks shapiro.test(resid(deaths_lm)) #List of residuals plot(density(resid(deaths_lm))) #A density plot #selection of predictors models_deaths <- regsubsets(dpm_N ~ Current.BCG.vaccination + TB_Z + Ob_N + temp_N + age_N + frac_N, data = DATA_COMPLETE, nvmax = 6) summary(models_deaths) res.sum <- summary(models_deaths) data.frame( Adj.R2 = which.max(res.sum$adjr2), CP = which.min(res.sum$cp), BIC = which.min(res.sum$bic) ) coef(models_deaths, 1:5) vcov(models_deaths, 5) deaths_lm_4 <- lm(dpm_N ~ Current.BCG.vaccination + TB_Z + age_N + frac_N, data = DATA_COMPLETE) summary(deaths_lm_4, corr=T) rsq.partial(deaths_lm_4, adj=TRUE) lm.beta(deaths_lm_4) shapiro.test(resid(deaths_lm_4)) #List of residuals ###strains### DATA_str <- CV19[, which (colnames(CV19) %in% c( "Country_code", #"BCG.TB", "cases_per_million_UP", "deaths_per_million_UP", #"Positive.cases.per.test", #"temp", "BCG.type","Current.BCG.vaccination", #"prevalence_overweight_2012", #"Fraction_urban_population", #"Population_ages_65_and_above_2018_percent", #"TB.incidents.per.10.thosand.mean_2000_2018", "Strains"#, #"TB.Incidents.100" ))] DATA_str <- subset(DATA_str, Strains=="G3+" | Strains=="G3-") #DATA_str <- subset(DATA_str, Current.BCG.vaccination=="Y") DATA_str_COMPLETE <- na.omit(DATA_str) DATA_str_COMPLETE$cpm_N <- predict(orderNorm(DATA_str_COMPLETE$cases_per_million_UP)) hist(DATA_str_COMPLETE$cpm_N) DATA_str_COMPLETE$cpm_Z <- scale(DATA_str_COMPLETE$cases_per_million_UP)[,1] hist(DATA_str_COMPLETE$cpm_Z) DATA_str_COMPLETE$dpm_N <- predict(orderNorm(DATA_str_COMPLETE$deaths_per_million_UP)) hist(DATA_str_COMPLETE$dpm_N) DATA_str_COMPLETE$dpm_Z <- scale(DATA_str_COMPLETE$deaths_per_million_UP)[,1] hist(DATA_str_COMPLETE$dpm_Z) #DATA_str_COMPLETE$TB_N <- predict(orderNorm(DATA_str_COMPLETE$TB.incidents.per.10.thosand.mean_2000_2018)) #hist(DATA_str_COMPLETE$TB_N) #DATA_str_COMPLETE$TB_Z <- scale(DATA_str_COMPLETE$TB.incidents.per.10.thosand.mean_2000_2018)[,1] kruskal.test(Strains~Current.BCG.vaccination,data=DATA_str_COMPLETE) chisq.test(DATA_str_COMPLETE$Strains,DATA_str_COMPLETE$Current.BCG.vaccination, correct=FALSE) t.test(cpm_N ~ Strains, data = DATA_str_COMPLETE) t.test(dpm_N ~ Strains, data = DATA_str_COMPLETE) ##### #PCA# cases_PCA_lm <- lm(cpm_N ~ Current.BCG.vaccination + TB_Z + RC1 + RC2 + RC3 , data = CV19) summary(cases_PCA_lm, corr=T) shapiro.test(resid(cases_PCA_lm)) #List of residuals plot(density(resid(cases_PCA_lm))) #A density plot qqnorm(rstandard(cases_PCA_lm)) qqline(rstandard(cases_PCA_lm)) rsq.partial(cases_PCA_lm, adj=TRUE) lm.beta(cases_PCA_lm) kruskal.test(RC1 ~ Current.BCG.vaccination, data = CV19) kruskal.test(RC2 ~ Current.BCG.vaccination, data = CV19) kruskal.test(RC3 ~ Current.BCG.vaccination, data = CV19) deaths_PCA_lm <- lm(dpm_N ~ Current.BCG.vaccination + TB_Z + RC1 + RC2 + RC3 , data = CV19) summary(deaths_PCA_lm, corr=T) shapiro.test(resid(deaths_PCA_lm)) #List of residuals plot(density(resid(deaths_PCA_lm))) #A density plot rsq.partial(deaths_PCA_lm, adj=TRUE) lm.beta(deaths_PCA_lm) ####PLOTS#### c = ggplot(DATA_COMPLETE, aes(x=Current.BCG.vaccination, y=cpm_N, fill=Current.BCG.vaccination)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) #c = c + geom_point() c = c + scale_fill_manual(name="BCG policy", breaks = c("N", "Y"), values=c("#00BA38","#F8766D")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("") c = c + ggtitle("Cases by BCG policy") c = c + theme(legend.position="bottom") #c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_cases_boxplot_Figure_1a.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_COMPLETE, aes(x=TB.Incidents.100, y=cpm_N, fill=TB.Incidents.100)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) #c = c + geom_point() c = c + scale_fill_manual(name="TB incidence 100", breaks = c("n", "y"), values=c("orange","purple")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("") c = c + ggtitle("Cases by TB incidence (100)") c = c + theme(legend.position="bottom") #c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_cases_boxplot_Figure_1e.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_COMPLETE, aes(x=Current.BCG.vaccination, y=dpm_N, fill=Current.BCG.vaccination)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) #c = c + geom_point() c = c + scale_fill_manual(name="BCG policy", breaks = c("N", "Y"), values=c("#00BA38","#F8766D")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("") c = c + ggtitle("Deaths by BCG policy") c = c + theme(legend.position="bottom") #c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_deaths_boxplot_Figure_1b.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_COMPLETE, aes(x=TB.Incidents.100, y=dpm_N, fill=TB.Incidents.100)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) #c = c + geom_point() c = c + scale_fill_manual(name="TB incidence 100", breaks = c("n", "y"), values=c("orange","purple")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("") c = c + ggtitle("Deaths by TB incidence (100)") c = c + theme(legend.position="bottom") #c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_deaths_boxplot_Figure_1f.png", plot=c, dpi=300, width=5, height=5) #plots Figures 2 d2 = ggplot(DATA_COMPLETE, aes(x=TB_Z, y=cpm_N, color=BCG.type)) +geom_point() d2 = d2 + stat_smooth(method="lm",size = 1, colour="black") d2 = d2 + scale_color_discrete( name="BCG type", labels = c("A: Currently administering", "B: Used to administer", "C: Never administered")) d2 = d2 + theme(legend.text=element_text(size=14)) d2 = d2 + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("TB incidents per 100K (scaled)") d2 = d2 + ggtitle("COVID-19 cases by TB incidents") d2 = d2 + geom_vline(xintercept=-0.25, linetype="dotted",size = 1) d2 = d2 + theme(text = element_text(size = 14)) d2 ggsave(file="scatterCOVID19_cases_TB_Figure_2a.png", plot=d2, dpi=300, width=7, height=5) d2 = ggplot(DATA_COMPLETE, aes(x=TB_Z, y=dpm_N, color=BCG.type)) +geom_point() d2 = d2 + stat_smooth(method="lm",size = 1, colour="black") d2 = d2 + scale_color_discrete( name="BCG type", labels = c("A: Currently administering", "B: Used to administer", "C: Never administered")) d2 = d2 + theme(legend.text=element_text(size=14)) d2 = d2 + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("TB incidents per 100K (scaled)") d2 = d2 + ggtitle("COVID-19 deaths by TB incidents") d2 = d2 + geom_vline(xintercept=-0.25, linetype="dotted",size = 1) d2 = d2 + theme(text = element_text(size = 14)) d2 ggsave(file="scatterCOVID19_deaths_TB_Figure_2d.png", plot=d2, dpi=300, width=7, height=5) DATA_COMPLETE_TB <- subset(DATA_COMPLETE, TB.Incidents.100=="n") a<-aggregate(TB_Z~Current.BCG.vaccination, DATA_COMPLETE_TB, FUN=mean) SEM<-(aggregate(TB_Z~Current.BCG.vaccination, DATA_COMPLETE_TB, FUN=sd)[,2])/sqrt(length(DATA_COMPLETE_TB[,1])) Tub<-cbind.data.frame(a,SEM) dat_temp <- subset(DATA_COMPLETE, temp>0 & temp<15) a<-aggregate(temp~Current.BCG.vaccination, dat_temp, FUN=mean) SEM<-(aggregate(temp~Current.BCG.vaccination, dat_temp, FUN=sd)[,2])#/sqrt(length(dat_temp[,1])) temp<-cbind.data.frame(a,SEM) c = ggplot(DATA_COMPLETE_TB, aes(x=Current.BCG.vaccination, y=cpm_N, fill=Current.BCG.vaccination)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) #c = c + geom_point() c = c + scale_fill_manual(name="BCG policy", breaks = c("N", "Y"), values=c("#00BA38","#F8766D")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("") c = c + ggtitle("Cases by BCG policy (filtered TB<100)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_cases_boxplot_Figure_2b.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_COMPLETE_TB, aes(x=Current.BCG.vaccination, y=dpm_N, fill=Current.BCG.vaccination)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) #c = c + geom_point() c = c + scale_fill_manual(name="BCG policy", breaks = c("N", "Y"), values=c("#00BA38","#F8766D")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("") c = c + ggtitle("Deaths by BCG policy (filtered TB<100)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_deaths_boxplot_Figure_2e.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_str_COMPLETE, aes(x=Strains, y=cpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("") c = c + ggtitle("COVID-19 cases by BCG vaccine Strains") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_cases_strains_Figure_1c.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_str_COMPLETE, aes(x=Strains, y=dpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("") c = c + ggtitle("COVID-19 deaths by BCG vaccine Strains") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_deaths_strains_Figure_1d.png", plot=c, dpi=300, width=5, height=5) DATA_str_fig <- CV19[, which (colnames(CV19) %in% c( "Country_code", #"BCG.TB", "cases_per_million_UP", "deaths_per_million_UP", #"Positive.cases.per.test", #"temp", "BCG.type","Current.BCG.vaccination", #"prevalence_overweight_2012", #"Fraction_urban_population", #"Population_ages_65_and_above_2018_percent", "TB.incidents.per.10.thosand.mean_2000_2018", "Strains", "TB.Incidents.100" ))] DATA_str_fig <- subset(DATA_str_fig, Strains=="G3+" | Strains=="G3-") DATA_str_fig <- subset(DATA_str_fig, TB.Incidents.100=="n") DATA_str_COMPLETE_fig <- na.omit(DATA_str_fig) DATA_str_COMPLETE_fig$cpm_N <- predict(orderNorm(DATA_str_COMPLETE_fig$cases_per_million_UP)) DATA_str_COMPLETE_fig$dpm_N <- predict(orderNorm(DATA_str_COMPLETE_fig$deaths_per_million_UP)) c = ggplot(DATA_str_COMPLETE_fig, aes(x=Strains, y=cpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("") c = c + ggtitle("BCG vaccine Strains (filtered TB<100)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_cases_strains_Figure_2c.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_str_COMPLETE_fig, aes(x=Strains, y=dpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("") c = c + ggtitle("BCG vaccine Strains (filtered TB<100)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_deaths_strains_Figure_2f.png", plot=c, dpi=300, width=5, height=5) DATA_str_fig65 <- CV19[, which (colnames(CV19) %in% c( "Country_code", #"BCG.TB", "cases_per_million_UP", "deaths_per_million_UP", #"Positive.cases.per.test", #"temp", "BCG.type","Current.BCG.vaccination", #"prevalence_overweight_2012", #"Fraction_urban_population", "Population_ages_65_and_above_2018_percent", #"TB.incidents.per.10.thosand.mean_2000_2018", "Strains"#, #"TB.Incidents.100" ))] DATA_str_fig65 <- subset(DATA_str_fig65, Strains=="G3+" | Strains=="G3-") DATA_str_fig65 <- subset(DATA_str_fig65, Population_ages_65_and_above_2018_percent>15) DATA_str_COMPLETE_fig65 <- na.omit(DATA_str_fig65) DATA_str_COMPLETE_fig65$cpm_N <- predict(orderNorm(DATA_str_COMPLETE_fig65$cases_per_million_UP)) DATA_str_COMPLETE_fig65$dpm_N <- predict(orderNorm(DATA_str_COMPLETE_fig65$deaths_per_million_UP)) c = ggplot(DATA_str_COMPLETE_fig65, aes(x=Strains, y=cpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("") c = c + ggtitle("BCG vaccine Strains (>65 years > 15%)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_cases_strains_age_Figure_3e.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_str_COMPLETE_fig65, aes(x=Strains, y=dpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("") c = c + ggtitle("BCG vaccine Strains (65 years > 15%)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_deaths_strains_age_Figure_3f.png", plot=c, dpi=300, width=5, height=5) DATA_str_fig_frac <- CV19[, which (colnames(CV19) %in% c( "Country_code", #"BCG.TB", "cases_per_million_UP", "deaths_per_million_UP", #"Positive.cases.per.test", #"temp", "BCG.type","Current.BCG.vaccination", #"prevalence_overweight_2012", "Fraction_urban_population", #"Population_ages_65_and_above_2018_percent", #"TB.incidents.per.10.thosand.mean_2000_2018", "Strains"#, #"TB.Incidents.100" ))] DATA_str_fig_frac <- subset(DATA_str_fig_frac, Strains=="G3+" | Strains=="G3-") DATA_str_fig_frac$Fraction_urban_population_100 <- DATA_str_fig_frac$Fraction_urban_population*100 DATA_str_fig_frac <- subset(DATA_str_fig_frac, Fraction_urban_population_100>75) DATA_str_COMPLETE_fig_frac <- na.omit(DATA_str_fig_frac) DATA_str_COMPLETE_fig_frac$cpm_N <- predict(orderNorm(DATA_str_COMPLETE_fig_frac$cases_per_million_UP)) DATA_str_COMPLETE_fig_frac$dpm_N <- predict(orderNorm(DATA_str_COMPLETE_fig_frac$deaths_per_million_UP)) c = ggplot(DATA_str_COMPLETE_fig_frac, aes(x=Strains, y=cpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 cases per 1M pop.") + xlab("") c = c + ggtitle("BCG vaccine Strains (Urbanization>75%)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_cases_strains_Figure_frac_4c.png", plot=c, dpi=300, width=5, height=5) c = ggplot(DATA_str_COMPLETE_fig_frac, aes(x=Strains, y=dpm_N, fill=Strains)) c = c + geom_boxplot() c = c + theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y = element_text( size = 14)) c = c + geom_jitter(shape=17, position=position_jitter(0.1)) c = c + scale_fill_manual(name="Strain", breaks = c("G3-", "G3+"), values=c("blue","yellow")) c = c + theme(legend.text=element_text(size=14)) c = c + ylab("Standardized COVID-19 deaths per 1M pop.") + xlab("") c = c + ggtitle("BCG vaccine Strains (Urbanization>75%)") c = c + theme(legend.position="bottom") c = c + stat_compare_means(method = "t.test", label.x = 1.3, label.y = 2.7, size=5) c = c + theme(text = element_text(size = 14)) c ggsave(file="COVID19_deaths_strains_frac_Figure_4f.png", plot=c, dpi=300, width=5, height=5)