## Call packages ## 

library(foreign)
library(dplyr)
library(ggplot2)
library(scales)
library(ggplot2)
library(palmerpenguins)
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(epitools)
library(meta)


## Set working directory ##

setwd("")

rm(list=ls(all=TRUE))

#### 1. Load and explore data #####

data<- read.xlsx("All Data Final - Objective 2_August_2022.xlsx", 1) 

options(scipen=999) ## turns off scientific notation
head(data)

dim(data) ## 115 344
names(data)

table(data$Country)
length( unique(data$Country)) ## 15 countries 

#### 2. Calculate early and late mortality rates by term types (AGA and LGA) #####

test<- 
  data %>% group_by(Country) %>% 
  dplyr::select(AGA_T_NBW_Livebirths,AGA_T_NBW_Early_neonatal,
                AGA_T_NBW_Late_neonatal, AGA_T_NBW_Neonatal,
                LGA_T_NBW_Livebirths, LGA_T_NBW_Early_neonatal, 
                LGA_T_NBW_Late_neonatal, LGA_T_NBW_Neonatal)%>%
  dplyr::summarise(sum_AGA_T_NBW_Livebirths=sum(AGA_T_NBW_Livebirths),
                   sum_AGA_T_NBW_Early_neonatal=sum(AGA_T_NBW_Early_neonatal), 
                   sum_AGA_T_NBW_Late_neonatal=sum(AGA_T_NBW_Late_neonatal), 
                   Rate_AGA_T_Early=(sum_AGA_T_NBW_Early_neonatal/sum_AGA_T_NBW_Livebirths)*1000,
                   Rate_AGA_T_Late=(sum_AGA_T_NBW_Late_neonatal/(sum_AGA_T_NBW_Livebirths-sum_AGA_T_NBW_Early_neonatal))*1000,
                   sum_LGA_T_NBW_Livebirths=sum(LGA_T_NBW_Livebirths),
                   sum_LGA_T_NBW_Early_neonatal=sum(LGA_T_NBW_Early_neonatal),
                   sum_LGA_T_NBW_Late_neonatal=sum(LGA_T_NBW_Late_neonatal),
                   Rate_LGA_T_Early=(sum_LGA_T_NBW_Early_neonatal/sum_LGA_T_NBW_Livebirths)*1000,
                   Rate_LGA_T_Late=(sum_LGA_T_NBW_Late_neonatal/(sum_LGA_T_NBW_Livebirths-sum_LGA_T_NBW_Early_neonatal))*1000,
                   Risk_Early=Rate_LGA_T_Early/Rate_AGA_T_Early,
                   Risk_Late=Rate_LGA_T_Late/Rate_AGA_T_Late)

test<-as.data.frame(test)

view(test)

#### 3. Calculate risk ratios for early mortality ####


test<- 
  data %>% group_by(Country) %>% 
  dplyr::select(AGA_T_NBW_Livebirths,AGA_T_NBW_Early_neonatal,
                AGA_T_NBW_Late_neonatal, AGA_T_NBW_Neonatal,
                LGA_T_NBW_Livebirths, LGA_T_NBW_Early_neonatal, 
                LGA_T_NBW_Late_neonatal, LGA_T_NBW_Neonatal) %>%
  dplyr::summarise(sum_LGA_T_NBW_Early_neonatal=sum(LGA_T_NBW_Early_neonatal),
                   sum_LGA_T_NBW_Livebirths=sum(LGA_T_NBW_Livebirths),
                   sum_AGA_T_NBW_Early_neonatal=sum(AGA_T_NBW_Early_neonatal),
                   sum_AGA_T_NBW_Livebirths=sum(AGA_T_NBW_Livebirths))


d<-as.data.frame(test)

y<-metabin(sum_LGA_T_NBW_Early_neonatal, sum_LGA_T_NBW_Livebirths, 
           sum_AGA_T_NBW_Early_neonatal,sum_AGA_T_NBW_Livebirths, studlab = Country, data=d, sm="RR", common = FALSE)

## Florest plot ##

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))

forest (y, leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


## random effect model 0.87 

#### 4. Calculate risk ratios for late mortality ####

test<- 
  data %>% group_by(Country) %>% 
  dplyr::select(AGA_T_NBW_Livebirths,AGA_T_NBW_Early_neonatal,
                AGA_T_NBW_Late_neonatal, AGA_T_NBW_Neonatal,
                LGA_T_NBW_Livebirths, LGA_T_NBW_Early_neonatal, 
                LGA_T_NBW_Late_neonatal, LGA_T_NBW_Neonatal) %>%
  dplyr::summarise(sum_LGA_T_NBW_Early_neonatal=sum(LGA_T_NBW_Early_neonatal),
                   sum_LGA_T_NBW_Livebirths=sum(LGA_T_NBW_Livebirths),
                   sum_AGA_T_NBW_Early_neonatal=sum(AGA_T_NBW_Early_neonatal),
                   sum_AGA_T_NBW_Livebirths=sum(AGA_T_NBW_Livebirths),
                   sum_LGA_T_NBW_Late_neonatal=sum(LGA_T_NBW_Late_neonatal),
                   sum_LGA_T_NBW_Late_livebirths=sum_LGA_T_NBW_Livebirths-sum_LGA_T_NBW_Early_neonatal,
                   sum_AGA_T_NBW_Late_neonatal=sum(AGA_T_NBW_Late_neonatal),
                   sum_AGA_T_NBW_Late_livebirths=sum_AGA_T_NBW_Livebirths-sum_AGA_T_NBW_Early_neonatal)


d<-as.data.frame(test)
d

y<-metabin(sum_LGA_T_NBW_Late_neonatal, sum_LGA_T_NBW_Late_livebirths, 
           sum_AGA_T_NBW_Late_neonatal, sum_AGA_T_NBW_Late_livebirths, 
           studlab = Country, data=d, sm="RR", common = FALSE)

## Florest plot ##

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


## Random effects model 0.81 

#### 5. Calculate risk ratios for neonatal mortality ####

test<- data %>% group_by(Country) %>% 
  dplyr::select(AGA_T_NBW_Livebirths,AGA_T_NBW_Early_neonatal,
                AGA_T_NBW_Late_neonatal, AGA_T_NBW_Neonatal,
                LGA_T_NBW_Livebirths, LGA_T_NBW_Early_neonatal, 
                LGA_T_NBW_Late_neonatal, LGA_T_NBW_Neonatal) %>%
  dplyr::summarise(sum_LGA_T_NBW_neonatal=sum(LGA_T_NBW_Neonatal),
                   sum_LGA_T_NBW_Livebirths=sum(LGA_T_NBW_Livebirths),
                   sum_AGA_T_NBW_neonatal=sum(AGA_T_NBW_Neonatal),
                   sum_AGA_T_NBW_Livebirths=sum(AGA_T_NBW_Livebirths))


d<-as.data.frame(test)
d

y<-metabin(sum_LGA_T_NBW_neonatal, sum_LGA_T_NBW_Livebirths, sum_AGA_T_NBW_neonatal, 
           sum_AGA_T_NBW_Livebirths, studlab = Country, data=d, sm="RR", common = FALSE)

## Florest plot ##

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
f<-forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


f


## Neonatal 
## random effect model 0.83 


#### 6. Calculate risk ratios for mortality, for macrosomia #####


test<- 
  data %>% group_by(Country) %>% 
  dplyr::select(B3_6000g_to_6500g_Livebirths, B3_6000g_to_6500g_Early_neonatal, B3_6000g_to_6500g_Late_neonatal, B3_6000g_to_6500g_Neonatal,
                B3_5500g_to_6000g_Livebirths, B3_5500g_to_6000g_Early_neonatal, B3_5500g_to_6000g_Late_neonatal, B3_5500g_to_6000g_Neonatal,
                B3_5000g_to_5500g_Livebirths, B3_5000g_to_5500g_Early_neonatal, B3_5000g_to_5500g_Late_neonatal, B3_5000g_to_5500g_Neonatal,
                B3_4500g_to_5000g_Livebirths, B3_4500g_to_5000g_Early_neonatal, B3_4500g_to_5000g_Late_neonatal, B3_4500g_to_5000g_Neonatal,
                B3_4000_to_4500g_Livebirths, B3_4000_to_4500g_Early_neonatal, B3_4000_to_4500g_Late_neonatal, B3_4000_to_4500g_Neonatal, 
                B3_3500_to_4000g_Livebirths, B3_3500_to_4000g_Early_neonatal, B3_3500_to_4000g_Late_neonatal, B3_3500_to_4000g_Neonatal, 
                B3_3000_to_3500g_Livebirths, B3_3000_to_3500g_Early_neonatal, B3_3000_to_3500g_Late_neonatal, B3_3000_to_3500g_Neonatal, 
                B3_2500_to_3000g_Livebirths, B3_2500_to_3000g_Early_neonatal,B3_2500_to_3000g_Late_neonatal, B3_2500_to_3000g_Neonatal) %>%
  
  dplyr::summarise(sum_B3_6000g_to_6500g_Livebirths=sum(B3_6000g_to_6500g_Livebirths), sum_B3_6000g_to_6500g_Early_neonatal=sum(B3_6000g_to_6500g_Early_neonatal),sum_B3_6000g_to_6500g_Late_neonatal=sum(B3_6000g_to_6500g_Late_neonatal), sum_B3_6000g_to_6500g_Neonatal=sum(B3_6000g_to_6500g_Neonatal),
                   sum_B3_5500g_to_6000g_Livebirths=sum(B3_5500g_to_6000g_Livebirths), sum_B3_5500g_to_6000g_Early_neonatal=sum(B3_5500g_to_6000g_Early_neonatal),sum_B3_5500g_to_6000g_Late_neonatal=sum(B3_5500g_to_6000g_Late_neonatal), sum_B3_5500g_to_6000g_Neonatal=sum(B3_5500g_to_6000g_Neonatal),
                   sum_B3_5000g_to_5500g_Livebirths=sum(B3_5000g_to_5500g_Livebirths), sum_B3_5000g_to_5500g_Early_neonatal=sum(B3_5000g_to_5500g_Early_neonatal),sum_B3_5000g_to_5500g_Late_neonatal=sum(B3_5000g_to_5500g_Late_neonatal), sum_B3_5000g_to_5500g_Neonatal=sum(B3_5000g_to_5500g_Neonatal),
                   sum_B3_4500g_to_5000g_Livebirths=sum(B3_4500g_to_5000g_Livebirths), sum_B3_4500g_to_5000g_Early_neonatal=sum(B3_4500g_to_5000g_Early_neonatal),sum_B3_4500g_to_5000g_Late_neonatal=sum(B3_4500g_to_5000g_Late_neonatal), sum_B3_4500g_to_5000g_Neonatal=sum(B3_4500g_to_5000g_Neonatal),
                   sum_B3_4000_to_4500g_Livebirths=sum(B3_4000_to_4500g_Livebirths), sum_B3_4000_to_4500g_Early_neonatal=sum(B3_4000_to_4500g_Early_neonatal),sum_B3_4000_to_4500g_Late_neonatal=sum(B3_4000_to_4500g_Late_neonatal), sum_B3_4000_to_4500g_Neonatal=sum(B3_4000_to_4500g_Neonatal),
                   sum_B3_2500_to_3000g_Livebirths=sum(B3_2500_to_3000g_Livebirths), sum_B3_2500_to_3000g_Early_neonatal=sum(B3_2500_to_3000g_Early_neonatal), sum_B3_3500_to_4000g_Late_neonatal=sum(B3_3500_to_4000g_Late_neonatal), sum_B3_3500_to_4000g_Neonatal=sum(B3_3500_to_4000g_Neonatal),
                   sum_B3_3000_to_3500g_Livebirths=sum(B3_3000_to_3500g_Livebirths), sum_B3_3000_to_3500g_Early_neonatal=sum(B3_3000_to_3500g_Early_neonatal), sum_B3_3000_to_3500g_Late_neonatal=sum(B3_3000_to_3500g_Late_neonatal), sum_B3_3000_to_3500g_Neonatal=sum(B3_3000_to_3500g_Neonatal),
                   sum_B3_3500_to_4000g_Livebirths=sum(B3_3500_to_4000g_Livebirths), sum_B3_3500_to_4000g_Early_neonatal=sum(B3_3500_to_4000g_Early_neonatal), sum_B3_2500_to_3000g_Late_neonatal=sum(B3_2500_to_3000g_Late_neonatal), sum_B3_2500_to_3000g_Neonatal=sum(B3_2500_to_3000g_Neonatal),
                   
                   sum_Livebirths_above_5000=sum_B3_6000g_to_6500g_Livebirths +sum_B3_5500g_to_6000g_Livebirths+ sum_B3_5000g_to_5500g_Livebirths,
                   sum_Livebirths_above_4500=sum_B3_6000g_to_6500g_Livebirths +sum_B3_5500g_to_6000g_Livebirths+ sum_B3_5000g_to_5500g_Livebirths+ sum_B3_4500g_to_5000g_Livebirths,
                   sum_Livebirths_above_4000=sum_B3_6000g_to_6500g_Livebirths +sum_B3_5500g_to_6000g_Livebirths+ sum_B3_5000g_to_5500g_Livebirths+ sum_B3_4500g_to_5000g_Livebirths + sum_B3_4000_to_4500g_Livebirths,
                   sum_Livebirths_2500_3999=sum_B3_2500_to_3000g_Livebirths+ sum_B3_3000_to_3500g_Livebirths+ sum_B3_3500_to_4000g_Livebirths,
                   
                   sum_Early_neonatal_above_5000=sum_B3_6000g_to_6500g_Early_neonatal+ sum_B3_5500g_to_6000g_Early_neonatal+ sum_B3_5000g_to_5500g_Early_neonatal,
                   sum_Early_neonatal_above_4500=sum_B3_6000g_to_6500g_Early_neonatal+ sum_B3_5500g_to_6000g_Early_neonatal+ sum_B3_5000g_to_5500g_Early_neonatal + sum_B3_4500g_to_5000g_Early_neonatal,
                   sum_Early_neonatal_above_4000=sum_B3_6000g_to_6500g_Early_neonatal+ sum_B3_5500g_to_6000g_Early_neonatal+ sum_B3_5000g_to_5500g_Early_neonatal + sum_B3_4500g_to_5000g_Early_neonatal + sum_B3_4000_to_4500g_Early_neonatal,
                   sum_Early_neonatal_2500_3999=sum_B3_2500_to_3000g_Early_neonatal + sum_B3_3000_to_3500g_Early_neonatal + sum_B3_3500_to_4000g_Early_neonatal,
                   
                   sum_Livebirths_Late_above_5000=sum_Livebirths_above_5000-sum_Early_neonatal_above_5000,
                   sum_Livebirths_Late_above_4500=sum_Livebirths_above_4500-sum_Early_neonatal_above_4500,
                   sum_Livebirths_Late_above_4000=sum_Livebirths_above_4000-sum_Early_neonatal_above_4000,
                   sum_Livebirths_Late_2500_3999=sum_Livebirths_2500_3999-sum_Early_neonatal_2500_3999,
                   
                   sum_Late_neonatal_above_5000=sum_B3_6000g_to_6500g_Late_neonatal+ sum_B3_5500g_to_6000g_Late_neonatal+ sum_B3_5000g_to_5500g_Late_neonatal,
                   sum_Late_neonatal_above_4500=sum_B3_6000g_to_6500g_Late_neonatal+ sum_B3_5500g_to_6000g_Late_neonatal+ sum_B3_5000g_to_5500g_Late_neonatal + sum_B3_4500g_to_5000g_Late_neonatal,
                   sum_Late_neonatal_above_4000=sum_B3_6000g_to_6500g_Late_neonatal+ sum_B3_5500g_to_6000g_Late_neonatal+ sum_B3_5000g_to_5500g_Late_neonatal + sum_B3_4500g_to_5000g_Late_neonatal + sum_B3_4000_to_4500g_Late_neonatal,
                   sum_Late_neonatal_2500_3999= sum_B3_2500_to_3000g_Late_neonatal + sum_B3_3000_to_3500g_Late_neonatal + sum_B3_3500_to_4000g_Late_neonatal, 
                   
                   sum_Neonatal_above_5000=sum_B3_6000g_to_6500g_Neonatal+ sum_B3_5500g_to_6000g_Neonatal+ sum_B3_5000g_to_5500g_Neonatal,
                   sum_Neonatal_above_4500=sum_B3_6000g_to_6500g_Neonatal+ sum_B3_5500g_to_6000g_Neonatal+ sum_B3_5000g_to_5500g_Neonatal+ sum_B3_4500g_to_5000g_Neonatal,
                   sum_Neonatal_above_4000=sum_B3_6000g_to_6500g_Neonatal+ sum_B3_5500g_to_6000g_Neonatal+ sum_B3_5000g_to_5500g_Neonatal+ sum_B3_4500g_to_5000g_Neonatal + sum_B3_4000_to_4500g_Neonatal,
                   sum_Neonatal_2500_3999=sum_B3_3500_to_4000g_Neonatal+ sum_B3_3000_to_3500g_Neonatal + sum_B3_2500_to_3000g_Neonatal)
                   
                   
              
            
d<-as.data.frame(test)
d


#### 7. Calculate risk ratios for babies above 4000 (vs normal) ####

## Early mortality, pooled RR 0.84 ##

y<-metabin(sum_Early_neonatal_above_4000 , sum_Livebirths_above_4000 , sum_Early_neonatal_2500_3999 , 
           sum_Livebirths_2500_3999, studlab = Country, data=d, sm="RR", common=FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


## Late neonatal, RR 0.74 ## 

y<-metabin(sum_Late_neonatal_above_4000  , sum_Livebirths_Late_above_4000 , sum_Late_neonatal_2500_3999  , 
           sum_Livebirths_Late_2500_3999, studlab = Country, data=d, sm="RR", common=FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)

## Neonatal 4000, RR 0.80 ## 

y<-metabin(sum_Neonatal_above_4000  , sum_Livebirths_above_4000 , sum_Neonatal_2500_3999  , 
           sum_Livebirths_2500_3999, studlab = Country, data=d, sm="RR", common=FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


#### 8. Calculate risk ratios for babies above 4500 (vs normal) ####

## Early mortality , RR 1.67 ##

y<-metabin(sum_Early_neonatal_above_4500 , sum_Livebirths_above_4500 , sum_Early_neonatal_2500_3999 , 
           sum_Livebirths_2500_3999, studlab = Country, data=d, sm="RR", common = FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


## Late mortality , RR  1.33 ##

y<-metabin(sum_Late_neonatal_above_4500  , sum_Livebirths_Late_above_4500 , sum_Late_neonatal_2500_3999  , 
           sum_Livebirths_Late_2500_3999, studlab = Country, data=d, sm="RR", common = FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)

## Neonatal mortality , RR 1.52 ##

y<-metabin(sum_Neonatal_above_4500 , sum_Livebirths_above_4500 , sum_Neonatal_2500_3999 , 
           sum_Livebirths_2500_3999, studlab = Country, data=d, sm="RR", common = FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


#### 9. Calculate risk ratios for babies above 5000 (vs normal) ####

## Early mortality , RR  5.96 ##

y<-metabin(sum_Early_neonatal_above_5000 , sum_Livebirths_above_5000 , sum_Early_neonatal_2500_3999 , 
           sum_Livebirths_2500_3999, studlab = Country, data=d, sm="RR", common = FALSE)
y

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)


## Late mortality, RR  4.24 ##

y<-metabin(sum_Late_neonatal_above_5000  , sum_Livebirths_Late_above_5000 , sum_Late_neonatal_2500_3999  , 
           sum_Livebirths_Late_2500_3999, studlab = Country, data=d, sm="RR", common = FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)

## Neonatal mortality, RR  4.54 ##

y<-metabin(sum_Neonatal_above_5000  , sum_Livebirths_above_5000 , sum_Neonatal_2500_3999  , 
           sum_Livebirths_2500_3999, studlab = Country, data=d, sm="RR", common = FALSE)

forest(y)
forest (y,leftcols = c("studlab", "effect"), rightcols = c("ci", "w.random"))
forest (y,leftcols = c("effect"), rightcols = c("ci", "w.random"), 
        col.square = "gray", col.diamond = "black", 
        print.tau2=F, print.Q=T, print.I2=T)



## End ##