#import library
library(haven)
library(tidyverse)
library(dplyr)
library(broom)


setwd("C:\\Users\\mike\\Documents\\r work\\IHS5")

food_security <- read_dta("HH_MOD_H.dta")
wellbeing_hh <- read_dta("HH_MOD_T.dta")
education <- read_dta("HH_MOD_C.dta")
view(education)

food_security2=food_security%>%
                        select(  "HHID","hh_h01","hh_h03b","hh_h03c","hh_h04","hh_h02a",
                        "hh_h02b")
wellbeing_hh2=wellbeing_hh%>%
                          select( "HHID","hh_t16","hh_t18")


education2=education%>% 
                    select("HHID","hh_c13","hh_c21","hh_c18","hh_c20")



#head(wellbeing_hh2)
#head(education2)
#head(food_security2)

#cleaning
wellbeing_clean= wellbeing_hh2%>% 
                    group_by(HHID)%>%
                     mutate( skipping_meal=case_when(
                                     hh_t16==2~1,
                                     hh_t16==1~0,
                                     TRUE~NA_real_ ),
                             went_hungry=case_when(
                                    hh_t18==2~1,
                                    hh_t18==1~0,
                                    TRUE~NA_real_
                                                   ),
                         wellbeing_score=skipping_meal + went_hungry)%>%       #combined scores
  
                                    select(HHID,skipping_meal,went_hungry,wellbeing_score)

                                                                                           

food_security_clean=food_security2%>% 
                       # filter(hh_h01==1)%>%
                            group_by(HHID)%>%
                              summarise(
                                     hdds=n() , #count food items consumed
                                     total_food_spent=sum(hh_h04, na.rm = "TRUE"),
                                     total_purchased=sum(hh_h02a, na.rm = "TRUE"),
                                     total_own_production=sum(hh_h02b, na.rm = "TRUE"),
                                     total_gifts=sum(hh_h03c,na.rm = "TRUE") )

summarise(food_security_clean)
nrow(food_security_clean)   # we have no groupings

education_clean= education2%>% 
                   group_by(HHID)%>%
                    mutate( 
                            attending=case_when(
                                                 hh_c13==1~1,
                                                 hh_c13==2~0,
                                                 TRUE~NA_real_),
                            ever_attended=case_when(
                                                 hh_c18==1~1,
                                                 hh_c18==2~0,
                                                  TRUE~NA_real_),
                            highest_grade=hh_c21, reason_not_attending=hh_c20)%>%
                                 select(HHID,attending,ever_attended,highest_grade,reason_not_attending)
                                                
    # View(education_clean)  

#merging
  data_merge= education_clean %>%
                left_join(wellbeing_clean, by ="HHID") %>%
                left_join(food_security_clean, by= "HHID")
   nrow(data_merge)
   View(data_merge)
   view(summarise(data_merge))
   
    #cleaning merged data
   data_merge_clean= data_merge%>%
                         filter(
                                 !is.na(attending),                 #filtering out missing data
                                 !is.na(ever_attended),
                                 !is.na(highest_grade),
                                 !is.na(reason_not_attending),
                                 !is.na(skipping_meal),
                                 !is.na(went_hungry),
                                 !is.na(hdds),
                                 !is.na(total_food_spent),
                                 !is.na(total_purchased),
                                 !is.na(total_own_production),
                                 !is.na(total_gifts)
                                                                   )
   view(data_merge_clean)
   
   
      
#tests
   #linear regression model
     model_f=lm(highest_grade ~skipping_meal+went_hungry,   #using selected variables
                  data = data_merge_clean
                 )
         summary(model_f)
         
#visualization
   result= barplot(c(3.19, 2.53), 
            names = c("Food Secure", "Skipped Meals"),
            col = c("steelblue", "coral"),
            ylim = c(0, 4),
            main = "Food Insecurity Reduces Grade Attainment",
            ylab = "Mean Highest Grade")
    text(c(0.7, 1.9), c(3.19, 2.53) + 0.15, c("3.19", "2.53"))
    text(1.3, 0.5, "Difference: -0.66 grades (p = 0.001)", cex = 0.9)   
      
  view(result)