* do-file for lecture 11 of VHM 802, Winter 2024
version 18 /* works also with versions 14-17 */
set more off
set scheme stcolor_alt
cd "r:\"

* Manly Example 4.1
import delimited r:\sparrow.csv, clear
encode survivorship, gen(surv)
foreach var of varlist total_length-l_keel_sternum {
  ttest `var', by(surv) unequal /* unequal variances are preferable, though not the version used in Manly */
  robvar `var', by(surv) /* Levene's test W50 for equal variances */
  }
hotelling total_length-l_keel_sternum, by(surv)
manova total_length-l_keel_sternum = surv /* same result as hotelling's T2 */
mvreg
matrix errorsscp=e(E)
matrix corr=corr(errorsscp)
matrix list corr
mvtest covariances total_length-l_keel_sternum, by(surv) /* Box's M-test */

* univariate two-sample permutation test
ttest total_length, by(surv) unequal /* to see the actual analysis */
permute surv meandiff=(r(mu_1)-r(mu_2)), saving(permdist, replace) reps(1000) seed(210327) nodrop nowarn: ttest total_length, by(surv)
preserve
use permdist, clear
histogram meandiff
sum meandiff, d
restore

* Manly Example 5.2
import delimited r:\skull.csv, clear
encode period, gen(Period)
discrim lda maximumbreadth-nasalheight, group(Period) notable
estat grdistances /* Mahalanobis distance matrix, with different group ordering than in Manly's Table 5.4 */
* semi-manual calculation of Penrose distances: note 10 pairs, 4 variables
matrix define penrose=(0,0,0,0,0,0,0,0,0,0)
foreach var of varlist maximumbreadth-nasalheight {
  quietly anova `var' Period
  scalar errorvar=e(rmse)^2
  quietly pwcompare Period
  matrix define penrose=penrose+vecdiag(r(b_vs)'*r(b_vs))/errorvar/4
}
matrix list penrose

* univariate one-way ANOVA permutation test
anova basibregmaticheight Period /* to see the actual analysis */
permute Period Fperiod=(e(F_1)), saving(permdist, replace) reps(1000) seed(210327) nodrop nowarn: anova basibregmaticheight Period