## Call or install relevant packages ## library(foreign) library(dplyr) library(ggplot2) library(scales) library(ggplot2) library(feathers) library (lubridate) library (tidyverse) library(plotly) library(readstata13) library(openxlsx) library(data.table) library(tidyverse) library(fs) library(readxl) library(stats) library(feathers) library(tibble) ## Set working directory ## setwd("C://Users//suare//OneDrive//Escritorio//Publications 2023//3. Mortality VN//3. Data Compass") rm(list=ls(all=TRUE)) #### Objective 1. Neonatal mortality risk by birth weight categories #### data<- read.xlsx("Mortality_National_Newborn_Types_Data.xlsx", 1 ) names(data) ## Calculate PAR ## data_test<- data %>% group_by(Country) %>% mutate(PAR_Below_1000g=((Prevalence_Below_1000g*(RR_Below_1000g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100), PAR_1000_1500g=((Prevalence_1000_1500g*(RR_1000_1500g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100), PAR_1500_2000g=((Prevalence_1500_2000g*(RR_1500_2000g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+ Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100), PAR_2000_2500g=((Prevalence_2000_2500g*(RR_2000_2500g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+ Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100), PAR_2500_4000g=((Prevalence_2500_4000g*(RR_2500_4000g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+ Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100), PAR_4000_4500g=((Prevalence_4000_4500g*(RR_4000_4500g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+ Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100), PAR_4500_5000g=((Prevalence_4500_5000g*(RR_4500_5000g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100), PAR_Above_5000g=((Prevalence_Above_5000g*(RR_Above_5000g-1))/ (Prevalence_Below_1000g*RR_Below_1000g+Prevalence_1000_1500g*RR_1000_1500g+ Prevalence_1500_2000g*RR_1500_2000g+ Prevalence_2000_2500g*RR_2000_2500g+ Prevalence_2500_4000g*RR_2500_4000g+Prevalence_4000_4500g*RR_4000_4500g+ Prevalence_4500_5000g*RR_4500_5000g+Prevalence_Above_5000g*RR_Above_5000g)*100)) names(data_test) #### select relevant columns for box plots #### data_2<- data_test %>% group_by(Country) %>% dplyr::select(c(RR_Below_1000g, RR_1000_1500g, RR_1500_2000g, RR_2000_2500g, RR_2500_4000g, RR_4000_4500g, RR_4500_5000g, RR_Above_5000g)) colnames(data_2)<- c("Country", "<1000 g", "1000 to 1500 g", "1500 to 2000 g", "2000 to 2500 g", "2500 to 4000 g", "4000 to 4500 g", "4500 to 5000 g", ">5000 g") data_3 <- melt(data_2, id.vars ="Country", variable.name = "Birthweight categories", value.name = "Mortality risk ratio", variable.factor = TRUE) names(data_3) ### g2 is ok ##### Fig2a<-ggplot(data_3, aes(x=`Birthweight categories`, y=`Mortality risk ratio`, color=`Birthweight categories`)) + geom_boxplot() + stat_summary(fun=median, geom="point", col="black")+ stat_summary(fun=median, geom="text", col="black", vjust=1.5, aes(label=paste("", round(..y..,digits=1 )))) Fig2a+scale_color_manual(values = c(">5000 g"="#9E8960", "4500 to 5000 g"="#DC8C3C", "4000 to 4500 g"="#FCC57C", "2500 to 4000 g"="#92AA70", "2000 to 2500 g"="#570951", "1500 to 2000 g"="#FC08C8", "1000 to 1500 g" = "#FF66FF", "<1000 g" ="#FFA3FF" ))+ theme(panel.background = element_blank()) + theme(axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'), axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid')) #### Objective 2. Neonatal mortality risk by gestational age#### data<- read.xlsx("Mortality_National_Newborn_Types_Data.xlsx", 2 ) names(data) data_test<- data %>% group_by(Country) %>% mutate(PAR_Below_28=((Prevalence_Below_28*(RR_Below_28-1))/ (Prevalence_37_42*RR_37_42+Prevalence_Below_28*RR_Below_28+ Prevalence_28_31*RR_28_31+Prevalence_32_33*RR_32_33+ Prevalence_34_36*RR_34_36+Prevalence_Above_42*RR_Above_42)*100), PAR_28_31=((Prevalence_28_31*(RR_28_31-1))/ (Prevalence_37_42*RR_37_42+Prevalence_Below_28*RR_Below_28+ Prevalence_28_31*RR_28_31+Prevalence_32_33*RR_32_33+ Prevalence_34_36*RR_34_36+Prevalence_Above_42*RR_Above_42)*100), PAR_32_33=((Prevalence_32_33*(RR_32_33-1))/ (Prevalence_37_42*RR_37_42+Prevalence_Below_28*RR_Below_28+ Prevalence_28_31*RR_28_31+Prevalence_32_33*RR_32_33+ Prevalence_34_36*RR_34_36+Prevalence_Above_42*RR_Above_42)*100), PAR_34_36=((Prevalence_34_36*(RR_34_36-1))/ (Prevalence_37_42*RR_37_42+Prevalence_Below_28*RR_Below_28+ Prevalence_28_31*RR_28_31+Prevalence_32_33*RR_32_33+ Prevalence_34_36*RR_34_36+Prevalence_Above_42*RR_Above_42)*100), PAR_37_42=((Prevalence_37_42*(RR_37_42-1))/ (Prevalence_37_42*RR_37_42+Prevalence_Below_28*RR_Below_28+ Prevalence_28_31*RR_28_31+Prevalence_32_33*RR_32_33+ Prevalence_34_36*RR_34_36+Prevalence_Above_42*RR_Above_42)*100), PAR_Above_42=((Prevalence_Above_42*(RR_Above_42-1))/ (Prevalence_37_42*RR_37_42+Prevalence_Below_28*RR_Below_28+ Prevalence_28_31*RR_28_31+Prevalence_32_33*RR_32_33+ Prevalence_34_36*RR_34_36+Prevalence_Above_42*RR_Above_42)*100)) ### select relevant columns for box plots ### data_2<- data_test %>% group_by(Country) %>% dplyr::select(c(RR_Below_28, RR_28_31, RR_32_33, RR_34_36, RR_37_42, RR_Above_42)) colnames(data_2)<- c("Country", "<28", "28-31", "32-33", "34-36", "37-42", ">42") ### to large ### data_3 <- melt(data_2, id.vars ="Country", variable.name = "GA categories", value.name = "Mortality risk ratio", variable.factor = TRUE) Fig2b<-ggplot(data_3, aes(x=`GA categories`, y=`Mortality risk ratio`, color=`GA categories`)) + geom_boxplot() + stat_summary(fun=median, geom="point", col="black")+ stat_summary(fun=median, geom="text", col="black", vjust=1.5, aes(label=paste("", round(..y..,digits=1 )))) Fig2b +scale_color_manual(values = c(">42"="#FCC57C", "37-42"="#92AA70", "34-36"="#570951", "32-33"="#FC08C8", "28-31" = "#FF66FF", "<28" ="#FFA3FF" ))+ theme(panel.background = element_blank()) + theme(axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'), axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid')) #### Objective 3. Neonatal mortality risk by newborn types #### data<- read.xlsx("Mortality_National_Newborn_Types_Data.xlsx", 3 ) data_test<- data %>% group_by(Country) %>% mutate(PAR_PT_SGA=((Prevalence_PT_SGA*(RR_PT_SGA-1))/ (Prevalence_T_AGA*RR_T_AGA+Prevalence_PT_SGA*RR_PT_SGA+ Prevalence_PT_AGA*RR_PT_AGA+Prevalence_PT_LGA*RR_PT_LGA+ Prevalence_T_SGA*RR_T_SGA+Prevalence_T_LGA*RR_T_LGA)*100), PAR_PT_AGA=((Prevalence_PT_AGA*(RR_PT_AGA-1))/ (Prevalence_T_AGA*RR_T_AGA+Prevalence_PT_SGA*RR_PT_SGA+ Prevalence_PT_AGA*RR_PT_AGA+Prevalence_PT_LGA*RR_PT_LGA+ Prevalence_T_SGA*RR_T_SGA+Prevalence_T_LGA*RR_T_LGA)*100), PAR_PT_LGA=((Prevalence_PT_LGA*(RR_PT_LGA-1))/ (Prevalence_T_AGA*RR_T_AGA+Prevalence_PT_SGA*RR_PT_SGA+ Prevalence_PT_AGA*RR_PT_AGA+Prevalence_PT_LGA*RR_PT_LGA+ Prevalence_T_SGA*RR_T_SGA+Prevalence_T_LGA*RR_T_LGA)*100), PAR_T_SGA=((Prevalence_T_SGA*(RR_T_SGA-1))/ (Prevalence_T_AGA*RR_T_AGA+Prevalence_PT_SGA*RR_PT_SGA+ Prevalence_PT_AGA*RR_PT_AGA+Prevalence_PT_LGA*RR_PT_LGA+ Prevalence_T_SGA*RR_T_SGA+Prevalence_T_LGA*RR_T_LGA)*100), PAR_T_LGA=((Prevalence_T_LGA*(RR_T_LGA-1))/ (Prevalence_T_AGA*RR_T_AGA+Prevalence_PT_SGA*RR_PT_SGA+ Prevalence_PT_AGA*RR_PT_AGA+Prevalence_PT_LGA*RR_PT_LGA+ Prevalence_T_SGA*RR_T_SGA+Prevalence_T_LGA*RR_T_LGA)*100), PAR_T_AGA=((Prevalence_T_AGA*(RR_T_AGA-1))/ (Prevalence_T_AGA*RR_T_AGA+Prevalence_PT_SGA*RR_PT_SGA+ Prevalence_PT_AGA*RR_PT_AGA+Prevalence_PT_LGA*RR_PT_LGA+ Prevalence_T_SGA*RR_T_SGA+Prevalence_T_LGA*RR_T_LGA)*100)) ### select relevant columns ### data_2<- data_test %>% group_by(Country) %>% dplyr::select(c(RR_PT_SGA, RR_PT_AGA, RR_PT_LGA, RR_T_SGA, RR_T_AGA, RR_T_LGA)) colnames(data_2)<- c("Country", "PT+SGA", "PT+AGA", "PT+LGA", "T+SGA", "T+AGA", "T+LGA") ### to large ### data_3 <- melt(data_2, id.vars ="Country", variable.name = "Newborn types", value.name = "Mortality risk ratio", variable.factor = TRUE) names(data_3) ### g2 is ok ##### Fig3<-ggplot(data_3, aes(x=`Newborn types`, y=`Mortality risk ratio`, color=`Newborn types`)) + geom_boxplot() + stat_summary(fun=median, geom="point", col="black")+ stat_summary(fun=median, geom="text", col="black", vjust=1.5, aes(label=paste("", round(..y..,digits=1 )))) Fig3 +scale_color_manual(values = c("T+LGA"="#9E8960", "T+AGA"="#92AA70", "T+SGA"="#E29F5C", "PT+LGA"="#570951", "PT+AGA"="#FC08C8", "PT+SGA" ="#FFA3FF" ))+ theme(panel.background = element_blank()) + theme(axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'), axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'))