/******************************************************************************
*Template from Elizabeth Pancotti
*Do file name: princ_meta_analysis_publish.do
*Project Title: The effect of principals' behaviors on schools and students: A meta-analysis
*Created by: David Liebowitz
*Created on: 11-1-18
*Last modified on: 4-16-19
*Last modified by: David Liebowitz
*Purpose: this .do file brings in the coded principal studies, assesses reliability of coding, ///
*			conducts meta-analysis, conducts meta-regression, tests for group heterogeneity and publication bias
*Inputs:
*	- meta_analytic_dataframe.dta

******************************************************************************/

clear all
set more off
set scheme s1mono
*Define directories
global projdir "C:\Users\daviddl\OneDrive - University of Oregon\Papers\Principal MetaAnalysis\"

global data "$projdir\Data"
global code "$projdir\Code"
global graph "$projdir\Graphs"
global table "$projdir\Tables"

cd "$projdir"

log using "$table\meta_an_log.smcl", replace

/*This portion commented out as not conducted on main dataset

*********************************************************************************************
*							Conducting IRR tests											*
*********************************************************************************************

*************GETTING DATA READY TO DO IRR******************

*Bring in Lorna's coding
use "$data\meta_analytic_dataframe_LP.dta", clear
gen case_id = _n


*Renaming variables to fit when renamed
ren cov_school_characteristics cov_school_char
ren	cov_student_charcteristics cov_stu_char 
ren cov_prin_characteristics cov_prin_char

*Generate count of number of outcomes per study
bys study_ID: egen cases = count(case_id)

*Generate flag for outcomes identified by rater in study
local outcomes stu_ach_assessment stu_ach_grades stu_ach_grad_rate stu_behavior ///
	teach_retention teach_well_being teach_practices prin_retention school_org_health
foreach x of local outcomes{
	egen `x'_count = tag(study_ID `x')
replace `x'_count=. if sample_size==.
	}

	
*Collapse to one row per study
collapse (mean) sample_size (max) ua* peer_review tech_report dissertation causal united_states ///
	school_elem school_secondary school_public school_private instructional_management internal_relations org_mangagement /// 
	admin_duties external_relations stu_ach_assessment_count stu_ach_grades_count stu_ach_grad_rate_count stu_behavior_count ///
	teach_retention_count teach_well_being_count teach_practices_count prin_retention_count school_org_health_count cov_school_prior_ach cov_school_char ///
	cov_student_prior_ach cov_stu_char cov_prin_char cases, by(study_ID)

gen rater=1

save "$data\Working\irr_LP.dta", replace

use "$data\meta_analytic_dataframe_DLIEBOWITZ.dta", clear
ren cov_school_characteristics cov_school_char
ren	cov_student_charcteristics cov_stu_char 
ren cov_prin_characteristics cov_prin_char

*Generate count of number of outcomes per study
bys study_ID: egen cases = count(case_id)

*Generate flag for outcomes identified by rater in study
local outcomes stu_ach_assessment stu_ach_grades stu_ach_grad_rate stu_behavior ///
	teach_retention teach_well_being teach_practices prin_retention school_org_health
foreach x of local outcomes{
	egen `x'_count = tag(study_ID `x')
replace `x'_count=. if sample_size==.
	}
	
*Collapse to one row per study
collapse (mean) sample_size (max) ua* peer_review tech_report dissertation causal united_states ///
	school_elem school_secondary school_public school_private instructional_management internal_relations org_mangagement /// 
	admin_duties external_relations stu_ach_assessment_count stu_ach_grades_count stu_ach_grad_rate_count stu_behavior_count teach_retention_count ///
	teach_well_being_count teach_practices_count prin_retention_count school_org_health_count cov_school_prior_ach cov_school_char ///
	cov_student_prior_ach cov_stu_char cov_prin_char cases, by(study_ID)

gen rater=2 
append using "$data\Working\irr_LP.dta"

*Fixing inconsistencies in study names
replace study_ID="clark_2009" if study_ID=="clark_2009 "
replace study_ID="cancio_2013" if study_ID=="canci_2013"
replace study_ID="jacob_2011" if study_ID=="jacob_2010"
replace study_ID="kruge_2007" if study_ID=="kruger_2007"
replace study_ID="kling_2011" if study_ID=="kling_2007"
replace study_ID="mazza_2014" if study_ID=="mazza_2013"

bys study_ID: egen study_in = total(sample_size), missing
drop if study_in==.

replace cov_school_prior_ach=. if sample_size==.
replace cov_school_char =. if sample_size==.
replace cov_student_prior_ach =. if sample_size==.
replace cov_stu_char =. if sample_size==.
replace cov_prin_char =. if sample_size==.


*Conduct Intraclass Correlation test for coded variables
local var sample_size ua_teacher ua_principal ua_student peer_review tech_report dissertation causal united_states ///
	school_elem school_secondary /*school_public school_private*/ instructional_management internal_relations org_mangagement /// 
	admin_duties external_relations stu_ach_assessment_count stu_ach_grades_count stu_ach_grad_rate_count stu_behavior_count teach_retention_count ///
	teach_well_being_count teach_practices_count prin_retention_count school_org_health_count cov_school_prior_ach cov_school_char ///
	cov_student_prior_ach cov_stu_char cov_prin_char
foreach y of local var {
	
	icc `y' study_ID rater, mixed
	gen icc_`y'=r(icc_avg)
	}
egen avg_icc=rowmean(icc*)

*Displaying cross-item average ICC
sum avg_icc, de
**EXCLUDING public/private where there is no "true" variation, ICC is .814 for 61 studies at least one of us coded
**If we include public/private it is .757

*/


*********************************************************************************************
*							Meta-analysis with master data									*
*********************************************************************************************


*********************************************
*		Prepare data						*
*********************************************

use "$data\meta_analytic_dataframe_master.dta", clear

*Generate single variable for predict
gen predict=.

replace predict=1 if instructional_management!=0
replace predict=2 if internal_relations!=0
replace predict=3 if org_mangagement!=0
replace predict=4 if admin_duties!=0
replace predict=5 if external_relations!=0
	lab def pred 1 "Instructional Management" 2 "Internal Relations" 3 "Organizational Management" ///
		4 "Administration" 5 "External Relations"
	lab val predict pred
tab predict, mi


*Gen single variable for outcome
gen outcome=.

replace outcome=1 if std_stu_ach_assessment!=.
replace outcome=2 if std_stu_ach_grades!=.
replace outcome=3 if std_stud_grad_rate!=.
replace outcome=4 if std_stu_behavior!=.
replace outcome=5 if std_teach_retention!=.
replace outcome=6 if std_teach_well_being!=.
replace outcome=7 if std_teach_practices!=.
replace outcome=8 if std_prin_retention!=.
replace outcome=9 if std_school_org_health!=.
	lab def outc 1 "Student Academic Achievement" 2 "Student Grades" 3 "Graduation rate" ///
		4 "Student Behavior/Attendance" 5 "Teacher Retention" 6 "Teacher Well-Being" 7 "Teaching Practices" ///
		8 "Principal Retention" 9 "School Organizational Health"
	lab val outcome outc
tab outcome, mi
	
ren instructional_management instr_mgmt
ren org_mangagement org_mgmt
ren cov_school_characteristics cov_school_chars
ren cov_student_charcteristics cov_stu_chars
ren cov_prin_characteristics cov_prin_chars
ren cov_teacher_characteristics cov_teach_chars
ren std_teach_retention std_teach_reten
ren std_stu_ach_assessment std_stu_ach_assess
	
*generate numeric ID variable
sort study_ID
egen num_id = group(study_ID)

ssc install estout
ssc install metan
ssc install robumeta
ssc install metatrim
*net install gr0033_1.pkg


gen var = std_standard_error*std_standard_error

label var std_stu_ach_assess "Student Achievement Effect Size" 
label var std_teach_well_being "Teacher Well-Being Effect Size" 
label var std_teach_practices "Teaching Practices Effect Size"
label var std_teach_reten "Teacher Retention Effect Size"
label var std_school_org_health "School Organizational Health Effect Size"

			
****************************************************************************
******					Descriptive work								****
****************************************************************************

split study_ID, p(_)
ren study_ID2 year
tab year
	destring year, replace
*Generating year ranges
gen year_band=.
	replace year_band=1 if year<=2005
	replace year_band=2 if year>2005 & year<=2010
	replace year_band=3 if year>2010 & year<=2015
	replace year_band=4 if year>2015 & year<=2020
		lab def years 1 "2001-2005" 2 "2006-2010" 3 "2011-2015" 4 "2016-2019"
		lab val year_band years
	tab year_band, mi

*Generate unit of analysis
gen unit_analysis=.
	replace unit_analysis=1 if ua_student==1
	replace unit_analysis=2 if ua_teacher==1
	replace unit_analysis=3 if ua_principal==1
	replace unit_analysis=4 if ua_student==0 & ua_teacher==0 & ua_principal==0
			gen ua_district=0
			replace ua_district =1 if unit_analysis==4
		lab def unit 1 "Student" 2 "Teacher" 3 "Principal" 4 "District"
		lab val unit_analysis unit
	tab unit_analysis, mi

*Generate pub type
gen pub_type=.
	replace pub_type=1 if peer_review==1
	replace pub_type=2 if tech_report==1
	replace pub_type=3 if dissertation==1
		lab def pub 1 "Peer-reviewed journal" 2 "Technical report" 3 "Dissertation/working paper"
		lab val pub_type pub
	tab pub_type, mi

*Label levels
lab def us 0 "International" 1 "United States"
lab def elem 0 "" 1 "Elementary"
lab def sec 0 "" 1 "Secondary"
lab def caus 0 "Observational" 1 "Causal (exogenous variation in treatment)"
	lab val united_states us
	lab val school_elem elem
	lab val school_secondary sec
	lab val causal caus

****************************************************************
****						TABLE 1							****
****					Study descriptives					****
****************************************************************

putexcel set "$table\descriptives_publish.xlsx", replace
*Publication year
putexcel A3=("Publication year")
tab year_band, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A4 = matrix(names)
	putexcel C4 = matrix(freq) D4 = "[" F4 = "]"   
	putexcel H4 = matrix(freq/r(N)) I4 = "["  K4 = "]", nformat(number_d2)
*Unit of analysis
putexcel A8=("Unit of analysis")
tab unit_analysis, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A9 = matrix(names)
	putexcel C9 = matrix(freq)  
	putexcel H9 = matrix(freq/r(N)), nformat(number_d2)
*Publication type
putexcel A13=("Publication type")
tab pub_type, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A14 = matrix(names)
	putexcel C14 = matrix(freq)  
	putexcel H14 = matrix(freq/r(N)), nformat(number_d2)
*Country of study
putexcel A17=("Country of study")
tab united_states, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A18 = matrix(names)
	putexcel C18 = matrix(freq)   
	putexcel H18 = matrix(freq/r(N)), nformat(number_d2)
*School level
putexcel A20=("School level")
tab school_elem, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A21 = matrix(names)
	putexcel C21 = matrix(freq)   
	putexcel H21 = matrix(freq/r(N)), nformat(number_d2)
tab school_secondary, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A23 = matrix(names)
	putexcel C23  = matrix(freq)  
	putexcel H23 = matrix(freq/r(N)), nformat(number_d2)
*Research design
putexcel A25=("Research design")
tab causal, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A26 = matrix(names)
	putexcel C26 = matrix(freq) 
	putexcel H26 = matrix(freq/r(N)), nformat(number_d2)
*Principal behavior
putexcel A28=("Principal behavior")
tab predict, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A29 = matrix(names)
	putexcel C29 = matrix(freq)  
	putexcel H29 = matrix(freq/r(N)), nformat(number_d2)
*Outcomes
putexcel A35=("Outcomes")
tab outcome, matcell(freq) matrow(names)
	matrix list freq
	matrix list names
	putexcel A36 = matrix(names)
	putexcel C36  = matrix(freq)
	putexcel H36 = matrix(freq/r(N)), nformat(number_d2)


/*Repeat for study-level*/  
tab year_band, gen(year_)
preserve
collapse (max) year_1 year_2 year_3 year_4 peer_review tech_report dissertation united_states school_elem school_secondary causal ///
	(sum) instr_mgmt internal_relations org_mgmt admin_duties external_relations ua_teacher ua_principal ua_student ua_district ///
	(count) std_stu_ach_asses std_stu_behavior std_teach_reten std_teach_well_being std_teach_practices std_prin_retention std_school_org_health, by(study_ID)
gen year_band=.
	replace year_band=1 if year_1==1
	replace year_band=2 if year_2==1
	replace year_band=3 if year_3==1
	replace year_band=4 if year_4==1
		tab year_band, mi
gen pub_type=.
	replace pub_type=1 if peer_review==1
	replace pub_type=2 if tech_report==1
	replace pub_type=3 if dissertation==1
		tab pub_type, mi
local outcome instr_mgmt internal_relations org_mgmt admin_duties external_relations ua_teacher ua_principal ///
	ua_student ua_district std_stu_ach_asses std_stu_behavior std_teach_reten std_teach_well_being ///
	std_teach_practices std_prin_retention std_school_org_health
	foreach o of local outcome{
		replace `o' = 1 if `o'!=0
		tab `o', mi
		}


*Publication year
tab year_band, matcell(freq)
	matrix list freq
	putexcel E4 = matrix(freq) 
	putexcel J4 = matrix(freq/r(N)), nformat(number_d2)
*Unit of analysis
tab ua_district, matcell(freq) 
	matrix list freq
	putexcel E11 = matrix(freq) 
	putexcel J11 = matrix(freq/r(N)), nformat(number_d2)
tab ua_principal, matcell(freq)
	matrix list freq
	putexcel E10 = matrix(freq) 
	putexcel J10 = matrix(freq/r(N)), nformat(number_d2)
tab ua_teacher, matcell(freq)
	matrix list freq
	putexcel E9 = matrix(freq) 
	putexcel J9 = matrix(freq/r(N)), nformat(number_d2)
tab ua_student, matcell(freq)
	matrix list freq
	putexcel E8 = matrix(freq) 
	putexcel J8 = matrix(freq/r(N)), nformat(number_d2)
*Publication type
tab pub_type, matcell(freq)
	matrix list freq
	putexcel E14 = matrix(freq) 
	putexcel J14 = matrix(freq/r(N)), nformat(number_d2)
*Country of study
tab united_states, matcell(freq)
	matrix list freq
	putexcel E18 = matrix(freq) 
	putexcel J18 = matrix(freq/r(N)), nformat(number_d2)
*School level
tab school_elem, matcell(freq)
	matrix list freq
	putexcel E21 = matrix(freq) 
	putexcel J21 = matrix(freq/r(N)), nformat(number_d2)
tab school_secondary, matcell(freq)
	matrix list freq
	putexcel E23 = matrix(freq) 
	putexcel J23 = matrix(freq/r(N)), nformat(number_d2)
*Research design
tab causal, matcell(freq)
	matrix list freq
	putexcel E26 = matrix(freq) 
	putexcel J26 = matrix(freq/r(N)), nformat(number_d2)
*Principal behavior
tab external_relations, matcell(freq)
	matrix list freq
	putexcel E32 = matrix(freq) 
	putexcel J32 = matrix(freq/r(N)), nformat(number_d2)
tab admin_duties, matcell(freq)
	matrix list freq
	putexcel E31 = matrix(freq) 
	putexcel J31 = matrix(freq/r(N)), nformat(number_d2)
tab org_mgmt, matcell(freq)
	matrix list freq
	putexcel E30 = matrix(freq) 
	putexcel J30 = matrix(freq/r(N)), nformat(number_d2)
tab internal_relations, matcell(freq)
	matrix list freq
	putexcel E29 = matrix(freq) 
	putexcel J29 = matrix(freq/r(N)), nformat(number_d2)
tab instr_mgmt, matcell(freq)
	matrix list freq
	putexcel E28 = matrix(freq) 
	putexcel J28 = matrix(freq/r(N)), nformat(number_d2)

*Outcomes
tab std_school_org_health, matcell(freq)
	matrix list freq
	putexcel E41 = matrix(freq) 
	putexcel J41 = matrix(freq/r(N)), nformat(number_d2) 
tab std_prin_retention, matcell(freq)
	matrix list freq
	putexcel E40 = matrix(freq) 
	putexcel J40= matrix(freq/r(N)), nformat(number_d2) 	
tab std_teach_practices, matcell(freq)
	matrix list freq
	putexcel E39 = matrix(freq) 
	putexcel J39= matrix(freq/r(N)), nformat(number_d2) 
tab std_teach_well_being, matcell(freq)
	matrix list freq
	putexcel E38 = matrix(freq) 
	putexcel J38= matrix(freq/r(N)), nformat(number_d2) 
tab std_teach_reten, matcell(freq)
	matrix list freq
	putexcel E37 = matrix(freq) 
	putexcel J37= matrix(freq/r(N)), nformat(number_d2)
tab std_stu_behavior, matcell(freq) 
	matrix list freq
	putexcel E36 = matrix(freq) 
	putexcel J36= matrix(freq/r(N)), nformat(number_d2)	
tab std_stu_ach_asses, matcell(freq)
	matrix list freq
	putexcel E35 = matrix(freq) 
	putexcel J35 = matrix(freq/r(N)), nformat(number_d2) 
	
restore



************************************
*	Create kdensity for outcomes	*
************************************

**Do this in two steps because of huge spread of school org health outcomes
local outcome std_stu_ach_assess std_teach_reten std_teach_well_being std_teach_practices
foreach y of local outcome{
	kdensity `y', ytitle("Density") legend(off) title("") yscale (r(0 2)) ylabel(0(1)2) xscale(r(-1 2)) xlabel(-1(1)2)
	graph save "$graph\kdensity_`y'", replace
	sum `y', de
	}

kdensity std_school_org_health, ytitle("Density") legend(off) title("") yscale(r(0 2)) ylabel(0(1)2) xscale(r(-2 4)) xlabel(-2(2)4) ///
	note("kernel = epanechnikov, bandwidth = 0.1405, axis length diff. to preserve distribution of outcomes", span size(vsmall))
	graph save "$graph\kdensity_std_school_org_health", replace
sum std_school_org_health, de	

	
*************************************
*		FIGURE 3					*
*************************************
	
cd "$graph"

graph combine kdensity_std_stu_ach_assess.gph kdensity_std_teach_well_being.gph ///
	kdensity_std_teach_practices.gph kdensity_std_school_org_health.gph
	graph export "$graph\kdensity_combined.emf", replace

	
*****************************************************
*		Create kdensity by outcome-predictor pair	*
*****************************************************

local predict 1 2 3
foreach x of local predict{

local outcome std_stu_ach_assess std_teach_well_being std_teach_practices std_school_org_health
foreach y of local outcome{
	kdensity `y' if predict==`x', ytitle("Density") legend(off) title("")
	graph save "$graph\kdensity_`x'_`y'", replace
	}
	}

local predict 4 5
foreach x of local predict{

local outcome std_stu_ach_assess std_school_org_health
foreach y of local outcome{
	kdensity `y' if predict==`x', ytitle("Density") legend(off) title("")
	graph save "$graph\kdensity_`x'_`y'", replace
	}
	}

*****************************
*		FIGURE A1			*
*****************************

cd "$graph"
graph combine kdensity_1_std_stu_ach_assess.gph kdensity_1_std_teach_well_being.gph kdensity_1_std_teach_practices.gph ///
	kdensity_1_std_school_org_health.gph kdensity_2_std_stu_ach_assess.gph kdensity_2_std_teach_well_being.gph kdensity_2_std_school_org_health.gph ///
	kdensity_3_std_stu_ach_assess.gph kdensity_3_std_school_org_health.gph kdensity_4_std_stu_ach_assess.gph kdensity_5_std_stu_ach_assess.gph ///
	, r(5) col(4) holes(7 10 11 14 15 16 18 19 20) 
	graph export "$graph\kdensity_predictor_combined.emf", replace

	
*********************************************************************************
*					Robust Variance Meta-Analysis								*
*********************************************************************************



*Empty models		
robumeta std_stu_ach_assess, study(num_id) variance(var) weighttype(random) rho(.90)			

*******************************************************
***			Create Table 2 - Main Analysis			***
*******************************************************

*Regression models
********************************************
****		Student achievement (1)		***
********************************************
*Testing for any variations with changes in rho, there is none!
/*forval i = 0.1(.1)0.9{
robumeta std_stu_ach_assess if instr_mgmt==1, study(num_id) variance(var) weighttype(random) rho(`i')
robumeta std_stu_ach_assess if internal_relations==1, study(num_id) variance(var) weighttype(random) rho(`i')
robumeta std_stu_ach_assess if org_mgmt==1, study(num_id) variance(var) weighttype(random) rho(`i')
	}
*/

eststo A_Instruct_Management: robumeta std_stu_ach_assess if instr_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_Internal_Relations: robumeta std_stu_ach_assess if internal_relations==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_Org_Management: robumeta std_stu_ach_assess if org_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_Admin_Duties: robumeta std_stu_ach_assess if admin_duties==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_External_Relations: robumeta std_stu_ach_assess if external_relations==1, study(num_id) variance(var) weighttype(random) rho(.9)


/**By unit of analysis
eststo Instructional_Management: robumeta std_stu_ach_assess ua_teacher ua_student if instr_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.9)
*eststo Internal_Relations: robumeta std_stu_ach_assess ua_teacher ua_student if internal_relations==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo Organizational_Management: robumeta std_stu_ach_assess ua_teacher ua_student if org_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo Admin_Duties: robumeta std_stu_ach_assess ua_teacher ua_student if admin_duties==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo External_Relations: robumeta std_stu_ach_assess ua_teacher ua_student if external_relations==1, study(num_id) variance(var) weighttype(random) rho(.9)
*/

********************************************
*****	Teacher retention (5)			****
********************************************
*eststo Instructional_Management: robumeta std_teach_reten if instr_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.9)

********************************************
*****	Teacher well-being (6)			****
********************************************
eststo B_Instruct_Management: robumeta std_teach_well_being if instr_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo B_Internal_Relations: robumeta std_teach_well_being if internal_relations==1, study(num_id) variance(var) weighttype(random) rho(.90)
tab std_teach_well_being study_ID if org_mgmt==1
tab std_teach_well_being study_ID if admin_duties==1
tab std_teach_well_being study_ID if external_relations==1

*********************************************
******	Teach practices (7)				*****
*********************************************
eststo C_Instruct_Management: robumeta std_teach_practices if instr_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.90)
tab std_teach_practices study_ID if internal_relations==1
tab std_teach_practices study_ID if org_mgmt==1
tab std_teach_practices study_ID if admin_duties==1
tab std_teach_practices study_ID if external_relations==1
			
*********************************************
******	School Org Health (9)			*****
*********************************************
eststo D_Instruct_Management: robumeta std_school_org_health if instr_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_Internal_Relations: robumeta std_school_org_health if internal_relations==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_Org_Management: robumeta std_school_org_health if org_mgmt==1, study(num_id) variance(var) weighttype(random) rho(.90)
*eststo D_External_Relations: robumeta std_school_org_health if external_relations==1, study(num_id) variance(var) weighttype(random) rho(.9)
tab std_school_org_health study_ID if admin_duties==1
tab std_school_org_health study_ID if external_relations==1

		estout * using "$table/main_results.txt", stats(N N_g, fmt(0) layout(@[@])) cells(b(star fmt(3)) se(par fmt(3))) replace
estimates clear




			
***************************************************************************
*		TABLES 3 and 4
*		Subgroup analysis 
*		-for predictors/outcomes for which we have sufficient coverage
***************************************************************************

***By level
local level school_elem school_secondary
	foreach x of local level{
eststo A_`x'_Instruct: robumeta std_stu_ach_assess if instr_mgmt==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_Internal: robumeta std_stu_ach_assess if internal_relations==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_Org: robumeta std_stu_ach_assess if org_mgmt==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_Admins: robumeta std_stu_ach_assess if admin_duties==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_External: robumeta std_stu_ach_assess if external_relations==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo B_`x'_Instruct: robumeta std_teach_well_being if instr_mgmt==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo B_`x'_Internal: robumeta std_teach_well_being if internal_relations==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo C_`x'_Instruct: robumeta std_teach_practices if instr_mgmt==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_`x'_Instruct: robumeta std_school_org_health if instr_mgmt==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_`x'_Internal: robumeta std_school_org_health if internal_relations==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_`x'_Org: robumeta std_school_org_health if org_mgmt==1 & `x'==1, study(num_id) variance(var) weighttype(random) rho(.90)
		estout * using "$table/level_results.txt", stats(N N_g, fmt(0) layout(@[@])) cells(b(star fmt(3)) se(par fmt(3))) replace
		}
estimates clear


***By country
local locale /*0*/ 1
	foreach x of local locale{
eststo A_`x'_Instruct: robumeta std_stu_ach_assess if instr_mgmt==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_Internal: robumeta std_stu_ach_assess if internal_relations==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_Org: robumeta std_stu_ach_assess if org_mgmt==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_Admins: robumeta std_stu_ach_assess if admin_duties==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_`x'_External: robumeta std_stu_ach_assess if external_relations==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.9)
eststo B_`x'_Instruct: robumeta std_teach_well_being if instr_mgmt==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.90)
eststo B_`x'_Internal: robumeta std_teach_well_being if internal_relations==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.90)
eststo C_`x'_Instruct: robumeta std_teach_practices if instr_mgmt==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_`x'_Instruct: robumeta std_school_org_health if instr_mgmt==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_`x'_Internal: robumeta std_school_org_health if internal_relations==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_`x'_Org: robumeta std_school_org_health if org_mgmt==1 & united_states==`x', study(num_id) variance(var) weighttype(random) rho(.90)
		estout * using "$table/locale_results.txt", stats(N N_g, fmt(0) layout(@[@])) cells(b(star fmt(3)) se(par fmt(3))) replace
		}
estimates clear	

/*****
******	Testing US-specific results in meta-regression framework
******
local predictor instr_mgmt internal_relations org_mgmt admin_duties external_relations
	foreach p of local predictor{
		gen usX`p' = united_states*`p'
	robumeta std_school_org_health usX`p' if `p'==1, study(num_id) variance(var) weighttype(random) rho(.9)
	}

****across all outcomes, there are no significant coefficients on interaction between US and predictors*/
	
***********************************************
*****	Within-study outcome correlations	***
***********************************************
bys study_ID: egen avg_std_ach=mean(std_stu_ach_assess)			
bys study_ID: egen avg_std_well=mean(std_teach_well_being)
bys study_ID: egen avg_std_tch_pract=mean(std_teach_practices)
bys study_ID: egen avg_std_schl_health=mean(std_school_org_health)

label var avg_std_ach "Student Achievement Effect Size"
label var avg_std_well "Teacher Well-Being Effect Size"
label var avg_std_tch_pract "Teacher Practices Effect Size"
label var avg_std_schl_health "School Organizational Health Effect Size"

/*
preserve
collapse (count) std_stu_ach_assess std_teach_reten std_teach_well_being std_teach_practices std_school_org_health, by(study_ID)
gen outcome=.
replace outcome=1 if std_stu_ach_assess!=0
replace outcome=2 if std_teach_reten!=0
replace outcome=3 if std_teach_well_being!=0
replace outcome=4 if std_teach_practices!=0
replace outcome=5 if std_school_org_health!=0
tab2 outcome outcome
restore
*/


************************************************************************
*							FIGURE 4								****
*		Student Achivement by Different Principal Behaviors			****
************************************************************************

local predict instr_mgmt internal_relations org_mgmt admin_duties external_relations
	foreach p of local predict{
	gen stu_achX`p'=.
	replace stu_achX`p' = `p' * std_stu_ach_assess if `p'!=0
	}
bys study_ID: egen avg_stu_ach_instr_mgmt=mean(stu_achXinstr_mgmt)
bys study_ID: egen avg_stu_ach_internal_rel=mean(stu_achXinternal_relations)
bys study_ID: egen avg_stu_ach_org_mgmt=mean(stu_achXorg_mgmt)
bys study_ID: egen avg_stu_ach_admin_duties=mean(stu_achXadmin_duties)
bys study_ID: egen avg_stu_ach_external_rel=mean(stu_achXexternal_relations)

*****
*PANEL A: Instructional Mgmt and Internal Relations on Student Achievement
*****
	
preserve
collapse avg_stu_ach_instr_mgmt avg_stu_ach_internal_rel avg_stu_ach_org_mgmt avg_stu_ach_admin_duties avg_stu_ach_external_rel var, by(study_ID)
	lab var avg_stu_ach_instr_mgmt "Instructional Management Effect Size"
	lab var avg_stu_ach_internal_rel "Internal Relations Effect Size"
corr avg_stu_ach_instr_mgmt avg_stu_ach_internal_rel
scatter avg_stu_ach_instr_mgmt avg_stu_ach_internal_rel [aw=var], ytitle("Instructional Management Effect Size") /*title("Leadership Effects on Student Achievement")*/  msymbol(circle_hollow) legend(off) ///
	|| lfit avg_stu_ach_instr_mgmt avg_stu_ach_internal_rel  [aw=var]
	graph export "$graph\instrmgmt_internrelXachieve.emf", replace
*scatter  avg_stu_ach_instr_mgmt avg_stu_ach_org_mgmt [aw=var], msymbol(circle_hollow) || lfit avg_stu_ach_instr_mgmt avg_stu_ach_org_mgmt [aw=var] 
restore


/*Only 4 studies...
*Teacher Well-Being by Different Principal Behaviors
local predict instr_mgmt internal_relations
	foreach p of local predict{
	gen well_beingX`p'=.
	replace well_beingX`p' = `p' * std_teach_well_being if `p'!=0
	}
bys study_ID: egen avg_wellbeing_instr_mgmt=mean(well_beingXinstr_mgmt)
bys study_ID: egen avg_wellbeing_internal_rel=mean(well_beingXinternal_relations)

preserve
collapse avg_wellbeing_instr_mgmt avg_wellbeing_internal_rel var, by(study_ID)
scatter avg_wellbeing_instr_mgmt avg_wellbeing_internal_rel [aw=var], msymbol(circle_hollow)
restore
*/

****
*PANEL B: Instructional Mgmt and Internal Relations on Organizational Health 
****

local predict instr_mgmt internal_relations
	foreach p of local predict{
	gen org_healthX`p'=.
	replace org_healthX`p' = `p' * std_school_org_health if `p'!=0
	}
bys study_ID: egen avg_orghealth_instr_mgmt=mean(org_healthXinstr_mgmt)
bys study_ID: egen avg_orghealth_internal_rel=mean(org_healthXinternal_relations)

preserve
collapse avg_orghealth_instr_mgmt avg_orghealth_internal_rel var, by(study_ID)
	lab var avg_orghealth_instr_mgmt "Instructional Management Effect Size"
	lab var avg_orghealth_internal_rel "Internal Relations Effect Size"
corr avg_orghealth_instr_mgmt avg_orghealth_internal_rel
	scatter avg_orghealth_instr_mgmt avg_orghealth_internal_rel [aw=var], ytitle("Instructional Management Effect Size") /*title("Leadership Effects on Organizational Health")*/ msymbol(circle_hollow) legend(off) ///
	|| lfit avg_orghealth_instr_mgmt avg_orghealth_internal_rel [aw=var]
		graph export "$graph\instrmgmt_internrelXorghealth.emf", replace
restore


*******************************************
***				FIGURE 5				***
*******************************************

*Achievement by School Health
preserve
collapse avg_std_ach avg_std_schl_health var, by(study_ID)
label var avg_std_ach "Student Achievement Effect Size"
label var avg_std_schl_health "School Organizational Health Effect Size"
corr avg_std_ach avg_std_schl_health
reg avg_std_ach avg_std_schl_health [aw=var]
scatter avg_std_ach avg_std_schl_health [aw=var], ytitle("Student Achievement Effect Size") msymbol(circle_hollow) legend(off) || lfit avg_std_ach avg_std_schl_health [aw=var]
	graph export "$graph\achieve_health.emf", replace
restore

/*Achievement by Teaching Practices
preserve
collapse avg_std_ach avg_std_tch_pract var, by(study_ID)
scatter avg_std_ach avg_std_tch_pract [aw=var], msymbol(circle_hollow) || lfit avg_std_ach avg_std_tch_pract [aw=var]
restore


/*Achievement by Teacher Well being
preserve
collapse avg_std_ach avg_std_well var, by(study_ID)
scatter avg_std_ach avg_std_well [aw=var], msymbol(circle_hollow) || lfit avg_std_ach avg_std_well [aw=var]
restore
*/

*Teacher Well-Being By Org Health
preserve
collapse avg_std_well avg_std_schl_health var, by(study_ID)
label var avg_std_well "Teacher Well-Being Effect Size"
label var avg_std_schl_health "School Organizational Health Effect Size"
scatter avg_std_well avg_std_schl_health [aw=var], ytitle("Teacher Well-Being Effect Size") msymbol(circle_hollow) legend(off) || lfit avg_std_well avg_std_schl_health [aw=var]
	graph export "$graph\wellbeing_health.emf", replace
restore

*Teacher Practices By Org Health
preserve
collapse avg_std_tch_pract avg_std_schl_health var, by(study_ID)
label var avg_std_tch_pract "Teacher Practices Effect Size"
label var avg_std_schl_health "School Organizational Health Effect Size"
scatter avg_std_tch_pract avg_std_schl_health [aw=var], ytitle("Teacher Practices Effect Size") msymbol(circle_hollow) legend(off) || lfit avg_std_tch_pract avg_std_schl_health [aw=var]
	graph export "$graph\teachpract_health.emf", replace
restore
*/

*******************************************************************************************
****		Comparing across principal behaviors in meta-regression	framework		*******
*******************************************************************************************

****************************
***		TABLE 5			****
****************************

*Conduct comparison of principal behaviors via metaregression
eststo Student_Achievement: robumeta std_stu_ach_assess internal_relations org_mgmt admin_duties external_relations, study(num_id) variance(var) weighttype(random) rho(.9)
eststo Teacher_Well_Being: robumeta std_teach_well_being internal_relations org_mgmt /*admin_duties*/ external_relations if admin_duties==0, study(num_id) variance(var) weighttype(random) rho(.9)
eststo Teaching_Practices: robumeta std_teach_practices internal_relations if external_relations==0, study(num_id) variance(var) weighttype(random) rho(.9)
eststo School_Org_Health: robumeta std_school_org_health internal_relations org_mgmt /*admin_duties*/ external_relations if admin_duties==0, study(num_id) variance(var) weighttype(random) rho(.9)
	estout * using "$table/metareg.txt", stats(N N_g, fmt(0) layout(@[@])) cells(b(star fmt(3)) se(par fmt(3))) replace
estimates clear



***********************************************************************************
***********			Sensitivity Analysis								***********
***********************************************************************************
*Funnel plots
confunnel std_stu_ach_assess std_standard_error, title(Student Achievement)
	graph export "$graph\achieve_funnel.emf", replace
confunnel std_teach_well_being std_standard_error, title(Teacher Well-Being)
	graph export "$graph\wellbeing_funnel.emf", replace
confunnel std_teach_practices std_standard_error, title(Teaching Practices)
	graph export "$graph\teachpractices_funnel.emf", replace
confunnel std_school_org_health std_standard_error, title(School Organizational Health)
	graph export "$graph\orghealth_funnel.emf", replace

*Trimmed
confunnel std_stu_ach_assess std_standard_error if std_stu_ach_assess<1		
			
			
			

****************************************************
******				TABLE 6						****
******		Sensitivity to top-bottom 5%		****
****************************************************

*Achievement
xtile pct_ach = std_stu_ach_assess, nq(99)
eststo A_Instructional_Management: robumeta std_stu_ach_assess if instr_mgmt==1 & pct_ach<95 & pct_ach>5, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_Internal_Relations: robumeta std_stu_ach_assess if internal_relations==1 & pct_ach<95 & pct_ach>5, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_Organizational_Management: robumeta std_stu_ach_assess if org_mgmt==1 & pct_ach<95 & pct_ach>5, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_Admin_Duties: robumeta std_stu_ach_assess if admin_duties==1 & pct_ach<95 & pct_ach>5, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_External_Relations: robumeta std_stu_ach_assess if external_relations==1 & pct_ach<95 & pct_ach>5, study(num_id) variance(var) weighttype(random) rho(.9)
			
*Well-Being			
xtile pct_well_being = std_teach_well_being, nq(99)
eststo B_Instructional_Management: robumeta std_teach_well_being if instr_mgmt==1 & pct_well_being<95 & pct_well_being>5, study(num_id) variance(var) weighttype(random) rho(.90)
eststo B_Internal_Relations: robumeta std_teach_well_being if internal_relations==1 & pct_well_being<95 & pct_well_being>5, study(num_id) variance(var) weighttype(random) rho(.90)

*Teaching Practices
xtile pct_pract = std_teach_practices, nq(99)
eststo C_Instructional_Management: robumeta std_teach_practices if instr_mgmt==1 & pct_pract<95 & pct_pract>5, study(num_id) variance(var) weighttype(random) rho(.90)

*Organizational Health
xtile pct_orghealth = std_school_org_health, nq(99)
eststo D_Instructional_Management: robumeta std_school_org_health if instr_mgmt==1 & pct_orghealth<95 & pct_orghealth>5, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_Internal_Relations: robumeta std_school_org_health if internal_relations==1 & pct_orghealth<95 & pct_orghealth>5, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_Organizational_Management: robumeta std_school_org_health if org_mgmt==1 & pct_orghealth<95 & pct_orghealth>5, study(num_id) variance(var) weighttype(random) rho(.90)
	
	estout * using "$table/sensitivity.txt", stats(N N_g, fmt(0) layout(@[@])) cells(b(star fmt(3)) se(par fmt(3))) replace
estimates clear

***************************************************
*****				TABLE 7						***
*****		Sensitivity to Trim-and-Fill		***
***************************************************	

preserve
collapse  std_stu_ach_assess std_teach_reten std_teach_well_being std_teach_practices std_school_org_health var (max) instr_mgmt internal_relations org_mgmt admin_duties external_relations , by(study_ID)		
*Student Achievement
metatrim std_stu_ach_assess var if instr_mgmt==1, var reffect funnel			
metatrim std_stu_ach_assess var if internal_relations==1, var reffect funnel				
metatrim std_stu_ach_assess var if org_mgmt==1, var reffect funnel
metatrim std_stu_ach_assess var if admin_duties==1, var reffect funnel
metatrim std_stu_ach_assess var if external_relations==1, var reffect funnel

/*Exclude these since no studies found "missing"
*Teacher Well-Being
metatrim std_teach_well_being var if instr_mgmt==1, var reffect funnel			
metatrim std_teach_well_being var if internal_relations==1, var reffect funnel			



*Teaching Practices
metatrim std_teach_practices var if instr_mgmt==1, var reffect funnel			

*School Organizational Health				
metatrim std_school_org_health var if instr_mgmt==1, var reffect funnel			
metatrim std_school_org_health var if internal_relations==1, var reffect funnel			
metatrim std_school_org_health var if org_mgmt==1, var reffect funnel			
*/
			
restore	
*Hand save these to trim_fill.xls file as metatrim does not save sufficient info in ereturn

*************************************************************************
******						Table A2								*****
******		Sensitivity to Bivariate and Partial Correlation		*****
*************************************************************************
gen bivariate=0
replace bivariate=1 if mean_difference==1 | correlation_coef==1

gen partial=0
replace partial=1 if regresion_coef==1 | standardized_slope==1

*Student Achievement
eststo A_Instructional_Management: robumeta std_stu_ach_assess if instr_mgmt==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A1_Instructional_Management: robumeta std_stu_ach_assess if instr_mgmt==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.9)			
eststo A_Internal_Relations: robumeta std_stu_ach_assess if internal_relations==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A1_Internal_Relations: robumeta std_stu_ach_assess if internal_relations==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.9)
*eststo A_Org_Management: robumeta std_stu_ach_assess if org_mgmt==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A1_Org_Management: robumeta std_stu_ach_assess if org_mgmt==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A_Admin_Duties: robumeta std_stu_ach_assess if admin_duties==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A1_Admin_Duties: robumeta std_stu_ach_assess if admin_duties==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.9)
*eststo A_External_Relations: robumeta std_stu_ach_assess if external_relations==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.9)
eststo A1_External_Relations: robumeta std_stu_ach_assess if external_relations==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.9)

*Teacher Well-Being
eststo B_Instruct_Management: robumeta std_teach_well_being if instr_mgmt==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo B1_Instruct_Management: robumeta std_teach_well_being if instr_mgmt==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo B_Internal_Relations: robumeta std_teach_well_being if internal_relations==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo B1_Internal_Relations: robumeta std_teach_well_being if internal_relations==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.90)

*Teaching Practices
eststo C_Instruct_Management: robumeta std_teach_practices if instr_mgmt==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo C1_Instruct_Management: robumeta std_teach_practices if instr_mgmt==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.90)

*School Organizational Health
eststo D_Instruct_Management: robumeta std_school_org_health if instr_mgmt==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D1_Instruct_Management: robumeta std_school_org_health if instr_mgmt==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D_Internal_Relations: robumeta std_school_org_health if internal_relations==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D1_Internal_Relations: robumeta std_school_org_health if internal_relations==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.90)
*eststo D_Org_Management: robumeta std_school_org_health if org_mgmt==1 & bivariate==1, study(num_id) variance(var) weighttype(random) rho(.90)
eststo D1_Org_Management: robumeta std_school_org_health if org_mgmt==1 & partial==1, study(num_id) variance(var) weighttype(random) rho(.90)
		
		estout * using "$table/bivariate_partial.txt", stats(N N_g, fmt(0) layout(@[@])) cells(b(star fmt(3)) se(par fmt(3))) replace
		
estimates clear

				

log close