*W4E SanQoL-5 paper analysis
*20 Oct '26

*--------------------------------------------------------
** Set directories and import
*--------------------------------------------------------

		if c(os) == "Windows" {
		global github "C:\Users\Blessings White\Documents\GitHub\W4E_trial"
	}
	else if c(os) == "MacOSX" {
		global github "/Users/`c(username)'/Documents/GitHub/W4E_trial"
	}

clear
cd "$github/2_input"
use w4e_sanqol5_data
cd "$github/3_output"

**------------------------------------------
* Outcomes to analyse
**------------------------------------------
	*1 - SanQoL-5 index (0-1)
	*2 - SanQoL-5 attribute scores (0-2)
	*3 - SanVAS (0-100)

**------------------------------------------
* Outcome 1 - SanQoL-5
**------------------------------------------

********* Means
		eststo clear
		 
		 forvalues i = 1(1)3 {
		 putexcel set means_sanqol_`i'.xlsx, sheet(sheet1) replace
		 mean sanqol if g1_name_of_ta == `i', over(baseline_endline)
		 matrix sanqol_`i'_mean = r(table)
		 putexcel A1 = `i'
		 putexcel A2 = matrix(sanqol_`i'_mean)
		} 

********* Main analysis

	*partially adjusted
	eststo: mixed sanqol ib2.g1_name_of_ta##i.baseline_endline comm_hhmed || village:, vce(robust)
	estat icc

	*fully adjusted
	eststo: mixed sanqol ib2.g1_name_of_ta##i.baseline_endline n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || village:, vce(robust)
	estat icc

********* CLTS+ relative to CLTS
	*partially adjusted
	eststo: mixed sanqol ib3.g1_name_of_ta##i.baseline_endline comm_hhmed || village:, vce(robust)
	estat icc

	*fully adjusted
	eststo: mixed sanqol ib3.g1_name_of_ta##i.baseline_endline n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || village:, vce(robust)
	estat icc

	*the two models give the same effect and standard error for SanQoL-5 

esttab using reg_sanqol5_main.csv, star(* 0.10 ** 0.05 *** 0.01) replace se(4) b(4) nogap label

********* with gender interaction

	codebook gender_hhid
	*0 = female, 1 = male

	*partially adjusted
	eststo: mixed sanqol ib2.g1_name_of_ta##i.baseline_endline##gender_hhid ///
	comm_hhmed || village:, vce(robust)
	estat icc

	*fully adjusted
	eststo: mixed sanqol ib2.g1_name_of_ta##i.baseline_endline##gender_hhid ///
	n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || ///
	village:, vce(robust)
	estat icc


**------------------------------------------
* Outcome 2 - SanQoL-5 attribute scores
**------------------------------------------

eststo clear

local sq_atts disg heal priv sham safe

foreach i of local sq_atts {
    putexcel set means_sq_`i'.xlsx, sheet("Results") replace
    mean sq_`i', over(g1_name_of_ta baseline_endline)
    matrix mean_results = r(table)
    putexcel A1 = "Variable: sq_`i'", bold
    putexcel A2 = matrix(mean_results)
}

*partially adjusted
foreach i of local sq_atts {
	eststo: mixed sq_`i' ib2.g1_name_of_ta##i.baseline_endline comm_hhmed || village:, vce(robust)
	estat icc
 }

*fully adjusted
foreach i of local sq_atts {
	eststo: mixed sq_`i' ib2.g1_name_of_ta##i.baseline_endline n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || village:, vce(robust)
	estat icc
 }
 
 
 *fully adjusted relative to CLTS
foreach i of local sq_atts {
	eststo: mixed sq_`i' ib3.g1_name_of_ta##i.baseline_endline n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || village:, vce(robust)
	estat icc
 }
 
 *fully adjusted gender interaction
foreach i of local sq_atts {
	eststo: mixed sq_`i' ib2.g1_name_of_ta##i.baseline_endline##gender_hhid n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || village:, vce(robust)
	estat icc
 }
 
esttab using reg_sq_atts_main.csv, star(* 0.10 ** 0.05 *** 0.01) replace se(4) b(4) nogap label
 
**------------------------------------------
* Outcome 3 - SanVAS
**------------------------------------------

eststo clear

 forvalues i = 1(1)3 {
 putexcel set means_vas_sanitation_`i'.xlsx, sheet(sheet1) replace
 mean vas_sanitation if g1_name_of_ta == `i', over(baseline_endline)
 matrix vas_sanitation_`i'_mean = r(table)
 putexcel A1 = `i'
 putexcel A2 = matrix(vas_sanitation_`i'_mean)
} 

eststo clear

bysort g1_name_of_ta: tab baseline_endline, sum(vas_sanitation)

*partially adjusted
	eststo: mixed vas_sanitation ib2.g1_name_of_ta##i.baseline_endline comm_hhmed || village:, vce(robust)
	estat icc


*fully adjusted
	eststo: mixed vas_sanitation ib2.g1_name_of_ta##i.baseline_endline n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || village:, vce(robust)
	estat icc

*CLTS+ relative to CLTS
	*partially adjusted
	eststo: mixed vas_sanitation ib3.g1_name_of_ta##i.baseline_endline comm_hhmed || village:, vce(robust)
	estat icc


	*fully adjusted
	eststo: mixed vas_sanitation ib3.g1_name_of_ta##i.baseline_endline n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || village:, vce(robust)
	estat icc
	


esttab using reg_sanvas_main.csv, star(* 0.10 ** 0.05 *** 0.01) replace se(4) b(4) nogap label

********* with gender interaction

	codebook gender_hhid
	*0 = female, 1 = male
	*replace gender_hhid = . if gender_hhid == 10 // for people where gender missing

	*partially adjusted
	eststo: mixed vas_sanitation ib2.g1_name_of_ta##i.baseline_endline##gender_hhid ///
	comm_hhmed || village:, vce(robust)
	estat icc

	*fully adjusted
	eststo: mixed vas_sanitation ib2.g1_name_of_ta##i.baseline_endline##gender_hhid ///
	n_hh_members disability_overall i.quintile water_source_onsite_onplot comm_hhmed || ///
	village:, vce(robust)
	estat icc