simultaneous Tukey HSD
check: instead of sorting, I use absolute value of pairwise differences in means. That’s irrelevant for the test, but maybe reporting actual differences would be better. CHANGED: meandiffs are with sign, studentized range uses abs
q_crit added for testing
TODO: error in variance calculation when nobs_all is scalar, missing 1/n