pacman::p_load(pacman, rio, raster, rgdal, stars, 
               ggplot2, maps,  sf, RColorBrewer,
               #rasterVis, 
               ggthemes, viridis, tidyverse,
               RColorBrewer, haven,
               extrafont, dplyr, mapdeck, echarts4r, DT,
               data.table,readr, plotly, htmlwidgets, mapview,
               leafpop, mapedit, reshape2, rasterVis, tidync, pastecs,
               foreign, haven, Rfast, forcats, gdata,
               biscale, cowplot, gtools, PerformanceAnalytics,
               fastDummies, ipumsr, stringr, rgeos, maptools, tmap,
               scales, sjlabelled, jtools, huxtable, texreg,
               robustbase, readxl, miceadds , estimatr, purrr) 


memory.limit(10000000)

########################################################
# IPUMS - DHS data 
########################################################
# IPUMS DHS - Women

ddi <- read_ipums_ddi("data/DHS Data Process/idhs_00008.xml")
md <- read_ipums_micro(ddi)

# IPUMS DHS - men

ddi <- read_ipums_ddi("data/DHS Data Process/idhs_00009.xml")
wd <- read_ipums_micro(ddi)
########################################################

###$$$$$$$$$

# Men 

md <- md %>%  mutate(vul_eth = ifelse( COUNTRY == 120 & (ETHNICITYMN_CM == 600 | ETHNICITYMN_CM == 800 | ETHNICITYMN_CM == 500 ) , 0 , 1))     # Cameroon - non-excluded

md <- md %>%  
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 4  & ETHNICITYMN_AF == 1 , 0)) %>%                # Afghanistan - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 140  & (ETHNICITYMN_CF == 7 | ETHNICITYMN_CF == 8) , 0)) %>%                # Central Africa - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 148  & (ETHNICITYMN_TD2 == 13 | ETHNICITYMN_TD2 == 7) , 0)) %>%                # Chad - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 180  & (ETHNICITYMN_CD == 7 | ETHNICITYMN_CD == 8) , 0)) %>%                # Congo, Democratic Republic of the - non-excluded ( from 1960 - 2015) - Lunda and Luba people 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 384  & (ETHNICITYMN_CI == 13 | ETHNICITYMN_CI == 19 | ETHNICITYMN_CI == 2) , 0)) %>%                # Cote d'Ivoire- non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 231  & (ETHNICITYMN_ET == 5) , 0)) %>%                # Ethiopia- non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 288  & (ETHNICITYMN_GH == 1 | ETHNICITYMN_GH == 2 | ETHNICITYMN_GH == 5 | ETHNICITYMN_GH == 11 | ETHNICITYMN_GH == 12 ) , 0)) %>%                # Ghana - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 324  & (ETHNICITYMN_GN == 1) , 0)) %>%                # Guinea- non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 356  & (ETHNICITYMN_IA == 10 | ETHNICITYMN_IA == 20) , 0)) %>%                # India - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 404  & (ETHNICITYMN_KE == 210 | ETHNICITYMN_KE == 230 | ETHNICITYMN_KE == 240 | ETHNICITYMN_KE == 250) , 0)) %>%                # Kenya - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 454  & (ETHNICITYMN_MW == 1) , 0)) %>%                # Malawi - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 466  & (ETHNICITYMN_ML != 72 | ETHNICITYMN_ML != 74 | ETHNICITYMN_ML != 300 ) , 0)) %>%          # Mali - Excluded  - Not identifying other ethnicities or not reported one as non-excluded
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 508  & (ETHNICITYMN_MZ == 602 | ETHNICITYMN_MZ == 610 | ETHNICITYMN_MZ == 612 ) , 0)) %>%                # Mozambique - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 566  & (ETHNICITYMN_NG1 == 2 | ETHNICITYMN_NG1 == 3) , 0)) %>%                # Nigeria - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 178  & (ETHNICITYMN_CG == 105 | ETHNICITYMN_CG == 136) , 0)) %>%                #  Congo Brazzaville - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 686  & (ETHNICITYMN_SN == 1 | ETHNICITYMN_SN == 2 | ETHNICITYMN_SN == 3 | ETHNICITYMN_SN == 4 | ETHNICITYMN_SN == 5) , 0)) %>%                #  Senegal - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 710  & (ETHNICITYMN_ZA == 3 ) , 0)) %>%                # South Africa - non-excluded  - we consider only whites as non-excluded populations
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 768  & (ETHNICITYMN_TG == 10 | ETHNICITYMN_TG == 20 ) , 0)) %>%                # Togo - non-excluded - Ebe and Kibre - the power has been shifted between the two largest ethnic groups - in the year of survey 2013 they shared the power
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 800  & (ETHNICITYMN_UG == 1 | ETHNICITYMN_UG == 3 | ETHNICITYMN_UG == 40 | ETHNICITYMN_UG ==  41 | 
                                                          ETHNICITYMN_UG == 45 | ETHNICITYMN_UG == 50 | ETHNICITYMN_UG == 52 | ETHNICITYMN_UG ==  53 |
                                                          ETHNICITYMN_UG == 58 | ETHNICITYMN_UG == 26) , 0)) %>%                # Uganda - non-excluded - 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 894  & (ETHNICITYMN_ZM == 3 | ETHNICITYMN_ZM == 15 | ETHNICITYMN_ZM == 37 | ETHNICITYMN_ZM == 57 |
                                                          ETHNICITYMN_ZM == 33 | ETHNICITYMN_ZM == 28) , 0)) %>%                # Zambia- non-excluded  
  
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 204  & ETHNICITYMN_BJ == 1 , 0)) %>%                # Benin - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 400  & ETHNICITYMN_JO == 1 , 0)) %>%                # Jordan - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 430  & (ETHNICITYMN_LR == 5 |  ETHNICITYMN_LR == 13 |  ETHNICITYMN_LR == 90) , 0))                # Liberia - non-excluded (after 1997 and state collapse in 1990) 


###$$$$$$$$$

# Women 

wd <- wd %>%  mutate(vul_eth = ifelse( COUNTRY == 120 & (ETHNICITYCM == 600 | ETHNICITYCM == 800 | ETHNICITYCM == 500 ) , 0 , 1))     # Cameroon - non-excluded

wd <- wd %>%  
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 4  & ETHNICITYAF == 1 , 0)) %>%                # Afghanistan - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 140  & (ETHNICITYCF == 7 | ETHNICITYCF == 8) , 0)) %>%                # Central Africa - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 148  & (ETHNICITYTD2 == 13 | ETHNICITYTD2 == 7) , 0)) %>%                # Chad - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 180  & (ETHNICITYCD == 7 | ETHNICITYCD == 8) , 0)) %>%                # Congo, Democratic Republic of the - non-excluded ( from 1960 - 2015) - Lunda and Luba people 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 384  & (ETHNICITYCI2011 == 14 | ETHNICITYCI2011 == 23 | ETHNICITYCI2011 == 2) , 0)) %>%                # Cote d'Ivoire- non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 231  & (ETHNICITYET == 5) , 0)) %>%                # Ethiopia- non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 288  & (ETHNICITYGH == 1 | ETHNICITYGH == 2 | ETHNICITYGH == 5 | ETHNICITYGH == 11 | ETHNICITYGH == 12 ) , 0)) %>%                # Ghana - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 324  & (ETHNICITYGN == 1) , 0)) %>%                # Guinea- non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 356  & (ETHNICITYIA == 10 | ETHNICITYIA == 20) , 0)) %>%                # India - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 404  & (ETHNICITYKE == 10 | ETHNICITYKE == 102 | ETHNICITYKE == 103 | ETHNICITYKE == 105) , 0)) %>%                # Kenya - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 454  & (ETHNICITYMW == 1) , 0)) %>%                # Malawi - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 466  & (ETHNICITYML != 132 | ETHNICITYML != 133 | ETHNICITYML != 300 ) , 0)) %>%            # Mali - Excluded  - Not identifying other ethnicities or not reported one as non-excluded
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 508  & (ETHNICITYMZ == 602 | ETHNICITYMZ == 610 | ETHNICITYMZ == 612 ) , 0)) %>%                # Mozambique - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 516  & (ETHNICITYNM == 2 | ETHNICITYNM == 4 ) , 0)) %>%                # Namibia - non-excluded - Only reported for women
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 566  & (ETHNICITYNG == 2 | ETHNICITYNG == 3) , 0)) %>%                # Nigeria - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 178  & (ETHNICITYCG == 105 | ETHNICITYCG == 136) , 0)) %>%                #  Congo Brazzaville - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 686  & (ETHNICITYSN == 1 | ETHNICITYSN == 2 | ETHNICITYSN == 3 | ETHNICITYSN == 4 | ETHNICITYSN == 5) , 0)) %>%                #  Senegal - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 710  & (ETHNICITYZA == 3 ) , 0)) %>%                # South Africa - non-excluded  - we consider only whites as non-excluded populations
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 768  & (ETHNICITYTG == 10 | ETHNICITYTG == 20 ) , 0)) %>%                # Togo - non-excluded - Ebe and Kibre - the power has been shifted between the two largest ethnic groups - in the year of survey 2013 they shared the power
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 800  & (ETHNICITYUG == 1 | ETHNICITYUG == 3 | ETHNICITYUG == 40 | ETHNICITYUG ==  41 | 
                                                          ETHNICITYUG == 45 | ETHNICITYUG == 50 | ETHNICITYUG == 52 | ETHNICITYUG ==  53 |
                                                          ETHNICITYUG == 58 | ETHNICITYUG == 26) , 0)) %>%                # Uganda - non-excluded -
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 894  & (ETHNICITYZM == 3 | ETHNICITYZM == 15 | ETHNICITYZM == 37 | ETHNICITYZM == 57 |
                                                          ETHNICITYZM == 33 | ETHNICITYZM == 28) , 0)) %>%                # Zambia- non-excluded  
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 204  & ETHNICITYBJ == 1 , 0)) %>%                # Benin - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 400  & ETHNICITYJO == 1 , 0)) %>%                # Jordan - non-excluded 
  mutate(vul_eth = replace(vul_eth,  COUNTRY == 430  & (ETHNICITYLR == 5 |  ETHNICITYLR == 13 |  ETHNICITYLR == 90) , 0))                 # Liberia - non-excluded (after 1997 and state collapse in 1990)



############################################### 
# Main file for EECs analysis (DHS data)
############################################### 

# combine men and women DHS data : match by HHID

md <- md %>% 
  select(-starts_with("GEO")) %>%  select(-starts_with("DHS_"))  %>%  select(-starts_with("ETHNIC")) %>%  select(-starts_with("TEMP")) %>% 
  select(-starts_with("TEMPMIN")) %>%  mutate(gender = "man")


wd <- wd %>% 
  select(-starts_with("GEO")) %>%  select(-starts_with("DHS_"))  %>%  select(-starts_with("ETHNIC")) %>%  select(-starts_with("TEMP")) %>% 
  select(-starts_with("TEMPMIN")) %>%  mutate(gender2 = "woman")

### No data for ethnicity

m <- merge(md, wd, by=c('COUNTRY' , 'IDHSHID'), all=TRUE)

m <- m %>% 
  mutate(gender = ifelse(gender == "man", "man" , "woman" )) 

m <- m %>% 
  mutate(vul_eth = ifelse(gender == "man", vul_eth.x , vul_eth.y )) %>% 
  mutate(ELECTRCHH = ifelse(gender == "man", ELECTRCHH.x , ELECTRCHH.y )) %>% 
  mutate(COOKFUEL = ifelse(gender == "man", COOKFUEL.x , COOKFUEL.y )) %>% 
  mutate(KITCHEN = ifelse(gender == "man", KITCHEN.x , KITCHEN.y )) %>% 
  mutate(COOKWHERE = ifelse(gender == "man", COOKWHERE.x , COOKWHERE.y )) %>% 
  mutate(BIKEHH = ifelse(gender == "man", BIKEHH.x , BIKEHH.y )) %>% 
  mutate(CARHH = ifelse(gender == "man", CARHH.x , CARHH.y )) %>% 
  mutate(MOTORCYCLHH = ifelse(gender == "man", MOTORCYCLHH.x , MOTORCYCLHH.y )) %>% 
  mutate(STOVE = ifelse(gender == "man", STOVE.x , STOVE.y )) %>% 
  mutate(urban = ifelse(gender == "man", URBANMN , URBAN )) %>% 
  select(-ends_with(".x")) %>% select(-ends_with(".y"))


m <- m %>% 
  mutate(WEALTH_Q = ifelse(gender == "man", WEALTHQMN , WEALTHQ )) %>%          # Household wealth index in quintiles
  mutate(WEALTH_S = ifelse(gender == "man", WEALTHSMN , WEALTHS ))  %>%             # Wealth index factor score (5 decimals)
  select(c("COUNTRY"  ,  "IDHSHID" , "WEALTH_Q" , "WEALTH_S" , "ELECTRCHH" ,  "COOKFUEL" , "KITCHEN" ,  "urban" , 
           "COOKWHERE"  , "BIKEHH"   ,  "CARHH"   ,   "MOTORCYCLHH"  ,  "STOVE" , "vul_eth" ))





##############################
# country incomes 
# country region
##############################

# source: https://datahelpdesk.worldbank.org/knowledgebase/articles/906519


# country regions 


m <- m %>% 
  mutate(regions = ifelse(  COUNTRY == 204 |  COUNTRY == 288  |  COUNTRY == 454 |  COUNTRY == 694 |  COUNTRY == 768 |  COUNTRY == 800 |  COUNTRY == 710 |
                              COUNTRY == 894 | COUNTRY == 716 | COUNTRY == 180 | COUNTRY == 178 | COUNTRY == 140 |  COUNTRY == 384 |  COUNTRY == 120 |
                              COUNTRY == 231 | COUNTRY == 324 | COUNTRY == 404 |  COUNTRY == 430 | COUNTRY == 466 |  COUNTRY == 508 |  COUNTRY == 566 |
                              COUNTRY == 516 |  COUNTRY == 686 |  COUNTRY == 148, "Africa" , 
                            ifelse(  COUNTRY == 156 |  COUNTRY == 242 | COUNTRY == 360 | COUNTRY == 496  |  COUNTRY == 458 |  COUNTRY == 608 , "EAP" ,
                                     ifelse(  COUNTRY == 50 | COUNTRY == 524 | COUNTRY == 4 | COUNTRY == 356   , "SAR",
                                              ifelse( COUNTRY == 400  , "MENA", 
                                                      ifelse(COUNTRY == 68 |  COUNTRY == 76 | COUNTRY == 152 | COUNTRY == 170 |
                                                               COUNTRY == 188 |  COUNTRY == 192 |  COUNTRY == 218 | COUNTRY == 320 | 
                                                               COUNTRY == 340  |  COUNTRY == 484  |  COUNTRY == 558  |  COUNTRY == 591  |
                                                               COUNTRY == 600 | COUNTRY == 222 | COUNTRY == 858  , "LAC",
                                                             ifelse( COUNTRY == 840 , "North America", NA)))))))

# Sub-Saharan Africa (Africa) - East Asia and Pacific (EAP) - Europe and Central Asia (ECA) - Latin America and the Caribbean (LAC) - South Asia Region (SAR) - Middle East and North Africa (MENA)


# country income 

m <- m %>% 
  mutate(income_class = ifelse( COUNTRY == 562 | COUNTRY == 646 | COUNTRY == 108 | COUNTRY == 454 | COUNTRY == 694 | COUNTRY == 768 | COUNTRY == 800 | COUNTRY == 894 | COUNTRY == 4  | COUNTRY == 180 | COUNTRY == 231 |
                                  COUNTRY == 324 | COUNTRY == 430 | COUNTRY == 466 | COUNTRY == 508 | COUNTRY == 148 , "LICs" , 
                                ifelse( COUNTRY == 834 |   COUNTRY == 748 |  COUNTRY == 426 | COUNTRY == 716 | COUNTRY == 50 | COUNTRY == 204 |  COUNTRY == 68   |  COUNTRY == 288 |  COUNTRY == 340 |  COUNTRY == 360 |  COUNTRY == 496 |
                                          COUNTRY == 558 | COUNTRY == 524 |  COUNTRY == 608 | COUNTRY == 222 | COUNTRY == 140  | COUNTRY == 178 | COUNTRY == 384 |
                                          COUNTRY == 120 | COUNTRY == 356 | COUNTRY == 404 | COUNTRY == 566 | COUNTRY == 686 | COUNTRY == 24, "LMICs" ,
                                        ifelse(  COUNTRY == 76 |  COUNTRY == 156 | COUNTRY == 170  | COUNTRY == 188 | COUNTRY == 192  | COUNTRY == 218 | 
                                                   COUNTRY == 242  | COUNTRY == 192 | COUNTRY == 320 | COUNTRY == 484 | COUNTRY == 600  | COUNTRY == 710 |
                                                   COUNTRY == 400 | COUNTRY == 516 , "UMICs" ,
                                                 ifelse(  COUNTRY == 152 | COUNTRY == 591 | COUNTRY == 840 | COUNTRY == 858   , "HICs", NA)))))



# putting China  in UMICs (However it was in LMICs in  2000 - the year data is reported)
# Classification is based on the last report of the WB

###

# UMICs only the South Africa and Jordan in our sample - only South africa has data on wealth DHS

###


##############################

main <- m %>%  distinct(COUNTRY  ,  IDHSHID , .keep_all= TRUE)

###=############################################ 

# Graphs 

###=############################################ 
# Cooking Fuels
###=############################################ 

#

unique(main$COOKFUEL)

main <- main %>% 
  mutate(Cooking_Fuel = ifelse(COOKFUEL == 240   | COOKFUEL == 410 | COOKFUEL == 410 | COOKFUEL == 411 | COOKFUEL == 500 | COOKFUEL == 510 |
                                 COOKFUEL == 520 | COOKFUEL == 540 |  COOKFUEL == 600 |  COOKFUEL == 700 |  COOKFUEL == 710 , "Dirtiest" , 
                               ifelse(  COOKFUEL == 200 | COOKFUEL == 210 , "medium pollutants" , 
                                        ifelse(COOKFUEL == 100 | COOKFUEL == 220 | COOKFUEL == 221 | COOKFUEL == 222 | COOKFUEL == 802  |
                                                 COOKFUEL == 804 | COOKFUEL == 995  | COOKFUEL == 230  | COOKFUEL == 300 , "Cleanest" ,  NA))))




main1 <-main[!is.na(main$Cooking_Fuel) , ]
main1 <-main1[!is.na(main1$WEALTH_Q) , ]

###

main2 <-main1 %>% 
  select(c("COUNTRY" , "IDHSHID"  ,  "WEALTH_Q"  ,  "ELECTRCHH" ,   "COOKFUEL" ,   "CARHH"  , "MOTORCYCLHH"   ,   "vul_eth"  ,  "income_class" , "Cooking_Fuel"))

write_dta(main2, "data/DHS_EEC.dta")

