* do-file for additional exercise 6.1 (VHM 802)
version 16 /* works also with versions 14-15 */
set more off
cd "r:\"

import delimited using hs06_1.csv, clear
* descriptive statistics and simple oneway ANOVA
tabstat iron, by(diet) stat(n mean sd min max)
graph box iron, over(diet)
oneway iron diet, tabulate
generate lniron=ln(iron)
tabstat lniron, by(diet) stat(n mean sd min max)
graph box lniron, over(diet)
oneway lniron diet, tabulate

xi: boxcox iron i.diet
generate invrootiron=-1/sqrt(iron)
tabstat invrootiron, by(diet) stat(n mean sd min max)
graph box invrootiron, over(diet)
oneway invrootiron diet, tabulate

encode diet, gen(Diet) /* anova allows only numerical predictors */
anova lniron Diet
* semi-manual calculation of SS and Scheffe's test, note 5=#groups
scalar mse=e(rmse)^2 /* MSE of last model */
lincom 1.Diet-2.Diet
scalar tval=r(estimate)/r(se)
scalar sstrt=e(mss)
di "SS: " tval^2*mse "  in %: " tval^2*mse/sstrt*100
di "Scheffe test:  F = " tval^2/(5-1) "  P = " Ftail(5-1,e(df_r),tval^2/(5-1))
lincom 1.Diet+2.Diet-2*3.Diet
scalar tval=r(estimate)/r(se)
di "SS: " tval^2*mse "  in %: " tval^2*mse/sstrt*100
di "Scheffe test:  F = " tval^2/(5-1) "  P = " Ftail(5-1,e(df_r),tval^2/(5-1))
lincom 2*1.Diet+2*2.Diet+2*3.Diet-3*4.Diet-3*5.Diet
scalar tval=r(estimate)/r(se)
di "SS: " tval^2*mse "  in %: " tval^2*mse/sstrt*100
di "Scheffe test:  F = " tval^2/(5-1) "  P = " Ftail(5-1,e(df_r),tval^2/(5-1))
lincom 4.Diet-5.Diet
scalar tval=r(estimate)/r(se)
di "SS: " tval^2*mse "  in %: " tval^2*mse/sstrt*100
di "Scheffe test:  F = " tval^2/(5-1) "  P = " Ftail(5-1,e(df_r),tval^2/(5-1))

* using matrix specification
test, showorder
matrix input mycon=(1,-1,0,0,0,0\1,1,-2,0,0,0\2,2,2,-3,-3,0\0,0,0,1,-1,0)
matrix checkorth=mycon*mycon'
matrix list checkorth /* diagonal matrix => orthogonal contrasts */
test, test(mycon) mtest      /* no correction */
test, test(mycon) mtest(bon) /* Bonferroni correction */
test, test(mycon) mtest(holm) /* Holm correction */