.Rhistory

legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank(),
plot.margin = margin(-3,-3,-3,-3, "cm")
) +
coord_polar() +
# Add labels on top of each bar
geom_text(data=label_data2, aes(x=id, y=tot-50, label= str_wrap(ind_value, width = 22), hjust=label_data2$hjust), color="black", fontface="bold",alpha=1, size=pillartext, angle= label_data2$angle, inherit.aes = FALSE )
ggsave(inc, file=paste0(dir_out, group,'/',i,"_score_graph_incl_v",version,".png"), width=10, height=10)
ggsave(inc, file=paste0(dir_out, group,'/',i,"_score_graph_incl_v",version,".pdf"), width=10, height=10)
}
##########
# 2. Efficiency  # Same set of indicators for V0 and V1 (still using version control in case the indicators change)
##########
#detach(package:plyr)
for(i in country){
# Load dataset with country scores generated in the previous section
country_score <- read.csv(file=paste0(dir_out, group,'/',i,"_score_graphs_v",version,".csv"), header=TRUE, sep=",")
#Split data for graphs
country_score_graph2 <- country_score %>% filter(category== '2.Efficiency')
data2<-country_score_graph2
data2$order <- seq(1, nrow(data2))
# Move order column to the front so it does not get rejected when transforming to long format (column necessary to keep order of dimensions for the graph)
data2 <- data2 %>%
filter(!is.na(score))%>%
select(order, everything())
# Transform data in a tidy format (long format)
data2 <- stats::reshape(data2,
direction = "long",
varying = c("score", "values"),
v.names = "value",
idvar = c("order", "indicator", "category", "raw_value"),
timevar = "observation",
times = 1:2)
# Set a number of 'empty bar' to add at the end of each group
empty_bar <- 2
nObsType <- nlevels(as.factor(data2$observation))
to_add <- data.frame( matrix(NA, empty_bar*nlevels(as.factor(data2$category))*nObsType*.5, ncol(data2)))
colnames(to_add) <- colnames(data2)
to_add$category <- rep(levels(as.factor(data2$category)), each=empty_bar*nObsType*.5 )
data2 <- rbind(data2, to_add)
grpid = function(x) match(x, unique(x))
data2 <- data2 %>% mutate(group_id = group_indices(., category) %>% grpid)
#To group inverse and no inverse value
data2<- data2 %>% mutate(group_order=case_when(
group_id == 1 ~ "A",
TRUE ~ NA_character_ # all else are NA
))
data2 <- data2[order(data2$indicator),]
data2$id <- rep( seq(1, nrow(data2)/nObsType) , each=nObsType)
#Create labels for each column and angles (no need to adjust. With this method labels are aligned automatically)
label_data2 <- data2 %>% dplyr::group_by(id, indicator, raw_value) %>% dplyr::summarize(tot=sum(value))
number_of_bar2 <- nrow(label_data2)
angle2 <- 90 - 360 * (label_data2$id-0.5) /number_of_bar2     # I substract 0.5 because the letter must have the angle of the center of the bars. Not extreme right(1) or extreme left (0)
label_data2$hjust <- ifelse(angle2 < -90, 1, 0)
label_data2$angle <- ifelse(angle2 < -90, angle2+180, angle2)
# prepare a data frame for base lines
base_data2 <- data2 %>%
group_by(group_order) %>%
dplyr::summarize(start=min(id), end=max(id)
- empty_bar
) %>%
rowwise() %>%
mutate(title=mean(c(start, end)))
# prepare a data frame for grid (scales)
grid_data2 <- base_data2
grid_data2$end <- grid_data2$end[ c( nrow(grid_data2), 1:nrow(grid_data2)-1)] + 1.4
grid_data2$start <- grid_data2$start - .4
#grid_data2 <- grid_data2[-1,]
#Paste categories and raw values
data2$inv  <- ifelse(grepl("\\) \\*",data2$indicator), 1, 0)
data2$inv[is.na(data2$indicator)] <- NA
data2$cat <- paste(data2$observation,data2$group_order,data2$inv)
label_data2$ind_value<-paste(label_data2$raw_value, "-", label_data2$indicator)
e <- ggplot(data2) +
# Add the stacked bar
geom_bar(aes(x=as.factor(id), y=value, fill=cat, colour=cat), stat="identity", alpha=0.5, size=0.2) +
#scale_fill_manual(values=c("white","gold1")) +
scale_fill_manual(values=c("white","white","darkmagenta","plum")) +
scale_colour_manual(values=c("black","black","black","black")) +
# Add a val=100/75/50/25 lines. I do it at the beginning to make sur barplots are OVER it.
#geom_segment(data=grid_data2, aes(x = end, y = 100, xend = start, yend = 100), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 80, xend = start, yend = 80), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 60, xend = start, yend = 60), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 40, xend = start, yend = 40), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 20, xend = start, yend = 20), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_label(x = 0, y = 0, label = i, hjust=0.5, fontface="bold", size=10) +
annotate("text", x = rep(max(data2$id),5), y = c(20, 40, 60, 80, 100), label = c("20", "40", "60", "80", "100") , color="grey", size=5 , angle=0, fontface="bold", hjust=.5) +
ylim(0 , 120) +
theme_minimal() +
theme(
legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank(),
plot.margin = margin(-3,-3,-3,-3, "cm")
) +
coord_polar() +
# Add labels on top of each bar
geom_text(data=label_data2, aes(x=id, y=tot-50, label= str_wrap(ind_value, width = 22), hjust=label_data2$hjust), color="black", fontface="bold",alpha=1, size=pillartext, angle= label_data2$angle, inherit.aes = FALSE )
ggsave(e, file=paste0(dir_out, group,'/',i,"_score_graph_eff_v",version,".png"), width=10, height=10)
ggsave(e, file=paste0(dir_out, group,'/',i,"_score_graph_eff_v",version,".pdf"), width=10, height=10)
}
##########
#   RESILIENCE #    # Same set of indicators for V0 and V1 (still using version control in case it indicators change)
##########
#detach(package:plyr)
for(i in country){
# Load dataset with country scores generated in the previous section
country_score <- read.csv(file=paste0(dir_out, group,'/',i,"_score_graphs_v",version,".csv"), header=TRUE, sep=",")
##### Resilience #######
#Split data for graphs
country_score_graph2 <- country_score %>% filter(category== '4.Resilience')
data2<-country_score_graph2
data2$order <- seq(1, nrow(data2))
# Move order column to the front so it does not get rejected when transforming to long format (column necessary to keep order of dimensions for the graph)
data2 <- data2 %>%
filter(!is.na(score))%>%
select(order, everything())
# Transform data in a tidy format (long format)
data2 <- stats::reshape(data2,
direction = "long",
varying = c("score", "values"),
v.names = "value",
idvar = c("order", "indicator", "category", "raw_value"),
timevar = "observation",
times = 1:2)
# Set a number of 'empty bar' to add at the end of each group
empty_bar <- 2
nObsType <- nlevels(as.factor(data2$observation))
to_add <- data.frame( matrix(NA, empty_bar*nlevels(as.factor(data2$category))*nObsType*.5, ncol(data2)))
colnames(to_add) <- colnames(data2)
to_add$category <- rep(levels(as.factor(data2$category)), each=empty_bar*nObsType*.5 )
data2 <- rbind(data2, to_add)
grpid = function(x) match(x, unique(x))
data2 <- data2 %>% mutate(group_id = group_indices(., category) %>% grpid)
#To group inverse and no inverse value
data2<- data2 %>% mutate(group_order=case_when(
group_id == 1 ~ "A",
TRUE ~ NA_character_ # all else are NA
))
data2 <- data2[order(data2$indicator),]
data2$id <- rep( seq(1, nrow(data2)/nObsType) , each=nObsType)
#Create labels for each column and angles (no need to adjust. With this method labels are aligned automatically)
label_data2 <- data2 %>% dplyr::group_by(id, indicator, raw_value) %>% dplyr::summarize(tot=sum(value))
number_of_bar2 <- nrow(label_data2)
angle2 <- 90 - 360 * (label_data2$id-0.5) /number_of_bar2     # I substract 0.5 because the letter must have the angle of the center of the bars. Not extreme right(1) or extreme left (0)
label_data2$hjust <- ifelse(angle2 < -90, 1, 0)
label_data2$angle <- ifelse(angle2 < -90, angle2+180, angle2)
# prepare a data frame for base lines
base_data2 <- data2 %>%
group_by(group_order) %>%
dplyr::summarize(start=min(id), end=max(id)
- empty_bar
) %>%
rowwise() %>%
mutate(title=mean(c(start, end)))
# prepare a data frame for grid (scales)
grid_data2 <- base_data2
grid_data2$end <- grid_data2$end[ c( nrow(grid_data2), 1:nrow(grid_data2)-1)] + 1.4
grid_data2$start <- grid_data2$start - .4
#grid_data2 <- grid_data2[-1,]
#Paste categories and raw values
data2$inv  <- ifelse(grepl("\\) \\*",data2$indicator), 1, 0)
data2$inv[is.na(data2$indicator)] <- NA
data2$cat <- paste(data2$observation,data2$group_order,data2$inv)
label_data2$ind_value<-paste(label_data2$raw_value, "-", label_data2$indicator)
r <- ggplot(data2) +
# Add the stacked bar
geom_bar(aes(x=as.factor(id), y=value, fill=cat, colour=cat), stat="identity", alpha=0.5, size=0.2) +
#scale_fill_manual(values=c("white","green3")) +
scale_fill_manual(values=c("white","white","darkcyan","lightblue2")) +
scale_colour_manual(values=c("black","black","black","black")) +
# Add a val=100/75/50/25 lines. I do it at the beginning to make sur barplots are OVER it.
#geom_segment(data=grid_data2, aes(x = end, y = 100, xend = start, yend = 100), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 80, xend = start, yend = 80), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 60, xend = start, yend = 60), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 40, xend = start, yend = 40), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_segment(data=grid_data2, aes(x = end, y = 20, xend = start, yend = 20), colour = "grey", alpha=1, size=0.3 , inherit.aes = FALSE ) +
geom_label(x = 0, y = 0, label = i, hjust=0.5, fontface="bold", size=10) +
annotate("text", x = rep(max(data2$id),5), y = c(20, 40, 60, 80, 100), label = c("20", "40", "60", "80", "100") , color="grey", size=5 , angle=0, fontface="bold", hjust=.5) +
ylim(0 , 120) +
theme_minimal() +
theme(
legend.position = "none",
axis.text = element_blank(),
axis.title = element_blank(),
panel.grid = element_blank(),
plot.margin = margin(-3,-3,-3,-3, "cm")
) +
coord_polar() +
# Add labels on top of each bar
geom_text(data=label_data2, aes(x=id, y=tot-50, label= str_wrap(ind_value, width = 22), hjust=label_data2$hjust), color="black", fontface="bold",alpha=1, size=pillartext, angle= label_data2$angle, inherit.aes = FALSE )
ggsave(r, file=paste0(dir_out, group,'/',i,"_score_graph_res_v",version,".png"), width=10, height=10)
ggsave(r, file=paste0(dir_out, group,'/',i,"_score_graph_res_v",version,".pdf"), width=10, height=10)
}
file.remove(paste0(dir_out, group,'/',i,"_percentil_v",version,".csv"))
file.remove(paste0(dir_out, group,'/',i,"_score_graphs_v",version,".csv"))
?percent_ranj
?percent_rank
dplyr::percent_rank; dplyr::minrank
dplyr::percent_rank
dplyr::minrank
knitr::opts_knit$set(root.dir = "C:/Users/####/OneDrive/")
library(dplyr)
dplyr::minrank
tidyr::separate()
dplyr::minrank
dplyr::percent_rank
knitr::opts_knit$set(root.dir = "C:/Users/####/OneDrive/")
library(WDI)
library(dplyr)
library(readxl)
library(xml2)
library(rvest)
library(tidyverse)
library(data.table)
library(ggplot2)
library(grid)
library(lattice)
library(survival)
library(Formula)
library(Hmisc)
#detach(package:plyr)
dir_int<-"./Integrated data/" #Folder with data generated in code Benchmarking_Data_V3
#dir_out<- "./Data/Benchmarking/SD_Profiles/Output/"
dir_out<-"./Output/"
#In your coumputer, this directory create  separate folders for income group - ex: ./My Files/Output/lmic -./My Files/Output/umic, etc, in oder to save countries info per income group
getwd()
#Income group (select same category)
group <- 'lmics' ### HERE CHOOSE INCOME GROUP WE WANT TO COMPARE TO FROM ABOVE VARIABES ###
country="MAR"  #use one country for tests
#Set files names to read according to comparison groups
comparison_percentile <- paste0('DB_SD_',group,'_score')
comparison_score <- paste0('DB_SD_',group,'_score')
### Version indicator
#choose 0, 1 or 20.
#v= 0 Refers to a set of indicators selected and vetted by the Chief Economist Office -Sustainable Development (work lead by Jason Daniel Russ)
# v=1 refers to a set of indicators selected by the MENA - Sustainable Development team focusing on Egypt results
#Each version would show a different set of indicators
#v=9 refers to a set of indicators selected by the Egyptian MoPED - Ministry of Planning and Economic Development
version <- 2 #  0,1,2 or 9 for MoPED indicators
if (version != 2){
perctext <- 3.5
pillartext <- 4.5
} else {
perctext <-1.5
pillartext <- 3
}
#Retreive database with all the variables generated in "Benchmarking_V3"
integrated3 <- read.csv(paste0(dir_int, "World_SD_values_V5.csv"))
# For each version add all indicators
base <- c("country","iso3","income","region")
v0 <- c("risk_wellbeing", "risk_to_assets", "food_sec", "pop_dis", "rain_exp", "EN.POP.SLUM.UR.ZS", "pov_1016", 'SL.TLF.CACT.FE.ZS',"SH.STA.SMSS.ZS", "SH.H2O.SMDW.ZS", "HD.HCI.OVRL", "SI.POV.GINI", "SE.SEC.ENRR", "SH.IMM.IDPT", "account.t.d.7", "LGBT_acc", "wbl", "trust_ppl", "ref_idp", "conf_gov", "NY.ADJ.SVNG.GN.ZS", "EN.ATM.PM25.MC.M3", "EG.FEC.RNEW.ZS", "ER.H2O.INTR.PC","wat_qua", "WWT.new" ,"forest_loss","agri_value", "agri_land", "gini_ghg", "gov_effect", "cont_corr", "lpi_score", "sdg9_intuse", "water_efficiency_ceo", "Tenure.insecurity" )
v0_per <- paste(v0, "_per", sep='')
v1 <- c("risk_wellbeing", "risk_to_assets", "food_sec", "pop_dis", "rain_exp", "EN.POP.SLUM.UR.ZS", "pov_1016", 'SL.TLF.CACT.FE.ZS',"SH.STA.SMSS.ZS", "SH.H2O.SMDW.ZS", "HD.HCI.OVRL", "SI.POV.GINI", "SE.SEC.ENRR", "SH.IMM.IDPT", "account.t.d.7", "LGBT_acc", "wbl", "trust_ppl", "ref_idp", "sec_aparat", "conf_gov", "EN.ATM.PM25.MC.M3", "EG.FEC.RNEW.ZS", "ER.H2O.INTR.PC","wat_qua", "WWT.new", "waste_rate", "waste_managed", "BHV.new", "FSH.new", "GIB.new", "EN.ATM.CO2E.PC" , "forest_loss", "agri_value", "agri_land", "gini_ghg", "gov_effect", "cont_corr", "lpi_score", "sdg9_intuse", "water_efficiency_FAO", "Tenure.insecurity")
v1_per <- paste(v1, "_per", sep='')
v2 <- c("pop_dis", "rain_exp", "pov_1016", 'SL.TLF.CACT.FE.ZS', "HD.HCI.OVRL", "SI.POV.GINI", "SE.SEC.ENRR", "SH.IMM.IDPT", "account.t.d.7", "LGBT_acc", "wbl", "ref_idp", "EN.ATM.PM25.MC.M3", "EG.FEC.RNEW.ZS","wat_qua", "WWT.new", "waste_rate", "waste_managed", "BHV.new", "FSH.new", "EN.ATM.CO2E.PC" , "forest_loss", "agri_value", "agri_land", "gini_ghg", "gov_effect", "cont_corr", "lpi_score", "water_efficiency_ceo", "SG.GEN.PARL.ZS",
"SH.MED.PHYS.ZS", "SH.UHC.SRVS.CV.XD", "acc_sani_SDG", "acc_water_SDG", "rai", "rule_law", "IT.CEL.SETS.P2", "IT.NET.USER.ZS", "rapid_trans", "displaced", "wat_stress_SDG", "glob_hunger_index", "terror", "pop_slums_SDG", "heal_01_NDG", "habi_01_NDG", "habi_02_NDG", "infr_02_NDG", "food_03_NDG", "infr_03_NDG", "heal_03_NDG", "SH.STA.AIRP.P5","SH.STA.WASH.P5","fao_withdraw","ecos_04_NDG","chloroph_SDG","energ_int","EG.ELC.ACCS.ZS","reg_qua","ghg_pc","pop_epid","patents","IS.SHP.GCNW.XQ","dai")
v2_per <- paste(v2, "_per", sep='')
v9 <- c("food_sec",   "EN.POP.SLUM.UR.ZS", "SH.STA.TRAF.P5", "gdpgr_1519" ,"u_1519","infl_1519","hepat_SDG", "pov_1016", "SL.TLF.CACT.FE.ZS","acc_sani_SDG","acc_water_SDG", "SI.POV.GINI", "EG.ELC.ACCS.ZS", "SI.POV.NAHC", "SG.GEN.PARL.ZS","wom_manage_SDG", "FX.OWN.TOTL.FE.ZS" , "per_si_allsi.cov_pop_tot", "SE.ADT.LITR.FE.ZS" , "WWT.new", "waste_rate", "waste_managed", "BHV.new", "EN.ATM.CO2E.PP.GD", "sdg9_intuse", "debt_gdp_IMF", "capstock_imf", "gov_spe_edu_SDG", "balance_imf", "taxrev_imf","you_neet_1317","smes_wom")
v9_per <- paste(v9, "_per", sep='')
# Subset the database based on the indicators used in each version
if (version == 0){
integrated3 <- subset(integrated3, select = c(base,v0,v0_per))
} else {
if (version == 1) {
integrated3 <- subset(integrated3, select = c(base,v1,v1_per))
} else {
if (version == 2) {
integrated3 <- subset(integrated3, select = c(base,v2,v2_per))
} else {
integrated3 <- subset(integrated3, select = c(base,v9,v9_per))
}}}
#V0 Resilience (same set v0 and v1)
if (version==0){
integrated3<-mutate(integrated3, Resilience_v0 = rowMeans(select(integrated3, c("risk_wellbeing_per", "risk_to_assets_per", "food_sec_per", "pop_dis_per", "rain_exp_per", "EN.POP.SLUM.UR.ZS_per")), na.rm = TRUE))
#V0 Inclusion Indicators not included compared with V1: "sec_aparat" (Security Apparatus)
integrated3<-mutate(integrated3, Inclusion_v0 = rowMeans(select(integrated3, c("pov_1016_per", 'SL.TLF.CACT.FE.ZS_per',"SH.STA.SMSS.ZS_per", "SH.H2O.SMDW.ZS_per", "HD.HCI.OVRL_per", "SI.POV.GINI_per", "SE.SEC.ENRR_per", "SH.IMM.IDPT_per", "account.t.d.7_per", "LGBT_acc_per", "wbl_per", "trust_ppl_per", "ref_idp_per", "conf_gov_per")), na.rm = TRUE))
#V0 Sustainability Indicators not included compared with V1: ,"waste_managed_per, BHV_new_per", "FSH_new_per", "GIB_new_per", "EN.ATM.CO2E.PC_per", "waste_rate_per"
#Note, forest loss is excluded when the country has very low levels of forest loss
integrated3<-mutate(integrated3, Sustainability_v0 = rowMeans(select(integrated3, c("NY.ADJ.SVNG.GN.ZS_per", "EN.ATM.PM25.MC.M3_per", "EG.FEC.RNEW.ZS_per", "ER.H2O.INTR.PC_per","wat_qua_per", "WWT.new_per" ,"forest_loss_per" )), na.rm = TRUE))
#V0 Efficiency use water_efficiency_ceo - GNI/(waterwithdral*10*9) (calculated in this code- aqua dataframe)
#Other indicators not included: "reg_qua_per"
integrated3<-mutate(integrated3, Efficiency_v0 = rowMeans(select(integrated3, c("agri_value_per", "agri_land_per", "gini_ghg_per", "gov_effect_per", "cont_corr_per", "lpi_score_per", "sdg9_intuse_per", "water_efficiency_ceo_per", "Tenure.insecurity_per")), na.rm = TRUE))
} else {
if (version == 1) {
#V1 Resilience
integrated3<-mutate(integrated3, Resilience_v1 = rowMeans(select(integrated3, c("risk_wellbeing_per", "risk_to_assets_per", "food_sec_per", "pop_dis_per", "rain_exp_per", "EN.POP.SLUM.UR.ZS_per")), na.rm = TRUE))
#V1 Inclusion
integrated3<-mutate(integrated3, Inclusion_v1 = rowMeans(select(integrated3, c("pov_1016_per", 'SL.TLF.CACT.FE.ZS_per',"SH.STA.SMSS.ZS_per", "SH.H2O.SMDW.ZS_per", "HD.HCI.OVRL_per", "SI.POV.GINI_per", "SE.SEC.ENRR_per", "SH.IMM.IDPT_per", "account.t.d.7_per", "LGBT_acc_per", "wbl_per", "trust_ppl_per", "ref_idp_per", "sec_aparat_per", "conf_gov_per")), na.rm = TRUE))
# V1 Sustainability Indicators not included compared with V0: "NY.ADJ.SVNG.GN.ZS_per"
#Note, forest loss is excluded when the country has very low levels of forest loss
integrated3<-mutate(integrated3, Sustainability_v1 = rowMeans(select(integrated3, c("EN.ATM.PM25.MC.M3_per", "EG.FEC.RNEW.ZS_per", "ER.H2O.INTR.PC_per","wat_qua_per", "WWT.new_per", "waste_rate_per", "waste_managed_per", "BHV.new_per", "FSH.new_per", "GIB.new_per", "EN.ATM.CO2E.PC_per" , "forest_loss_per")), na.rm = TRUE))
#v1 Efficiency use water_efficiency_FAO_per, value added per water withdrawn, expressed in USD/m3 from FAO (details also in this code)
#Other indicators not included: "reg_qua_per"
integrated3<-mutate(integrated3, Efficiency_v1 = rowMeans(select(integrated3, c("agri_value_per", "agri_land_per", "gini_ghg_per", "gov_effect_per", "cont_corr_per", "lpi_score_per", "sdg9_intuse_per", "water_efficiency_FAO_per", "Tenure.insecurity_per")), na.rm = TRUE))
} else {
if (version == 2) {
#V2 Resilience
integrated3<-mutate(integrated3, Resilience_v2 = rowMeans(select(integrated3, c( "pop_dis_per", "rain_exp_per", "pop_epid_per", "pop_slums_SDG_per", "displaced_per" , "terror_per" , "heal_01_NDG_per" , "habi_01_NDG_per",  "habi_02_NDG_per",   "infr_02_NDG_per","infr_03_NDG_per" ,  "food_03_NDG_per",  "SH.MED.PHYS.ZS_per", "heal_03_NDG_per" ,"glob_hunger_index_per","wat_stress_SDG_per" )), na.rm = TRUE))
#V2 Inclusion
integrated3<-mutate(integrated3, Inclusion_v2 = rowMeans(select(integrated3, c("pov_1016_per", 'SL.TLF.CACT.FE.ZS_per', "HD.HCI.OVRL_per", "SI.POV.GINI_per", "SE.SEC.ENRR_per", "SH.IMM.IDPT_per", "account.t.d.7_per", "LGBT_acc_per", "wbl_per",  "ref_idp_per", "SG.GEN.PARL.ZS_per","rai_per" ,"acc_sani_SDG_per", "acc_water_SDG_per","SH.UHC.SRVS.CV.XD_per", "EG.ELC.ACCS.ZS_per" )), na.rm = TRUE))
# V2 Sustainability
integrated3<-mutate(integrated3, Sustainability_v2 = rowMeans(select(integrated3, c("EN.ATM.PM25.MC.M3_per", "EG.FEC.RNEW.ZS_per", "wat_qua_per", "WWT.new_per", "waste_rate_per", "waste_managed_per", "BHV.new_per", "FSH.new_per", "ghg_pc_per", "EN.ATM.CO2E.PC_per" , "forest_loss_per", "SH.STA.AIRP.P5_per" , "SH.STA.WASH.P5_per", "ecos_04_NDG_per", "chloroph_SDG_per", "fao_withdraw_per" )), na.rm = TRUE))
#v2 Efficiency
integrated3<-mutate(integrated3, Efficiency_v2 = rowMeans(select(integrated3, c("agri_value_per", "agri_land_per", "gini_ghg_per", "gov_effect_per", "cont_corr_per", "lpi_score_per",  "water_efficiency_ceo_per","rule_law_per" , "IT.CEL.SETS.P2_per", "IT.NET.USER.ZS_per","rapid_trans_per" ,"energ_int_per" , "reg_qua_per","patents_per","IS.SHP.GCNW.XQ_per","dai_per")), na.rm = TRUE))
} else {
#v9 (MoPED) Resilience
integrated3<-mutate(integrated3, Resilience_v9 = rowMeans(select(integrated3, c( "food_sec_per",   "EN.POP.SLUM.UR.ZS_per", "SH.STA.TRAF.P5_per", "gdpgr_1519_per","u_1519_per","infl_1519_per","hepat_SDG_per")), na.rm = TRUE))
#MoPED Inclusion
integrated3<-mutate(integrated3, Inclusion_v9 = rowMeans(select(integrated3, c("pov_1016_per", 'SL.TLF.CACT.FE.ZS_per',"acc_sani_SDG_per","acc_water_SDG_per" , "SI.POV.GINI_per", "EG.ELC.ACCS.ZS_per", "SI.POV.NAHC_per", "SG.GEN.PARL.ZS_per","wom_manage_SDG_per", "FX.OWN.TOTL.FE.ZS_per" , "per_si_allsi.cov_pop_tot_per", "SE.ADT.LITR.FE.ZS_per","you_neet_1317_per","smes_wom_per" )), na.rm = TRUE))
# MoPED Sustainability
integrated3<-mutate(integrated3, Sustainability_v9 = rowMeans(select(integrated3, c("WWT.new_per", "waste_rate_per", "waste_managed_per", "BHV.new_per", "EN.ATM.CO2E.PP.GD_per")), na.rm = TRUE))
# MoPED Efficiency
integrated3<-mutate(integrated3, Efficiency_v9 = rowMeans(select(integrated3, c("sdg9_intuse_per", "debt_gdp_IMF_per", "capstock_imf_per", "gov_spe_edu_SDG_per", "balance_imf_per", "taxrev_imf_per")), na.rm = TRUE))
}}}
#Save all the database with raw values and percentile data (worldwide data)
write.csv(integrated3, paste0(dir_int, "World_SD_values_V5_index_v",version,".csv"), row.names = FALSE)
# Load country dataset to create list of countries per income group that we will used in the loop
country_estimates_global <- read.csv(paste0(dir_int,"World_SD_values_V5_index_v",version,".csv"))
# Lower MICS
country_estimates <- country_estimates_global %>% filter(income=='Lower middle income')
lmics <- unique(country_estimates$iso3)
lmics <- as.character(lmics)
# Low income
country_estimates <- country_estimates_global %>% filter(income=='Low income')
lics <- unique(country_estimates$iso3)
lics <- as.character(lics)
# Upper MICS
country_estimates <- country_estimates_global %>% filter(income=='Upper middle income')
umics <- unique(country_estimates$iso3)
umics <- as.character(umics)
# High Income
country_estimates <- country_estimates_global %>% filter(income=='High income')
hics <- unique(country_estimates$iso3)
hics <- as.character(hics)
country_all <- read.csv(paste0(dir_int,"World_SD_values_V5_index_v",version,".csv"))
View(country_all)
density <- country_all >%> filter('income'=="Lower middle income")
density <- country_all >%> filter(income=="Lower middle income")
density <- country_all %>% filter(income=="Lower middle income")
View(density)
d <- density(density$ref_idp)
d <- density(density$ref_idp,na.rm=T)
View(d)
plot(d, main="Kernel Density of Miles Per Gallon")
d <- density(density$displaced,na.rm=T)
plot(d, main="Kernel Density of Miles Per Gallon")
knitr::opts_knit$set(root.dir = "C:/Users/####/OneDrive/")
ddd <- country_all %>% filter(income=="Lower middle income")
ggplot(ddd, aes(x=displaced)) +
geom_density() +
geom_vline(data=cdat, aes(xintercept=displaced.mean,  colour=cond),
linetype="dashed", size=1)
cddd <- ddply(ddd, summarise, rating.mean=mean(rating))
mean_var <- ddd %>%
pull(displaced) %>%
mean() %>%
signif(6)
mean_var <- ddd %>%
pull(displaced) %>%
mean(na.rm=T) %>%
signif(6)
med_var <- ddd %>%
pull(displaced) %>%
median(na.rm=T) %>%
signif(6)
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(data=mean_var, aes(xintercept=displaced.mean,  colour="tomato2"),
linetype="dashed", size=1) +
geom_vline(data=med_var, aes(xintercept=displaced.mean,  colour="tomato1"),
linetype="dashed", size=1)
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1)
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(aes(x=mean_var+50, label=paste0("Mean\n",mean_var), y=0.003,colour="tomato2"))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(aes(x=mean_var+1000, label=paste0("Mean\n",mean_var), y=0.003,colour="tomato2"))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(aes(x=mean_var+50, label=paste0("Mean\n",mean_var), y=0.003,colour="tomato2"))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=10,aes(x=mean_var+50, label=paste0("Mean\n",mean_var), y=0.003,colour="tomato2")) +
geom_text(size=10,aes(x=median_var+50, label=paste0("Median\n",mean_var), y=0.003,colour="tomato2"))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=10,aes(x=mean_var+50, label=paste0("Mean\n",mean_var), y=0.003,colour="tomato2")) +
geom_text(size=10,aes(x=med_var+50, label=paste0("Median\n",med_var), y=0.003,colour="tomato2"))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=1,aes(x=mean_var+50, label=paste0("Mean\n",mean_var), y=0.003,colour="tomato2")) +
geom_text(size=1,aes(x=med_var+50, label=paste0("Median\n",med_var), y=0.003,colour="tomato2"))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=3,aes(x=mean_var+50, label=paste0("Mean\n",mean_var), y=0.003,colour="tomato2")) +
geom_text(size=3,aes(x=med_var+50, label=paste0("Median\n",med_var), y=0.003,colour="tomato2"))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=3,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=3,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004))
ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=3,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=3,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
theme_minimal()
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=3,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=3,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
theme_minimal()
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.pdf"), width=10, height=10)
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
theme_minimal()
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.pdf"), width=10, height=10)
View(ddd)
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
theme_minimal()+
ylab("IDPs due to natural disasters (cases per 100,000")
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.pdf"), width=10, height=10)
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
ylab("IDPs due to natural disasters (cases per 100,000") +
theme_minimal()
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.pdf"), width=10, height=10)
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
xlab("IDPs due to natural disasters (cases per 100,000") +
theme_minimal()
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
xlab("IDPs due to natural disasters (cases per 100,000") +
theme_minimal()
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.pdf"), width=10, height=10)
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
xlab("IDPs due to natural disasters (cases per 100,000)") +
theme_minimal()
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.pdf"), width=10, height=10)
View(ddd)
View(d)
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
xlab("IDPs due to natural disasters (cases per 100,000)") +
theme_minimal()
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.png"), width=10, height=10)
r <- ggplot(ddd, aes(x=displaced)) +
geom_density(alpha=.5, fill="tomato3") +
geom_vline(xintercept=mean_var,  colour="tomato2",linetype="dashed", size=1) +
geom_vline(xintercept=med_var,  colour="tomato1",linetype="dashed", size=1) +
geom_text(size=4,aes(x=mean_var+80, label=paste0("Mean\n",mean_var), y=0.003)) +
geom_text(size=4,aes(x=med_var+80, label=paste0("Median\n",med_var), y=0.004)) +
xlab("IDPs due to natural disasters (cases per 100,000)") +
theme_minimal()
ggsave(r, file=paste0(dir_out,'/',"MAR_idpdistasters.png"), width=16, height=8)