#The NORC Ghana/CDI survey analyses: create the R data files to use for significance tests

#Remove any history
rm(list=ls())

#Load the packages required for this analysis
#Please note that the "pacman" package must be installed in R before loading these packages
pacman::p_load(survey, readstata13, readxl)

#Set the working directory to where the data set is saved
#Either forward or backward slashes may be required
#CHANGE DIRECTORY HERE TO YOUR DIRECTION PLUS "/Datasets and Codes"
setwd("C:/Directory/Datasets and Codes")

#Read in the data sets
Full_data_set_child = read.dta13("child_survey_2008_2013_2018_final.dta")
Full_data_set_HH = read.dta13("HHH_all_years_final.dta")
Full_data_set_child_2018 = read.dta13("all_surveys_2018_final.dta")

#The size of the data set
dim(Full_data_set_child)
dim(Full_data_set_HH)
dim(Full_data_set_child_2018)

#View the first few entries
Full_data_set_child[1:10,1:25]
head(Full_data_set_HH)
Full_data_set_child_2018[1:10,1:25]

#Some data cleaning
x = which(Full_data_set_child[,"hhid"] == 7305)
Full_data_set_child[x,]
Full_data_set_child[x,"num_children_roster"] = Full_data_set_child[x,"FPC_Ind_within_HH"] = 3
x = which(Full_data_set_child[,"hhid"] == 10881)
Full_data_set_child[x,]
Full_data_set_child[x,"num_children_roster"] = Full_data_set_child[x,"FPC_Ind_within_HH"] = 4

#Some frequency counts
summary(Full_data_set_child)
table(Full_data_set_child[,"Strata"], useNA="always")
table(Full_data_set_child[,"Year"], useNA="always")
table(Full_data_set_child[,"ghana"], useNA="always")
table(Full_data_set_child[,"cdi"], useNA="always")
table(Full_data_set_child[,"weight"], useNA="always")
table(Full_data_set_child[,"cocoa_weight"], useNA="always")
table(Full_data_set_child[,"EA_2013"], useNA="always")
table(Full_data_set_child[,"EA_2018"], useNA="always")
table(Full_data_set_child[,"cocoa_hh"], useNA="always")
table(Full_data_set_child[,"cocoachild_f"], useNA="always")
table(Full_data_set_child[,"agchild_f"], useNA="always")
table(Full_data_set_child[,"child_labor_12m"], useNA="always")
table(Full_data_set_child[,"hazardous_12m"], useNA="always")
table(Full_data_set_child[,"cl_ag"], useNA="always")
table(Full_data_set_child[,"haz_ag"], useNA="always")
table(Full_data_set_child[,"cl_cocoa"], useNA="always")
table(Full_data_set_child[,"haz_cocoa"], useNA="always")
table(Full_data_set_child[,"num_children_roster"], useNA="always")
table(Full_data_set_child[,"hhid"], useNA="always")

##########################################################################
##########################################################################
#Ghana 2008 study - Child
##########################################################################
##########################################################################

#Subset out the corresponding entries
Ghana_2008_data_set = subset(Full_data_set_child, ghana==1 & Year=="2008")
nrow(Ghana_2008_data_set)

#Add a row of id indicators to declare the design
C = ncol(Ghana_2008_data_set)
Ghana_2008_data_set[,C+1] = 1:nrow(Ghana_2008_data_set)
colnames(Ghana_2008_data_set)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(Ghana_2008_data_set)
Ghana_2008_data_set[,C+1] = 1
colnames(Ghana_2008_data_set)[C+1] = "ones"

#Declare the sampling design
d_unequal_Ghana_2008 = svydesign(ids=~ids, data=Ghana_2008_data_set, weights=~weight)

#Total number of children
svytotal(~ones, d_unequal_Ghana_2008)

#Survey means 
svymean(~cocoachild_f, d_unequal_Ghana_2008)
svymean(~cl_cocoa, d_unequal_Ghana_2008)
svymean(~haz_cocoa, d_unequal_Ghana_2008)

#Survey totals 
svytotal(~cocoachild_f, d_unequal_Ghana_2008)
svytotal(~cl_cocoa, d_unequal_Ghana_2008)
svytotal(~haz_cocoa, d_unequal_Ghana_2008)

##########################################################################
##########################################################################
#Ghana 2008 study - Household
##########################################################################
##########################################################################

#Subset out the corresponding entries
Ghana_2008_data_set_HH = subset(Full_data_set_HH, ghana==1 & Year=="2008")
nrow(Ghana_2008_data_set_HH)

#Add a row of id indicators to declare the design
C = ncol(Ghana_2008_data_set_HH)
Ghana_2008_data_set_HH[,C+1] = 1:nrow(Ghana_2008_data_set_HH)
colnames(Ghana_2008_data_set_HH)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(Ghana_2008_data_set_HH)
Ghana_2008_data_set_HH[,C+1] = 1
colnames(Ghana_2008_data_set_HH)[C+1] = "ones"

#Declare the sampling design
d_unequal_Ghana_2008_HH = svydesign(ids=~ids, data=Ghana_2008_data_set_HH, weights=~weight)

#Total number of households
svytotal(~ones, d_unequal_Ghana_2008_HH)

##########################################################################
##########################################################################
#CDI 2008 study - Child
##########################################################################
##########################################################################

CDI_2008_data_set = subset(Full_data_set_child, cdi==1 & Year=="2008")
nrow(CDI_2008_data_set)

#Add a row of id indicators to declare the design
C = ncol(CDI_2008_data_set)
CDI_2008_data_set[,C+1] = 1:nrow(CDI_2008_data_set)
colnames(CDI_2008_data_set)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(CDI_2008_data_set)
CDI_2008_data_set[,C+1] = 1
colnames(CDI_2008_data_set)[C+1] = "ones"

#Declare the design
d_unequal_CDI_2008 = svydesign(ids=~ids, data=CDI_2008_data_set, weights=~weight)

#Total number of children
svytotal(~ones, d_unequal_CDI_2008)

#Survey means
svymean(~cocoachild_f, d_unequal_CDI_2008)
svymean(~cl_cocoa, d_unequal_CDI_2008)
svymean(~haz_cocoa, d_unequal_CDI_2008)

#Survey totals
svytotal(~cocoachild_f, d_unequal_CDI_2008)
svytotal(~cl_cocoa, d_unequal_CDI_2008)
svytotal(~haz_cocoa, d_unequal_CDI_2008)

##########################################################################
##########################################################################
#CDI 2008 study - Household
##########################################################################
##########################################################################

CDI_2008_data_set_HH = subset(Full_data_set_HH, cdi==1 & Year=="2008")
nrow(CDI_2008_data_set_HH)

#Add a row of id indicators to declare the design
C = ncol(CDI_2008_data_set_HH)
CDI_2008_data_set_HH[,C+1] = 1:nrow(CDI_2008_data_set_HH)
colnames(CDI_2008_data_set_HH)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(CDI_2008_data_set_HH)
CDI_2008_data_set_HH[,C+1] = 1
colnames(CDI_2008_data_set_HH)[C+1] = "ones"

#Declare the design
d_unequal_CDI_2008_HH = svydesign(ids=~ids, data=CDI_2008_data_set_HH, weights=~weight)

#Total number of households
svytotal(~ones, d_unequal_CDI_2008_HH)

##########################################################################
##########################################################################
#Ghana and CDI 2008 studies combined - Child
##########################################################################
##########################################################################

Ghana_CDI_2008_data_set = rbind(Ghana_2008_data_set, CDI_2008_data_set)
nrow(Ghana_CDI_2008_data_set)

#Declare the sampling design
d_unequal_Ghana_CDI_2008 = svydesign(ids=~ids, data=Ghana_CDI_2008_data_set, weights=~weight)

#Total number of children
svytotal(~ones, d_unequal_Ghana_CDI_2008)

#Survey means 
svymean(~cocoachild_f, d_unequal_Ghana_CDI_2008)
svymean(~cl_cocoa, d_unequal_Ghana_CDI_2008)
svymean(~haz_cocoa, d_unequal_Ghana_CDI_2008)

#Survey totals
svytotal(~cocoachild_f, d_unequal_Ghana_CDI_2008)
svytotal(~cl_cocoa, d_unequal_Ghana_CDI_2008)
svytotal(~haz_cocoa, d_unequal_Ghana_CDI_2008)

##########################################################################
##########################################################################
#Ghana and CDI 2008 studies combined - Household
##########################################################################
##########################################################################

Ghana_CDI_2008_data_set_HH = rbind(Ghana_2008_data_set_HH, CDI_2008_data_set_HH)
nrow(Ghana_CDI_2008_data_set_HH)

#Declare the sampling design
d_unequal_Ghana_CDI_2008_HH = svydesign(ids=~ids, data=Ghana_CDI_2008_data_set_HH, weights=~weight)

#Total number of households
svytotal(~ones, d_unequal_Ghana_CDI_2008_HH)

##########################################################################
##########################################################################
#Ghana 2013 study - Child
##########################################################################
##########################################################################

#Subset out the corresponding entries
Ghana_2013_data_set = subset(Full_data_set_child, ghana==1 & Year=="2013" & !is.na(cocoa_weight))
nrow(Ghana_2013_data_set)

#Add a row of id indicators to declare the design
C = ncol(Ghana_2013_data_set)
Ghana_2013_data_set[,C+1] = 1:nrow(Ghana_2013_data_set)
colnames(Ghana_2013_data_set)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(Ghana_2013_data_set)
Ghana_2013_data_set[,C+1] = 1
colnames(Ghana_2013_data_set)[C+1] = "ones"

#Declare the sampling design
d_unequal_Ghana_2013 = svydesign(ids=~ids, data=Ghana_2013_data_set, weights=~cocoa_weight)

#Total number of children
svytotal(~ones, d_unequal_Ghana_2013)

#Survey means 
svymean(~cocoachild_f, d_unequal_Ghana_2013)
svymean(~cl_cocoa, d_unequal_Ghana_2013)
svymean(~haz_cocoa, d_unequal_Ghana_2013)

#Survey totals
svytotal(~cocoachild_f, d_unequal_Ghana_2013)
svytotal(~cl_cocoa, d_unequal_Ghana_2013)
svytotal(~haz_cocoa, d_unequal_Ghana_2013)

##########################################################################
##########################################################################
#Ghana 2013 study - Household
##########################################################################
##########################################################################

#Subset out the corresponding entries
Ghana_2013_data_set_HH = subset(Full_data_set_HH, ghana==1 & Year=="2013" & noncocoa==0)
nrow(Ghana_2013_data_set_HH)

#Add a row of id indicators to declare the design
C = ncol(Ghana_2013_data_set_HH)
Ghana_2013_data_set_HH[,C+1] = 1:nrow(Ghana_2013_data_set_HH)
colnames(Ghana_2013_data_set_HH)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(Ghana_2013_data_set_HH)
Ghana_2013_data_set_HH[,C+1] = 1
colnames(Ghana_2013_data_set_HH)[C+1] = "ones"

#Declare the sampling design
d_unequal_Ghana_2013_HH = svydesign(ids=~ids, data=Ghana_2013_data_set_HH, weights=~cocoa_weight)

##########################################################################
##########################################################################
#CDI 2013 study - Child
##########################################################################
##########################################################################

CDI_2013_data_set = subset(Full_data_set_child, cdi==1 & Year=="2013" & !is.na(cocoa_weight))
nrow(CDI_2013_data_set)

#Add a row of id indicators to declare the design
C = ncol(CDI_2013_data_set)
CDI_2013_data_set[,C+1] = 1:nrow(CDI_2013_data_set)
colnames(CDI_2013_data_set)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(CDI_2013_data_set)
CDI_2013_data_set[,C+1] = 1
colnames(CDI_2013_data_set)[C+1] = "ones"

#Declare the design
d_unequal_CDI_2013 = svydesign(ids=~ids, data=CDI_2013_data_set, weights=~cocoa_weight)

#Total number of children
svytotal(~ones, d_unequal_CDI_2013)

#Survey means 
svymean(~cocoachild_f, d_unequal_CDI_2013)
svymean(~cl_cocoa, d_unequal_CDI_2013)
svymean(~haz_cocoa, d_unequal_CDI_2013)

#Survey totals 
svytotal(~cocoachild_f, d_unequal_CDI_2013)
svytotal(~cl_cocoa, d_unequal_CDI_2013)
svytotal(~haz_cocoa, d_unequal_CDI_2013)

##########################################################################
##########################################################################
#CDI 2013 study - Household
##########################################################################
##########################################################################

CDI_2013_data_set_HH = subset(Full_data_set_HH, cdi==1 & Year=="2013" & noncocoa==0)
nrow(CDI_2013_data_set_HH)

#Add a row of id indicators to declare the design
C = ncol(CDI_2013_data_set_HH)
CDI_2013_data_set_HH[,C+1] = 1:nrow(CDI_2013_data_set_HH)
colnames(CDI_2013_data_set_HH)[C+1] = "ids"

#Add a row of 1s to estimate total population size 
C = ncol(CDI_2013_data_set_HH)
CDI_2013_data_set_HH[,C+1] = 1
colnames(CDI_2013_data_set_HH)[C+1] = "ones"

#Declare the design
d_unequal_CDI_2013_HH = svydesign(ids=~ids, data=CDI_2013_data_set_HH, weights=~cocoa_weight)

#Total number of children
svytotal(~ones, d_unequal_CDI_2013_HH)

##########################################################################
##########################################################################
#Ghana and CDI 2013 studies combined - Child
##########################################################################
##########################################################################

Ghana_CDI_2013_data_set = rbind(Ghana_2013_data_set, CDI_2013_data_set)
nrow(Ghana_CDI_2013_data_set)

#Declare the sampling design
d_unequal_Ghana_CDI_2013 = svydesign(ids=~ids, data=Ghana_CDI_2013_data_set, weights=~cocoa_weight)

#Total number of children
svytotal(~ones, d_unequal_Ghana_CDI_2013)

#Survey means 
svymean(~cocoachild_f, d_unequal_Ghana_CDI_2013)
svymean(~cl_cocoa, d_unequal_Ghana_CDI_2013)
svymean(~haz_cocoa, d_unequal_Ghana_CDI_2013)

#Survey totals
svytotal(~cocoachild_f, d_unequal_Ghana_CDI_2013)
svytotal(~cl_cocoa, d_unequal_Ghana_CDI_2013)
svytotal(~haz_cocoa, d_unequal_Ghana_CDI_2013)

##########################################################################
##########################################################################
#Ghana and CDI 2013 studies combined - Household
##########################################################################
##########################################################################

Ghana_CDI_2013_data_set_HH = rbind(Ghana_2013_data_set_HH, CDI_2013_data_set_HH)
nrow(Ghana_CDI_2013_data_set_HH)

#Declare the sampling design
d_unequal_Ghana_CDI_2013_HH = svydesign(ids=~ids, data=Ghana_CDI_2013_data_set_HH, weights=~cocoa_weight)

#Total number of children
svytotal(~ones, d_unequal_Ghana_CDI_2013_HH)

##########################################################################
##########################################################################
#Ghana 2018 study - Child
##########################################################################
##########################################################################

#The subset of the data
Ghana_2018_data_set = subset(Full_data_set_child_2018, ghana==1) 
nrow(Ghana_2018_data_set)

#Add a row of id indicators for declaring the design; note that these are a crossing of EA, HH, and individual in HH
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = 1:nrow(Ghana_2018_data_set)
colnames(Ghana_2018_data_set)[C+1] = "ids"
Ghana_2018_data_set[,"ids"] = 1:nrow(Ghana_2018_data_set)

#Add a row of 1s to estimate total population size 
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = 1
colnames(Ghana_2018_data_set)[C+1] = "ones"

#Read in the data set with the EA and HH counts
Ghana_2018_Stratum_EA_HH_counts = read_excel("Ghana-CDI_Survey_Weight_Count.xlsx", sheet = "Ghana_Sample_inf")

Ghana_2018_Stratum_EA_HH_counts = data.frame(Ghana_2018_Stratum_EA_HH_counts)
head(Ghana_2018_Stratum_EA_HH_counts)

#Read in the roster data set
Roster_data_set = read.dta13("Roster_2008_2018.dta")
Roster_data_set = subset(Roster_data_set, ghana==1 & Year==2018)
Roster_data_set = as.data.frame(Roster_data_set)
head(Roster_data_set)

length(which(!(unique(Ghana_2018_data_set[,"hhid"]) %in% Roster_data_set[,"hh_id"])))

x = which(!(Roster_data_set[,"hh_id"] %in% Ghana_2018_data_set[,"hhid"]))
Roster_data_set[x,"hh_id"]
length(x)
x = which(Roster_data_set[x,"children"] != 0)
Roster_data_set[x,"hh_id"]
length(x)

for(k in unique(Ghana_2018_data_set[,"hhid"]))
	{
	x = which(Ghana_2018_data_set[,"hhid"] == k)
	if(length(unique(Ghana_2018_data_set[x,"hhid"])) > 1)
		cat("hhid", k, "needs review", '\n')
	}



z = numeric()
for(k in unique(Ghana_2018_data_set[,"hhid"]))
	{
	x = which(Ghana_2018_data_set[,"hhid"] == k)	
	y = which(Roster_data_set[,"hh_id"] == k)	
	if(unique(Ghana_2018_data_set[x,"num_children_roster"]) != Roster_data_set[y,"children"])
		{
		cat("hhid", k, "child count does not coincide", '\n')
		z = c(z, k)
		
		Ghana_2018_data_set[x,"num_children_roster"] = Roster_data_set[y,"children"]		
		}
	}

#Let's ensure the child count is greater than the cocoa labor count
which(Ghana_2018_data_set[,"cl_cocoa"] > Ghana_2018_data_set[,"num_children_roster"])

#Here is how sampling was done:
#1) Select EAs at random within each stratum
#2) Select eligible HHs at random within each EA (note that some HHs did not have children present at the time of the survey)
#3) Select as many children as possible to interview in each selected household

##########################################################################
#The expert group weighting approach for the Ghana 2018 study
##########################################################################

#We will need to evaluate the weights manually for each observation
#First level, EA-level weight (selecting EAs within each stratum)
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = NA
colnames(Ghana_2018_data_set)[C+1] = "EA_weight"

x = which(Ghana_2018_data_set[,"Strata"] == 1)
Ghana_2018_data_set[x,"EA_weight"] = 1875/40
x = which(Ghana_2018_data_set[,"Strata"] == 2)
Ghana_2018_data_set[x,"EA_weight"] = 3659/23
x = which(Ghana_2018_data_set[,"Strata"] == 3)
Ghana_2018_data_set[x,"EA_weight"] = 1063/10

#Second level, HH-level (selecting HHs within each EA)	
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = NA
colnames(Ghana_2018_data_set)[C+1] = "HH_weight"

for(k in unique(Ghana_2018_data_set[,"EA_2018"]))
	{
	x = which(Ghana_2018_data_set[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(Ghana_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		Ghana_2018_data_set[x,"HH_weight"] = 
			Ghana_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]/
				Ghana_2018_Stratum_EA_HH_counts[y,"Eligible.AG.HH.interviewed"]
		}
	}

#Third level, individual-level (averaging over number of eligible individuals in selected HHs in EA against number interviewed in EA)
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = NA
colnames(Ghana_2018_data_set)[C+1] = "Ind_weight"

for(k in unique(Ghana_2018_data_set[,"EA_2018"]))
	{
	x = which(Ghana_2018_data_set[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(Ghana_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		Ghana_2018_data_set[x,"Ind_weight"] = 
			Ghana_2018_Stratum_EA_HH_counts[y,"Total.children.in.the.roster.of.sampled.Eligible.AG.HH"]/
				Ghana_2018_Stratum_EA_HH_counts[y,"Total.children.interviewed"]
		}
	}

#The final "EG" weights
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = NA
colnames(Ghana_2018_data_set)[C+1] = "EG_weight"
Ghana_2018_data_set[,"EG_weight"] = 
	Ghana_2018_data_set[,"EA_weight"] * Ghana_2018_data_set[,"HH_weight"] * Ghana_2018_data_set[,"Ind_weight"]

#We need finite population correction factors corresponding to each level to apply 

#The first finite population correct factor corresponds to selecting EAs within each stratum
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = NA
colnames(Ghana_2018_data_set)[C+1] = "FPC_EA_within_stratum"

x = which(Ghana_2018_data_set[,"Strata"] == 1)
Ghana_2018_data_set[x,"FPC_EA_within_stratum"] = 1875
x = which(Ghana_2018_data_set[,"Strata"] == 2)
Ghana_2018_data_set[x,"FPC_EA_within_stratum"] = 3659
x = which(Ghana_2018_data_set[,"Strata"] == 3)
Ghana_2018_data_set[x,"FPC_EA_within_stratum"] = 1063

#The second finite population correction factor corresponds to selecting eligible HHs within each EA
#However, note that some do not respond, so we will use the actual number of HH surveyed here
#To do this correctly one would have to come up with a set of dummy entries, the number of HH which do not respond within each EA and
#then code these with NAs; there are 84 of these 
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = NA
colnames(Ghana_2018_data_set)[C+1] = "FPC_HH_within_EA"

for(k in unique(Full_data_set_child[,"EA_2018"]))
	{
	x = which(Ghana_2018_data_set[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(Ghana_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		Ghana_2018_data_set[x,"FPC_HH_within_EA"] = 
			Ghana_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]
		}
	}

#Note that we do not actually make use of the third finite population correction factor since we model selection of individuals within
#each EA as a simple random sample among the pool of individuals available within the households selected in the EA
d_EG_Ghana_2018 <- svydesign(ids=~EA_2018+hhid, strata=~Strata, 
								fpc=~FPC_EA_within_stratum+FPC_HH_within_EA, weights=~EG_weight,
									data = Ghana_2018_data_set, nest = FALSE, check.strata = FALSE)

#A histogram of the EG weights
hist(weights(d_EG_Ghana_2018), breaks=100, col="lightblue", 
	main="Histogram of Expert Group Suggested Sampling Weights", xlab="Sampling Weight")

#Total number of children
svytotal(~ones, d_EG_Ghana_2018)

#The estimates of the means 
svymean(~cocoachild_f, d_EG_Ghana_2018)
svymean(~cl_cocoa, d_EG_Ghana_2018)
svymean(~haz_cocoa, d_EG_Ghana_2018)

#The estimates of the totals 
svytotal(~cocoachild_f, d_EG_Ghana_2018)
svytotal(~cl_cocoa, d_EG_Ghana_2018)
svytotal(~haz_cocoa, d_EG_Ghana_2018)

##The Rao-Wu Bootstrap Routine to properly calculate standard errors
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_EG_Ghana_2018 <- as.svrepdesign(design=d_EG_Ghana_2018, type="subbootstrap", replicates=2000)
d_boot_EG_Ghana_2018

d_boot_EG_Ghana_2018$type
d_boot_EG_Ghana_2018$rho
d_boot_EG_Ghana_2018$scale
d_boot_EG_Ghana_2018$rscales

#Total number of children
svytotal(~ones, d_boot_EG_Ghana_2018)

#The bootstrap estimates of the means
svymean(~cocoachild_f, d_boot_EG_Ghana_2018)
svymean(~cl_cocoa, d_boot_EG_Ghana_2018)
svymean(~haz_cocoa, d_boot_EG_Ghana_2018)

#The bootstrap estimates of the totals
svytotal(~cocoachild_f, d_boot_EG_Ghana_2018)
svytotal(~cl_cocoa, d_boot_EG_Ghana_2018)
svytotal(~haz_cocoa, d_boot_EG_Ghana_2018)




##########################################################################
#The sampling design exact weighting approach for the Ghana 2018 study
##########################################################################

#We need finite population correction factors corresponding to each level
#We already have the first two
#The third finite population correction factor corresponds to selecting individuals within each HH
C = ncol(Ghana_2018_data_set)
Ghana_2018_data_set[,C+1] = NA
colnames(Ghana_2018_data_set)[C+1] = "FPC_Ind_within_HH"
Ghana_2018_data_set[,"FPC_Ind_within_HH"] = Ghana_2018_data_set[,"num_children_roster"]

d_multi_Ghana_2018 <- svydesign(ids=~EA_2018+hhid+ids, strata=~Strata, 
								fpc=~FPC_EA_within_stratum+FPC_HH_within_EA+FPC_Ind_within_HH, 
									data = Ghana_2018_data_set, nest = FALSE, check.strata = FALSE)

#A histogram of the design weights
hist(weights(d_multi_Ghana_2018), breaks=80, col="lightblue", 
	main="Histogram of Sampling Design Weights", xlab="Sampling Weight")

#There may be a (small) issue with lonely PSUs
options(survey.lonely.psu="adjust")

#Total number of children
svytotal(~ones, d_multi_Ghana_2018)

#The estimates of the means 
svymean(~cocoachild_f, d_multi_Ghana_2018)
svymean(~cl_cocoa, d_multi_Ghana_2018)
svymean(~haz_cocoa, d_multi_Ghana_2018)

#The estimates of the totals
svytotal(~cocoachild_f, d_multi_Ghana_2018)
svytotal(~cl_cocoa, d_multi_Ghana_2018)
svytotal(~haz_cocoa, d_multi_Ghana_2018)

##The Rao-Wu Bootstrap Routine
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_multi_Ghana_2018 <- as.svrepdesign(design=d_multi_Ghana_2018, type="subbootstrap", replicates=2000)
d_boot_multi_Ghana_2018

d_boot_multi_Ghana_2018$type
d_boot_multi_Ghana_2018$rho
d_boot_multi_Ghana_2018$scale
d_boot_multi_Ghana_2018$rscales

#Total number of children
svytotal(~ones, d_boot_multi_Ghana_2018)

#The bootstrap estimates of the mean
svymean(~cocoachild_f, d_boot_multi_Ghana_2018)
svymean(~cl_cocoa, d_boot_multi_Ghana_2018)
svymean(~haz_cocoa, d_boot_multi_Ghana_2018)

#The bootstrap estimates of the total
svytotal(~cocoachild_f, d_boot_multi_Ghana_2018)
svytotal(~cl_cocoa, d_boot_multi_Ghana_2018)
svytotal(~haz_cocoa, d_boot_multi_Ghana_2018)

##########################################################################
#Hypothesis Testing for the Ghana 2018 study
##########################################################################

#H_0: The means are the same
#H_A: The mean in 2018 is larger than that in 2008
#This entails a one-tailed test. For a two-tailed test, just double the p-value.

#For child working in cocoa, here is all the output we need
svymean(~cocoachild_f, d_unequal_Ghana_2008)
svymean(~cocoachild_f, d_boot_EG_Ghana_2018)
svymean(~cocoachild_f, d_boot_multi_Ghana_2018)

Z = (0.59945 - 0.46155) / (sqrt(2)*0.0290)
1-pnorm(Z)

Z = (0.60155 - 0.46155) / (sqrt(2)*0.0298)
1-pnorm(Z)

#For child labor in cocoa, here is all the output we need
svymean(~cl_cocoa, d_unequal_Ghana_2008)
svymean(~cl_cocoa, d_boot_EG_Ghana_2018)
svymean(~cl_cocoa, d_boot_multi_Ghana_2018)

Z = (0.54932 - 0.43861) / (sqrt(2)*0.0257)
1-pnorm(Z)

Z = (0.55071 - 0.43861) / (sqrt(2)*0.0266)
1-pnorm(Z)

#For hazardous child labor in cocoa, here is all the output we need
svymean(~haz_cocoa, d_unequal_Ghana_2008)
svymean(~haz_cocoa, d_boot_EG_Ghana_2018)
svymean(~haz_cocoa, d_boot_multi_Ghana_2018)

Z = (0.51177 - 0.43085) / (sqrt(2)*0.0233)
1-pnorm(Z)

Z = (0.51175 - 0.43085) / (sqrt(2)*0.0237)
1-pnorm(Z)

##########################################################################
##########################################################################
#Ghana 2018 study - Household
##########################################################################
##########################################################################

#The subset of the data
Ghana_2018_data_set_HH = subset(Full_data_set_HH, ghana==1 & Year=="2018") 
nrow(Ghana_2018_data_set_HH)

#Add a row of id indicators for declaring the design; note that these are a crossing of EA, HH, and individual in HH
C = ncol(Ghana_2018_data_set_HH)
Ghana_2018_data_set_HH[,C+1] = 1:nrow(Ghana_2018_data_set_HH)
colnames(Ghana_2018_data_set_HH)[C+1] = "ids"
Ghana_2018_data_set_HH[,"ids"] = 1:nrow(Ghana_2018_data_set_HH)

#Add a row of 1s to estimate total population size 
C = ncol(Ghana_2018_data_set_HH)
Ghana_2018_data_set_HH[,C+1] = 1
colnames(Ghana_2018_data_set_HH)[C+1] = "ones"

#Read in the data set with the EA and HH counts
Ghana_2018_Stratum_EA_HH_counts = read_excel("Ghana-CDI_Survey_Weight_Count.xlsx", sheet = "Ghana_Sample_inf")

Ghana_2018_Stratum_EA_HH_counts = data.frame(Ghana_2018_Stratum_EA_HH_counts)
head(Ghana_2018_Stratum_EA_HH_counts)

#Here is how sampling was done:
#1) Select EAs at random within each stratum
#2) Select eligible HHs at random within each EA (note that some HHs did not have children present at the time of the survey)
#3) Select as many children as possible to interview in each selected household

##########################################################################
#The expert group weighting approach for the Ghana 2018 study
##########################################################################

#We will need to evaluate the weights manually for each observation
#First level, EA-level weight (selecting EAs within each stratum)
C = ncol(Ghana_2018_data_set_HH)
Ghana_2018_data_set_HH[,C+1] = NA
colnames(Ghana_2018_data_set_HH)[C+1] = "EA_weight"

x = which(Ghana_2018_data_set_HH[,"strata"] == 1)
Ghana_2018_data_set_HH[x,"EA_weight"] = 1875/40
x = which(Ghana_2018_data_set_HH[,"strata"] == 2)
Ghana_2018_data_set_HH[x,"EA_weight"] = 3659/23
x = which(Ghana_2018_data_set_HH[,"strata"] == 3)
Ghana_2018_data_set_HH[x,"EA_weight"] = 1063/10

#Second level, HH-level (selecting HHs within each EA)	
C = ncol(Ghana_2018_data_set_HH)
Ghana_2018_data_set_HH[,C+1] = NA
colnames(Ghana_2018_data_set_HH)[C+1] = "HH_weight"

for(k in unique(Ghana_2018_data_set_HH[,"EA_2018"]))
	{
	x = which(Ghana_2018_data_set_HH[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(Ghana_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		Ghana_2018_data_set_HH[x,"HH_weight"] = 
			Ghana_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]/
				Ghana_2018_Stratum_EA_HH_counts[y,"Eligible.AG.HH.interviewed"]
		}
	}

#The final "EG" weights
C = ncol(Ghana_2018_data_set_HH)
Ghana_2018_data_set_HH[,C+1] = NA
colnames(Ghana_2018_data_set_HH)[C+1] = "EG_weight"
Ghana_2018_data_set_HH[,"EG_weight"] = 
	Ghana_2018_data_set_HH[,"EA_weight"] * Ghana_2018_data_set_HH[,"HH_weight"] 

#We need finite population correction factors corresponding to each level to apply 
#The first finite population correct factor corresponds to selecting EAs within each stratum
C = ncol(Ghana_2018_data_set_HH)
Ghana_2018_data_set_HH[,C+1] = NA
colnames(Ghana_2018_data_set_HH)[C+1] = "FPC_EA_within_stratum"

x = which(Ghana_2018_data_set_HH[,"strata"] == 1)
Ghana_2018_data_set_HH[x,"FPC_EA_within_stratum"] = 1875
x = which(Ghana_2018_data_set_HH[,"strata"] == 2)
Ghana_2018_data_set_HH[x,"FPC_EA_within_stratum"] = 3659
x = which(Ghana_2018_data_set_HH[,"strata"] == 3)
Ghana_2018_data_set_HH[x,"FPC_EA_within_stratum"] = 1063

#The second finite population correction factor corresponds to selecting eligible HHs within each EA
#However, note that some do not respond, so we will use the actual number of HH surveyed here
#To do this correctly one would have to come up with a set of dummy entries, the number of HH which do not respond within each EA and
#then code these with NAs; there are 84 of these 
C = ncol(Ghana_2018_data_set_HH)
Ghana_2018_data_set_HH[,C+1] = NA
colnames(Ghana_2018_data_set_HH)[C+1] = "FPC_HH_within_EA"

for(k in unique(Full_data_set_HH[,"EA_2018"]))
	{
	x = which(Ghana_2018_data_set_HH[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(Ghana_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		Ghana_2018_data_set_HH[x,"FPC_HH_within_EA"] = 
			Ghana_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]
		}
	}
		
#The design
d_EG_Ghana_2018_HH <- svydesign(ids=~EA_2018, strata=~strata, 
								fpc=~FPC_EA_within_stratum, weights=~EG_weight,
									data = Ghana_2018_data_set_HH, nest = FALSE, check.strata = FALSE)

#A histogram of the EG weights
hist(weights(d_EG_Ghana_2018_HH), breaks=100, col="lightblue", 
	main="Histogram of Expert Group Suggested Sampling Weights", xlab="Sampling Weight")

#Total number of children
svytotal(~ones, d_EG_Ghana_2018_HH)

##The Rao-Wu Bootstrap Routine to properly calculate standard errors
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_EG_Ghana_2018_HH <- as.svrepdesign(design=d_EG_Ghana_2018_HH, type="subbootstrap", replicates=2000)
d_boot_EG_Ghana_2018_HH

d_boot_EG_Ghana_2018_HH$type
d_boot_EG_Ghana_2018_HH$rho
d_boot_EG_Ghana_2018_HH$scale
d_boot_EG_Ghana_2018_HH$rscales

#Total number of children
svytotal(~ones, d_boot_EG_Ghana_2018_HH)

##########################################################################
##########################################################################
#CDI 2018 study - Child
##########################################################################
##########################################################################

#The subset of the data
CDI_2018_data_set = subset(Full_data_set_child_2018, cdi==1) 
nrow(CDI_2018_data_set)

#Add a row of id indicators for declaring the design; note that these are a crossing of EA, HH, and individual in HH
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = 1:nrow(CDI_2018_data_set)
colnames(CDI_2018_data_set)[C+1] = "ids"
CDI_2018_data_set[,"ids"] = 1:nrow(CDI_2018_data_set)

#Add a row of 1s to estimate total population size 
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = 1
colnames(CDI_2018_data_set)[C+1] = "ones"

#Read in the data set with the EA and HH counts
CDI_2018_Stratum_EA_HH_counts = read_excel("Ghana-CDI_Survey_Weight_Count.xlsx", sheet = "CDI_Sample_inf")

CDI_2018_Stratum_EA_HH_counts = data.frame(CDI_2018_Stratum_EA_HH_counts)
head(CDI_2018_Stratum_EA_HH_counts)

#Read in the roster data set
Roster_data_set = read.dta13("Roster_2008_2018.dta")
Roster_data_set = subset(Roster_data_set, cdi==1 & Year==2018)
Roster_data_set = as.data.frame(Roster_data_set)
head(Roster_data_set)

length(which(!(unique(CDI_2018_data_set[,"hhid"]) %in% Roster_data_set[,"hh_id"])))
 
x = which(!(Roster_data_set[,"hh_id"] %in% CDI_2018_data_set[,"hhid"]))
Roster_data_set[x,"hh_id"]
length(x)
x = which(Roster_data_set[x,"children"] != 0)
Roster_data_set[x,"hh_id"]
length(x)

#Does every household have a unique "num_child_roster" entry?
for(k in unique(CDI_2018_data_set[,"hhid"]))
	{
	x = which(CDI_2018_data_set[,"hhid"] == k)
	if(length(unique(CDI_2018_data_set[x,"hhid"])) > 1)
		cat("hhid", k, "needs review", '\n')
	}

z = numeric()
for(k in unique(CDI_2018_data_set[,"hhid"]))
	{
	x = which(CDI_2018_data_set[,"hhid"] == k)	
	y = which(Roster_data_set[,"hh_id"] == k)	
	if(unique(CDI_2018_data_set[x,"num_children_roster"]) != Roster_data_set[y,"children"])
		{
		cat("hhid", k, "child count does not coincide", '\n')
		z = c(z, k)
		
		CDI_2018_data_set[x,"num_children_roster"] = Roster_data_set[y,"children"]		
		}
	}

#Let's ensure the child count is greater than the cocoa labor count
which(CDI_2018_data_set[,"cl_cocoa"] > CDI_2018_data_set[,"num_children_roster"])

#Here is how sampling was done:
#1) Select EAs at random within each stratum
#2) Select eligible HHs at random within each EA (note that some HHs did not have children present at the time of the survey)
#3) Select as many children as possible to interview in each selected household

##########################################################################
#The expert group weighting approach for the CDI 2018 study
##########################################################################

#We will need to evaluate the weights manually for each observation

#First level, EA-level weight (selecting EAs within each stratum)
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = NA
colnames(CDI_2018_data_set)[C+1] = "EA_weight"

x = which(CDI_2018_data_set[,"Strata"] == 1)
CDI_2018_data_set[x,"EA_weight"] = 3409/40
x = which(CDI_2018_data_set[,"Strata"] == 2)
CDI_2018_data_set[x,"EA_weight"] = 3348/25
x = which(CDI_2018_data_set[,"Strata"] == 3)
CDI_2018_data_set[x,"EA_weight"] = 2803/10


#Second level, HH-level (selecting HHs within each EA)	
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = NA
colnames(CDI_2018_data_set)[C+1] = "HH_weight"

for(k in unique(CDI_2018_data_set[,"EA_2018"]))
	{
	x = which(CDI_2018_data_set[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(CDI_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		CDI_2018_data_set[x,"HH_weight"] = 
			CDI_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]/
				CDI_2018_Stratum_EA_HH_counts[y,"Eligible.AG.HH.interviewed"]
		}
	}

#Third level, individual-level (averaging over number of eligible individuals in selected HHs in EA against number interviewed in EA)
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = NA
colnames(CDI_2018_data_set)[C+1] = "Ind_weight"

for(k in unique(CDI_2018_data_set[,"EA_2018"]))
	{
	x = which(CDI_2018_data_set[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(CDI_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		CDI_2018_data_set[x,"Ind_weight"] = 
			CDI_2018_Stratum_EA_HH_counts[y,"Total.children.in.the.roster.of.sampled.Eligible.AG.HH"]/
				CDI_2018_Stratum_EA_HH_counts[y,"Total.children.interviewed"]
		}
	}

#The final "EG" weights
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = NA
colnames(CDI_2018_data_set)[C+1] = "EG_weight"
CDI_2018_data_set[,"EG_weight"] = 
	CDI_2018_data_set[,"EA_weight"] * CDI_2018_data_set[,"HH_weight"] * CDI_2018_data_set[,"Ind_weight"]

#We need finite population correction factors corresponding to each level to apply 
#The first finite population correct factor corresponds to selecting EAs within each stratum
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = NA
colnames(CDI_2018_data_set)[C+1] = "FPC_EA_within_stratum"

x = which(CDI_2018_data_set[,"Strata"] == 1)
CDI_2018_data_set[x,"FPC_EA_within_stratum"] = 3409
x = which(CDI_2018_data_set[,"Strata"] == 2)
CDI_2018_data_set[x,"FPC_EA_within_stratum"] = 3348
x = which(CDI_2018_data_set[,"Strata"] == 3)
CDI_2018_data_set[x,"FPC_EA_within_stratum"] = 2803

#The second finite population correction factor corresponds to selecting eligible HHs within each EA
#However, note that some do not respond, so we will use the actual number of HH surveyed here
#To do this correctly one would have to come up with a set of dummy entries, the number of HH which do not respond within each EA and
#then code these with NAs; there are 84 of these 
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = NA
colnames(CDI_2018_data_set)[C+1] = "FPC_HH_within_EA"

for(k in unique(Full_data_set_child[,"EA_2018"]))
	{
	x = which(CDI_2018_data_set[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(CDI_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		CDI_2018_data_set[x,"FPC_HH_within_EA"] = 
			CDI_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]
		}
	}
		
#Note that we do not actually make use of the third finite population correction factor since we model selection of individuals within
#each EA as a simple random sample among the pool of individuals available within the households selected in the EA

d_EG_CDI_2018 <- svydesign(ids=~EA_2018+hhid, strata=~Strata, 
								fpc=~FPC_EA_within_stratum+FPC_HH_within_EA, weights=~EG_weight,
									data = CDI_2018_data_set, nest = FALSE, check.strata = FALSE)

#A histogram of the EG weights
hist(weights(d_EG_CDI_2018), breaks=80, col="lightblue", 
	main="Histogram of Expert Group Suggested Sampling Weights", xlab="Sampling Weight")

#Total number of children
svytotal(~ones, d_EG_CDI_2018)

#The estimates of the means 
svymean(~cocoachild_f, d_EG_CDI_2018)
svymean(~cl_cocoa, d_EG_CDI_2018)
svymean(~haz_cocoa, d_EG_CDI_2018)

#The estimates of the standard errors 
svytotal(~cocoachild_f, d_EG_CDI_2018)
svytotal(~cl_cocoa, d_EG_CDI_2018)
svytotal(~haz_cocoa, d_EG_CDI_2018)

##The Rao-Wu Bootstrap Routine to properly calculate standard errors
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_EG_CDI_2018 <- as.svrepdesign(design=d_EG_CDI_2018, type="subbootstrap", replicates=2000)
d_boot_EG_CDI_2018

d_boot_EG_CDI_2018$type
d_boot_EG_CDI_2018$rho
d_boot_EG_CDI_2018$scale
d_boot_EG_CDI_2018$rscales

#The bootstrap estimates
svymean(~cocoachild_f, d_boot_EG_CDI_2018)
svymean(~cl_cocoa, d_boot_EG_CDI_2018)
svymean(~haz_cocoa, d_boot_EG_CDI_2018)

#The bootstrap estimates
svytotal(~cocoachild_f, d_boot_EG_CDI_2018)
svytotal(~cl_cocoa, d_boot_EG_CDI_2018)
svytotal(~haz_cocoa, d_boot_EG_CDI_2018)

##########################################################################
#The sampling design exact weighting approach for the CDI 2018 study
##########################################################################

#We need finite population correction factors corresponding to each level
#We already have the first two
#The third finite population correction factor corresponds to selecting individuals within each HH
C = ncol(CDI_2018_data_set)
CDI_2018_data_set[,C+1] = NA
colnames(CDI_2018_data_set)[C+1] = "FPC_Ind_within_HH"
CDI_2018_data_set[,"FPC_Ind_within_HH"] = CDI_2018_data_set[,"num_children_roster"]

d_multi_CDI_2018 <- svydesign(ids=~EA_2018+hhid+ids, strata=~Strata, 
								fpc=~FPC_EA_within_stratum+FPC_HH_within_EA+FPC_Ind_within_HH, 
									data = CDI_2018_data_set, nest = FALSE, check.strata = FALSE)

#A histogram of the design weights
hist(weights(d_multi_CDI_2018), breaks=80, col="lightblue", 
	main="Histogram of Sampling Design Weights", xlab="Sampling Weight")

#There may be a (small) issue with lonely PSUs
options(survey.lonely.psu="adjust")

#Total number of children
svytotal(~ones, d_multi_CDI_2018)

#The estimates of the means 
svymean(~cocoachild_f, d_multi_CDI_2018)
svymean(~cl_cocoa, d_multi_CDI_2018)
svymean(~haz_cocoa, d_multi_CDI_2018)

#The estimates of the totals
svytotal(~cocoachild_f, d_multi_CDI_2018)
svytotal(~cl_cocoa, d_multi_CDI_2018)
svytotal(~haz_cocoa, d_multi_CDI_2018)

##The Rao-Wu Bootstrap Routine
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_multi_CDI_2018 <- as.svrepdesign(design=d_multi_CDI_2018, type="subbootstrap", replicates=2000)
d_boot_multi_CDI_2018

d_boot_multi_CDI_2018$type
d_boot_multi_CDI_2018$rho
d_boot_multi_CDI_2018$scale
d_boot_multi_CDI_2018$rscales

#Total number of children
svytotal(~ones, d_boot_multi_CDI_2018)

#The bootstrap estimates
svymean(~cocoachild_f, d_boot_multi_CDI_2018)
svymean(~cl_cocoa, d_boot_multi_CDI_2018)
svymean(~haz_cocoa, d_boot_multi_CDI_2018)

#The bootstrap estimates
svytotal(~cocoachild_f, d_boot_multi_CDI_2018)
svytotal(~cl_cocoa, d_boot_multi_CDI_2018)
svytotal(~haz_cocoa, d_boot_multi_CDI_2018)

##########################################################################
#Hypothesis Testing for the CDI 2018 study
##########################################################################

#H_0: The means are the same
#H_A: The mean in 2018 is larger than that in 2008
#This entails a one-tailed test. For a two-tailed test, just double the p-value.
#For child working in cocoa, here is all the output we need
svymean(~cocoachild_f, d_unequal_CDI_2008)
svymean(~cocoachild_f, d_boot_EG_CDI_2018)
svymean(~cocoachild_f, d_boot_multi_CDI_2018)

Z = (0.39949 - 0.23099) / (sqrt(2)*0.0239)
1-pnorm(Z)

Z = (0.4055 - 0.23099) / (sqrt(2)*0.0227)
1-pnorm(Z)

#For child labor in cocoa, here is all the output we need
svymean(~cl_cocoa, d_unequal_CDI_2008)
svymean(~cl_cocoa, d_boot_EG_CDI_2018)
svymean(~cl_cocoa, d_boot_multi_CDI_2018)

Z = (0.37966 - 0.23013) / (sqrt(2)*0.0227)
1-pnorm(Z)

Z = (0.38701 - 0.23013) / (sqrt(2)*0.0217)
1-pnorm(Z)

#For hazardous child labor in cocoa, here is all the output we need
svymean(~haz_cocoa, d_unequal_CDI_2008)
svymean(~haz_cocoa, d_boot_EG_CDI_2018)
svymean(~haz_cocoa, d_boot_multi_CDI_2018)

Z = (0.36746 - 0.22687) / (sqrt(2)*0.0219)
1-pnorm(Z)

Z = (0.37512 - 0.22687) / (sqrt(2)*0.0209)
1-pnorm(Z)

##########################################################################
##########################################################################
#CDI 2018 study - Household
##########################################################################
##########################################################################

#The subset of the data
CDI_2018_data_set_HH = subset(Full_data_set_HH, cdi==1 & Year=="2018") 
nrow(CDI_2018_data_set_HH)

#Add a row of id indicators for declaring the design; note that these are a crossing of EA, HH, and individual in HH
C = ncol(CDI_2018_data_set_HH)
CDI_2018_data_set_HH[,C+1] = 1:nrow(CDI_2018_data_set_HH)
colnames(CDI_2018_data_set_HH)[C+1] = "ids"
CDI_2018_data_set_HH[,"ids"] = 1:nrow(CDI_2018_data_set_HH)

#Add a row of 1s to estimate total population size 
C = ncol(CDI_2018_data_set_HH)
CDI_2018_data_set_HH[,C+1] = 1
colnames(CDI_2018_data_set_HH)[C+1] = "ones"

#Read in the data set with the EA and HH counts
CDI_2018_Stratum_EA_HH_counts = read_excel("Ghana-CDI_Survey_Weight_Count.xlsx", sheet = "CDI_Sample_inf")

CDI_2018_Stratum_EA_HH_counts = data.frame(CDI_2018_Stratum_EA_HH_counts)
head(CDI_2018_Stratum_EA_HH_counts)

#Here is how sampling was done:
#1) Select EAs at random within each stratum
#2) Select eligible HHs at random within each EA (note that some HHs did not have children present at the time of the survey)
#3) Select as many children as possible to interview in each selected household

##########################################################################
#The expert group weighting approach for the CDI 2018_HH study
##########################################################################

#We will need to evaluate the weights manually for each observation
#First level, EA-level weight (selecting EAs within each stratum)
C = ncol(CDI_2018_data_set_HH)
CDI_2018_data_set_HH[,C+1] = NA
colnames(CDI_2018_data_set_HH)[C+1] = "EA_weight"

x = which(CDI_2018_data_set_HH[,"strata"] == 1)
CDI_2018_data_set_HH[x,"EA_weight"] = 3409/40
x = which(CDI_2018_data_set_HH[,"strata"] == 2)
CDI_2018_data_set_HH[x,"EA_weight"] = 3348/25
x = which(CDI_2018_data_set_HH[,"strata"] == 3)
CDI_2018_data_set_HH[x,"EA_weight"] = 2803/10

#Second level, HH-level (selecting HHs within each EA)	
C = ncol(CDI_2018_data_set_HH)
CDI_2018_data_set_HH[,C+1] = NA
colnames(CDI_2018_data_set_HH)[C+1] = "HH_weight"

for(k in unique(CDI_2018_data_set_HH[,"EA_2018"]))
	{
	x = which(CDI_2018_data_set_HH[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(CDI_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		CDI_2018_data_set_HH[x,"HH_weight"] = 
			CDI_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]/
				CDI_2018_Stratum_EA_HH_counts[y,"Eligible.AG.HH.interviewed"]
		}
	}

#The final "EG" weights
C = ncol(CDI_2018_data_set_HH)
CDI_2018_data_set_HH[,C+1] = NA
colnames(CDI_2018_data_set_HH)[C+1] = "EG_weight"
CDI_2018_data_set_HH[,"EG_weight"] = 
	CDI_2018_data_set_HH[,"EA_weight"] * CDI_2018_data_set_HH[,"HH_weight"]

#We need finite population correction factors corresponding to each level to apply 
#The first finite population correct factor corresponds to selecting EAs within each stratum
C = ncol(CDI_2018_data_set_HH)
CDI_2018_data_set_HH[,C+1] = NA
colnames(CDI_2018_data_set_HH)[C+1] = "FPC_EA_within_stratum"

x = which(CDI_2018_data_set_HH[,"strata"] == 1)
CDI_2018_data_set_HH[x,"FPC_EA_within_stratum"] = 3409
x = which(CDI_2018_data_set_HH[,"strata"] == 2)
CDI_2018_data_set_HH[x,"FPC_EA_within_stratum"] = 3348
x = which(CDI_2018_data_set_HH[,"strata"] == 3)
CDI_2018_data_set_HH[x,"FPC_EA_within_stratum"] = 2803

#The second finite population correction factor corresponds to selecting eligible HHs within each EA
#However, note that some do not respond, so we will use the actual number of HH surveyed here
#To do this correctly one would have to come up with a set of dummy entries, the number of HH which do not respond within each EA and
#then code these with NAs; there are 84 of these (
C = ncol(CDI_2018_data_set_HH)
CDI_2018_data_set_HH[,C+1] = NA
colnames(CDI_2018_data_set_HH)[C+1] = "FPC_HH_within_EA"

for(k in unique(Full_data_set_HH[,"EA_2018"]))
	{
	x = which(CDI_2018_data_set_HH[,"EA_2018"] == k)
	if(length(x) > 0)
		{
		y = which(CDI_2018_Stratum_EA_HH_counts[,"FULL_EA_CODE"] == k)
		CDI_2018_data_set_HH[x,"FPC_HH_within_EA"] = 
			CDI_2018_Stratum_EA_HH_counts[y,"Total.Eligible.AG.HH.listed"]
		}
	}
		
#The design
d_EG_CDI_2018_HH <- svydesign(ids=~EA_2018, strata=~strata, 
								fpc=~FPC_EA_within_stratum, weights=~EG_weight,
									data = CDI_2018_data_set_HH, nest = FALSE, check.strata = FALSE)

#A histogram of the EG weights
hist(weights(d_EG_CDI_2018_HH), breaks=80, col="lightblue", 
	main="Histogram of Expert Group Suggested Sampling Weights", xlab="Sampling Weight")

#Total number of households
svytotal(~ones, d_EG_CDI_2018_HH)

##The Rao-Wu Bootstrap Routine to properly calculate standard errors
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_EG_CDI_2018_HH <- as.svrepdesign(design=d_EG_CDI_2018_HH, type="subbootstrap", replicates=2000)
d_boot_EG_CDI_2018_HH

d_boot_EG_CDI_2018_HH$type
d_boot_EG_CDI_2018_HH$rho
d_boot_EG_CDI_2018_HH$scale
d_boot_EG_CDI_2018_HH$rscales

#Total number of households
svytotal(~ones, d_boot_EG_CDI_2018_HH)

##########################################################################
##########################################################################
#Combined Ghana 2018 and CDI 2018 studies - Child
##########################################################################
##########################################################################

x = which(Ghana_2018_data_set[,"Strata"] == 1)
Ghana_2018_data_set[x,"Strata"] = "Ghana_1"
x = which(Ghana_2018_data_set[,"Strata"] == 2)
Ghana_2018_data_set[x,"Strata"] = "Ghana_2"
x = which(Ghana_2018_data_set[,"Strata"] == 3)
Ghana_2018_data_set[x,"Strata"] = "Ghana_3"

x = which(CDI_2018_data_set[,"Strata"] == 1)
CDI_2018_data_set[x,"Strata"] = "CDI_1"
x = which(CDI_2018_data_set[,"Strata"] == 2)
CDI_2018_data_set[x,"Strata"] = "CDI_2"
x = which(CDI_2018_data_set[,"Strata"] == 3)
CDI_2018_data_set[x,"Strata"] = "CDI_3"

Ghana_CDI_2018_data_set = rbind(Ghana_2018_data_set, CDI_2018_data_set)
nrow(Ghana_CDI_2018_data_set)

##########################################################################
#The expert group weighting approach for the combining Ghana 2018 and CDI 2018 studies
##########################################################################

d_EG_Ghana_CDI_2018 <- svydesign(ids=~EA_2018+hhid, strata=~Strata, 
								fpc=~FPC_EA_within_stratum+FPC_HH_within_EA, weights=~EG_weight,
									data = Ghana_CDI_2018_data_set, nest = FALSE, check.strata = FALSE)

d_EG_Ghana_CDI_2018

#A histogram of the EG weights
hist(weights(d_EG_Ghana_CDI_2018), breaks=80, col="lightblue", 
	main="Histogram of Expert Group Suggested Sampling Weights", xlab="Sampling Weight")

#The total number of children
svytotal(~ones, d_EG_Ghana_CDI_2018)

#The estimates of the means 
svymean(~cocoachild_f, d_EG_Ghana_CDI_2018)
svymean(~cl_cocoa, d_EG_Ghana_CDI_2018)
svymean(~haz_cocoa, d_EG_Ghana_CDI_2018)

#The estimates of the standard errors
svytotal(~cocoachild_f, d_EG_Ghana_CDI_2018)
svytotal(~cl_cocoa, d_EG_Ghana_CDI_2018)
svytotal(~haz_cocoa, d_EG_Ghana_CDI_2018)

##The Rao-Wu Bootstrap Routine to properly calculate standard errors
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_EG_Ghana_CDI_2018 <- as.svrepdesign(design=d_EG_Ghana_CDI_2018, type="subbootstrap", replicates=2000)
d_boot_EG_Ghana_CDI_2018

d_boot_EG_Ghana_CDI_2018$type
d_boot_EG_Ghana_CDI_2018$rho
d_boot_EG_Ghana_CDI_2018$scale
d_boot_EG_Ghana_CDI_2018$rscales

#The total number of children
svytotal(~ones, d_boot_EG_Ghana_CDI_2018)

#The bootstrap estimates
svymean(~cocoachild_f, d_boot_EG_Ghana_CDI_2018)
svymean(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)
svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)

#The bootstrap estimates
svytotal(~cocoachild_f, d_boot_EG_Ghana_CDI_2018)
svytotal(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)
svytotal(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)

##########################################################################
#The sampling design exact weighting approach for the Ghana 2018 and CDI 2018 studies
##########################################################################

d_multi_Ghana_CDI_2018 <- svydesign(ids=~EA_2018+hhid+ids, strata=~Strata, 
								fpc=~FPC_EA_within_stratum+FPC_HH_within_EA+FPC_Ind_within_HH, 
									data = Ghana_CDI_2018_data_set, nest = FALSE, check.strata = FALSE)

#A histogram of the design weights
hist(weights(d_multi_Ghana_CDI_2018), breaks=80, col="lightblue", 
	main="Histogram of Sampling Design Weights", xlab="Sampling Weight")

#There may be a (small) issue with lonely PSUs
options(survey.lonely.psu="adjust")

#The total number of children
svytotal(~ones, d_multi_Ghana_CDI_2018)

#The estimates of the means (note that the standard errors are valid here)
svymean(~cocoachild_f, d_multi_Ghana_CDI_2018)
svymean(~cl_cocoa, d_multi_Ghana_CDI_2018)
svymean(~haz_cocoa, d_multi_Ghana_CDI_2018)

#The estimates of the totals (note that the standard errors are valid here)
svytotal(~cocoachild_f, d_multi_Ghana_CDI_2018)
svytotal(~cl_cocoa, d_multi_Ghana_CDI_2018)
svytotal(~haz_cocoa, d_multi_Ghana_CDI_2018)

##The Rao-Wu Bootstrap Routine
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_multi_Ghana_CDI_2018 <- as.svrepdesign(design=d_multi_Ghana_CDI_2018, type="subbootstrap", replicates=2000)
d_boot_multi_Ghana_CDI_2018

d_boot_multi_Ghana_CDI_2018$type
d_boot_multi_Ghana_CDI_2018$rho
d_boot_multi_Ghana_CDI_2018$scale
d_boot_multi_Ghana_CDI_2018$rscales

#The total number of children
svytotal(~ones, d_boot_multi_Ghana_CDI_2018)

#The bootstrap estimates
svymean(~cocoachild_f, d_boot_multi_Ghana_CDI_2018)
svymean(~cl_cocoa, d_boot_multi_Ghana_CDI_2018)
svymean(~haz_cocoa, d_boot_multi_Ghana_CDI_2018)

#The bootstrap estimates
svytotal(~cocoachild_f, d_boot_multi_Ghana_CDI_2018)
svytotal(~cl_cocoa, d_boot_multi_Ghana_CDI_2018)
svytotal(~haz_cocoa, d_boot_multi_Ghana_CDI_2018)

##########################################################################
#Hypothesis Testing for the Ghana and CDI 2018 studies combined
##########################################################################

#H_0: The means are the same
#H_A: The mean in 2018 is larger than that in 2008
#This entails a one-tailed test. For a two-tailed test, just double the p-value.
#For child working in cocoa, here is all the output we need
svymean(~cocoachild_f, d_unequal_Ghana_CDI_2008)
svymean(~cocoachild_f, d_boot_EG_Ghana_CDI_2018)
svymean(~cocoachild_f, d_boot_multi_Ghana_CDI_2018)

Z = (0.47967 - 0.31822) / (sqrt(2)*0.0185)
1-pnorm(Z)

Z = (0.48332 - 0.31822) / (sqrt(2)*0.0183)
1-pnorm(Z)

#For child labor in cocoa, here is all the output we need
svymean(~cl_cocoa, d_unequal_Ghana_CDI_2008)
svymean(~cl_cocoa, d_boot_EG_Ghana_CDI_2018)
svymean(~cl_cocoa, d_boot_multi_Ghana_CDI_2018)

Z = (0.44769 - 0.30901) / (sqrt(2)*0.017)
1-pnorm(Z)

Z = (0.45199 - 0.30901) / (sqrt(2)*0.0169)
1-pnorm(Z)

#For hazardous child labor in cocoa, here is all the output we need
svymean(~haz_cocoa, d_unequal_Ghana_CDI_2008)
svymean(~haz_cocoa, d_boot_EG_Ghana_CDI_2018)
svymean(~haz_cocoa, d_boot_multi_Ghana_CDI_2018)

Z = (0.42532 - 0.30404) / (sqrt(2)*0.0159)
1-pnorm(Z)

Z = (0.42935 - 0.30404) / (sqrt(2)*0.0155)
1-pnorm(Z)

##########################################################################
##########################################################################
#Combined Ghana 2018 and CDI 2018 studies - Household
##########################################################################
##########################################################################

x = which(Ghana_2018_data_set_HH[,"strata"] == 1)
Ghana_2018_data_set_HH[x,"strata"] = "Ghana_1"
x = which(Ghana_2018_data_set_HH[,"strata"] == 2)
Ghana_2018_data_set_HH[x,"strata"] = "Ghana_2"
x = which(Ghana_2018_data_set_HH[,"strata"] == 3)
Ghana_2018_data_set_HH[x,"strata"] = "Ghana_3"

x = which(CDI_2018_data_set_HH[,"strata"] == 1)
CDI_2018_data_set_HH[x,"strata"] = "CDI_1"
x = which(CDI_2018_data_set_HH[,"strata"] == 2)
CDI_2018_data_set_HH[x,"strata"] = "CDI_2"
x = which(CDI_2018_data_set_HH[,"strata"] == 3)
CDI_2018_data_set_HH[x,"strata"] = "CDI_3"

Ghana_CDI_2018_data_set_HH = rbind(Ghana_2018_data_set_HH, CDI_2018_data_set_HH)

nrow(Ghana_CDI_2018_data_set_HH)

##########################################################################
#The expert group weighting approach for the combining Ghana 2018 and CDI 2018 studies
##########################################################################

d_EG_Ghana_CDI_2018_HH <- svydesign(ids=~EA_2018, strata=~strata, 
								fpc=~FPC_EA_within_stratum, weights=~EG_weight,
									data = Ghana_CDI_2018_data_set_HH, nest = FALSE, check.strata = FALSE)

d_EG_Ghana_CDI_2018_HH

#A histogram of the EG weights
hist(weights(d_EG_Ghana_CDI_2018_HH), breaks=80, col="lightblue", 
	main="Histogram of Expert Group Suggested Sampling Weights", xlab="Sampling Weight")

#The total number of children
svytotal(~ones, d_EG_Ghana_CDI_2018_HH)

##The Rao-Wu Bootstrap Routine to properly calculate standard errors
#Set the seed to recreate all the (bootstrap) output that has been reported
set.seed(1)
d_boot_EG_Ghana_CDI_2018_HH <- as.svrepdesign(design=d_EG_Ghana_CDI_2018_HH, type="subbootstrap", replicates=2000)
d_boot_EG_Ghana_CDI_2018_HH

d_boot_EG_Ghana_CDI_2018_HH$type
d_boot_EG_Ghana_CDI_2018_HH$rho
d_boot_EG_Ghana_CDI_2018_HH$scale
d_boot_EG_Ghana_CDI_2018_HH$rscales

#The total number of children
svytotal(~ones, d_boot_EG_Ghana_CDI_2018_HH)

##############For Table 7 in the report
###A quick modification to the 2013 data sets and designs to produce point estimates
temp_CDI_2013_HH = subset(Full_data_set_HH, cdi==1 & Year=="2013")
C = ncol(temp_CDI_2013_HH)
temp_CDI_2013_HH[,C+1] = 1:nrow(temp_CDI_2013_HH)
colnames(temp_CDI_2013_HH)[C+1] = "ids"
d_temp_CDI_2013_HH = svydesign(ids=~ids, data=temp_CDI_2013_HH, weights=~weight)

temp_Ghana_2013_HH = subset(Full_data_set_HH, ghana==1 & Year=="2013")
C = ncol(temp_Ghana_2013_HH)
temp_Ghana_2013_HH[,C+1] = 1:nrow(temp_Ghana_2013_HH)
colnames(temp_Ghana_2013_HH)[C+1] = "ids"
d_temp_Ghana_2013_HH = svydesign(ids=~ids, data=temp_Ghana_2013_HH, weights=~weight)

temp_CDI_Ghana_2013_HH = rbind(temp_CDI_2013_HH, temp_Ghana_2013_HH)
d_temp_CDI_Ghana_2013_HH = svydesign(ids=~ids, data=temp_CDI_Ghana_2013_HH, weights=~weight)

#Set the working directory to where the tables will be saved
save.image(file="Updated Tables.RData")