\section*{An example in Splus and SAS}

\subsection*{Multiple Comparison Methods}

Recall that if we want a (family) of tests at the $\alpha$ level, then
we form confidence intervals of the form
$$
         \hat{D} \pm Z \times \hat{s}(\hat{D}) $$

All that differs in the various tests is the multiplier, $Z$.  The
following table summarized the various multipliers.

\begin{tabular}{lll}
Test & Critical Statistic & Remarks\\
Single Comparison & $t[1 - \alpha/2; \nu]$ & $\nu$ df for Error \\
Tukey all pairs & $T = {1\over{\sqrt{2}} } q[1 - \alpha; r; \nu]$ & $r$ groups\\
Scheffe's method & $ S^2 = (r-1)F[1 - \alpha; r-1; \nu]$ \\
Bonferroni method &$ B = t[1 - \alpha/{2g}; \nu]$ & $g$ comparisons \\
\end{tabular}


\subsection*{A Randomized Block Design}

The following data are found on page ~28, of {\em Planning of
  Experiments\/} by D.~R~.Cox, Wiley, 1958, and were previously
described in Cochran and Cox's (a different Cox) book.

The treatments are five elves of potash (26, 54, 7, 108, 144 lb 
$K_2 O$ per acre applied to a cotton crop.  On observation analyzed
was a measure of single-fiber strength in arbitrary units, obtained as
an average of a number of test on the cotton from each plot.  There
were three blocks each containing five plots.  We have no details of
the arrangement of plots within blocks.   A summary of the data (which
ignores the order in which the data were collected) is

\begin{tabular}{llllll}
       & $T_1$ & $T_2$ & $T_3$ & $T_4$ & $T_5$ \\
Block 1& 7.62 & 8.14 & 7.76 & 7.17 & 7.46 \\
Block 2& 8.00 & 8.15 & 7.73 & 7.57 & 7.68 \\
Block 3& 7.93 & 7.87 & 7.74 & 7.80 & 7.21
\end{tabular}

We are interested in knowing if the treatments differ.  We are
prepared to live with block to block variations, but don't really care
about them.  

I entered the data into a computer file as follows,
\begin{verbatim}
Yield    Tmt    Block
7.62     T1     B1
8.14     T2     B1
7.76     T3     B1
7.17     T4     B1
7.46     T5     B1
8.00     T1     B2
8.15     T2     B2
7.73     T3     B2
7.57     T4     B2
7.68     T5     B2
7.93     T1     B3
7.87     T2     B3
7.74     T3     B3
7.80     T4     B3
7.21     T5     B3

\end{verbatim}

and then read the data into Splus
\begin{verbatim}
> coxdat <- read.table("coxdata.dat", head=T)
> coxdat
   Yield Tmt Block 
 1  7.62  T1    B1
 2  8.14  T2    B1
......
\end{verbatim}

and now, to a simple analysis (assume I've presented all the plots---I
did look at them and Treatment 3 is a little less variable than one
might expect.).

\psfig{figure=coxdat.eps,angle=-90,width=4in}

\begin{verbatim}
> mymod <- aov(Yield ~ Block + Tmt)      

> anova(mymod)
Analysis of Variance Table

Response: Yield

Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
Block      2   0.09712 0.048560 1.111594 0.3749855
Tmt        4   0.73244 0.183110 4.191599 0.0403682
Residuals  8   0.34948 0.043685                   
> 
\end{verbatim}

It is quite clear that there is a reasonable treatment effect, and an
ignorable Block effect (more on that in a moment), but which treatment
is the best?   In the absence of any particular interesting contrasts
we might want to test all pairwise contrasts to see which treatments
are better than the others.  The plot shows us that the highest
treatment level doesn't seem to generate the highest yield.


\subsubsection*{The multiple comparison procedures}

With 5 treatments (and hence 10 possible pairwise comparisons) we
can't do simple t-tests and the simple Bonferroni adjusts would be too
severe, so lets try both the Tukey test and the Scheffe test.  To get
started we need to calculate the treatment means.  (We could do this
with a fancy use of contrasts, but brute force is just as good.)
This is yet another one of those ``cool Splus tricks''---you should
look up the {\tt tapply} and {\tt apply} functions.
\begin{verbatim}
> mymeans<-tapply(Yield, Tmt, mean)
> mymeans
   T1       T2       T3       T4   T5 
 7.85 8.053333 7.743333 7.513333 7.45
\end{verbatim}
The rest of the calculation is essentially by hand.  Lets compare
level T2 (the biggest) with level T5 (the smallest).  We expect this
to be significant.
\begin{verbatim}
> mymeans[2] - mymeans[5]
        T2 
 0.6033333
\end{verbatim}
A 90\% confidence interval for this difference, using Tukey's
studentized range intervals, is 
$$
         \hat{D} \pm Z \times \hat{s}(\hat{D}) $$
where $$D = .60 $$, 
$$ \hat{s}(\hat{D}) = \sqrt{MSE \times ({{1}\over{3}} + {{1}\over{3}})}
= \sqrt{0.043685 * 2/3} = 0.17 $$ and
$$Z = {1\over{\sqrt{2}} } q[0.90; 5; 8] = 2.95$$

Hence the interval is
$$ 0.6 \pm 2.95 \times 0.17 = 0.6 \pm 0.50 $$.

From the calculation it is clear that no pairwise difference that is
less than $0.5$ will be significant, which rules out everything else
expect T2 against T4 or T5.

If we had done the test just as a pre-known contrast, using the
t-test, the confidence 90\% confidence region would have been
$$ 0.6 \pm 1.86 \times 0.17 = 0.6 \pm 0.31 $$.  Thus, in addition,
the comparison of T2 with T3 and T1 with T4 or T5 would all have been
significant. (But we would have overstated things).

How about as an F-test?  For Scheffe's method we need, 
$$ S^2 = (r-1)F[1 - \alpha; r-1; \nu]$$.  If we are only going to do the
pairwise differences, then $\hat{s}(\hat{D}$ is the same as before.
Now
\begin{verbatim}
> qf(.90, 4, 8)
[1] 2.806426
\end{verbatim}
so, $ S = \sqrt{(r-1)F[1 - \alpha; r-1; \nu]}$ is
\begin{verbatim}
> sqrt(4*qf(.90, 4, 8))
[1] 3.350478
\end{verbatim}
Now the multiplier is 3.35, so the confidence interval is
$$ 0.6 \pm 3.35 \times 0.17 = 0.6 \pm 0.57 $$

So with this method we only get a significant difference between T2
and T5.  

\subsubsection*{An analysis in SAS}
Here is a SAS program to do the same thing we did in Splus.
{\small
\begin{verbatim}
options linesize=78;
data cox;
        input yield tmt $ block $ ;
        cards;
7.62     T1     B1
8.14     T2     B1
7.76     T3     B1
7.17     T4     B1
7.46     T5     B1
8.00     T1     B2
8.15     T2     B2
7.73     T3     B2
7.57     T4     B2
7.68     T5     B2
7.93     T1     B3
7.87     T2     B3
7.74     T3     B3
7.80     T4     B3
7.21     T5     B3
;
proc anova;
        class tmt block;
        model yield = tmt block;
        means tmt/bon tukey scheffe alpha=0.10;

\end{verbatim}}
(I will explain what the various lines mean in class---you should
annotate this).    The log file from that SAS run was
{\small
\begin{verbatim}
......
NOTE: AUTOEXEC processing beginning; file is /usr/local/lib/sas/autoexec.sas.

NOTE: SAS initialization used:
      real time           0.19 seconds
      cpu time            0.18 seconds


NOTE: AUTOEXEC processing completed.

1          options linesize=78;
2          data cox;
3               input yield tmt $ block $ ;
4               cards;

NOTE: The data set WORK.COX has 15 observations and 3 variables.
NOTE: DATA statement used:
      real time           0.28 seconds
      cpu time            0.06 seconds


20         ;
21         proc anova;
22              class tmt block;
23              model yield = tmt block;
24              means tmt/bon tukey scheffe alpha=0.10;
NOTE: The PROCEDURE ANOVA printed pages 1-5.
NOTE: PROCEDURE ANOVA used:
      real time           0.54 seconds
      cpu time            0.09 seconds


NOTE: The SAS System used:
      real time           1.06 seconds
      cpu time            0.37 seconds
.....
\end{verbatim}}
and the susbstantive output follows.  Notice that the significant
differences we found are confirmed with SAS.   
{\small
\begin{verbatim}
                                The SAS System                               1
                                             09:19 Wednesday, January 31, 1996

                        Analysis of Variance Procedure
                           Class Level Information

                      Class    Levels    Values

                      TMT           5    T1 T2 T3 T4 T5

                      BLOCK         3    B1 B2 B3


                   Number of observations in data set = 15


                                The SAS System                               2
                                             09:19 Wednesday, January 31, 1996

                        Analysis of Variance Procedure

Dependent Variable: YIELD   
                                    Sum of           Mean
Source                  DF         Squares         Square   F Value     Pr > F

Model                    6      0.82956000     0.13826000      3.16     0.0677

Error                    8      0.34948000     0.04368500

Corrected Total         14      1.17904000

                  R-Square            C.V.       Root MSE           YIELD Mean

                  0.703589        2.706677       0.209010             7.722000


Source                  DF        Anova SS    Mean Square   F Value     Pr > F

TMT                      4      0.73244000     0.18311000      4.19     0.0404
BLOCK                    2      0.09712000     0.04856000      1.11     0.3750



                                The SAS System                               3
                                             09:19 Wednesday, January 31, 1996

                        Analysis of Variance Procedure

           Tukey's Studentized Range (HSD) Test for variable: YIELD

         NOTE: This test controls the type I experimentwise error rate, but 
               generally has a higher type II error rate than REGWQ.

                       Alpha= 0.1  df= 8  MSE= 0.043685
                  Critical Value of Studentized Range= 4.169
                    Minimum Significant Difference= 0.503

         Means with the same letter are not significantly different.

                 Tukey Grouping              Mean      N  TMT
                                      
                              A            8.0533      3  T2
                              A       
                      B       A            7.8500      3  T1
                      B       A       
                      B       A            7.7433      3  T3
                      B               
                      B                    7.5133      3  T4
                      B               
                      B                    7.4500      3  T5



                                The SAS System                               4
                                             09:19 Wednesday, January 31, 1996

                        Analysis of Variance Procedure

                Bonferroni (Dunn) T tests for variable: YIELD

         NOTE: This test controls the type I experimentwise error rate, but 
               generally has a higher type II error rate than REGWQ.

                       Alpha= 0.1  df= 8  MSE= 0.043685
                          Critical Value of T= 3.36
                    Minimum Significant Difference= 0.5726

         Means with the same letter are not significantly different.

                  Bon Grouping              Mean      N  TMT
                                     
                             A            8.0533      3  T2
                             A       
                     B       A            7.8500      3  T1
                     B       A       
                     B       A            7.7433      3  T3
                     B       A       
                     B       A            7.5133      3  T4
                     B               
                     B                    7.4500      3  T5



                                The SAS System                               5
                                             09:19 Wednesday, January 31, 1996

                        Analysis of Variance Procedure

                      Scheffe's test for variable: YIELD

         NOTE: This test controls the type I experimentwise error rate but 
               generally has a higher type II error rate than REGWF for all 
               pairwise comparisons

                       Alpha= 0.1  df= 8  MSE= 0.043685
                         Critical Value of F= 2.80643
                    Minimum Significant Difference= 0.5718

         Means with the same letter are not significantly different.

                Scheffe Grouping              Mean      N  TMT
                                       
                               A            8.0533      3  T2
                               A       
                       B       A            7.8500      3  T1
                       B       A       
                       B       A            7.7433      3  T3
                       B       A       
                       B       A            7.5133      3  T4
                       B               
                       B                    7.4500      3  T5

\end{verbatim}}
\section*{Randomized Block Experiments}
The following data are drawn from Box, Hunter, and Hunter, {\em
  Statistics for Experimenters}.  
4 different methods (treatments) are used in the manufacture of
penicillin.  5 different blends of the raw material are available for
the process.  (We want to use different blends to ensure that that
whatever method we use is good over a wide variety of situations).  
The data are \\
\begin{tabular}[c]{l r r r r}
block & \multicolumn{4}{c}{treatment} \\
      & A & B & C & D \\ \hline
blend 1 & $89^1$ & $88^3$ & $97^2$ & $94^4$ \\
blend 2 & $84^4$ & $77^2$ & $92^3$ & $79^1$ \\
blend 3 & $81^2$ & $87^1$ & $87^4$ & $85^3$ \\
blend 4 & $87^1$ & $92^3$ & $89^2$ & $84^2$ \\
blend 5 & $79^3$ & $81^4$ & $80^1$ & $88^2$
\end{tabular}

The superscripts represent the order in which the experiments were
performed.  We are interested in the treatment effects---we don't
really care about the block (blend) effects.

We need to consider both models for the data, and ways of looking at
the data.  First, lets think about how to enter the data.  I
re-expressed the data in the following format
\begin{verbatim}
Yld Blend Tmt Order
89 blend1  A    1
88 blend1  B    3
97 blend1  C    2
94 blend1  D    4
.....
\end{verbatim}

and then tried to do some simple plots within S,
\begin{verbatim}
> corn <- read.table("bhh-corn.dat", header=T)
> corn
   Yld  Blend Tmt Order 
 1  89 blend1   A     1
 2  88 blend1   B     3
 3  97 blend1   C     2
 4  94 blend1   D     4
 . . .
\end{verbatim}
\newpage
\begin{verbatim}
> plot.factor(corn)
> interaction.plot(Blend, Tmt, Yld)
> interaction.plot(Tmt, Blend, Yld)   
\end{verbatim}
These plots appear below.

\psfig{figure=rblock1.eps,angle=-90,width=5in}

\psfig{figure=rblock2.eps,angle=-90,width=5in}


What is an appropriate model for these data?

How about
$$
 y_{ij} = \mu + \beta_i + \tau_j + \epsilon_{ij} $$
where i identifies the blend and j identifies the tmt.

Our next step is to set up an ANOVA.   Here is a hint of what we might
do
\begin{verbatim}
> attach(corn, pos=1)
> ls()
[1] "Blend" "Order" "Tmt"   "Yld"  
> mymod <- aov(Yld ~ Blend + Tmt)
> anova(mymod)
Analysis of Variance Table

Response: Yld

Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
Blend      4       264 66.00000 3.504425 0.0407462
Tmt        3        70 23.33333 1.238938 0.3386581
Residuals 12       226 18.83333                   
\end{verbatim}
and
\begin{verbatim}
> anova( aov(Yld ~ Blend + Tmt + Order) )
Analysis of Variance Table

Response: Yld

Terms added sequentially (first to last)
          Df Sum of Sq  Mean Sq  F Value     Pr(F) 
Blend      4  264.0000 66.00000 3.531415 0.0435472
Tmt        3   70.0000 23.33333 1.248480 0.3392314
Order      1   20.4167 20.41667 1.092420 0.3183604
Residuals 11  205.5833 18.68939                  
\end{verbatim}

\section*{Things not to forget}
Remember that we are doing linear analyses just like regression, so
many of the ``tricks'' you learnt, like transformations of the response
variable and diagnostics about influence, are appropriate here.  The
careful data analyst will take all these into account with real data.


\section*{More examples}
The following are some gedanken examples and comments about
multifactor experiments, followed by some data.
\begin{description}
\item[Example 1] A company investigated the effect of selling price
  and type of promotion campaign on sales of one if its products.
  Three selling prices (59, 60, and 64 cents) were studies, as were
  two types of promotional campaigns (radio, and newspaper).  Consider
  the prices to be factor A and the campaign to be factor B.  The
  experiment was done using each combination of A and B (that is there
  were 6 runs of the experiment).  That isn't strictly necessary, but
  does make the analysis a little easier.  Each combination was tried
  in 2 communities, for a total of 12 experimental runs.

  This is an {\bf experiment} because the experimenters exercised
  control of the assignment of factors (A and B) to treatment units.
  This is an example of a completely randomized design, with two
  factors.
\item[Example 2] A steel company studied the effects of carbon content
  and tempering temperature on the strength of steel.  Carbon content
  was investigated at a high and low level.  Tempering temperature was
  also studied at two levels, high and low.  Thus we have a $2 \times
  2$ design with 4 possible treatment combinations.

  The 4 treatments were assigned to 12 production batches in a
  randomized fashion.  This is again an experiment, and again a
  completely randomized design.

\item[Example 3]  A researcher studied the effects of family income
  (under \$10,000, 10,000--29,999, 30,000--49,999, over \$50,000) and
  stage in life family life cycle (4 stages from just married to just
  retired) on appliance purchases.  Here we have $ 4 \times 4 = 16 $
  possible treatments.

  The researcher selected 20 families with the required income and
  life-cycle characteristics for each treatment combination, yielding
  320 families.

  This is an observational study---the data were collected without
  assigning the factor levels to experimental unit.  This data would
  generally be analyzed as though it were a completely randomized
  design.  (Lots of things can go wrong.  One obvious problem is that
  the chances of being in a family go down as one gets older, \ldots)
  
\end{description}

If one only has one observation per treatment combination, then one
can extract ``main effects'' for the factors and ``interactions'', but
there is nothing left to estimate random variability (error).  With
multiple observations per cell, one can estimate main effects and
interactions an still have something left with which to estimate
error.

One does not have to do the complete replication of the experiment
(all factors at all levels).  For example, suppose we wished to
investigate an industrial process that had three factors with 6, 5,
and 4 levels.  That would yield $6 \times 5 \times 4 = 120$
combinations.  It is probably not possible to do 120 runs of the
process.  However, by carefully selecting which combinations to run,
one can still recover the main effects, while sacrificing the ability
to estimate some interactions.  

We will mainly focus on completely replicated experiments, unless
someone wants to see something else.  

\section*{Some more data}
Here is some more data.  Again, data on blood pressure.  There are two
types of treatment to consider,  Biofeedback (Yes or No) and Drugs
(three types, A, B, and C).  An experiment is performed and there are
5 replications per treatment combination.  The raw data is

\begin{tabular}{l c c c}
      &Drug A & Drug B & Drug C \\ \hline
Bio Y & 170   & 186    & 180 \\
      & 175   & 194    & 187 \\
      & 165   & 201    & 199 \\
      & 180   & 215    & 170 \\ 
      & 160   & 219    & 204 \\ \hline
Bio N & 173   & 189    & 202 \\
      & 194   & 194    & 228 \\
      & 197   & 217    & 190 \\
      & 190   & 206    & 206 \\
      & 176   & 199    & 224 \\ \hline
\end{tabular} \\
\ \\
The easiest way to lay out this data is in the following form
\begin{verbatim}
Resp    Biofeedback    Drug
170        Y            A
175        Y            A
165        Y            A
180        Y            A
160        Y            A
186        Y            B
194        Y            B
201        Y            B
215        Y            B
219        Y            B
180        Y            C
187        Y            C
199        Y            C
170        Y            C
204        Y            C
173        N            A
194        N            A
197        N            A
......
\end{verbatim}

There are many questions to ask of the data, starting with
\begin{itemize}
\item What are the main effects of the two treatments, and do they
  matter?
\item Are there any interactions between the treatments?
\end{itemize}

A good place to start is with an interaction plot of the full data
\begin{verbatim}
> blood2 <- read.table("blood2.dat", header=T)
> blood2
   Resp Biofeedback Drug 
 1  170           Y    A
 . . . .
> plot.factor(blood2)
> attach(blood2, pos=1)
> ls()
[1] "Biofeedback" "Drug"        "Resp"       
> interaction.plot(Drug, Biofeedback, Resp)
> interaction.plot(Biofeedback, Drug, Resp)      
\end{verbatim}

\psfig{figure=rblock3.eps,angle=-90,width=5in}

The interaction plots are very interesting. It certainly looks like
there might be an interaction between the two factors.  Lets try a
very simple analysis to confirm/deny this.
\begin{verbatim}
> ls()
[1] "Biofeedback" "Drug"        "Resp"       
> mymod <- aov( Resp ~ Biofeedback * Drug ) 
> anova(mymod)
Analysis of Variance Table

Response: Resp

Terms added sequentially (first to last)
                 Df Sum of Sq Mean Sq  F Value     Pr(F) 
Biofeedback       1      1080 1080.00  6.93419 0.0145635
Drug              2      3420 1710.00 10.97913 0.0004113
Biofeedback:Drug  2       780  390.00  2.50401 0.1028767
Residuals        24      3738  155.75                   
> 
\end{verbatim}

The ANOVA table suggests that both the Bio and Drug treatments are
effective (why?) and that there is borderline evidence for the
interaction(why?).

What to do now?
\begin{itemize}
\item Do some more analyses in class
\item Wonder about what component of the interaction is important?
\end{itemize}

\section*{The Plan}
\begin{itemize}
\item What does interaction mean?
\item GLIM notation for main-effects, interactions (Venables and
  Ripley, p.~153)
\item Tests for Main Effects and Interactions %Section 18.6
\item Contrasts %Section 19.2 ...
\item Unequal sample sizes in 2-way ANOVA %(Chpt 20)
\item Multifactor Studies %(Chpt 22)
\item Analysis of Covariance %(Chpt 23)
\item Some Standard Designs (Latin Square and Split Plot).
\end{itemize}
See also Venables and Ripley p.~163 et.~seq.\ for a discussion of
ANOVA models.

% Local Variables: 
% mode: latex
% TeX-master: t
% TeX-master: t
% End: 


