**********************************************************************
**********************************************************************
**** ****
**** Chapter 13 -- Imperfect Hierarchies ****
**** Stata do file to reproduce examples ****
**** in Chapter 13 of ****
**** Snijders, Tom A.B., and Bosker, Roel J. ****
**** Multilevel Analysis: ****
**** An Introduction to Basic and Advanced Multilevel Modeling,****
**** second edition ****
**** London etc.: Sage Publishers, 2012 ****
**** ****
**** ****
**********************************************************************
**********************************************************************
** Crossed random effects **
* The dataset from the book is not available, so we are using a dataset from
* Rabe-Hesketh and Skrondal (2008)
* Multilevel and Longitudinal Modeling Using Stata, Second Edition *
* See also chapter 8.4 for this particular example *
capture log close
log using "J:\Multilevel\Chp13_Crossed_Random_Effects.log", replace
clear all
set more off
*We get the data online
use http://www.stata-press.com/data/mlmus2/fife.dta, clear
/*The dependent var is attainment score at age 16 for pupils that attended various combinations
of primary and secondary schools.
pid - primary school identifier
sid - secondary school identifier
vrq - verbal reasoning score
sex - sex */
sum pid sid vrq sex
center vrq
*We explore the special nesting structure first.
*The tag() function lets us define unique primary/secondary school combinations:
egen pick_comb=tag(pid sid)
*We count the unique sid combinations for each primary school.
*To how many different secondary schools did a primary school sent the pupils?
egen numsid=total(pick_comb), by(pid)
sort pid sid
list pid sid numsid if pick_comb & pid <10 in 1/100
*Kids who went to primary school 1 ended up in 3 different secondary schools: 1, 9, 18.
*Now we want to know for each primary school, up to how many secondary schools were sent to.
egen pick_pid=tag(pid)
tab numsid if pick_pid
*90% of primary schools send the kids to 3 different secondary schools.
*First we try a normal 2-level specification (nested in primary schools)
*The trick is to pretend that there exists a new level, in which all observations are nested. (_all)
xtmixed attain || pid: , mle var
est store mod1
*Crossed effects
xtmixed attain || _all: R.sid || pid:, mle var
est store mod2
esttab mod1 mod2, se wide nostar ///
transform(ln*: exp(2*@) exp(2*@)) ///
eqlabels(, none)
*Very few changes. A little bit of variation is now distributed to secondary schools.
*Different ICC for different primary/secondary school combination.
*Same primary, but different secondary schools:
dis 1.12/(.35+1.12+8.11)
*.11691023
*Same secondary, but not the same primary school:
dis .35/(.35+1.12+8.11)
*.03653445
*Same primary AND secondary school:
dis (.35+1.12)/(.35+1.12+8.11)
*.15344468
*Crossed effects, fixed explanatory variables. Verbal aptitude and gender.
xtmixed attain c_vrq sex || _all: R.sid || pid:, mle var
*Let's try to include a random slope (c_vrq) on both levels.
*It's a bit complicated. For primary schools, we just include it after pid: .
* But for the secondary school level, we have to generate interaction terms of c_vrq and all secondary schools.
*Create dummies for each secondary school
tab sid, gen(id_sid)
unab idvar: id_sid*
foreach v of local idvar{
gen inter`v' = c_vrq*`v'
}
*Now we can include a random slope on both levels
*We restrict the covariance between all interaction dummies to 0.
*We also could specify random slopes for different variables on each level.
xtmixed attain c_vrq sex || _all: R.sid || ///
_all:inter*, cov(identity) nocons || pid: c_vrq, var mle
/*No significant random slopes for c_vrq. Maybe there are gender differences?
drop inter*
unab idvar: id_sid*
foreach v of local idvar{
gen inter`v' = sex*`v'
}
xtmixed attain c_vrq sex || _all: R.sid || ///
_all:inter*, cov(identity) nocons || pid: sex, var mle
*Mean Verbal apptitude score for primary schools
bysort pid: egen mvrq=mean(c_vrq)
xtmixed attain c_vrq c.mvrq##i.sex || _all: R.sid || ///
_all:inter*, cov(identity) nocons || pid: sex, var mle
*/
*Clearly, there are more things to discover, such as interactions of random intercepts.
*Specific combinations of primary/secondary schools might matter. But that would be too much now...
capture log close