*set up the grouping variable as numeric 
generate ms_gp = 0 if msmoke == 0
replace ms_gp = 1 if msmoke > 0 & msmoke <= 20
replace ms_gp = 2 if msmoke > 20

graph box weight,over(ms_gp) b1title("Packs of cigarettes mother smoked")

*do one-way analysis using anova command
anova weight ms_gp
anova,regress /* get contrasts vs. last group*/
lincom(_b[ms_gp[1]]-_b[ms_gp[2]]) /*non- vs. light-smoke contrast*/

*get same analysis using xi command (parameter constraints differ)
xi:regress weight i.ms_gp
testparm _Ims_gp*
lincom(_b[_Ims_gp_1]-_b[_Ims_gp_2])

*fit pre-pregnancy weight as covariate using anova command
anova weight ms_gp mweight,cont(mweight)
anova,regress /*get contrasts vs. last group*/
lincom(_b[ms_gp[1]]-_b[ms_gp[2]]) /*non- vs. light-smoke contrast*/

*repeat using xi (get different contrasts)
xi: regress weight i.ms_gp mweight
testparm _Ims_gp*
lincom(_b[_Ims_gp_1]-_b[_Ims_gp_2])

*get fitted lines plotted on scatter plots
predict yhat
twoway (scatter weight mweight, mlabel(ms_gp)) ///
       (line yhat mweight if ms_gp==0) (line yhat mweight if ms_gp==1) ///
       (line yhat mweight if ms_gp ==2), legend(off)


*fit unequal slopes
anova weight ms_gp mweight ms_gp*mweight,cont(mweight)
drop yhat
predict yhat
twoway (scatter weight mweight, mlabel(ms_gp)) ///
       (line yhat mweight if ms_gp==0) (line yhat mweight if ms_gp==1) ///
       (line yhat mweight if ms_gp ==2), legend(off)
anova,regress

*do 3 separate linear regressions to see same as above
regress weight mweight if ms_gp==0
regress weight mweight if ms_gp==1
regress weight mweight if ms_gp==2

*fit mother's age, height, and gest length in addition to mother's age
anova weight ms_gp mweight mage mheight gest,cat(ms_gp)
anova,regress
lincom(_b[ms_gp[1]]-_b[ms_gp[2]]) /*non- vs. light-smoke contrast*/

*Do simulated examples
*Load data outside of do file
graph box y,over(group)
anova y group /* ignore covariate */
anova y group x,cont(x) /* include covariate */
predict yhat
twoway (scatter y x, mlabel(group)) ///
       (line yhat x if group==1) (line yhat x if group==2) ///
       (line yhat x if group==3), legend(off)

tabstat y,by(group) stat(mean semean)
anova y group x, cont(x)
adjust x, by(group) se