\section*{Loglinear Models}

Notice that for both of the $2 \times 2$ table examples the fitted
values were of the form
$$
\hat{m_{i,j}} = \frac{x_{i,+} \cdot x_{+,j}}{x_{+,+}} $$

Now, take logarithms of this expression and notice that
$$
\log(\hat{m_{i,j}}) = - \log(x_{+,+}) + \log(x_{i,+}) +\log(x_{+,j})
$$

In other words, the fitted values are {\em linear\/} functions of the
margins on the logarithmic scale.

Now lets consider the model for the fitted values which says,
$$
\log({m_{i,j}}) = u + u_{1(i)} + u_{2(j)} + u_{12(i,j)} $$

with the constraints that
$\sum_i{u_{1(i)}} = 0$, and
$\sum_j{u_{2(j)}} = 0$, and
$\sum_i{u_{12(i,j)}} = \sum_j{u_{12(i,j)}} = 0$.

These models are the analogue of ANOVA models and are ``natural''
models for Poisson/Multinomial data.  In particular, the notions of
independence and conditional independence can be easily considered
using these models.  

Another way of viewing the loglinear model is to say that it
represents the logarithm of the expected values ($m_{i,j}$) in terms
of main effects and interactions---which is just what the ANOVA models
did.

A third, and most general, way of describing the models is to say that
the loglinear models posit that
$$ \log(m) \in {\cal M} $$
where ${\cal M}$ is a linear space which describes the hypothesis of
interest.  That is about the same as saying that an ANOVA model is
just a linear model, and just as in the ANOVA situation it ignores the
usefulness of various forms for the linear space ${\cal M}$.  

The first formulation, in terms of the u-term parameterization, is one
useful way of describing the linear space ${\cal M}$. 

To make this a little more concrete consider a $3 \times 3$ table,
\begin{center}
  \begin{tabular}{| l | l | l |}
    \hline
  $x_{1,1}$ & $x_{1,2}$ & $x_{1,3}$ \\ \hline
  $x_{2,1}$ & $x_{2,2}$ & $x_{2,3}$ \\ \hline
  $x_{3,1}$ & $x_{3,2}$ & $x_{3,3}$ \\ \hline
  \end{tabular}
\end{center}
where the two categorizing variables are both responses. (as in the
lizard example).   Lets assume that the expected values of the
response (under some model to be determined) are

\begin{center}
$m = $  \begin{tabular}{| l | l | l |}
    \hline
  $m_{1,1}$ & $m_{1,2}$ & $m_{1,3}$ \\ \hline
  $m_{2,1}$ & $m_{2,2}$ & $m_{2,3}$ \\ \hline
  $m_{3,1}$ & $m_{3,2}$ & $m_{3,3}$ \\ \hline
  \end{tabular}
\end{center}

The model which asserts independence of the two categories is
$$
ln(m) \in {\rm span}\left\{
\begin{array}{| l | l | l |}
\hline
1 & 1 & 1 \\ \hline
1 & 1 & 1 \\ \hline
1 & 1 & 1 \\ \hline
\multicolumn{3}{c}{u}
\end{array} , 
\begin{array}{| l | l | l |}
\hline
1 & 1 & 1 \\ \hline
0 & 0 & 0 \\ \hline
0 & 0 & 0 \\ \hline
\multicolumn{3}{c}{u_{1(1)}}
\end{array} , 
\begin{array}{| l | l | l |}
\hline
0 & 0 & 0 \\ \hline
1 & 1 & 1 \\ \hline
0 & 0 & 0 \\ \hline
\multicolumn{3}{c}{u_{1(2)}}
\end{array} , 
\begin{array}{| l | l | l |}
\hline
0 & 0 & 0 \\ \hline
0 & 0 & 0 \\ \hline
1 & 1 & 1 \\ \hline
\multicolumn{3}{c}{u_{1(3)}}
\end{array} , 
\begin{array}{| l | l | l |}
\hline
1 & 0 & 0 \\ \hline
1 & 0 & 0 \\ \hline
1 & 0 & 0 \\ \hline
\multicolumn{3}{c}{u_{2(1)}}
\end{array} , 
\begin{array}{| l | l | l |}
\hline
0& 1  & 0 \\ \hline
0& 1  & 0 \\ \hline
0& 1  & 0 \\ \hline
\multicolumn{3}{c}{u_{2(2)}}
\end{array} , 
\begin{array}{| l | l | l |}
\hline
0  & 0& 1 \\ \hline
0  & 0& 1 \\ \hline
0  & 0& 1 \\ \hline
\multicolumn{3}{c}{u_{2(3)}}
\end{array}
\right\}
$$

That is, any set of fitted values which can be written as a loglinear
combination of the above vectors satisfies independence.  One can also
see from this discussion where the constraints of our particular
parameterization come from.  (I will probably write out a few more of
these gory spanning vectors, but not too many.  They are just too hard
for me to type!).

In general, we will oscillate between talking about abstract loglinear
models, specific parameterizations, and the compromise model formulae
of a package like S-Plus.  Mostly (hopefully all the time, but I'll
hedge a little), the model embodies some underlying notion and one
should be careful to talk about the underlying model and not
over-interpret the specific parameter estimates.  

Now that we have written a model in this general (loglinear) form, an
obvious question arises. 
\begin{itemize}
\item Why not model $ m \in {\cal M}$, a linear model?
\end{itemize}
or on some other scale?

The main reason is that if the loglinear model is written as a model
for probabilities (which come directly from the fitted values), then
we can be sure that the estimated probabilities will lie between 0 and
1.  Without a lot of extra constraints and complications, linear
models can not guarantee this.   A more (or maybe less) compelling
reason is that under the loglinear model the likelihood function
factorizes in a natural way.


\section*{Bigger Two-way Tables}

Consider the following cross-classification of 1019 school children
into a teachers rating of schoolwork and an ``independent'' rating of
homework conditions.  (Independent in the non-technical sense.)  We
would like to ask whether the two responses are independent of each
other.
\begin{center}
  \begin{tabular}{l | r r r r r | r}
    &\multicolumn{5}{c}{Homework Conditions} \\
    Rating & A & B & C & D & E \\\hline
    A & 141 & 67 & 114 & 79 & 39 & 440 \\
    B & 131 & 66 & 143 & 72 & 35 & 447 \\
    C& 36 & 14 & 38 & 28 & 16 & 132 \\ \hline
     & 308 & 147 & 295 & 179 & 90 & 1019 \\
  \end{tabular}
\end{center}
Under independence the fitted values are generated by multiplying the
row margins but the column margins and dividing by the grand total.
This doesn't have to be obvious, but it is true.

Looking at the data we see that as homework conditions get worse (A to
E) more people get a bad teachers rating.  Is this significant, or am
I imagining it?  The table of fitted values, under independence is.
\begin{center}
  \begin{tabular}{| r r r r r |}
    \hline
    133 & 63.5 & 121.4 & 77.3 & 38.9 \\
    135.1 & 64.5 & 129.4 & 78.5 & 39.5 \\
    39.9 & 19.0 & 38.2 & 23.2 & 11.6 \\\hline
  \end{tabular}
\end{center}
with ${\rm X}^2 = 8.21 $ on $(3-1) \times (5-1) = 8$ degrees of
freedom.
The value of ${\rm X}^2$ tells us that there is not strong reason to
reject the assumption of independence.  What about looking at the
residuals (generally a good thing to do), but which residuals? (there
are at least two commonly used residuals \ldots).

The observations $ x_{i,j} - \hat{m_{i,j}}$ have different variances
(which is approximately $\hat{m_{i,j}}$) so it makes sense to look at
the Pearson standardized residuals
$$
r_{i,j} = \frac{ x_{i,j} - \hat{m_{i,j}}}{\sqrt{\hat{m_{i,j}}}} $$

For these data these numbers are
\begin{center}
  \begin{tabular}{| r r r r r |}
    \hline
    0.69 & 0.44 & -0.67 & 0.19 & 0.02 \\
    -0.35 & 0.19 & 1.20 & -0.73 & -0.72 \\
    -0.62 & -1.15 & 0.0 & 1.0 & 1.29 \\ \hline
  \end{tabular}
\end{center}

None of these residuals is very large (compare to ${\rm N}(0,1)$ nor
does there appear to be any strong pattern to them.  It seems that it
is pretty safe to assume that the homework conditions and teachers
ratings are independent.  


\paragraph{Informal Exercise}
\begin{enumerate}
\item Verify that the table of fitted values can be written as 
  $\log({m_{i,j}}) = \mu + \mu_{1(i)} + \mu_{2(j)} + \mu_{12(i,j)} $
  with the constraints that
  $\sum_i{\mu_{1(i)}} = 0$, and
  $\sum_j{\mu_{2(j)}} = 0$.
\item Use the residuals to calculate ${\rm X}^2$.
\end{enumerate}


\paragraph{Example}
The following example is drawn from Problem 2.2 in the Xmas-Tree book.
The data below, from the National Crime Survey, represent
repeat-victimization information (i.e., successive pairs of
victimization for households in the survey).  If, over th course of the
period in question, a household is victimized {\em m\/} times, then these
{\em m\/} victimizations will lead to $m - 1$ entries in the table,
one for each successive pair of victimizations.  Here A is assault; B
is burglary; HL is household larceny; MV is motor vehicle theft; PL is
personal larceny; PP/PS is pocket picking and purse snatching; Ra is
rape; and Ro is robbery.

\begin{center}
  \begin{tabular}{l | r r r r r r r r | r}
    & \multicolumn{8}{c}{Second Victimization} & \\
    First & Ra & A & Ro & PP/PS & PL & B & HL & MV & Totals\\\hline
    Ra&    26&    50&    11&     6&    82&    39&    48&    11& 273 \\     
    A&    65&  2997&   238&    85&  2553&  1083&  1349&   216& 8586 \\     
    Ro&    12&   279&   197&    36&   459&   197&   221&    47& 1448 \\    
    PP/PS&     3&   102&    40&    61&   243&   115&   101&    38& 703 \\  
    PL&    75&  2628&   413&   229& 12137&  2658&  3689&   687& 22526 \\   
    B&    52&  1117&   191&   102&  2649&  3210&  1973&   301& 9595 \\     
    HL&    42&  1251&   206&   117&  3757&  1962&  4646&   391& 12372 \\   
    MV&     3&   221&    51&    24&   678&   301&   367&   269& 1914 \\\hline
    Total & 278 & 8645 & 1347 & 660 & 22558 & 9565 & 12394 & 1960 & \\     
  \end{tabular}
\end{center}

There are many things to ask of these data, including:
\begin{itemize}
\item Test these data for homogeneity of proportions
\item Examine a table of the square roots of the contribution from
  each cell to ${\rm X}^2$, i.e. 
$({\rm Observed} - {\rm Expected})/\sqrt{{\rm Expected}}$  Which cell
stands out?
\item Provide an interpretation to your findings above.
\item Why might you not have bothered to do the cull calculation of
  ${\rm X}^2$ after you found the contribution of the $(1,1)$ cell?
\end{itemize}

Here is a cursory set of answers to these questions.  A quick
examination of the data shows that the row and colun margins are very
similar and the diagonal cells of teh table (corresponding to repeat
victimizations of the same type) clearly dominate the table. 

Thus prior to formal modelling we have good reason to expect that the
model of homogeneity of row proportions will {\em not\/} be
appropriate.  

Under the model of homogeneity of row proportions, with data $ x_{i,j}$, and estimated expected values
$$
\hat{m_{i,j}} = ( x_{i,+} \cdot x_{+,j}) / x_{+,+} $$
then the fitted values are as follows.
\begin{verbatim}
Reiss <- fac.design(c(8,8),factor.names=c("Second","First"))
Reiss$Resp <- scan("reiss.dat")
> Reiss[1:10,]
   Second First Resp 
 1  Scnd1 Frst1   26
 2  Scnd2 Frst1   50
 3  Scnd3 Frst1   11
 4  Scnd4 Frst1    6
 5  Scnd5 Frst1   82
 6  Scnd6 Frst1   39
 7  Scnd7 Frst1   48
 8  Scnd8 Frst1   11
 9  Scnd1 Frst2   65
10  Scnd2 Frst2 2997
\end{verbatim}
Now, fit the model and look at the Analysis of Deviance statistics,
\begin{verbatim}
> mymod<-glm(Resp ~ First + Second, family=poisson, data=Reiss)  
> anova(mymod)
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev 
  NULL                    63   127032.9
 First  7 58883.15        56    68149.7
Second  7 59387.44        49     8762.3
\end{verbatim}
What do the fitted values look like?

\begin{verbatim}
> t(matrix(round(fitted(mymod),digits=3),8,8))                            
        [,1]     [,2]    [,3]    [,4]     [,5]     [,6]     [,7]    [,8] 
[1,]   1.322   41.111   6.406   3.139  107.275   45.487   58.940   9.321
[2,]  41.579 1292.978 201.462  98.712 3373.857 1430.576 1853.692 293.145
[3,]   7.012  218.056  33.976  16.647  568.990  241.262  312.619  49.438
[4,]   3.404  105.866  16.495   8.082  276.243  117.132  151.776  24.002
[5,] 109.036 3390.716 528.316 258.863 8847.631 3751.555 4861.137 768.745
[6,]  46.465 1444.924 225.137 110.312 3770.342 1598.693 2071.532 327.594
[7,]  59.913 1863.117 290.297 142.239 4861.560 2061.389 2671.078 422.407
[8,]   9.269  288.232  44.910  22.005  752.104  318.906  413.227  65.348
> 
\end{verbatim}
Notice the ``trick'' to get the data to print out in a nice tabular
format.

The goodness-of-fit statistics for the model of homogeneity of
proportions is ${\rm X}^2 = 11,189$ and ${\rm G}^2=8,762$, both of
with should be asymptotically distributed as a $\chi^2_{49}$ random
variable.  The likelihood ratio statistic comes directly from the
output, 
\begin{verbatim}
> anova(mymod)
Analysis of Deviance Table

Poisson model

Response: Resp

Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev 
  NULL                    63   127032.9
 First  7 58883.15        56    68149.7
Second  7 59387.44        49     8762.3
> 
\end{verbatim}

The Pearson statistic requires a custom function.  I wrote a little
function, called Pearson, to print the statistics and degrees of
freedom.
\begin{verbatim}
> Pearson(mymod)
Degrees of Freedom: 64 Total; 49 Residual
Residual Pearson X^2 11189.88 
\end{verbatim}

The Pearson function itself if,
\begin{verbatim}
> Pearson
function(mymod)
{
# This is not completely bullet proof.
        nobs <- length(mymod$residuals)
        rdf <- mymod$df.resid
        if(is.null(rdf))
                rdf <- nobs - rank
        cat("Degrees of Freedom:", nobs, "Total;", rdf, "Residual\n")
        cat("Residual Pearson X^2", format(sum((mymod$y - fitted(mymod))^2/
                fitted(mymod))), "\n")
        invisible(mymod)
}
\end{verbatim}

It also makes a lot of sense to examine the residuals from this
model.  There are many different residuals one could consider, but a
natural candidate is the so-called pearson residual which is the
component of Pearson's ${\rm X}^2$,
$$
({\rm Observed} - {\rm Expected}) / \sqrt{{\rm Expected}} $$

For these data, the numbers are
\begin{verbatim}

> t(matrix(round( residuals(mymod,type="pearson"), dig=3), 8,8))   
        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8] 
[1,]  21.457   1.386   1.815   1.615  -2.440  -0.962  -1.425   0.550
[2,]   3.632  47.388   2.574  -1.380 -14.132  -9.189 -11.722  -4.506
[3,]   1.883   4.127  27.967   4.743  -4.611  -2.850  -5.182  -0.347
[4,]  -0.219  -0.376   5.787  18.613  -2.000  -0.197  -4.121   2.857
[5,]  -3.259 -13.098  -5.017  -1.856  34.970 -17.854 -16.812  -2.948
[6,]   0.812  -8.627  -2.275  -0.791 -18.262  40.299  -2.165  -1.469
[7,]  -2.314 -14.181  -4.947  -2.116 -15.842  -2.189  38.213  -1.528
[8,]  -2.059  -3.960   0.909   0.425  -2.702  -1.003  -2.274  25.192
\end{verbatim}

The contribution to the overall ${\rm X}^2$ from just the first
component of the residuals is $(21.457)^2 = 460.40 $, which is much
bigger than any upper bound of a  $\chi^2_{49}$.  


% 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: 


