use http://www.ats.ucla.edu/stat/stata/examples/mlm_ma_hox/popular.dta, clear

* EDA
hist popular, freq


* Random intercept model 
xtreg popular, re i(school) mle 
xtmixed popular || school:, mle
gllamm popular, i(school) adapt 

* get xtmixed to give you variances, not sd
xtmixed popular || school:, mle var

* EDA of relationship between popular and texp

* ignore clustering
lowess popular texp, bw(0.5) jitter(3)

* the average school score and texp
sort school
by school: egen mscore_sch=mean(popular)
lowess mscore_sch texp, ytitle("School Average Popularity score") xtitle("Teacher Experience") bw(0.5)

* create an indicator for 'girl' out of the sex variable
rename sex girl

* EDA of relationship between popular and pupil gender

sort school girl
by school girl:egen mscore_sch_girl=mean(popular)
sort girl
graph box mscore_sch_girl,over(girl) ytitle("School Average Score")

* interaction between girl and texp? 
lowess mscore_sch_girl texp, by(girl)

* Random intercept with the fixed effects of texp and girl
xtmixed popular texp girl || school:, mle


* between-school heterogeneity in the gender effect?
sort school girl
by school: gen meanpop_boys = mscore_sch_girl[1] 
gen boy = 1-girl
sort school boy
by school: gen meanpop_girls = mscore_sch_girl[1] 

gen genddiff = meanpop_girls - meanpop_boys 
replace genddiff =. if pupil!=1
hist genddiff, norm freq xtitle(School-specific difference between female and male average scores)


*Random coeffecient model, both intercept and girl are random effects 
xtmixed popular texp girl || school: girl, cov(unstr) mle 

*Random coeffecient model with cross-level interaction
gen girlXtexp = girl*texp
xtmixed popular texp girl girlXtexp || school: girl, cov(unstr) mle