* do-file for Additional multivariate exercise 12, VHM 802, Winter 2023
version 17 /* works also with versions 14-16 */
set more off
cd "r:\"

* Manly Example 4.3
import delimited r:\skull.csv, clear
encode period, gen(Period)
foreach var of varlist maximumbreadth-nasalheight {
  oneway `var' Period, tab bon  
  }
manova maximumbreadth-nasalheight=Period
matrix errorsscp=e(E)
matrix list errorsscp
matrix errorcorr=corr(errorsscp)
matrix list errorcorr
mvtest covariances maximumbreadth-nasalheight, by(Period)
mvreg maximumbreadth-nasalheight=i.Period
predict res, eq(#1) residual
margins Period
margins Period, predict(equation(#2))
pwcompare Period, equation(#2) pv mcomp(bon) /* equation 1 is the default */
lincom [#2]1.Period - [#2]2.Period /* equation 1 is the default */

* comparison with univariate analysis
anova maximumbreadth Period
regress
predict res1, res
margins Period

* time as a quantitative predictor
gen time=-4000
replace time=-3300 if Period==3
replace time=-1850 if Period==1
replace time=-200 if Period==4
replace time=150 if Period==5
generate t=time/1000 /* useful to rescale time to avoid too large values, in particular for higher orders */
manova maximumbreadth-nasalheight=c.t
mvreg
manova maximumbreadth-nasalheight=c.t##c.t
mvreg
manova maximumbreadth-nasalheight=c.t##c.t##c.t
mvreg