///////KINDERGARTEN STUDY///////

*** CODE TO REPRODUCE FOREST PLOT ***
// Adapted from code by Dr Baptiste Leurent: https://figshare.com/articles/journal_contribution/Subgroup_forest_plot_in_Stata/686060

// One dataset needed:
// ~/WISE_KG_study_pupil_dataset_WIDE.dta (pupil-level wide dataset - one row per pupil)

//---------------------------------------------//
** 1) Create table with results **
//---------------------------------------------//

use WISE_KG_study_pupil_dataset_WIDE, clear
keep if eligible==1 & completed_1==1 // restrict analysis to pupils eligible by age (aged 3–6) and enrolled in the study (successfully completed call 1)

keep days_resp days_diar days_absent total_days school_days arm strata pupil_gender pupil_age subcity school_ID pupil_ID

gen overall=1  // Create dummy variable to have overall effect on plot	

local varlist1 overall pupil_gender
local varlist2 days_absent days_diar days_resp

* 1.1 Create regression models for each subgroup and save results *

//foreach out in `varlist2' {
	foreach var in `varlist1' {
		levelsof `var', local(levels_`var') 
		foreach l of local levels_`var' {
			meglm days_absent i.arm i.strata if `var'==`l' || school_ID:, family(binomial school_days) or // Write here analysis model wanted
			matrix res_days_absent_`var'_`l' = r(table)
			local n_days_absent_`var'_`l' = `e(N)'
			meglm days_diar i.arm i.strata if `var'==`l' || school_ID:, family(binomial total_days) or
			matrix res_days_diar_`var'_`l' = r(table)
			local n_days_diar_`var'_`l' = `e(N)'
			meglm days_resp i.arm i.strata if `var'==`l' || school_ID:, family(binomial total_days) or
			matrix res_days_resp_`var'_`l' = r(table)
			local n_days_resp_`var'_`l' = `e(N)'
			}
		}
//	}		
		
* 1.2 Store results in a dataset *
		
clear
set obs 50
gen name=""
gen wgt=.
gen dif=.
gen low=.
gen up=.
gen outcome=.
gen group=.
		
local line=1  // Starting value for line number to copy results
local outcome=1 // Starting value for outcome numbering

foreach out in `varlist2' {
	local group=1 // Starting value for subgroup numbering
	foreach var in `varlist1' {
		foreach l of local levels_`var' {
			qui: replace name = "`out'_`var'_`l'"   if _n==`line'
			qui: replace wgt=`n_`out'_`var'_`l''   if _n==`line'										
			qui: replace dif=res_`out'_`var'_`l'[1,2]  if _n==`line' // Check this correponds to wanted results in saved results matrix
			qui: replace low=res_`out'_`var'_`l'[5,2]  if _n==`line'   
			qui: replace up=res_`out'_`var'_`l'[6,2]  if _n==`line'   
			qui: replace group=`group'   if _n==`line'  
			qui: replace outcome=`outcome'   if _n==`line'  
			qui: local line=`line'+1				
			}
		local group=`group'+1
		}
	local outcome=`outcome'+1
	}

drop if name=="" // Drop empty rows
list, noobs sepby(group)

gen pun1=" ("
gen pun2=", "
gen pun3=")"
egen OR_CI=concat(dif pun1 low pun2 up pun3), format(%8.2f) // Create data label

drop pun*
		
//---------------------------------------------//
** 2) Forest plot **
//---------------------------------------------//
	
* 2.1 Position on vertical axis for each subgroup *

gen y=-1 if _n==1 
replace y=y[_n-1]-1-0.8*(outcome[_n]!=outcome[_n-1]) if _n>1 // Go to next line + 0.8 space when changing group

* Store subgroups label 
	
	*Subgroups
	
gen grplb="" 
local line = 1
foreach group in "Overall" "Girls" "Boys" "Overall" "Girls" "Boys" "Overall" "Girls" "Boys" {	
	qui: replace grplb="`group'" if _n==`line'
	local line=`line'+1	
	}
	
	*Variable (=set of subgroups)	
			
gen varlb=""
gen yvar=.  //Position on y axis for variable name
local group = 0
foreach title in "Caregiver-reported absence" "Caregiver-reported diarrhoea" "Caregiver-reported respiratory illness" {
	local group=`group'+1
	dis "`title'  " `group'
	qui: bysort outcome: replace varlb="`title'" if outcome==`group' & _n==1
	qui: bysort outcome: replace yvar=y+0.8 if outcome==`group' & _n==1 // Place variable name 0.8 above 1st subgroup
	}		

gen xvar=0.015  //Horizontal position of variable label
gen xgrp=0.018  //Horizontal position of subgroups label
gen xCI=4.5 //Horizontal position of OR and CI label
local ref1 = 1 // No-effect line

gen CIlb="aOR (95% CI)" if group==1 & outcome==1
gen yCI=0.4 if group==1 & outcome==1

* 2.2 Plot *

graph twoway  (pci 0.5 1 -11 1, lpattern(dash))  ///
	(scatter  y dif [aweight=wgt], msymbol(square) msize(0.8)   color(black))  (rspike low up y, horizontal color(black))  ///	
	(scatter  yvar xvar, msymbol(i) mlabel(varlb) mlabpos(3) mlabc(black) mlabsize(medsmall) ) ///
	(scatter  yCI xCI, msymbol(i) mlabel(CIlb) mlabpos(3) mlabc(black) mlabsize(medsmall) ) ///
	(scatter  y xgrp, msymbol(i) mlabel(grplb) mlabpos(3) mlabc(black) mlabsize(medsmall) ) ///
	(scatter  y xCI, msymbol(i) mlabel(OR_CI) mlabpos(3) mlabc(black) mlabsize(medsmall) ) ///
	, yscale(range(-11 1)) yscale(off) ylabel(none) ///
	xscale(log range(0.05 15)) xlabel(0.25 0.5 1 2 4)  xtitle("") ///
	legend(rows(2) holes(3) region(lcolor(black)) order(2 3 1) size(10pt)  position(11) ring(1) lab(2 "Odds Ratio") lab(3 "95% CI") lab(1 "No effect") ) ///
	text(2 1 " Favours Intervention      Favours Control      ")
// 						title("Effect on primary outcomes by gender", size(medium)) xsize(10) ysize(8)