*************************************************;
* Multiple comparisons:
*
* Suppose that a one-way ANOVA asserts that treatment means do
*   vary significantly.
*
* A natural question is, which PAIRS of treatments are different?
*   Can you conclude that ANY pair of treatments have statistically
*   different means? If so, which pairs?
*
* If the ANOVA is significant, is the significance of the ANOVA due
*   to one treatment being different from all the others, with no
*   significant differences among the others, or are all pairs of
*   treatment means significantly different?
*
* Note that this question is clouded by `multiple comparison'
*   considerations: 
*
* With K=5 treatments, there are 10 (=K*(K-1)/2) different treatment
*   pairs. Even if all of the true treatment means are the same,
*   some of these 10 treatment pairs may be significantly different
*   purely by chance. For example, if you toss a coin ten times that
*   has Prob(Heads)=0.05, then you expect to get one or more heads
*   about half the time, not one time in 20.
*
* With K=7 treatment groups, there are 7(6)/2=21 pairwise comparisons,
*   and you might expect at least one pair to be significantly
*   different at level P=0.05, again purely by chance.
*
* SAS has several different procedures for inferring which pairs of
*   treatments have different means that are sensitive the problem
*   of multiple comparisons. These are discussed below.
*
* In the first data step, the command `substr(xx,1,3)' for the text
*   variable `xx' isolates the first three letters of `xx'. Thus
*   substr(xx,1,3)='Trt' if and only if xx begins with `Trt'. In
*   general, `substr(xx,a,b)' returns the substring of `xx' of
*   length b or less that begins with the a-th character of `xx'.
*   Thus substr('abcde',3,2)='cd' and substr('defgh',4,1)='g'.
*
*************************************************;

title 'MULTIPLE COMPARISONS AMONG 5 TREATMENTS - YOURNAME';
options ls=75 ps=60 pageno=1 nocenter;

data multcomp;
     retain treat;
     input xx$ @@;
     if substr(xx,1,3)='Trt' then treat=xx;
     else do; yy = input(xx,12.0);  output;  end;
datalines;
  Trt1  183  191  197  198  211  219  220  225  246
  Trt2  171  182  186  187  197  199  207  237  237
  Trt3  148  153  157  175  178  193  193  204  208
  Trt4  150  160  166  168  179  183  184  189  203
  Trt5  118  150  159  161  165  168  179  189  193
      ;

options ps=40;     * Shorten page height for a nicer chart;

proc chart;
     title2 "VISUALIZE THE DATA BEFORE PROCEEDING:";
     var yy / subgroup=treat;
     run;

options ps=60;     * Return page height to previous value;

proc means;
     title2 'LIST TREATMENT MEANS AND STANDARD DEVIATIONS:';
     class treat;  var yy;
     run;



title2 'WHICH TREATMENT DIFFERENCES ARE SIGNIFICANTLY DIFFERENT?';
title3 'MULTIPLE COMPARISON TESTS:';

*************************************************;
* Various methods of determining which pairs of treatments
*   are significantly different:
*
* lsd - All pairwise t-tests using the one-way ANOVA MSE error
*         estimate instead of pairwise error estimates
*     - No multiple comparison correction
*     - These are `naive' pairwise comparisons, but are slightly better
*      than the usual (equal-variance) t-tests since they estimate
*      the error variance using all treatment data instead of just
*      data from the two treatments.
*
* tukey - `Tukey's studentized range statistic'
*     - Exact MC-correct P-value for the biggest mean difference,
*         conditional on that it was the biggest difference.
*     - Conservative for all other pairwise differences
*     - The Tukey test for the maximum pairwise difference is a
*         valid test for H_0:All treatment means are the same,
*         if you get bored with the one-way ANOVA test.
*
* duncan, snk - Duncan, (Studentized) Newman-Keuls
*     - Same as Tukey for biggest mean difference
*     - Corrections if not the biggest observed difference
*     - Duncan is the most commonly used, with Snk a close second
*
* scheffe - exact MC-correct P-value for the biggest CONTRAST
*        (see below for a discussion of CONTRASTS)
*     - Valid but usually VERY conservative for pairwise differences
*
* Bonferroni - this is LSD for alpha = alpha/(no.of pairs)
*     - Widely-used general multiple-comparison procedure that is
*        valid for all multiple-comparison situations
*     - Easy to implement
*     - Exact if Pvalues are small and tests are independent
*     - More conservative than Tukey or Duncan
*
* All these procedures assume normally-distributed errors with
*   error variances the same across treatment groups.
*
* Tukey,duncan,snk,scheffe assume equal treatment-group sizes,
*   and use approximations for non-equal group sizes.
* Lsd,Bonferroni work fine as is for unequal treatment sizes.
*
* In the `proc glm' procedure below, 
*
*   (i) `means treat / .....'  tells SAS to carry out listed
*     multiple-comparison procedures (like duncan, tukey, etc.)
*
*   (ii) `output out=newdat r=resid' creates a new dataset `newdat'
*     and writes the residuals to `newdat' along with the current
*     data. The residuals will be used to test the normality of the
*     errors.
*
* If the option `out=newdat' was left out, the residuals would be
*   written to the current dataset (`multcomp'). However, the extra 
*   column is `fragile' and can be lost in some situations.
*************************************************;


*************************************************;
* TESTING CONTRASTS AMONG THE TREATMENT MEANS: 
*
* A CONTRAST is defined as a linear combination of treatment means
*   with overall mean zero, as in
*       C(mu) = mu_1 - (mu_2 + mu_3)/2   or
*       C(mu) = 3mu_1 + 8mu_2 - 5mu_4 - 6mu_5
*
* A CONTRAST TEST is a test of the hypothesis
*
*        H_0:  C(mu)=0
*
*   For example, a test of E(X_1) = (1/2)(E(X_2)+E(X_3)).
*
* For any contrast C(mu), one can derive a simple t-test for H_0:
*   C(mu)=0 using the MSE estimator for sigma^2 from the ANOVA table.
*
* The last three lines of the `proc glm' block below tell SAS to
*   carry out contrast tests for three different contrasts. The
*   contrast itself is listed after the `treat' statement. Thus
*
*        treat -2 1 1 0 0
*
*   tells SAS to test  C(mu) = -2mu_1 + mu_2 + mu_3.  Here
*   H_0:C(mu)=0  is the same as   E(X_1) = (1/2)(E(X_2)+E(X_3)).
*
* The text just after each `contrast' command gives a row heading
*   for the output to help identify the results of that particular
*   contrast test in the output. SAS does not try to parse this
*   row heading.
*************************************************;

proc glm;
     class treat;
     model yy = treat;
     means treat / lsd tukey duncan snk scheffe;
     output out=newdat r=resid;
     *************************************************;
     * The Bonferroni correction for alpha=0.05 for k=5 and
     *   k(k-1)/2 = 10 pairs is alpha = 0.05/10 = 0.005
     *************************************************;
     means treat / lsd alpha=0.005;
     *************************************************;
     * Tests for three contrasts (see above for description)
     * The text is a title for the contrast and the
     *   numbers are the coefficients.
     *************************************************;
     contrast 'C1 with av of C3 C4'  treat -2 0 1 1 0;
     contrast 'Av of C3 C4 with C5'  treat 0 0 1 1 -2;
     contrast 'C2 with av C3 C4 C5'  treat 0 3 -1 -1 -1;
     run;


*************************************************;
* Check variability of residuals by treatment group, preferably
*   with a plot with a plotting symbol for treatment group.
*
* A problem here is that we can't use the variable `treat' as the
*   plotting symbol: The `plot' command uses the first letter
*   from a plotting variable, which would give the same symbol `T'
*   for all 5 treatments.
*
* One way to handle this would be to include a plotting variable
*   in the original `multcomp' data step. This would be copied to
*   the current dataset (`newdat') in the `output' statement in
*   `proc glm'.
*
* Another way is to reopen the current dataset and add a new
*   plotting-symbol variable just before it is used.
*
* The following code creates a new data set without using a
*   `datalines' block. The syntax is the same as in a data step
*   with a `datalines' block, with the records being read from
*   the previous dataset instead of a `datalines' block. SAS allows
*   you to use the same name for both datasets: SAS presumably
*   writes a new dataset, renames the new dataset with the same
*   name as the old, and then deletes the old dataset.
*
* `substr(treat,4,1)' creates a string of length one
*   beginning at the 4th character, which is what we need.
*************************************************;

data newdat;
     set newdat;
     trtsym = substr(treat,4,1);  * 1,2,3,4,5 for the treatments;
     run;


options ps=40;     * Shorten page height for a nicer chart;

proc plot;
     title2 'ARE TREATMENT-GROUP VARIANCES THE SAME?';
     title3 'PLOT OF RESIDUALS:';
     plot resid*treat = trtsym;
     run;

options ps=60;     * Return page height to previous value;

*************************************************;
* The treatment variances seem to be the same across treatment
*    groups.
* Let's also check for normality of the residuals for all treatments
*    together:
*************************************************;

proc univariate normal plot;
     title2 'PROC UNIVARIATE: Testing normality';
     title3 'LOOK FOR W STATISTIC';
     var resid;
     run;