* Stata syntax for Chapter 6 analyses (Version 17).

* Import Autism data.
insheet using http://www-personal.umich.edu/~bwest/autism.csv, clear

gen age_2 = age - 2
gen age_2sq = age_2 * age_2

* Model 6.1 (REML)
mixed vsae ib3.sicdegp age_2 age_2sq ///
   ib3.sicdegp#c.age_2 ib3.sicdegp#c.age_2sq ///
   || childid: age_2 age_2sq, covariance(unstruct) ///
   variance reml dfmethod(kroger)
   
* Model 6.2
mixed vsae ib3.sicdegp age_2 age_2sq ///
   ib3.sicdegp#c.age_2 ib3.sicdegp#c.age_2sq ///
   || childid: age_2 age_2sq, covariance(unstruct) ///
   variance reml dfmethod(kroger) noconst
   
estat ic

* Model 6.2A
mixed vsae ib3.sicdegp age_2 age_2sq ///
   ib3.sicdegp#c.age_2 ib3.sicdegp#c.age_2sq ///
   || childid: age_2, covariance(unstruct) ////
   variance reml dfmethod(kroger) noconst

* Mixture-based likelihood ratio test
di (-2 * -2349.6013) - (-2 * -2307.6378)
di 0.5 * chi2tail(1,83.927) + 0.5 * chi2tail(2,83.927)

* Model 6.2 (ML)
mixed vsae ib3.sicdegp age_2 age_2sq ///
   ib3.sicdegp#c.age_2 ib3.sicdegp#c.age_2sq ///
   || childid: age_2 age_2sq, covariance(unstruct) ///
   variance mle noconst
   
est store model6_2_ml

* Model 6.3 (ML)
mixed vsae ib3.sicdegp age_2 age_2sq ///
   ib3.sicdegp#c.age_2 ///
   || childid: age_2 age_2sq, covariance(unstruct) ///
   variance mle noconst
   
est store model6_3_ml

lrtest model6_2_ml model6_3_ml

* Model 6.3 (REML)
mixed vsae ib3.sicdegp age_2 age_2sq ///
   ib3.sicdegp#c.age_2 ///
   || childid: age_2 age_2sq, covariance(unstruct) ///
   variance reml dfmethod(kroger) noconst
estat ic

mixed vsae ib3.sicdegp age_2 ///
   c.age_2#c.age_2 ib3.sicdegp#c.age_2 ///
   || childid: age_2 age_2sq, covariance(unstruct) ///
   variance reml dfmethod(kroger) noconst

* Visualize interaction
margins, at(age_2=(0(1)11)) by(sicdegp)
marginsplot