#This script produces the significance tests for the Final Cocoa Labor Report

#Call up a plotting window
#X11()

#Source the file that will produce all the weights and design with bootstrap replications 
#source("Create R datasets.R")

#Set the working directory to where the tables will be saved
#CHANGE DIRECTORY HERE
setwd("C:\Directory")
#save.image(file="Updated Tables.RData")
load("Updated Tables.RData")

#Load the packages required
pacman::p_load(openxlsx, survey)

#Here are the sampling designs we need for updating the tables
d_unequal_Ghana_CDI_2008
d_unequal_CDI_2008
d_unequal_Ghana_2008

d_unequal_Ghana_CDI_2013
d_unequal_CDI_2013
d_unequal_Ghana_2013

d_boot_EG_Ghana_CDI_2018
d_boot_EG_CDI_2018
d_boot_EG_Ghana_2018

################################
#Adding Gender2, Gender1, all
########################
# Adding Gender column
d_boot_EG_Ghana_2018$variables['Gender2'] <- (as.numeric(d_boot_EG_Ghana_2018$variables['gender'] == 'Male'))
d_boot_EG_Ghana_2018$variables['Gender1'] <- (as.numeric(d_boot_EG_Ghana_2018$variables['gender'] == 'Female'))

d_boot_EG_CDI_2018$variables['Gender2'] <- (as.numeric(d_boot_EG_CDI_2018$variables['gender'] == 'Male'))
d_boot_EG_CDI_2018$variables['Gender1'] <- (as.numeric(d_boot_EG_CDI_2018$variables['gender'] == 'Female'))

d_boot_EG_Ghana_CDI_2018$variables['Gender2'] <- (as.numeric(d_boot_EG_Ghana_CDI_2018$variables['gender'] == 'Male'))
d_boot_EG_Ghana_CDI_2018$variables['Gender1'] <- (as.numeric(d_boot_EG_Ghana_CDI_2018$variables['gender'] == 'Female'))

# Adding all age groups column
d_boot_EG_Ghana_2018$variables['all'] <- (as.numeric(d_boot_EG_Ghana_2018$variables['age_group1'] == 1 | 
                                                       d_boot_EG_Ghana_2018$variables['age_group2'] == 1 |
                                                       d_boot_EG_Ghana_2018$variables['age_group3'] == 1))

d_boot_EG_CDI_2018$variables['all'] <- (as.numeric(d_boot_EG_CDI_2018$variables['age_group1'] == 1 | 
                                                     d_boot_EG_CDI_2018$variables['age_group2'] == 1 |
                                                     d_boot_EG_CDI_2018$variables['age_group3'] == 1))

d_boot_EG_Ghana_CDI_2018$variables['all'] <- (as.numeric(d_boot_EG_Ghana_CDI_2018$variables['age_group1'] == 1 | 
                                                           d_boot_EG_Ghana_CDI_2018$variables['age_group2'] == 1 |
                                                           d_boot_EG_Ghana_CDI_2018$variables['age_group3'] == 1))

d_EG_Ghana_2018$variables['Gender2'] <- (as.numeric(d_EG_Ghana_2018$variables['gender'] == 'Male'))
d_EG_Ghana_2018$variables['Gender1'] <- (as.numeric(d_EG_Ghana_2018$variables['gender'] == 'Female'))

d_EG_CDI_2018$variables['Gender2'] <- (as.numeric(d_EG_CDI_2018$variables['gender'] == 'Male'))
d_EG_CDI_2018$variables['Gender1'] <- (as.numeric(d_EG_CDI_2018$variables['gender'] == 'Female'))

d_EG_Ghana_CDI_2018$variables['Gender2'] <- (as.numeric(d_EG_Ghana_CDI_2018$variables['gender'] == 'Male'))
d_EG_Ghana_CDI_2018$variables['Gender1'] <- (as.numeric(d_EG_Ghana_CDI_2018$variables['gender'] == 'Female'))

########################################################
# Table Values
##############################################

# Table Template ----------------------------------------------------------

# vars is a list of variable names as strings
cat('\n Table #')
vars <- c()

# Loop over variable names
for(x in vars){
  cat('\n \n Variable: ', x)
  
  # Get percentages (the second argument of subset is the variable(s) that
  # we filter by to get percentages)
  
  # 2008
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, ),na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008,  ),na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, ),na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2008, ),na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2008,  ),na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2008, ),na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

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

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

# Tables 1-5 do not require significance tests 

# Table 6 -----------------------------------------------------------------

d_boot_EG_Ghana_CDI_2018
d_boot_EG_CDI_2018
d_boot_EG_Ghana_2018

# Population counts
# Work 7 day
# TOTAL
svytotal(~workchild_7d, d_boot_EG_Ghana_CDI_2018)
svymean(~workchild_7d, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~workchild_7d, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~workchild_7d, d_boot_EG_CDI_2018)
svymean(~workchild_7d, d_boot_EG_CDI_2018)
confint(svymean(~workchild_7d, d_boot_EG_CDI_2018))

# Ghana
svytotal(~workchild_7d, d_boot_EG_Ghana_2018)
svymean(~workchild_7d, d_boot_EG_Ghana_2018)
confint(svymean(~workchild_7d, d_boot_EG_Ghana_2018))

# Work 12 month
# TOTAL
svytotal(~workchild_f, d_boot_EG_Ghana_CDI_2018)
svymean(~workchild_f, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~workchild_f, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~workchild_f, d_boot_EG_CDI_2018)
svymean(~workchild_f, d_boot_EG_CDI_2018)
confint(svymean(~workchild_f, d_boot_EG_CDI_2018))

# Ghana
svytotal(~workchild_f, d_boot_EG_Ghana_2018)
svymean(~workchild_f, d_boot_EG_Ghana_2018)
confint(svymean(~workchild_f, d_boot_EG_Ghana_2018))

# Work child labor
# TOTAL
svytotal(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)
svymean(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~cl_cocoa, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~cl_cocoa, d_boot_EG_CDI_2018)
svymean(~cl_cocoa, d_boot_EG_CDI_2018)
confint(svymean(~cl_cocoa, d_boot_EG_CDI_2018))

# Ghana
svytotal(~cl_cocoa, d_boot_EG_Ghana_2018)
svymean(~cl_cocoa, d_boot_EG_Ghana_2018)
confint(svymean(~cl_cocoa, d_boot_EG_Ghana_2018))

# Work hazardous child labor
# TOTAL
svytotal(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)
svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~haz_cocoa, d_boot_EG_CDI_2018)
svymean(~haz_cocoa, d_boot_EG_CDI_2018)
confint(svymean(~haz_cocoa, d_boot_EG_CDI_2018))

# Ghana
svytotal(~haz_cocoa, d_boot_EG_Ghana_2018)
svymean(~haz_cocoa, d_boot_EG_Ghana_2018)
confint(svymean(~haz_cocoa, d_boot_EG_Ghana_2018))

# Work V1
# TOTAL
svytotal(~V1a, d_boot_EG_Ghana_CDI_2018)
svymean(~V1a, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~V1a, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~V1a, d_boot_EG_CDI_2018)
svymean(~V1a, d_boot_EG_CDI_2018)
confint(svymean(~V1a, d_boot_EG_CDI_2018))

# Ghana
svytotal(~V1a, d_boot_EG_Ghana_2018)
svymean(~V1a, d_boot_EG_Ghana_2018)
confint(svymean(~V1a, d_boot_EG_Ghana_2018))

# Work V2
# TOTAL
svytotal(~V2a, d_boot_EG_Ghana_CDI_2018)
svymean(~V2a, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~V2a, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~V2a, d_boot_EG_CDI_2018)
svymean(~V2a, d_boot_EG_CDI_2018)
confint(svymean(~V2a, d_boot_EG_CDI_2018))

# Ghana
svytotal(~V2a, d_boot_EG_Ghana_2018)
svymean(~V2a, d_boot_EG_Ghana_2018)
confint(svymean(~V2a, d_boot_EG_Ghana_2018))

# Work V3
# TOTAL
svytotal(~V3a, d_boot_EG_Ghana_CDI_2018)
svymean(~V3a, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~V3a, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~V3a, d_boot_EG_CDI_2018)
svymean(~V3a, d_boot_EG_CDI_2018)
confint(svymean(~V3a, d_boot_EG_CDI_2018))

# Ghana
svytotal(~V3a, d_boot_EG_Ghana_2018)
svymean(~V3a, d_boot_EG_Ghana_2018)
confint(svymean(~V3a, d_boot_EG_Ghana_2018))

# Work V4
# TOTAL
svytotal(~V4a, d_boot_EG_Ghana_CDI_2018)
svymean(~V4a, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~V4a, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~V4a, d_boot_EG_CDI_2018)
svymean(~V4a, d_boot_EG_CDI_2018)
confint(svymean(~V4a, d_boot_EG_CDI_2018))

# Ghana
svytotal(~V4a, d_boot_EG_Ghana_2018)
svymean(~V4a, d_boot_EG_Ghana_2018)
confint(svymean(~V4a, d_boot_EG_Ghana_2018))

# Work V5
# TOTAL
svytotal(~V5a, d_boot_EG_Ghana_CDI_2018)
svymean(~V5a, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~V5a, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~V5a, d_boot_EG_CDI_2018)
svymean(~V5a, d_boot_EG_CDI_2018)
confint(svymean(~V5a, d_boot_EG_CDI_2018))

# Ghana
svytotal(~V5a, d_boot_EG_Ghana_2018)
svymean(~V5a, d_boot_EG_Ghana_2018)
confint(svymean(~V5a, d_boot_EG_Ghana_2018))

# Work V6
# TOTAL
svytotal(~V6a, d_boot_EG_Ghana_CDI_2018)
svymean(~V6a, d_boot_EG_Ghana_CDI_2018)
confint(svymean(~V6a, d_boot_EG_Ghana_CDI_2018))

# CDI
svytotal(~V6a, d_boot_EG_CDI_2018)
svymean(~V6a, d_boot_EG_CDI_2018)
confint(svymean(~V6a, d_boot_EG_CDI_2018))

# Ghana
svytotal(~V6a, d_boot_EG_Ghana_2018)
svymean(~V6a, d_boot_EG_Ghana_2018)
confint(svymean(~V6a, d_boot_EG_Ghana_2018))

# Table 7 -----------------------------------------------------------------

###Total
#2008/18
svymean(~(I(noncocoa==1)), d_unequal_Ghana_CDI_2008_HH)
svymean(~(I(noncocoa==1)), d_boot_EG_Ghana_CDI_2018_HH)

Z = abs(svymean(~noncocoa, d_boot_EG_Ghana_CDI_2018_HH)[1] - svymean(~noncocoa, d_unequal_Ghana_CDI_2008_HH)[1]) / 
  (sqrt(2)*SE(svymean(~noncocoa, d_boot_EG_Ghana_CDI_2018_HH)))
1 - pnorm(Z)

###CDI
#2008/18
svymean(~(I(noncocoa==1)), d_unequal_CDI_2008_HH)
svymean(~(I(noncocoa==1)), d_boot_EG_CDI_2018_HH)

Z = abs(svymean(~noncocoa, d_boot_EG_CDI_2018_HH)[1] - svymean(~noncocoa, d_unequal_CDI_2008_HH)[1]) / 
  (sqrt(2)*SE(svymean(~noncocoa, d_boot_EG_CDI_2018_HH)))
1 - pnorm(Z)

###Ghana
#2008/18
svymean(~(I(noncocoa==1)), d_unequal_Ghana_2008_HH)
svymean(~(I(noncocoa==1)), d_boot_EG_Ghana_2018_HH)

Z = abs(svymean(~noncocoa, d_boot_EG_Ghana_2018_HH)[1] - svymean(~noncocoa, d_unequal_Ghana_2008_HH)[1]) / 
  (sqrt(2)*SE(svymean(~noncocoa, d_boot_EG_Ghana_2018_HH)))
1 - pnorm(Z)

# Table 8 -----------------------------------------------------------------

#####fert
###Total
#2013/18
svyby(~fert, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean)
svyby(~fert, ~noncocoa, d_EG_Ghana_CDI_2018_HH, svymean)
Z = abs(svyby(~fert, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean)[1,2] - svyby(~fert, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~fert, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean))[1])
1 - pnorm(Z)

###CDI
#2013/18
svyby(~fert, ~noncocoa, d_unequal_CDI_2013_HH, svymean)
svyby(~fert, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean)
Z = abs(svyby(~fert, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean)[1,2] - svyby(~fert, ~noncocoa, d_unequal_CDI_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~fert, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean))[1])
1 - pnorm(Z)

###Ghana
#2013/18
svyby(~fert, ~noncocoa, d_unequal_Ghana_2013_HH, svymean)
svyby(~fert, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean)
Z = abs(svyby(~fert, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean)[1,2] - svyby(~fert, ~noncocoa, d_unequal_Ghana_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~fert, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean))[1])
1 - pnorm(Z)

#####pest
###Total
#2013/18
svyby(~pest, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean)
svyby(~pest, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean)
Z = abs(svyby(~pest, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean)[1,2] - svyby(~pest, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~pest, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean))[1])
1 - pnorm(Z)

###CDI
#2013/18
svyby(~pest, ~noncocoa, d_unequal_CDI_2013_HH, svymean)
svyby(~pest, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean)
Z = abs(svyby(~pest, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean)[1,2] - svyby(~pest, ~noncocoa, d_unequal_CDI_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~pest, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean))[1])
1 - pnorm(Z)

###Ghana
#2013/18
svyby(~pest, ~noncocoa, d_unequal_Ghana_2013_HH, svymean)
svyby(~pest, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean)
Z = abs(svyby(~pest, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean)[1,2] - svyby(~pest, ~noncocoa, d_unequal_Ghana_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~pest, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean))[1])
1 - pnorm(Z)

#####herb
###Total
#2013/18
svyby(~herb, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean)
svyby(~herb, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean)
Z = abs(svyby(~herb, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean)[1,2] - svyby(~herb, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~herb, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean))[1])
1 - pnorm(Z)

###CDI
#2013/18
svyby(~herb, ~noncocoa, d_unequal_CDI_2013_HH, svymean)
svyby(~herb, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean)
Z = abs(svyby(~herb, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean)[1,2] - svyby(~herb, ~noncocoa, d_unequal_CDI_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~herb, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean))[1])
1 - pnorm(Z)

###Ghana
#2013/18
svyby(~herb, ~noncocoa, d_unequal_Ghana_2013_HH, svymean)
svyby(~herb, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean)
Z = abs(svyby(~herb, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean)[1,2] - svyby(~herb, ~noncocoa, d_unequal_Ghana_2013_HH, svymean)[1,2]) / 
  (sqrt(2)*SE(svyby(~herb, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean))[1])
1 - pnorm(Z)

# Table 9 -----------------------------------------------------------------
cat('\n Table 9')
filter_vars <- c('all','Gender2', 'Gender1', 'age_group1', 'age_group2', 'age_group3')
for(y in filter_vars){
  cat('\n Variable: ', y)
  
  #2008
  G1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_2008, get(y)==1), na.rm=TRUE)
  CDI1 <-  svymean(~SchoolAttendance, subset(d_unequal_CDI_2008, get(y)==1), na.rm=TRUE)
  TOT1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_CDI_2008, get(y)==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~SchoolAttendance, subset(d_EG_Ghana_2018, get(y)==1), na.rm=TRUE)
  CDI2 <- svymean(~SchoolAttendance, subset(d_EG_CDI_2018, get(y)==1), na.rm=TRUE)
  TOT2 <- svymean(~SchoolAttendance, subset(d_EG_Ghana_CDI_2018, get(y)==1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
  
}

# Table 10 -----------------------------------------------------------------
cat('\n Table 10')
vars <- c('workchild_f', 'workchild_7d')
for(x in vars){
  cat('\n Variable: ', x)
  
  #2008
  G1 <-  svymean(~get(x), d_unequal_Ghana_2008, na.rm=TRUE)
  CDI1 <-  svymean(~get(x), d_unequal_CDI_2008, na.rm=TRUE)
  TOT1 <-  svymean(~get(x), d_unequal_Ghana_CDI_2008, na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
  
}

# Table 11 -----------------------------------------------------------------
cat('\n Table 11')
vars <- c('agchild_f', 'agchild_7d')
for(x in vars){
  cat('\n Variable: ', x)
  
  #2008
  G1 <-  svymean(~get(x), d_unequal_Ghana_2008, na.rm=TRUE)
  CDI1 <-  svymean(~get(x), d_unequal_CDI_2008, na.rm=TRUE)
  TOT1 <-  svymean(~get(x), d_unequal_Ghana_CDI_2008, na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 12 ----------------------------------------------------------------
cat('\n Table 12')
vars <- c('agchild_f', 'cl_ag')
for(x in vars){
  cat('\n Variable: ', x)
  #2008
  G1 <-  svymean(~get(x), d_unequal_Ghana_2008, na.rm=TRUE)
  CDI1 <-  svymean(~get(x), d_unequal_CDI_2008, na.rm=TRUE)
  TOT1 <-  svymean(~get(x), d_unequal_Ghana_CDI_2008, na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}


# Table 13 ----------------------------------------------------------------
cat('\n Table 13')
vars <- c('cocoachild_f', 'cocoachild_7d')
for(x in vars){
  cat('\n Variable: ', x)
  #2008
  G1 <-  svymean(~get(x), d_unequal_Ghana_2008, na.rm=TRUE)
  CDI1 <-  svymean(~get(x), d_unequal_CDI_2008, na.rm=TRUE)
  TOT1 <-  svymean(~get(x), d_unequal_Ghana_CDI_2008, na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 14  ----------------------------------------------------------------
cat('\n Table 14')
filter_vars <- c('age_group1', 'age_group2', 'age_group3', 'all')
for(y in filter_vars){
  
  # Percentages
  cat('\n Variable: ', y)
  cat(' \n Percentages')
  #2008
  G1 <- svymean(~vio_max_cocoa, subset(d_unequal_Ghana_2008, get(y)==1), na.rm=TRUE)
  CDI1 <- svymean(~vio_max_cocoa, subset(d_unequal_CDI_2008, get(y)==1), na.rm=TRUE)
  TOT1 <- svymean(~vio_max_cocoa, subset(d_unequal_Ghana_CDI_2008, get(y)==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_Ghana_2018, get(y)== 1), na.rm=TRUE)
  CDI2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_CDI_2018, get(y)== 1, na.rm=TRUE))
  TOT2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_Ghana_CDI_2018, get(y)== 1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
  
  # Hours Worked
  cat(' \n Average hours worked')
  
  #2008
  G1 <- svymean(~TotalHours, subset(d_unequal_Ghana_2008, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  CDI1 <- svymean(~TotalHours, subset(d_unequal_CDI_2008, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  TOT1 <- svymean(~TotalHours, subset(d_unequal_Ghana_CDI_2008, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~TotalHours, subset(d_boot_EG_Ghana_2018, cocoachild_f==1 & get(y)== 1), na.rm=TRUE)
  CDI2 <- svymean(~TotalHours, subset(d_boot_EG_CDI_2018, cocoachild_f==1 & get(y)== 1, na.rm=TRUE))
  TOT2 <- svymean(~TotalHours, subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f==1 & get(y)== 1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
  
}

# Table 15 ----------------------------------------------------------------
cat('\n Table 15')
vars <- c('land_preparation', 'planting_activities', 'farm_maintenance',
          'harvest_activities', 'postharv_activities')
for(x in vars){
  cat('\n Variable: ', x)
  #2008
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, cocoachild_f==1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008, cocoachild_f==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, cocoachild_f==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoachild_f== 1), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoachild_f== 1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f== 1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 16 ----------------------------------------------------------------
#All 
svymean(~cocoachild_f, d_unequal_Ghana_CDI_2008)
svymean(~cocoachild_f, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~cocoachild_f, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~cocoachild_f, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~cocoachild_f, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~cl_cocoa, d_unequal_Ghana_CDI_2008)
svymean(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~cl_cocoa, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~haz_cocoa, d_unequal_Ghana_CDI_2008)
svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~haz_cocoa, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

#CDI
svymean(~cocoachild_f, d_unequal_CDI_2008)
svymean(~cocoachild_f, d_boot_EG_CDI_2018)
Z = abs(svymean(~cocoachild_f, d_boot_EG_CDI_2018)[1] - svymean(~cocoachild_f, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~cocoachild_f, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~cl_cocoa, d_unequal_CDI_2008)
svymean(~cl_cocoa, d_boot_EG_CDI_2018)
Z = abs(svymean(~cl_cocoa, d_boot_EG_CDI_2018)[1] - svymean(~cl_cocoa, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~cl_cocoa, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~haz_cocoa, d_unequal_CDI_2008)
svymean(~haz_cocoa, d_boot_EG_CDI_2018)
Z = abs(svymean(~haz_cocoa, d_boot_EG_CDI_2018)[1] - svymean(~haz_cocoa, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~haz_cocoa, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

#Ghana
svymean(~cocoachild_f, d_unequal_Ghana_2008)
svymean(~cocoachild_f, d_boot_EG_Ghana_2018)
Z = abs(svymean(~cocoachild_f, d_boot_EG_Ghana_2018)[1] - svymean(~cocoachild_f, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~cocoachild_f, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~cl_cocoa, d_unequal_Ghana_2008)
svymean(~cl_cocoa, d_boot_EG_Ghana_2018)
Z = abs(svymean(~cl_cocoa, d_boot_EG_Ghana_2018)[1] - svymean(~cl_cocoa, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~cl_cocoa, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~haz_cocoa, d_unequal_Ghana_2008)
svymean(~haz_cocoa, d_boot_EG_Ghana_2018)
Z = abs(svymean(~haz_cocoa, d_boot_EG_Ghana_2018)[1] - svymean(~haz_cocoa, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~haz_cocoa, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

# Table 17 ----------------------------------------------------------------

##Total
svymean(~V1a, d_unequal_Ghana_CDI_2008)
svymean(~V1a, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~V1a, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~V1a, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V1a, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~V2a, d_unequal_Ghana_CDI_2008)
svymean(~V2a, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~V2a, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~V2a, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V2a, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~V3a, d_unequal_Ghana_CDI_2008)
svymean(~V3a, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~V3a, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~V3a, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V3a, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~V4a, d_unequal_Ghana_CDI_2008)
svymean(~V4a, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~V4a, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~V4a, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V4a, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~V5a, d_unequal_Ghana_CDI_2008)
svymean(~V5a, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~V5a, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~V5a, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V5a, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~V6a, d_unequal_Ghana_CDI_2008)
svymean(~V6a, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~V6a, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~V6a, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V6a, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~haz_cocoa, d_unequal_Ghana_CDI_2008)
svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~haz_cocoa, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

svymean(~Via, d_unequal_Ghana_CDI_2008)
svymean(~Via, d_boot_EG_Ghana_CDI_2018)
Z = abs(svymean(~Via, d_boot_EG_Ghana_CDI_2018)[1] - svymean(~Via, d_unequal_Ghana_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~Via, d_boot_EG_Ghana_CDI_2018)))
1 - pnorm(Z)

##CDI
svymean(~V1a, d_unequal_CDI_2008)
svymean(~V1a, d_boot_EG_CDI_2018)
Z = abs(svymean(~V1a, d_boot_EG_CDI_2018)[1] - svymean(~V1a, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V1a, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~V2a, d_unequal_CDI_2008)
svymean(~V2a, d_boot_EG_CDI_2018)
Z = abs(svymean(~V2a, d_boot_EG_CDI_2018)[1] - svymean(~V2a, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V2a, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~V3a, d_unequal_CDI_2008)
svymean(~V3a, d_boot_EG_CDI_2018)
Z = abs(svymean(~V3a, d_boot_EG_CDI_2018)[1] - svymean(~V3a, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V3a, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~V4a, d_unequal_CDI_2008)
svymean(~V4a, d_boot_EG_CDI_2018)
Z = abs(svymean(~V4a, d_boot_EG_CDI_2018)[1] - svymean(~V4a, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V4a, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~V5a, d_unequal_CDI_2008)
svymean(~V5a, d_boot_EG_CDI_2018)
Z = abs(svymean(~V5a, d_boot_EG_CDI_2018)[1] - svymean(~V5a, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V5a, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~V6a, d_unequal_CDI_2008)
svymean(~V6a, d_boot_EG_CDI_2018)
Z = abs(svymean(~V6a, d_boot_EG_CDI_2018)[1] - svymean(~V6a, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V6a, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~haz_cocoa, d_unequal_CDI_2008)
svymean(~haz_cocoa, d_boot_EG_CDI_2018)
Z = abs(svymean(~haz_cocoa, d_boot_EG_CDI_2018)[1] - svymean(~haz_cocoa, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~haz_cocoa, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

svymean(~Via, d_unequal_CDI_2008)
svymean(~Via, d_boot_EG_CDI_2018)
Z = abs(svymean(~Via, d_boot_EG_CDI_2018)[1] - svymean(~Via, d_unequal_CDI_2008)[1]) / 
  (sqrt(2)*SE(svymean(~Via, d_boot_EG_CDI_2018)))
1 - pnorm(Z)

##Ghana
svymean(~V1a, d_unequal_Ghana_2008)
svymean(~V1a, d_boot_EG_Ghana_2018)
Z = abs(svymean(~V1a, d_boot_EG_Ghana_2018)[1] - svymean(~V1a, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V1a, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~V2a, d_unequal_Ghana_2008)
svymean(~V2a, d_boot_EG_Ghana_2018)
Z = abs(svymean(~V2a, d_boot_EG_Ghana_2018)[1] - svymean(~V2a, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V2a, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~V3a, d_unequal_Ghana_2008)
svymean(~V3a, d_boot_EG_Ghana_2018)
Z = abs(svymean(~V3a, d_boot_EG_Ghana_2018)[1] - svymean(~V3a, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V3a, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~V4a, d_unequal_Ghana_2008)
svymean(~V4a, d_boot_EG_Ghana_2018)
Z = abs(svymean(~V4a, d_boot_EG_Ghana_2018)[1] - svymean(~V4a, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V4a, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~V5a, d_unequal_Ghana_2008)
svymean(~V5a, d_boot_EG_Ghana_2018)
Z = abs(svymean(~V5a, d_boot_EG_Ghana_2018)[1] - svymean(~V5a, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V5a, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~V6a, d_unequal_Ghana_2008)
svymean(~V6a, d_boot_EG_Ghana_2018)
Z = abs(svymean(~V6a, d_boot_EG_Ghana_2018)[1] - svymean(~V6a, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~V6a, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~haz_cocoa, d_unequal_Ghana_2008)
svymean(~haz_cocoa, d_boot_EG_Ghana_2018)
Z = abs(svymean(~haz_cocoa, d_boot_EG_Ghana_2018)[1] - svymean(~haz_cocoa, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~haz_cocoa, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

svymean(~Via, d_unequal_Ghana_2008)
svymean(~Via, d_boot_EG_Ghana_2018)
Z = abs(svymean(~Via, d_boot_EG_Ghana_2018)[1] - svymean(~Via, d_unequal_Ghana_2008)[1]) / 
  (sqrt(2)*SE(svymean(~Via, d_boot_EG_Ghana_2018)))
1 - pnorm(Z)

# Table 18 ----------------------------------------------------------------
cat('\n Table 18')
vars <- c('V_1', 'V_2', 'V_3', 'V_4', 'V_5', 'V_6', 'V_7', 'Vi')
for(x in vars){
  cat('\n Variable: ', x)
  
  #2008
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, cocoachild_f==1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008, cocoachild_f==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, cocoachild_f==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoachild_f== 1), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoachild_f== 1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f== 1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 19 ----------------------------------------------------------------
cat('\n Table 19')
vars <- c('V3_12m', 'spray_cocoa', 'present_cocoa', 'enter_cocoa', 'water_cocoa', 'agrochem_cocoa')
for(x in vars){
  cat('\n Variable: ', x)
  #2008
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, cocoachild_f==1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008, cocoachild_f==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, cocoachild_f==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoachild_f== 1), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoachild_f== 1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f== 1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 20  ---------------------------------------------------------------
cat('\n Table 20')

# Population

# In Cocoa Production
Gpop1 <- svytotal(~ones, subset(d_boot_EG_Ghana_2018, agchild_f==1 & cocoa_hh==1), na.rm=TRUE)
CDIpop1 <- svytotal(~ones, subset(d_boot_EG_CDI_2018, agchild_f==1 & cocoa_hh==1), na.rm=TRUE)
TOTpop1 <- svytotal(~ones, subset(d_boot_EG_Ghana_CDI_2018, agchild_f==1 & cocoa_hh==1), na.rm=TRUE)

# Non-Cocoa Agriculture
Gpop2 <- svytotal(~ones, subset(d_boot_EG_Ghana_2018, agchild_f==1 & cocoa_hh==0), na.rm=TRUE)
CDIpop2 <- svytotal(~ones, subset(d_boot_EG_CDI_2018, agchild_f==1 & cocoa_hh==0), na.rm=TRUE)
TOTpop2 <- svytotal(~ones, subset(d_boot_EG_Ghana_CDI_2018, agchild_f==1 & cocoa_hh==0), na.rm=TRUE)


cat('\n Ghana cocoa population: ', Gpop1, '\n CDI cocoa population: ', CDIpop1, 
    '\n Total cocoa population: ', TOTpop1)
cat('\n Ghana noncocoa population: ', Gpop2, '\n CDI noncocoa population: ', CDIpop2,
    '\n Total noncocoa population: ', TOTpop2)

# Percentages
vars <- c('F2_1_ch_m', 'F2_2_ch_m', 'F2_3_ch_m', 'F2_4_ch_m', 'F2_5_ch_m',
          'F2_6_ch_m', 'F2_7_ch_m', 'F2_8_ch_m', 'F2_9_ch_m')
for(x in vars){
  cat('\n Variable: ', x)
  
  #In Cocoa Production
  G1 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, agchild_f==1 & cocoa_hh==1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, agchild_f==1 & cocoa_hh==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, agchild_f==1 & cocoa_hh==1), na.rm=TRUE)
  
  #Not in cocoa production
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, agchild_f==1 & cocoa_hh==0), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, agchild_f==1 & cocoa_hh==0), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, agchild_f==1 & cocoa_hh==0), na.rm=TRUE)
  
  cat('\n Ghana cocoa: ', G1, '\n CDI cocoa: ', CDI1, '\n Total cocoa: ', TOT1)
  cat('\n Ghana noncocoa: ', G2, '\n CDI noncocoa: ', CDI2, '\n Total noncocoa: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 21  ----------------------------------------------------------------
cat('\n Table 21')

# Population

# Hazardous work
cat('\n Hazardous')
Gpop1 <- svytotal(~ones, subset(d_boot_EG_Ghana_2018,  agchild_f==1 & hazardous_12m==1), na.rm=TRUE)
CDIpop1 <- svytotal(~ones, subset(d_boot_EG_CDI_2018,  agchild_f==1 & hazardous_12m==1), na.rm=TRUE)
TOTpop1 <- svytotal(~ones, subset(d_boot_EG_Ghana_CDI_2018,  agchild_f==1 & hazardous_12m==1), na.rm=TRUE)

# Not Hazardous
cat('\n not hazardous')
Gpop2 <- svytotal(~ones, subset(d_boot_EG_Ghana_2018,  agchild_f==1 & hazardous_12m==0), na.rm=TRUE)
CDIpop2 <- svytotal(~ones, subset(d_boot_EG_CDI_2018,  agchild_f==1 & hazardous_12m==0), na.rm=TRUE)
TOTpop2 <- svytotal(~ones, subset(d_boot_EG_Ghana_CDI_2018, agchild_f==1 & hazardous_12m==0), na.rm=TRUE)


cat('\n Ghana hazardous population: ', Gpop1, '\n CDI hazardous population: ', CDIpop1, 
    '\n Total hazardous population: ', TOTpop1)
cat('\n Ghana not hazardous population: ', Gpop2, '\n CDI not hazardouspopulation: ', CDIpop2,
    '\n Total not hazardous population: ', TOTpop2)

# Percentages
vars <- c('F2_1_ch_m', 'F2_2_ch_m', 'F2_3_ch_m', 'F2_4_ch_m', 'F2_5_ch_m',
          'F2_6_ch_m', 'F2_7_ch_m', 'F2_8_ch_m', 'F2_9_ch_m')
for(x in vars){
  cat('\n Variable: ', x)
  
  #In Hazardous work
  G1 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, agchild_f==1 & hazardous_12m==1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, agchild_f==1 & hazardous_12m==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, agchild_f==1 & hazardous_12m==1), na.rm=TRUE)
  
  #Not in Hazardous work
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, agchild_f==1 & hazardous_12m==0), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, agchild_f==1 & hazardous_12m==0), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, agchild_f==1 & hazardous_12m==0), na.rm=TRUE)
  
  cat('\n Ghana hazardous: ', G1, '\n CDI hazardous: ', CDI1, '\n Total hazardous: ', TOT1)
  cat('\n Ghana non-hazardous: ', G2, '\n CDI non-hazardous: ', CDI2, '\n Total non-hazardous: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 22 ----------------------------------------------------------------
cat('\n Table 22')

# Population
Gpop1 <- svytotal(~ones, subset(d_boot_EG_Ghana_2018, cocoachild_f == 1), na.rm=TRUE)
CDIpop1 <- svytotal(~ones, subset(d_boot_EG_CDI_2018, cocoachild_f == 1), na.rm=TRUE)
TOTpop1 <- svytotal(~ones, subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f == 1), na.rm=TRUE)

cat('\n Ghana cocoa population: ', Gpop1, '\n CDI cocoa population: ', CDIpop1, 
    '\n Total cocoa population: ', TOTpop1)

# Percentages
vars <- c('H8_3_ch', 'H8_1_ch', 'H8_2_ch', 'H8_9_ch', 'H8_5_ch', 'H8_4_ch', 
          'H8_6_ch', 'H8_7_ch', 'H8_8_ch')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2018
  G1 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoachild_f == 1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoachild_f == 1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f == 1), na.rm=TRUE)
  
  cat('\n Ghana: ', G1, '\n CDI: ', CDI1, '\n Total: ', TOT1)
  
}

# Table 23  ----------------------------------------------------------------
cat('\n Table 23')
vars <- c('noncocoachild_f', 'nonagchild_f')
for(x in vars){
  cat('\n Variable: ', x)
  
  #2008
  G1 <- svymean(~get(x), d_unequal_Ghana_2008, na.rm=TRUE)
  CDI1 <- svymean(~get(x), d_unequal_CDI_2008, na.rm=TRUE)
  TOT1 <- svymean(~get(x), d_unequal_Ghana_CDI_2008, na.rm=TRUE)
  
  #2018
  G2 <- svymean(~get(x), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 24 ----------------------------------------------------------------

cat('\n Table 24')
vars <- c('cocoachild_f','cocoachild_7d')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2013
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1), na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1), na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 25 ---------------------------------------------------------------
cat('\n Table 25')
vars <- c('cocoachild_f','cl_cocoa', 'haz_cocoa')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2013
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1), na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1), na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 26 ----------------------------------------------------------------
cat('\n Table 26')
vars <- c('V1a', 'V2a', 'V3a', 'V4a', 'V5a', 'V6a', 'haz_cocoa', 'Via')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2013
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1), na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1), na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 27 ----------------------------------------------------------------

cat('\n Table 27')
vars <- c('V3_12m', 'spray_cocoa', 'present_cocoa', 'enter_cocoa', 'water_cocoa', 'agrochem_cocoa')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2013
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 28 ----------------------------------------------------------------
cat('\n Table 28')
filter_vars <- c('all','Gender2', 'Gender1', 'age_group1', 'age_group2', 'age_group3')
for(y in filter_vars){
  cat('\n Variable: ', y)
  
  #2008
  G1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_2008, get(y)==1), na.rm=TRUE)
  CDI1 <-  svymean(~SchoolAttendance, subset(d_unequal_CDI_2008, get(y)==1), na.rm=TRUE)
  TOT1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_CDI_2008, get(y)==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~SchoolAttendance, subset(d_boot_EG_Ghana_2018, get(y)==1), na.rm=TRUE)
  CDI2 <- svymean(~SchoolAttendance, subset(d_boot_EG_CDI_2018, get(y)==1), na.rm=TRUE)
  TOT2 <- svymean(~SchoolAttendance, subset(d_boot_EG_Ghana_CDI_2018, get(y)==1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
  
}

# Table 29 ----------------------------------------------------------------
cat('\n Table 29')
filter_vars <- c('all','Gender2', 'Gender1', 'age_group1', 'age_group2', 'age_group3')
for(y in filter_vars){
  cat('\n Variable: ', y)
  #2008
  G1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_2008, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  CDI1 <-  svymean(~SchoolAttendance, subset(d_unequal_CDI_2008, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  TOT1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_CDI_2008, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  
  #2018
  G2 <- svymean(~SchoolAttendance, subset(d_boot_EG_Ghana_2018, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  CDI2 <- svymean(~SchoolAttendance, subset(d_boot_EG_CDI_2018, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  TOT2 <- svymean(~SchoolAttendance, subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f==1 & get(y)==1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 30  ----------------------------------------------------------------
cat('\n Table 30')
age_vars <- c('age_group1', 'age_group2', 'age_group3')
vars <- c('cl_cocoa', 'haz_cocoa')

for(y in age_vars){
  for (x in vars){
    cat('\n Variable: ', x)
    cat('\n filter variable: ', y)
    
    #2008
    G1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_2008, get(x)==1 & get(y)==1), na.rm=TRUE)
    CDI1 <-  svymean(~SchoolAttendance, subset(d_unequal_CDI_2008, get(x)==1 & get(y)==1), na.rm=TRUE)
    TOT1 <-  svymean(~SchoolAttendance, subset(d_unequal_Ghana_CDI_2008, get(x)==1 & get(y)==1), na.rm=TRUE)
    
    #2018
    G2 <- svymean(~SchoolAttendance, subset(d_boot_EG_Ghana_2018, get(x)==1 & get(y)==1), na.rm=TRUE)
    CDI2 <- svymean(~SchoolAttendance, subset(d_boot_EG_CDI_2018, get(x)==1 & get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~SchoolAttendance, subset(d_boot_EG_Ghana_CDI_2018, get(x)==1 & get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval)
  }
}

# Table 31 ----------------------------------------------------------------
cat('\n Table 31')

vars <- c('read', 'write', 'math')
for(x in vars){
  cat('\n Variable: ', x)
  
  # Cocoa Household
  G1 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018,  cocoa_hh==1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_boot_EG_CDI_2018,  cocoa_hh==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018,  cocoa_hh==1), na.rm=TRUE)
  
  # Population
  Gpop1 <- svytotal(~get(x), subset(d_boot_EG_Ghana_2018,cocoa_hh==1), na.rm=TRUE)
  CDIpop1 <- svytotal(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  TOTpop1 <- svytotal(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1), na.rm=TRUE)
  
  # Non Cocoa Household
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==0), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==0), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==0), na.rm=TRUE)
  
  # Population
  Gpop2 <- svytotal(~get(x), subset(d_boot_EG_Ghana_2018,cocoa_hh==0), na.rm=TRUE)
  CDIpop2 <- svytotal(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==0), na.rm=TRUE)
  TOTpop2 <- svytotal(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==0), na.rm=TRUE)
  
  cat('\n Ghana cocoa: ', G1, '\n CDI cocoa: ', CDI1, '\n Total cocoa: ', TOT1)
  cat('\n Ghana cocoa population: ', Gpop1, '\n CDI cocoa population: ', CDIpop1, 
      '\n Total cocoa population: ', TOTpop1)
  
  cat('\n Ghana noncocoa: ', G2, '\n CDI noncocoa: ', CDI2, '\n Total noncocoa: ', TOT2)
  cat('\n Ghana non-cocoa population: ', Gpop2, '\n CDI non-cocoa population: ', CDIpop2, 
      '\n Total non-cocoa population: ', TOTpop2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 32 ----------------------------------------------------------------
cat('\n Table 32')
vars <- c('NoSchool', 'InterferenceSchool', 'V7V8')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2008
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008,  cocoachild_f==1 & (age_range_ch==1 | age_range_ch==2)), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008,  cocoachild_f==1 & (age_range_ch==1 | age_range_ch==2)), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008,  cocoachild_f==1 & (age_range_ch==1 | age_range_ch==2)), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoachild_f==1 & (age_range_ch==1 | age_range_ch==2)), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoachild_f==1 & (age_range_ch==1 | age_range_ch==2)), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoachild_f==1 & (age_range_ch==1 | age_range_ch==2)), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 33  ----------------------------------------------------------------
cat('\n Table 33')
filter_vars <- c('stratum1', 'stratum2', 'stratum3', 'all')
vars <- c('cl_cocoa', 'haz_cocoa')
for(y in filter_vars){
  for (x in vars){
    cat('\n Variable: ', x)
    #2008
    G1 <-  svymean(~get(x), subset(d_unequal_Ghana_2008, get(y)==1), na.rm=TRUE)
    CDI1 <-  svymean(~get(x), subset(d_unequal_CDI_2008, get(y)==1), na.rm=TRUE)
    TOT1 <-  svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, get(y)==1), na.rm=TRUE)
    
    #2018
    G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, get(y)==1), na.rm=TRUE)
    CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval)
  }
}

# Table 34 ---------------------------------------------------

#####fert_exp_ton_p10
###Total
#2013/18
svyby(~fert_exp_ton_p10, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean, na.rm=TRUE)
svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~fert_exp_ton_p10, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

###CDI
#2013/18
svyby(~fert_exp_ton_p10, ~noncocoa, d_unequal_CDI_2013_HH, svymean, na.rm=TRUE)
svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~fert_exp_ton_p10, ~noncocoa, d_unequal_CDI_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

###Ghana
#2013/18
svyby(~fert_exp_ton_p10, ~noncocoa, d_unequal_Ghana_2013_HH, svymean, na.rm=TRUE)
svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~fert_exp_ton_p10, ~noncocoa, d_unequal_Ghana_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~fert_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

#####pest_exp_ton_p10
###Total
#2013/18
svyby(~pest_exp_ton_p10, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean, na.rm=TRUE)
svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~pest_exp_ton_p10, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

###CDI
#2013/18
svyby(~pest_exp_ton_p10, ~noncocoa, d_unequal_CDI_2013_HH, svymean, na.rm=TRUE)
svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~pest_exp_ton_p10, ~noncocoa, d_unequal_CDI_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

###Ghana
#2013/18
svyby(~pest_exp_ton_p10, ~noncocoa, d_unequal_Ghana_2013_HH, svymean, na.rm=TRUE)
svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~pest_exp_ton_p10, ~noncocoa, d_unequal_Ghana_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~pest_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

#####herb_exp_ton_p10
###Total
#2013/18
svyby(~herb_exp_ton_p10, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean, na.rm=TRUE)
svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~herb_exp_ton_p10, ~noncocoa, d_unequal_Ghana_CDI_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_CDI_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

###CDI
#2013/18
svyby(~herb_exp_ton_p10, ~noncocoa, d_unequal_CDI_2013_HH, svymean, na.rm=TRUE)
svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~herb_exp_ton_p10, ~noncocoa, d_unequal_CDI_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_CDI_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

###Ghana
#2013/18
svyby(~herb_exp_ton_p10, ~noncocoa, d_unequal_Ghana_2013_HH, svymean, na.rm=TRUE)
svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE)
Z = abs(svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE)[1,2] - svyby(~herb_exp_ton_p10, ~noncocoa, d_unequal_Ghana_2013_HH, svymean, na.rm=TRUE)[1,2]) / 
  (sqrt(2)*SE(svyby(~herb_exp_ton_p10, ~noncocoa, d_boot_EG_Ghana_2018_HH, svymean, na.rm=TRUE))[1])
1 - pnorm(Z)

# Table 35 from roster survey, so does not need to be included

# Tables 36-46 do not require significance tests

# Table 47 ----------------------------------------------------------------
cat('\n Table 47')
ages <- c('all', 'age_group1', 'age_group2', 'age_group3')
for(x in ages){
  cat('\n Variable: ', x)
  
  # Attending School
  TOT1 <-svymean(~child_labor_12m, subset(d_boot_EG_Ghana_CDI_2018, SchoolAttendance==1 & get(x)==1), na.rm=TRUE)
  
  # Not Attending School
  TOT2 <-svymean(~child_labor_12m, subset(d_boot_EG_Ghana_CDI_2018, SchoolAttendance==0 & get(x)==1), na.rm=TRUE)
  
  cat( '\n Attending School: ', TOT1)
  cat( '\n Not Attending School: ', TOT2)
  
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Total p value: ', totpval)
  
}

# Table 48 ----------------------------------------------------------------
cat('\n Table 48')
gend_vars <- c('Gender2', 'Gender1')
vars <- c('age_group1', 'age_group2', 'age_group3', 'all')
for(y in gend_vars){
  for (x in vars){
    cat('\n Percentages Variable: ', x)
    cat('\n Gender Variable: ', y)
    
    #Percentages    
    #2008
    G1 <-  svymean(~vio_max_cocoa, subset(d_unequal_Ghana_2008, get(x)==1 & get(y)==1), na.rm=TRUE)
    CDI1 <-  svymean(~vio_max_cocoa, subset(d_unequal_CDI_2008, get(x)==1 & get(y)==1), na.rm=TRUE)
    TOT1 <-  svymean(~vio_max_cocoa, subset(d_unequal_Ghana_CDI_2008, get(x)==1 & get(y)==1), na.rm=TRUE)
    
    #2018
    G2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_Ghana_2018, get(x)==1 & get(y)==1), na.rm=TRUE)
    CDI2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_CDI_2018, get(x)==1 & get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_Ghana_CDI_2018, get(x)==1 & get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval)
    
    #Total Hours
    cat('\n Total Hours variable: ', x)
    
    #2008
    G1 <-  svymean(~TotalHours, subset(d_unequal_Ghana_2008, get(x)==1 & get(y)==1 & cocoachild_f == 1), na.rm=TRUE)
    CDI1 <-  svymean(~TotalHours, subset(d_unequal_CDI_2008, get(x)==1 & get(y)==1 & cocoachild_f == 1), na.rm=TRUE)
    TOT1 <-  svymean(~TotalHours, subset(d_unequal_Ghana_CDI_2008, get(x)==1 & get(y)==1 & cocoachild_f == 1), na.rm=TRUE)
    
    #2018
    G2 <- svymean(~TotalHours, subset(d_boot_EG_Ghana_2018, get(x)==1 & get(y)==1 & cocoachild_f == 1), na.rm=TRUE)
    CDI2 <- svymean(~TotalHours, subset(d_boot_EG_CDI_2018, get(x)==1 & get(y)==1 & cocoachild_f == 1), na.rm=TRUE)
    TOT2 <- svymean(~TotalHours, subset(d_boot_EG_Ghana_CDI_2018, get(x)==1 & get(y)==1 & cocoachild_f == 1), na.rm=TRUE)
    
    cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval)
  }
}

# Table 49  ----------------------------------------------------------------
cat('\n Table 49')
vars <- c('light_cocoa', 'legal_cocoa')
for(x in vars){
  cat('\n Variable: ', x)
  
  # 2008
  G1 <- svymean(~get(x), d_unequal_Ghana_2008, na.rm=TRUE)
  CDI1 <- svymean(~get(x), d_unequal_CDI_2008, na.rm=TRUE)
  TOT1 <- svymean(~get(x), d_unequal_Ghana_CDI_2008, na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 50 ----------------------------------------------------------------
cat('\n Table 50')
vars <- c('Cocoa12Mo_LandClearing_m', 'Cocoa12Mo_FellingChopping_m', 'Cocoa12Mo_Burning_m',
          'Cocoa12Mo_Stumping_m', 'Cocoa12Mo_PreparationSeedlings_m', 'Cocoa12Mo_HolingSeedlings_m',
          'Cocoa12Mo_SowingStake_m', 'Cocoa12Mo_Weeding_m', 'Cocoa12Mo_SprayingInsecticides_m',
          'Cocoa12Mo_Fertilizer_m', 'Cocoa12Mo_Fungicides_m', 'Cocoa12Mo_WaterCarrying_m',
          'Cocoa12Mo_SanitationPruning_m', 'Cocoa12Mo_MistletonControl_m', 
          'Cocoa12Mo_PluckingCocoa_m', 'Cocoa12Mo_GatheringPods_m', 'Cocoa12Mo_PodBreaking_m',
          'Cocoa12Mo_CartingBeans_m', 'Cocoa12Mo_DryingBeans_m', 'Cocoa12Mo_CartingBeansShed_m')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2008
  G1 <- svymean(~get(x), d_unequal_Ghana_2008, na.rm=TRUE)
  CDI1 <- svymean(~get(x), d_unequal_CDI_2008, na.rm=TRUE)
  TOT1 <- svymean(~get(x), d_unequal_Ghana_CDI_2008, na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 51  ----------------------------------------------------------------
cat('\n Table 51')
filter_vars <- c('cl_cocoa', 'haz_cocoa')
vars <- c('Gender2', 'Gender1', 'age_group1', 'age_group2', 'age_group3')

for(y in filter_vars){
  for(x in vars){
    cat('\n Variable: ', x)
    cat('\n Filter variable: ', y)
    
    # 2008
    G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, get(y)==1), na.rm=TRUE)
    CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008, get(y)==1), na.rm=TRUE)
    TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, get(y)==1), na.rm=TRUE)
    CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval)
  }
}

# Table 52 ---------------------------------------------------------------
cat('\n Table 52')
filter_vars <- c('age_group1', 'age_group2', 'age_group3')
vars <- c('V1a', 'V2a', 'V3a', 'V4a', 'V5a', 'V6a', 'haz_cocoa', 'Via')

for(y in filter_vars){
  for(x in vars){
    cat('\n Variable: ', x)
    cat('\n Filter variable: ', y)
    
    # 2008
    G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, get(y)==1), na.rm=TRUE)
    CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008, get(y)==1), na.rm=TRUE)
    TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, get(y)==1), na.rm=TRUE)
    CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval)
  }
}

# Table 53 ----------------------------------------------------------------
cat('\n Table 53')
filter_vars <- c('Gender2', 'Gender1')
vars <- c('V1a', 'V2a', 'V3a', 'V4a', 'V5a', 'V6a', 'haz_cocoa', 'Via')

for(y in filter_vars){
  for(x in vars){
    cat('\n Variable: ', x)
    cat('\n Filter variable: ', y)
    
    # 2008
    G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, get(y)==1), na.rm=TRUE)
    CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008, get(y)==1), na.rm=TRUE)
    TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, get(y)==1), na.rm=TRUE)
    CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval)
  }
}

# Table 54 ----------------------------------------------------------------
cat('\n Table 54')
vars <- c( 'age_group1', 'age_group2', 'age_group3', 'Gender2', 'Gender1')
for(x in vars){
  cat('\n Variable: ', x)
  
  # 2008
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2008, V3_cocoa == 1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2008, V3_cocoa == 1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2008, V3_cocoa == 1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, V3_cocoa == 1), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, V3_cocoa == 1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, V3_cocoa == 1), na.rm=TRUE)
  
  cat('\n Ghana 2008: ', G1, '\n CDI 2008: ', CDI1, '\n Total 2008: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 55 ----------------------------------------------------------------
cat('\n Table 55')
vars <- c('child_labor_12m', 'hazardous_12m')
for(x in vars){
  cat('\n Variable: ', x)
  
  #Percentages
  # Cocoa
  G1 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh == 1), na.rm=TRUE)
  CDI1 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh == 1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh == 1), na.rm=TRUE)
  
  # non-cocoa
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh == 0), na.rm=TRUE)
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh == 0), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh == 0), na.rm=TRUE)
  
  #Population
  # Cocoa 
  Gpop1 <- svytotal(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh == 1), na.rm=TRUE)
  CDIpop1 <- svytotal(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh == 1), na.rm=TRUE)
  TOTpop1 <- svytotal(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh == 1), na.rm=TRUE)
  
  # Non-Cocoa 
  Gpop2 <- svytotal(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh == 0), na.rm=TRUE)
  CDIpop2 <- svytotal(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh == 0), na.rm=TRUE)
  TOTpop2 <- svytotal(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh == 0), na.rm=TRUE)
  
  cat('\n Ghana cocoa: ', G1, '\n CDI cocoa: ', CDI1, '\n Total cocoa: ', TOT1)
  cat('\n Ghana non-cocoa: ', G2, '\n CDI non-cocoa: ', CDI2, '\n Total non-cocoa: ', TOT2)
  cat('\n Ghana cocoa population: ', Gpop1, '\n CDI cocoa population: ', CDIpop1, 
      '\n Total cocoa population: ', TOTpop1)
  cat('\n Ghana noncocoa population: ', Gpop2, '\n CDI noncocoa population: ', CDIpop2,
      '\n Total noncocoa population: ', TOTpop2)
  
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 56 ----------------------------------------------------------------
cat('\n Table 56')
filter_vars <- c('age_group1', 'age_group2', 'age_group3', 'all')

for(y in filter_vars){
  cat('\n Variable: ', y)
  
  # 2013
  G1 <- svymean(~vio_max_cocoa, subset(d_unequal_Ghana_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
  CDI1 <- svymean(~vio_max_cocoa, subset(d_unequal_CDI_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
  CDI2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_CDI_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval  )
  
  #Average Hours Worked
  cat('\n Hours Worked')
  
  G1 <- svymean(~TotalHours, subset(d_unequal_Ghana_2013, cocoachild_f==1 & cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
  CDI1 <- svymean(~TotalHours, subset(d_unequal_CDI_2013,cocoachild_f==1 & cocoa_hh==1 & get(y)==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~TotalHours, subset(d_boot_EG_Ghana_2018, cocoachild_f==1 & cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
  CDI2 <- svymean(~TotalHours, subset(d_boot_EG_CDI_2018, cocoachild_f==1 & cocoa_hh==1 & get(y)==1), na.rm=TRUE)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval  )
}

# Table 57 ----------------------------------------------------------------
cat('\n Table 57')

gend_vars <- c('Gender2', 'Gender1')
filter_vars <- c('age_group1', 'age_group2', 'age_group3', 'all')

for(x in gend_vars){
  for(y in filter_vars){
    cat('\n Variable: ', x)
    cat('\n filter by: ', y)
    
    # 2013
    G1 <- svymean(~vio_max_cocoa, subset(d_unequal_Ghana_2013, cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE) 
    CDI1 <- svymean(~vio_max_cocoa, subset(d_unequal_CDI_2013, cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE) 
    CDI2 <- svymean(~vio_max_cocoa, subset(d_boot_EG_CDI_2018, cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2)
    
    # Hypothesis Testing
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval  )
    
    # Hours Worked
    cat('\n Hours worked')
    #2008
    G1 <- svymean(~TotalHours, subset(d_unequal_Ghana_2013, cocoachild_f == 1 & cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE) 
    CDI1 <- svymean(~TotalHours, subset(d_unequal_CDI_2013, cocoachild_f == 1 & cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~TotalHours, subset(d_boot_EG_Ghana_2018, cocoachild_f == 1 & cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE) 
    CDI2 <- svymean(~TotalHours, subset(d_boot_EG_CDI_2018, cocoachild_f == 1 & cocoa_hh==1 & get(x)==1 & get(y)==1), na.rm=TRUE)
    
    # Hypothesis Testing
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval  )
    
  }
}

# Table 58 ----------------------------------------------------------------
cat('\n Table 58')  

#Sub-table 1
vars <-c('Gender2', 'Gender1','age_group1','age_group2','age_group3')
x <-'age_group1'
for(x in vars){
  cat('\n Variable: ', x)
  
  # 2013
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1 & cocoachild_f == 1 & cl_cocoa==1), na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1 & cocoachild_f==1 & cl_cocoa==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1 & cocoachild_f==1 & cl_cocoa==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & cocoachild_f==1 & cl_cocoa==1), na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1 & cocoachild_f==1 & cl_cocoa==1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1 & cocoachild_f==1 & cl_cocoa==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  #Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval  )
}

#Sub-table2
vars <-c('Gender2', 'Gender1','age_group1','age_group2','age_group3')

for(x in vars){
  cat('\n \n Variable: ', x)
  
  # 2013
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1 & cocoachild_f == 1 & haz_cocoa==1), na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1 & cocoachild_f==1 & haz_cocoa==1), na.rm=TRUE)
  TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1 & cocoachild_f==1 & haz_cocoa==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & cocoachild_f==1 & haz_cocoa==1), na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1 & cocoachild_f==1 & haz_cocoa==1), na.rm=TRUE)
  TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1 & cocoachild_f==1 & haz_cocoa==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval  )
}

# Table 59 ----------------------------------------------------------------
cat('\n Table 59')
filter_vars <- c('age_group1', 'age_group2', 'age_group3')
vars <- c('V1a', 'V2a', 'V3a', 'V4a', 'V5a', 'V6a', 'haz_cocoa', 'Via')

for(y in filter_vars){
  for(x in vars){
    cat('\n Variable: ', x)
    cat('\n Filter variable: ', y)
    
    # 2013
    G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
    CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
    CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Hypothesis Testing
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval  )
  }
}

# Table 60 ----------------------------------------------------------------
cat('\n Table 60')
filter_vars <- c('Gender2', 'Gender1')
vars <- c('V1a', 'V2a', 'V3a', 'V4a', 'V5a', 'V6a', 'haz_cocoa', 'Via')

for(y in filter_vars){
  for(x in vars){
    cat('\n Variable: ', x)
    cat('\n Filter variable: ', y)
    
    # 2013
    G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
    CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE) 
    CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1 & get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Hypothesis Testing
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval  )
  }
}

# Table 61 ----------------------------------------------------------------
cat('\n Table 61')  
vars <-c('V_7', 'V_6', 'V_5', 'V_4', 'V_3', 'V_2', 'V_1', 'Vi')

for(x in vars){
  cat('\n Variable: ', x)
  
  # 2013
  G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1 & cocoachild_f == 1), na.rm=TRUE) 
  CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE)
  
  # 2018
  G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018,  cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE) 
  CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018,  cocoa_hh==1 & cocoachild_f==1), na.rm=TRUE)
  
  cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1)
  cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2)
  
  # Hypothesis Testing
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval  )
}

# Table 62 ----------------------------------------------------------------
cat('\n Table 62')
filter_vars <- c('spray_cocoa', 'present_cocoa', 'enter_cocoa',
                 'water_cocoa', 'agrochem_cocoa')
vars <- c('age_group1', 'age_group2', 'age_group3', 'Gender2', 'Gender1')

for(y in filter_vars){
  for(x in vars){
    cat('\n Variable: ', x)
    cat('\n Filter variable: ', y)
    
    # 2013
    G1 <- svymean(~get(x), subset(d_unequal_Ghana_2013, cocoa_hh==1 & V3_cocoa == 1 & get(y)==1), na.rm=TRUE) 
    CDI1 <- svymean(~get(x), subset(d_unequal_CDI_2013, cocoa_hh==1 & V3_cocoa == 1 & get(y)==1), na.rm=TRUE)
    TOT1 <- svymean(~get(x), subset(d_unequal_Ghana_CDI_2013, cocoa_hh==1 & V3_cocoa == 1 & get(y)==1), na.rm=TRUE)
    
    # 2018
    G2 <- svymean(~get(x), subset(d_boot_EG_Ghana_2018, cocoa_hh==1 & V3_cocoa == 1 & get(y)==1), na.rm=TRUE) 
    CDI2 <- svymean(~get(x), subset(d_boot_EG_CDI_2018, cocoa_hh==1 & V3_cocoa == 1 & get(y)==1), na.rm=TRUE)
    TOT2 <- svymean(~get(x), subset(d_boot_EG_Ghana_CDI_2018, cocoa_hh==1 & V3_cocoa == 1 & get(y)==1), na.rm=TRUE)
    
    cat('\n Ghana 2013: ', G1, '\n CDI 2013: ', CDI1, '\n Total 2013: ', TOT1)
    cat('\n Ghana 2018: ', G2, '\n CDI 2018: ', CDI2, '\n Total 2018: ', TOT2)
    
    # Hypothesis Testing
    
    # Ghana Difference
    Z = abs(G1[1] - G2[1]) / 
      (sqrt(2)*SE(G2))
    gpval <- 1 - pnorm(Z)
    
    # CDI Difference
    Z = abs(CDI1[1] - CDI2[1]) / 
      (sqrt(2)*SE(CDI2))
    cdipval <- 1 - pnorm(Z)
    
    # Total Difference
    Z = abs(TOT1[1] - TOT2[1]) / 
      (sqrt(2)*SE(TOT2))
    totpval <- 1 - pnorm(Z)
    
    cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
        '\n Total p value: ', totpval  )
  }
}

# Table 63 has no significance tests

# Table 64 ----------------------------------------------------------------
cat('\n Table 64')
vars <- list(c('agchild_7d' ,'agchild_12m'), c('cl_ag1_7d', 'cl_ag1_12m'), c('haz_ag1_7d', 'haz_ag1_12m'))

for(x in vars){
  cat('\n Variable: ', x)
  
  
  # 7 days population
  Gpop1 <- svytotal(~get(x[[1]]), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDIpop1 <- svytotal(~get(x[[1]]), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOTpop1 <- svytotal(~get(x[[1]]), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  # 12 months population
  Gpop2 <- svytotal(~get(x[[2]]), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDIpop2 <- svytotal(~get(x[[2]]), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOTpop2 <- svytotal(~get(x[[2]]), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  # All Children population
  GFULLPOP <- Gpop1 + Gpop2
  CDIFULLPOP <- CDIpop1 + CDIpop2
  TOTFULLPOP <- TOTpop1 + TOTpop2
  
  # 7 days
  G1 <- svymean(~get(x[[1]]), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI1 <- svymean(~get(x[[1]]), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT1 <- svymean(~get(x[[1]]), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  # 12 months
  G2 <- svymean(~get(x[[2]]), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x[[2]]), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x[[2]]), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 7 days: ', G1, '\n CDI 7 days: ', CDI1, '\n Total 7 days: ', TOT1)
  cat('\n Ghana 12 months: ', G2, '\n CDI 12 months: ', CDI2, '\n Total 12 months: ', TOT2)
  
  cat('\n Ghana 7 days population: ', Gpop1, '\n CDI 7 days population: ', CDIpop1, 
      '\n Total 7 days population: ', TOTpop1)
  cat('\n Ghana 12 months population: ', Gpop2, '\n CDI 12 months population: ', CDIpop2,
      '\n Total 12 months population: ', TOTpop2)
  cat('\n Ghana All population: ', GFULLPOP, '\n CDI All population: ', CDIFULLPOP,
      '\n Total All population: ', TOTFULLPOP)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Table 65 ----------------------------------------------------------------
cat('\n Table 65')

# Population
Gpop1 <- svytotal(~ones, d_boot_EG_Ghana_2018, na.rm=TRUE)
CDIpop1 <- svytotal(~ones, d_boot_EG_CDI_2018, na.rm=TRUE)
TOTpop1 <- svytotal(~ones, d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)

cat('\n Ghana population: ', Gpop1, '\n CDI population: ', CDIpop1, 
    '\n Total population: ', TOTpop1)

vars <- list(c('V1_7d',  'V1_12m'), c('V2_7d', 'V2_12m'), c('V3_7d','V3_12m'),
             c('V4_7d', 'V4_12m'), c('V5_7d', 'V5_12m'), c('V6_7d','V6_12m'),
             c('hazardous_7d', 'hazardous_12m'), c('Vi_7d', 'Vi'))

for(x in vars){
  cat('\n Variable: ', x)
  
  # 7 days
  G1 <- svymean(~get(x[[1]]), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI1 <- svymean(~get(x[[1]]), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT1 <- svymean(~get(x[[1]]), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  # 12 months
  G2 <- svymean(~get(x[[2]]), d_boot_EG_Ghana_2018, na.rm=TRUE)
  CDI2 <- svymean(~get(x[[2]]), d_boot_EG_CDI_2018, na.rm=TRUE)
  TOT2 <- svymean(~get(x[[2]]), d_boot_EG_Ghana_CDI_2018, na.rm=TRUE)
  
  cat('\n Ghana 7 days: ', G1, '\n CDI 7 days: ', CDI1, '\n Total 7 days: ', TOT1)
  cat('\n Ghana 12 months: ', G2, '\n CDI 12 months: ', CDI2, '\n Total 12 months: ', TOT2)
  
  # Ghana Difference
  Z = abs(G1[1] - G2[1]) / 
    (sqrt(2)*SE(G2))
  gpval <- 1 - pnorm(Z)
  
  # CDI Difference
  Z = abs(CDI1[1] - CDI2[1]) / 
    (sqrt(2)*SE(CDI2))
  cdipval <- 1 - pnorm(Z)
  
  # Total Difference
  Z = abs(TOT1[1] - TOT2[1]) / 
    (sqrt(2)*SE(TOT2))
  totpval <- 1 - pnorm(Z)
  
  cat('\n Ghana p value: ', gpval, '\n CDI p value: ', cdipval, 
      '\n Total p value: ', totpval)
}

# Tables 66-69 don't have significance tests

# Table 70 ----------------------------------------------------------------

###CDI
#Education of Mother
svyby(~cl_cocoa, ~mother_educ_sec_mo_ro, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~mother_educ_sec_mo_ro, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~mother_educ_sec_mo_ro, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~mother_educ_sec_mo_ro, d_boot_EG_CDI_2018)

#Education of Father
svyby(~cl_cocoa, ~father_educ_sec_mo_ro, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~father_educ_sec_mo_ro, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~father_educ_sec_mo_ro, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~father_educ_sec_mo_ro, d_boot_EG_CDI_2018)

#Income Group
svyby(~cl_cocoa, ~income_quintile_45_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~income_quintile_45_hhh, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~income_quintile_45_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~income_quintile_45_hhh, d_boot_EG_CDI_2018)

#Area Under Cocoa Cultivation
svyby(~cl_cocoa, ~cult_cocoa_5_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~cult_cocoa_5_hhh, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~cult_cocoa_5_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~cult_cocoa_5_hhh, d_boot_EG_CDI_2018)

#Household Heads' Awareness of Child Labor
svyby(~cl_cocoa, ~aware_cl_reg_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~aware_cl_reg_hhh, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~aware_cl_reg_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~aware_cl_reg_hhh, d_boot_EG_CDI_2018)

#Migration: Head ever migrated to and from another village
svyby(~cl_cocoa, ~head_migrated_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~head_migrated_hhh, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~head_migrated_hhh, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~head_migrated_hhh, d_boot_EG_CDI_2018)

###Ghana
#Education of Mother
svyby(~cl_cocoa, ~mother_educ_sec_mo_ro, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~mother_educ_sec_mo_ro, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~mother_educ_sec_mo_ro, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~mother_educ_sec_mo_ro, d_boot_EG_Ghana_2018)

#Education of Father
svyby(~cl_cocoa, ~father_educ_sec_mo_ro, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~father_educ_sec_mo_ro, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~father_educ_sec_mo_ro, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~father_educ_sec_mo_ro, d_boot_EG_Ghana_2018)

#Income Group
svyby(~cl_cocoa, ~income_quintile_45_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~income_quintile_45_hhh, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~income_quintile_45_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~income_quintile_45_hhh, d_boot_EG_Ghana_2018)

#Area Under Cocoa Cultivation
svyby(~cl_cocoa, ~cult_cocoa_5_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~cult_cocoa_5_hhh, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~cult_cocoa_5_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~cult_cocoa_5_hhh, d_boot_EG_Ghana_2018)

#Household Heads' Awareness of Child Labor
svyby(~cl_cocoa, ~aware_cl_reg_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~aware_cl_reg_hhh, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~aware_cl_reg_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~aware_cl_reg_hhh, d_boot_EG_Ghana_2018)

#Migration: Head ever migrated to and from another village
svyby(~cl_cocoa, ~head_migrated_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~head_migrated_hhh, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~head_migrated_hhh, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~head_migrated_hhh, d_boot_EG_Ghana_2018)

# Table 71 ----------------------------------------------------------------

###CDI
#Primary School Location
svyby(~cl_cocoa, ~dist_pri_sch_1km_cl, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~dist_pri_sch_1km_cl, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~dist_pri_sch_1km_cl, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~dist_pri_sch_1km_cl, d_boot_EG_CDI_2018)

#Junior High School (JHS) Location
svyby(~cl_cocoa, ~dist_jhs_sch_2km_cl, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~dist_jhs_sch_2km_cl, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~dist_jhs_sch_2km_cl, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~dist_jhs_sch_2km_cl, d_boot_EG_CDI_2018)

#Senior High School (SHS) Location
svyby(~cl_cocoa, ~dist_shs_sch_5km_cl, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~dist_shs_sch_5km_cl, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~dist_shs_sch_5km_cl, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~dist_shs_sch_5km_cl, d_boot_EG_CDI_2018)

#Improved Road
svyby(~cl_cocoa, ~roads_imp_cl, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~roads_imp_cl, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~roads_imp_cl, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~roads_imp_cl, d_boot_EG_CDI_2018)

#Grid Electricity Available
svyby(~cl_cocoa, ~grid_elec_comm_cl_dum, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~grid_elec_comm_cl_dum, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~grid_elec_comm_cl_dum, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~grid_elec_comm_cl_dum, d_boot_EG_CDI_2018)

#Micro-finance Institution
svyby(~cl_cocoa, ~num_mfi1_cl_dum, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~num_mfi1_cl_dum, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~num_mfi1_cl_dum, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~num_mfi1_cl_dum, d_boot_EG_CDI_2018)

###Ghana
#Primary School Location
svyby(~cl_cocoa, ~dist_pri_sch_1km_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~dist_pri_sch_1km_cl, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~dist_pri_sch_1km_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~dist_pri_sch_1km_cl, d_boot_EG_Ghana_2018)

#Junior High School (JHS) Location
svyby(~cl_cocoa, ~dist_jhs_sch_2km_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~dist_jhs_sch_2km_cl, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~dist_jhs_sch_2km_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~dist_jhs_sch_2km_cl, d_boot_EG_Ghana_2018)

#Senior High School (SHS) Location
svyby(~cl_cocoa, ~dist_shs_sch_5km_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~dist_shs_sch_5km_cl, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~dist_shs_sch_5km_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~dist_shs_sch_5km_cl, d_boot_EG_Ghana_2018)

#Improved Road
svyby(~cl_cocoa, ~roads_imp_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~roads_imp_cl, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~roads_imp_cl, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~roads_imp_cl, d_boot_EG_Ghana_2018)

#Grid Electricity Available
svyby(~cl_cocoa, ~grid_elec_comm_cl_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~grid_elec_comm_cl_dum, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~grid_elec_comm_cl_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~grid_elec_comm_cl_dum, d_boot_EG_Ghana_2018)

#Micro-finance Institution
svyby(~cl_cocoa, ~num_mfi1_cl_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~num_mfi1_cl_dum, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~num_mfi1_cl_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~num_mfi1_cl_dum, d_boot_EG_Ghana_2018)

# Table 72 ----------------------------------------------------------------

###CDI
#Toilet in School
svyby(~cl_cocoa, ~toilets_inside_sch_dum, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~toilets_inside_sch_dum, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~toilets_inside_sch_dum, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~toilets_inside_sch_dum, d_boot_EG_CDI_2018)

#Piped Drinking Water in School
svyby(~cl_cocoa, ~drinking_water_pipe_sch_dum, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~drinking_water_pipe_sch_dum, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~drinking_water_pipe_sch_dum, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~drinking_water_pipe_sch_dum, d_boot_EG_CDI_2018)

#School has Canteen/Kitchen
svyby(~cl_cocoa, ~kitchen_sch_dum, d_boot_EG_CDI_2018, svymean)
svyttest(cl_cocoa~kitchen_sch_dum, d_boot_EG_CDI_2018)

svyby(~haz_cocoa, ~kitchen_sch_dum, d_boot_EG_CDI_2018, svymean)
svyttest(haz_cocoa~kitchen_sch_dum, d_boot_EG_CDI_2018)

###Ghana
#Toilet in School
svyby(~cl_cocoa, ~toilets_inside_sch_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~toilets_inside_sch_dum, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~toilets_inside_sch_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~toilets_inside_sch_dum, d_boot_EG_Ghana_2018)

#Piped Drinking Water in School
svyby(~cl_cocoa, ~drinking_water_pipe_sch_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~drinking_water_pipe_sch_dum, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~drinking_water_pipe_sch_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~drinking_water_pipe_sch_dum, d_boot_EG_Ghana_2018)

#School has Canteen/Kitchen
svyby(~cl_cocoa, ~kitchen_sch_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(cl_cocoa~kitchen_sch_dum, d_boot_EG_Ghana_2018)

svyby(~haz_cocoa, ~kitchen_sch_dum, d_boot_EG_Ghana_2018, svymean)
svyttest(haz_cocoa~kitchen_sch_dum, d_boot_EG_Ghana_2018)

# Tables 73-85 are regression tables and do not need to be included