\section*{A Categorical Data Example}
Lets do the Red Dye number 2 example from Fienberg in Splus.  You might try
to do the same thing in SAS ( and ask questions about PROC CATMOD when
you get stuck).
Here is the problem and data.  

The following table contains one version of the results of an
experiment in which the food additive Red Dye No.~2 was fed to four
different groups of rats at varying dosages. Some of the rats died
during the experiment (i.e., between 0 and 131 weeks).  The remainder
were sacrificed at the end of 131 weeks.  All rats were examined for
tumors.   Try to understand the data, and in particular what are the
explanatory and response variables, and find some good models to
describe the data.  
\begin{center}
  \begin{tabular}{l l | r r | r r}
\hline
\multicolumn{2}{r|}{Age at Death} & \multicolumn{2}{c|}{0--131
  weeks}&\multicolumn{2}{c}{Terminal sacrifice}\\\hline
\multicolumn{2}{r|}{Dosage Group}  & Low & High & Low & High \\ \hline
Cancer & Present & 4 & 7 & 0 & 7 \\
 & Absent & 26 & 16 & 14 & 14 \\ \hline
 & Totals & 30 & 23 & 14 & 21 \\
  \end{tabular}
\end{center}

This is a $2 \times 2 \time 2$ table, with $3$ categorizing
variables. 
It seems reasonable to regard the dosage level as an explanatory
variable with age at death and presence or absence of chance as both
response variables.  One could argue that the age at death variable
might be an explanatory variable, but the number of observations in
this group is not fixed (in advance) so for most purposes it makes
sense to regard it as a response variable.   

The distinction of response or explanatory variable does {\em not\/}
change how we read in or manipulate the data, but does change the
class of models we might consider.  

In the following, I choose to enter the data the fancy way, using the
{\tt fac.design} command.  Had I been lazy, I would have just
generated the design by hand.
\begin{verbatim}
> reddye <- fac.design(c(2,2,2),c("Dose","Age","Cancer"))
> reddye$Resp <- c(4,7,0,7,26,16,14,14) 
> reddye
   Dose  Age Cancer Resp 
1 Dose1 Age1  Cncr1    4
2 Dose2 Age1  Cncr1    7
3 Dose1 Age2  Cncr1    0
4 Dose2 Age2  Cncr1    7
5 Dose1 Age1  Cncr2   26
6 Dose2 Age1  Cncr2   16
7 Dose1 Age2  Cncr2   14
8 Dose2 Age2  Cncr2   14
\end{verbatim}
Now lets try to explore some models for these data.  The command we
need to use it {\tt glm} and the appropriate family is {\tt Poisson}.

The simplest model we can fit includes just the {\tt dose} term.  If
we do not include {\tt dose} then we risk fitting models which have
fitted numbers of animals at the various dose levels different from
the oberserved numbers---and that is not a good idea.

Thus to fit the simplest model (in this case as dose is explanatory),
\begin{verbatim}
> anova(mymod <- glm( Resp ~ Dose, family=poisson, data=reddye))
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
     Df      Deviance Resid. Df Resid. Dev 
NULL                          7   49.47734
Dose  1 -4.214286e-09         6   49.47734
\end{verbatim}
and the corresponding Pearson X2 is
\begin{verbatim}
> Pearson(mymod)
Degrees of Freedom: 8 Total; 6 Residual
Residual Pearson X^2 42.72691 
\end{verbatim}

There is little to be said here.  The dosage group is not a predictor
of the raw response.  We shouldn't be surprised by this.

The next step is to fit various models to the data.  Here goes with an
unstructured look at things.  We'll see more structured approaches
shortly. 

\begin{verbatim}
> anova(mymod <- glm( Resp ~ Dose + Age, family=poisson, data=reddye))
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev 
NULL                     7   49.47734
Dose  1 0.000000         6   49.47734
 Age  1 3.707932         5   45.76941

> anova(mymod <- glm( Resp ~ Dose*Age, family=poisson, data=reddye))  
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
         Df Deviance Resid. Df Resid. Dev 
    NULL                     7   49.47734
    Dose  1 0.000000         6   49.47734
     Age  1 3.707932         5   45.76941
Dose:Age  1 2.336721         4   43.43269
> 
\end{verbatim}
That is pretty boring.  We have learnt that {\tt Age} might be
appropriate, but we have no strategy.  Instead, let us  try to
actually struggle our way through this, like we would an ANOVA.

\begin{verbatim}
> anova(mymod <- glm( Resp ~ Dose*Age + Age*Cancer + Dose*Cancer, family=poi>
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev 
       NULL                     7   49.47734
       Dose  1  0.00000         6   49.47734
        Age  1  3.70793         5   45.76941
     Cancer  1 32.82534         4   12.94407
   Dose:Age  1  2.33672         3   10.60735
 Age:Cancer  1  0.00739         2   10.59995
Dose:Cancer  1  7.60034         1    2.99961
\end{verbatim}
Even though things are {\bf NOT} orthogonal, it looks like one might
wish to remove the {\tt Age} times {\tt Cancer} interaction (leaving the main
effects)
\begin{verbatim}
> anova(mymod <- update(mymod, . ~ . - Age:Cancer))
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev 
       NULL                     7   49.47734
       Dose  1  0.00000         6   49.47734
        Age  1  3.70793         5   45.76941
     Cancer  1 32.82534         4   12.94407
   Dose:Age  1  2.33672         3   10.60735
Dose:Cancer  1  7.31735         2    3.29000
\end{verbatim}

Now things are not clear, one may---or may not---remove the {\tt Dose} times
{\tt Age} interaction.  The upper 95\% point of $\chi^2_1$ is
\begin{verbatim}
> qchisq(.95,1)
[1] 3.841459
\end{verbatim}
so maybe we should just remove it.

\begin{verbatim}
> anova(mymod <- update(mymod, . ~ . - Age:Dose))  
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev 
       NULL                     7   49.47734
       Dose  1  0.00000         6   49.47734
        Age  1  3.70793         5   45.76941
     Cancer  1 32.82534         4   12.94407
Dose:Cancer  1  7.31735         3    5.62672

> Pearson(mymod)                                    
Degrees of Freedom: 8 Total; 3 Residual
Residual Pearson X^2 4.228423 

\end{verbatim}

Now it is clear we can't proceed very much.  Either of these last two
models seems adequate and both are comparatively easy to interpret.

Oh yeah?  How does one interpret these models, and what does it all
mean.  The final model (the one I would probably end up with) has an
overall goodness of fit statistic of (likelihood ratio) 5.62 or
(Pearson) 4.22, on 3 degrees of freedom.  The upper 95\% point of the
appropriate chi-square is
\begin{verbatim}
> qchisq(.95,3)                                    
[1] 7.814728
\end{verbatim}

So, on an overall basis, the model fits very well.  The model itself
says that {\tt Age} is indepedent of everything, but that {\tt Cancer}
does depend on {\tt Dose}.  We could look at the actual parameter
estimates, but would not learn very much.  

However, we are not yet done.  We certainly want to look at some plots
and diagnostics, even though we only have 8 data points.  We {\em
  think\/} we have a good model, do the diagnostics confuse the issue.

Start with a plot of the model/residuals.
\begin{verbatim}
> par(mfrow=c(2,3))
> plot(mymod)      
\end{verbatim}

\psfig{figure=cat3.eps,angle=-90,width=6in}

This doesn't show anything at all untoward.  


\section*{Another example, Graduate Admissions.}

Consider the following data which gives admissions data for the six
largest graduate programs at U.C. Berkeley for 1973.  
\begin{center}
  \begin{tabular}{c | r r | r r}
\hline
 & \multicolumn{2}{c|}{Men} & \multicolumn{2}{c}{Women} \\
 & \multicolumn{1}{c}{Number} & \multicolumn{1}{c|}{Number not} &
 \multicolumn{1}{c}{Number} & \multicolumn{1}{c}{Number not} \\
Major &\multicolumn{1}{c}{admitted} & \multicolumn{1}{c|}{admitted} &
\multicolumn{1}{c}{admitted} & \multicolumn{1}{c}{admitted} \\ \hline
A& 512& 313&  89&  19\\
B& 353& 207&  17&   8\\
C& 120& 205& 202& 391\\
D& 138& 279& 131& 244\\
E&  53& 138&  94& 299\\
F&  22& 351&  24& 317\\\hline
  \end{tabular}
\end{center}

Is there an issue of sex bias in admissions, or not?  

One can argue about what is a response or explanatory variable, but to
me it makes most sense to consider both {\tt Major} and {\rm Sex} as
explanatory variables and consider just {\tt Admit} as a response.

Start by reading in the data, then think about how you want to answer
the question.
\begin{verbatim}
> berk <- fac.design( c(2,2,6), fac=list(Admit=c("Yes","No"),   
+ Sex=c("Men","Women"), Major=c("A","B","C","D","E","F")))   
> berk$Resp <- scan("/tmp/berk")                                
> berk
   Admit   Sex Major Resp 
 1   Yes   Men     A  512
 2    No   Men     A  313
 3   Yes Women     A   89
 4    No Women     A   19
 5   Yes   Men     B  353
 6    No   Men     B  207
 7   Yes Women     B   17
 8    No Women     B    8
 9   Yes   Men     C  120
10    No   Men     C  205
11   Yes Women     C  202
12    No Women     C  391
13   Yes   Men     D  138
14    No   Men     D  279
15   Yes Women     D  131
16    No Women     D  244
17   Yes   Men     E   53
18    No   Men     E  138
19   Yes Women     E   94
20    No Women     E  299
21   Yes   Men     F   22
22    No   Men     F  351
23   Yes Women     F   24
24    No Women     F  317
\end{verbatim}

The first thing one might want to do is to compare the admission rates
for men and women.  One might start by comparing the percent of women
who are admitted and the percent of men admitted.
\begin{verbatim}
> sum(Resp[Sex=="Men"])           
[1] 2691
> sum(Resp[Sex=="Men" & Admit=="Yes"])
[1] 1198
> 
> sum(Resp[Sex=="Women"])               
[1] 1835
> sum(Resp[Sex=="Women" & Admit=="Yes"])
[1] 557
> 
> 1198/2691
[1] 0.4451877
> 557/1835
[1] 0.3035422
\end{verbatim}

That is quite revealing.  Only 30\% of the women are admitted while
almost 45\% of the men are.  Is it time to go and call the lawyers and
argue discrimination?   How could we have been mislead by this?  

The next logical step is to look at the admission rates for the
different programs, which are
\begin{verbatim}
> applied <- Resp[Admit=="Yes"] + Resp[Admit=="No"]
> applied
   1   3   5  7   9  11  13  15  17  19  21  23 
 825 108 560 25 325 593 417 375 191 393 373 341
> 
> round(Resp[Admit=="Yes"]/applied, dig=3)
     1     3    5    7     9    11    13    15    17    19    21   23 
 0.621 0.824 0.63 0.68 0.369 0.341 0.331 0.349 0.277 0.239 0.059 0.07
\end{verbatim}
Here the first two numbers are the admission rates for Men and Women
in Program A, the next two for Program B, etc.   There is much less
difference here.  This is easier to read below,
\begin{verbatim}
> dodo <- round(Resp[Admit=="Yes"]/applied, dig=3)
> dodo[(0:5)*2 + 1]
     1    5     9    13    17    21 
 0.621 0.63 0.369 0.331 0.277 0.059
> dodo[(0:5)*2 + 2]
     3    7    11    15    19   23 
 0.824 0.68 0.341 0.349 0.239 0.07
\end{verbatim}

Now there appears to be much less difference between Men and Women,
what is going on?   Time to call off the lawyers, or have we made a
mistake.  Perhaps the best way to proceed from here (after we decide
we think we know what is going on), is to try to fit models.

We have two explanatory variables, {\tt Sex} and {\tt Major}, so the
minimal model must include at least the interaction between these two.

Thus our starting model is
\begin{verbatim}
> anova(mymod <- glm(Resp ~ Sex*Major, family=poisson)) 
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
          Df Deviance Resid. Df Resid. Dev 
     NULL                    23   2650.095
      Sex  1  162.874        22   2487.221
    Major  5  159.522        17   2327.699
Sex:Major  5 1220.615        12   1107.084
\end{verbatim}

and we need to ask if there is any effect due to the {\tt Admit}
variable.  But lets be careful here.  What we are really asking about
is the interaction between admissions rates and sex.  We {\em know},
or at least should suspect that admissions rates will vary by major,
so maybe the null (that is uninteresting) model should be

\begin{verbatim}
> anova(mymod <- glm(Resp ~ Sex*Major + Admit*Major, family=poisson))        
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev 
       NULL                    23   2650.095
        Sex  1  162.874        22   2487.221
      Major  5  159.522        17   2327.699
      Admit  1  230.028        16   2097.671
  Sex:Major  5 1220.615        11    877.056
Admit:Major  5  855.321         6     21.736
\end{verbatim}

This is a reasonable model.  On an overall level, it does not
completely explain that data, {\tt  pchisq(21,6) $=$ 0.9981654}, but
it is a great improvement over the null model.  Now to answer the real
question.
\begin{verbatim}
> anova(mymod <- glm(Resp ~ Sex*Major + Admit*Major + Admit*Sex              
+ ,family=poisson))
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev 
       NULL                    23   2650.095
        Sex  1  162.874        22   2487.221
      Major  5  159.522        17   2327.699
      Admit  1  230.028        16   2097.671
  Sex:Major  5 1220.615        11    877.056
Admit:Major  5  855.321         6     21.736
  Admit:Sex  1    1.531         5     20.204
\end{verbatim}

Here we see that the improvement in goodness of fit by adding the
interaction between admit and sex is only 1.53, with a p-value of
0.78.  This can hardly be considered strong (or any evidence).  

So what is going on.  If it wasn't already obvious, the men are
applying to programs that have higher admission rates, so their
overall admission rate is higher.  

What about diagnostics for this model?  Our resulting model only has 5
degrees of freedom left, so we don't have a lot of leeway.  

\psfig{figure=cat3-2.eps,angle=-90,width=6in}

The one really big residual corresponds to the small number of women
not admitted to program A.  The model predicts that more women should
not be admitted.  With this discussion you also see some of the
problems with residuals in this case.  Although we have 24 residuals
we know that there are {\bf lots} of constraints imposed by the model,
so the residuals are far from independent from each other. 


% Local Variables: 
% mode: latex
% eval: (LaTeX-math-mode t)
% eval: (read-abbrev-file "~/.abbrev_defs")
% abbrev-mode: t
% TeX-master: t
% TeX-master: t
% End: 
