***********************************************************;
* A ``two-way layout'' for reading comprehension scores for
*   elementary-school children:
*
* We compare three different teaching methods, with possible
*   sex-dependent effects.
* (That is, boys may do better on some methods and girls may do
*   better on other methods.)
*
* The model is called a two-factor or two-way ANOVA with 5
*   observations per cell. That is, we have test results (reading
*   comprehension scores) for 5 children for each combination of
*   sex and reading method.
*
* Some useful general terminology for ANOVAs: Class-valued
*   covariates are called `FACTORS'. There are two factors here,
*   reading method (with three values) and sex (with two values).
*   A one-way ANOVA has one factor.
*
* The values of a factor (corresponding to `treatment groups' for
*   a one-way ANOVA) are called the `LEVELS' of that factor.
*
* In this example, the first factor (teaching method) has three
*   levels, which we will call (unimaginatively)
*   Method1, Method2, and Method3
*
* The second factor (sex) has two levels, Male and Female.
*
* Which (if any) of the two factors cause significant variation
*   in the reading scores?
*
* In principle, the simplest two-factor ANOVA model is:
*
*       Yij = mu + C(Sex_i) + D(Method_j) + error         (1)
*
* where Yij is a test score, mu is an overal mean, C_F and C_M are
*   constants for Sex effects, and D_1, D_2, D_3 are constants for
*   the three reading methods.  With 5 observations per cell, the
*   statistical model for the data is
*
*       Yijk = mu + C(Sex_i) + D(Method_j) + error(ijk)   (2)
*
*   where 1 le k le 5, since there were 5 observations for each
*   combination (C,D). Here error(ijk) are independent normally
*   distributed random variables with the same unknown variance.
*
* However, this model may not fit the data:
*
* Specifically, (1) implies that if girls do better than boys for
*   Method1, then girls will also do better than boys on Method2 and
*   Method3 as well, and by exactly the same amount.
*
* Specifically, if C_F > C_M, then (1) implies that, on the average,
*   girls will do better than boys by the amount  C_F - C_M  on each
*   of the three reading methods.
*
* If girls better than boys on Method_1 and Method_3, but do worse
*   on Method_3, then the model (1) cannot be correct. Similarly,
*   if girls do worse than boys on Methods 1 and 3 but the same on
*   Method 2, then also (1) cannot be correct.
*
* In fact, if girls do better than boys on Method_1, do the same as
*   boys on Method_2, and do worse on Method_3, then girls and boys
*   may do as well on the average. If you tried to fit the purely
*   additive model (2) to the data, then you might conclude
*   incorrectly that Reading Method had no effect, when in fact
*   there may be strong effects of both Reading Method and Sex.
*
* ANOVAs in which the effect of factors is not purely additive are
*   called models with INTERACTIONS.  The average overall effect of
*   each factor individually, in effect averaging over non-additive
*   interactions with other factors, are called the MAIN EFFECTS.
*
* Often the INTERACTION is more interesting than the effect of the
*   factors individually.  In this case, it might be known in
*   advance that the three reading methods do about the same overall
*   and that males and females also do about the same overall on the
*   three reading methods. Then the main interest would be in testing
*   for variability in sex differences among the tests --- perhaps
*   one method should be used for girls and the other for boys? This
*   is exactly equivalent to testing for an interaction between the
*   two factors.
*
* The most general model for two factors is a one-way ANOVA for all 
*   six possible combinations of the two factor levels, which is
*
*    Yij = G(Sex_i,Method_j) + error(ij)                (3)
*
* This has a separate constant for each combination of sex and
*   reading method, and covers main effects as well as interactions.
*   However, the constants in (3) do not give any intuition about
*   the relative importance of main effects and interactions.
*
* One can show using linear algebra that the model (3) (which has
*   6 free parameters constants for Sex=M,F and Method=1,2,3) is
*   exactly equivalent to the model
*
*    Yij = mu + C(Sex_i) + D(Method_j)
*
*              + H(Sex_i,Method_j) + error(ij)           (4)
*
* subject to the four constraints
*
*   Sum C(Sex_i)=0,         (Main Effect for Sex)
*
*   Sum D(Method_j)=0,      (Main Effect for Reading Method)
*
*   Sum_i H(Sex_i,Method_j) = 0  for all j  (Interactions)
*
*   Sum_j H(Sex_i,Method_j) = 0  for all i
*
* Under these constraints, the constants C_i, D_j, H_ij in (4) are
*   exactly determined by the constants G_ij in (3) and vice versa.
*
* In general, if factor C has r levels (here r=2) and factor D has
*   s levels (here s=3), then the constraints on the Main Effects
*   say that there are r-1 free constants for the C's and s-1 free
*   constants for the D's. The constraints on the interaction terms
*   H(Sex_i,Method_j) means that the interaction terms have 
*   (r-1)(s-1) free constants. (Think of an r by s matrix with one
*   row and one column determined by the constraints.)  This means
*   that the total number of free constants on the right-hand side
*   of (4) is
*
*     1 + (r-1) + (s-1) + (r-1)(s-1) = rs
*
*   This is the same number of free constants as in (3). This in
*   fact could be the first step of a proof that the models (2)
*   and (4) are equivalent.
*
* Statistically, estimates of the terms  C(Sex_i),  D(Method_j), and
*   H(Sex_i,Method_j)  in (4) lead to independent sum-of-squares
*   contributions that are independent of the error sum of squares.
*   This leads to ANOVA tables with three separate F-tests : for the
*   Sex main effect, the Reading Method main effect, and the Sex by
*   Reading Method Interaction effect. (See more discussion of these
*   F-tests below.)
*
* One can check graphically if there is a significant interaction
*   by plotting estimates of the constants G(Sex,Method) in (3) and
*   seeing if the estimates appear to be additive in i and j.  This
*   is called an ``interaction plot'' and is a very useful device
*   for seeing qualitatively if there is a (statistical) interaction
*   (meaning a lack of additivity) between the two factors.
*
* Some more useful terminology for ANOVAs:
*
*   An EFFECT is a term in a model like (1) or (2) or (3) or (4).
*   Effects can be either MAIN EFFECTS or else INTERACTION EFFECTS .
*   Except in simple models, there are more effects than factors.
*
*   A CELL is a particular combination of factor levels for all
*   factors, for example Method_1 with Sex=Female. Thus model (3)
*   can be described as a one-way ANOVA on the CELLS of a two-factor
*   ANOVA.
*
*   A FULL FACTORIAL MODEL is a model like (4) that has enough
*   effects to be equivalent to the one-way ANOVA on cells. For
*   two factors, a full factorial model has two main effects and
*   one interaction.
*
*   A BALANCED ANOVA is a model with the same number of observations
*   for each cell. This usually guarantees that the effects in the
*   model are `orthogonal'. The latter means that the Type I and
*   Type III tables in SAS output will be identical.
*
*   SAS can handle unbalanced ANOVAs as easily as balanced ANOVAs,
*   but the nonorthogonality means that the statistical significance
*   of the effects in the model depend on the order in which they
*   are considered by the model. (That is, Type I P-values for main
*   effects and interactions are likely to be different than
*   Type III P-values.) For a complicated ANOVA, this may mean that
*   the output will be very difficult to uderstand for a
*   statistician, let alone a lay person.
*
* DEGREES OF FREEDOM: The full factorial model for two factors is
*
*    Yij = mu + C(Sex_i) + D(Method_j)
*
*            + H(Sex_i,Method_j) + error(ij)               (5)
*
*   Suppose that the first factor (`Sex') has r levels (here r=2)
*   and the second factor (`Method') has s factors (here s=3).
*   The constraints
*
*      Sum C(Sex_i)=0,  Sum D(Method_j)=0,
*
*   suggest that the main effects SEX and METHOD have r-1 and s-1
*   degrees of freedom, respectively, exactly as for a one-way ANOVA.
*   Specifically, there are r-1 free constants among the C(Sex_i)
*   and s-1 free constants among the D(Method_j). The two
*   constraints for the interaction
*
*      Sum_i H(Sex_i,Method_j) = 0  for all j
*      Sum_j H(Sex_i,Method_j) = 0  for all i
*
*   suggest that the interaction effect SEX*METHOD has (r-1)(s-1)
*   degrees of freedom. Thus the total number of degrees of freedom
*   in (5) should be
*
*   1(mu) + (r-1)(SEX) + (s-1)(METHOD) + (r-1)(s-1)(SEX*METHOD) = rs
*
*   This is consistent with the one-way ANOVA (3) on cells.
*
*   These suggestions are in fact correct, and identify the NUMERATOR
*   degrees of freedom for F-tests of the three effects. The
*   DENOMINATOR degrees of freedom is the same as that of the Error
*   Sum of Squares (SSE). For a full factorial model, SSE is exactly
*   the same as the SSE for a one-way ANOVA on the cells, which 
*   means  krs-rs = (k-1)rs  degrees of freedom in this case.
*
*
* The example data below has observations (reading comprehension
*   scores) for 3 reading method levels and 2 sexes, with 5
*   observations for each cell.  We will use SAS to display
*   an interaction plot as well as do an ANOVA analysis.
*
* The data is from the text ``Applied Statistics and the SAS
*   programming language'', R. Cody and J. Smith, 2004, 5th ed,
*   Prentice Hall, New Jersey, page 210.
* 
* SAS NOTE: Note the use of `substr()' to identify what input words
*   are `Method' labels:  `Substr()' has the syntax
*
*        `substr(string,start,length)'.
* 
*   Thus, if zz=`Method1', then substr(zz,1,6)='Method' and
*   substr(zz,3,4)='thod'.
***********************************************************;

title 'TWO-WAY ANOVA WITH POSSIBLE INTERACTIONS: Reading and Sex - YOURNAME';
options ls=75 ps=60 pageno=1 nocenter;

data rdtest;
     retain Sex Method;
     input xx$ @@;
     if xx='Male' or xx='Female' then Sex=xx;
     else if substr(xx,1,4)='Meth' then  Method=xx;
     else do; yy=input(xx,12.0); output; end;
     drop xx;

datalines;
     Male     Method1  700  850  820  640  920
              Method2  480  460  500  570  580
              Method3  500  550  480  600  610
     Female   Method1  900  880  899  780  899
              Method2  590  540  560  570  555
              Method3  520  660  525  610  645
            ;

proc print;
     title2 'THE DATA AS SAS SEES IT';
     run;

proc glm;
     title2 'PROC GLM FOR MAIN EFFECTS only';
     title3 'This COMBINES INTERACTIONS with the error terms.';
     title3 'This is a BAD IDEA if there are significant interactions';
     classes Method Sex;
     model yy=Method Sex;
     run;


**********************************************;
* INTERACTION PLOT:
* How can we see if there is an interaction graphically?
*
* An interaction plot for two factors involves three values:
*     (i) The estimated mean value for a cell or group of cells,
*    (ii) The level (value) of the first factor
*   (iii) The level (value) of the second factor.
*
* A graph has only two axes, but we can let
*    Y-axis:  Mean value of cell
*    X-axis:  Levels of first factor
*    Plotting symbol:  Level of second factor
*
* If the Y-axis DIFFERENCES of the second factor in this group are
*   NOT THE SAME for all values of the first factor (that is, if the
*   vertical pattern varies with X), and if these differences are
*   statistically significant, THEN THERE IS AN INTERACTION.
*   In some cases, the interaction is strong enough so that lines
*   connecting points with the same plotting symbol cross. This is
*   the same saying that the Y-axis differences are of different sign.
*
* By convention, the factor with the LARGER NUMBER OF LEVELS is put
*   on the X axis. This is because reading from left to right is
*   far easier than comparing vertical distance by eye.
*
* By definition, `interaction plots' are based on CELL MEANS. If all
*   observations were included in the interaction plot, then there
*   could be too much clutter to see the interaction.
*
* This means that we must first create a NEW DATASET that has
*   CELL MEANS instead of individual values. Such a new dataset can
*   be created using `proc means' with an `output' statement (see
*   below). We can then use `proc plot' on the new dataset of cell
*   means. Since the new dataset is now the SAS default dataset, we
*   will have to refer to the original dataset explicitly if we want
*   to use it again.
*
* SAS notes: (1) The argument `nway' in `proc means' means that SAS
*   should write mean values for cell means only. Otherwise, SAS
*   will also write entries to the cell-mean dataset for row,
*   column, and overall means.
*     In the proc means output dataset, only the cell means have
*   nonmissing values for both variables. Thus, the only difference
*   that leaving out the argument `nway' makes is to avoid error
*   messages from SAS about missing values in the interaction plot.
*     (2) The command `noprint' means no printed output: Here we use
*   `proc means' just to write a dataset of cell means, and don't
*   want SAS to write extra clutter in the form of tables of mean
*   values to the output file.
*
**********************************************;

proc means nway noprint;
     class Method Sex;
     var yy;
     * The following writes a new SAS dataset called `Iplotmeans';
     * The name of the column with cell means will be `CellMean' ;
     * If you say `mean=', the cell-mean column name defaults to `yy';
     output out=Iplotmeans  mean=CellMean;
     run;


**********************************************;
* The cell-means dataset `Iplotmeans' is now the SAS default, but we
*   refer to it explicitly for clarity.
* The following displays Iplotmeans, but not the interaction plot.
* Note that there are only 6 records in `Iplotmeans' as opposed to
*   30 records in `rdtest'.
**********************************************;

proc print data=Iplotmeans;
     title2 'INTERACTION PLOT for Method versus Sex';
     title3 '(Full factorial two-way model)';
     title4 'CELL MEANS FOR INTERACTION PLOT ARE';
     run;

**********************************************;
* The INTERACTION PLOT is below.
* If M-F differences on the Y-axis are not the same for all values
*   on the X-axis, then there is an interaction.  Can you see an
*   interaction in this case?
* The page size (ps) command restricts the vertical range of the
*   plot in lineprinter plots to make it easier to read.
**********************************************;

options ps=35;

proc plot data=Iplotmeans;
     title4 'INTERACTION PLOT FOR IPLOTMEANS';
     plot CellMean*Method = Sex;
     run;

options ps=60;


**********************************************;
* In contrast, let's see what `proc means' does without the `nway'
*   argument. Note that an additional 6 records are written for
*   row, column, and overall means, each with missing values for
*   row or column or both.
**********************************************;

proc means noprint data=rdtest;
     class Method Sex;
     var yy;
     output out=Fullmeans  mean=CellMean;
     run;


proc print data=Fullmeans;
     title4 'IN CONTRAST, THE UNRESTRICTED PROC MEAN OUTPUT IS';
     run;

**********************************************;
* Now do a two-way ANOVA with the INTERACTION term.
*
* The output will now have an additional line in the Type I and
*   Type II ANOVA tables for the interaction.
*
* Is the interaction effect significant?
* Are the Type I and Type III tables the same?
*
* Don't forget to refer back to the original dataset, since
*   the default SAS dataset is now Fullmeans.
*
**********************************************;

proc glm  data=rdtest;
     title2 'TWO-WAY ANOVA with interaction:';
     title3 '(FULL FACTORIAL two-way model)';
     classes Method Sex;
     model yy=Method Sex Method*Sex;
     * Let's also test for Main-effect differences:';
     means Method Sex / duncan;
     run;


**********************************************;
* Alternatively, we could do a one-way ANOVA on all possible
*   (Method,Sex) combinations.  This is the model (3) in the
*   introduction above.
*
* How does the model test compare with the model test for the
*   full factorial model?
* How significantly do the 6 cell means vary?
*
* The first step is to create a new class variable that
*   contains all (Method,Sex) combinations.
*
* In SAS, the command `cellname = method||sex'' does this by 
*   creating 6 distinct values values for `cellname' by
*   concatenating the `method' level name with the `sex' level
*   names.
*
* Reopen `rdtest' to add a new column:
**********************************************;

data rdtest;
     set rdtest;
     RdCell = Method || Sex;
     run;

proc glm;
     title2 'ONE-WAY ANOVA FOR JOINT (METHOD,SEX) CELLS';
     title3 '(Note that Model P-value and F-statistic are';
     title4 'identical with the full factorial model)';
     class RdCell;
     model yy=RdCell;
     means RdCell / duncan;
     run;