% !Rnw weave = Sweave
%\VignetteEngine{Sweave}
%\VignetteIndexEntry{Tutorial: Working with categorical data with R and the vcd package}
%\VignetteDepends{vcd,gmodels,ca}
%\VignetteKeywords{contingency tables, mosaic plots, sieve plots, categorical data, independence, conditional independence, R}
%\VignettePackage{vcdExtra}

\documentclass[10pt,twoside]{article}
\usepackage{Sweave}
\usepackage{bm}
\usepackage[toc]{multitoc}   % for table of contents

% from Z.cls
\usepackage[authoryear,round,longnamesfirst]{natbib}
\bibpunct{(}{)}{;}{a}{}{,}
\bibliographystyle{jss}

\usepackage{hyperref}
\usepackage{color}
%% colors
\definecolor{Red}{rgb}{0.7,0,0}
\definecolor{Blue}{rgb}{0,0,0.8}
\hypersetup{%
  hyperindex = {true},
  colorlinks = {true},
%  linktocpage = {true},
  plainpages = {false},
  linkcolor = {Blue},
  citecolor = {Blue},
  urlcolor = {Red},
  pdfstartview = {Fit},
  pdfpagemode = {UseOutlines},
  pdfview = {XYZ null null null}
}
%\AtBeginDocument{
%  \hypersetup{%
%    pdfauthor = {Michael Friendly},
%    pdftitle = {Tutorial: Working with categorical data with R and the vcd package},
%    pdfkeywords = {contingency tables, mosaic plots, sieve plots, categorical data, independence, conditional independence, R}
%  }
%}


% math stuff
\newcommand*{\given}{\ensuremath{\, | \,}}
\renewcommand*{\vec}[1]{\ensuremath{\bm{#1}}}
\newcommand{\mat}[1]{\ensuremath{\bm{#1}}}
\newcommand{\trans}{\ensuremath{^\mathsf{T}}}
\newcommand{\diag}[1]{\ensuremath{\mathrm{diag} (#1)}}
\def\binom#1#2{{#1 \choose #2}}%
\newcommand{\implies}{ \ensuremath{\mapsto} }

\newenvironment{equation*}{\displaymath}{\enddisplaymath}%

\newcommand{\tabref}[1]{Table~\ref{#1}}
\newcommand{\figref}[1]{Figure~\ref{#1}}
\newcommand{\secref}[1]{Section~\ref{#1}}
\newcommand{\loglin}{loglinear }

%\usepackage{thumbpdf}

% page dimensions
\addtolength{\hoffset}{-1.5cm}
\addtolength{\textwidth}{3cm}
\addtolength{\voffset}{-1cm}
\addtolength{\textheight}{2cm}

% Vignette examples
\newcommand*{\Example}{\fbox{\textbf{\emph{Example}}:} }

% R stuff
\newcommand{\var}[1]{\textit{\texttt{#1}}}
\newcommand{\data}[1]{\texttt{#1}}
\newcommand{\class}[1]{\textsf{"#1"}}
%% \code without `-' ligatures
\def\nohyphenation{\hyphenchar\font=-1 \aftergroup\restorehyphenation}
\def\restorehyphenation{\hyphenchar\font=`-}
{\catcode`\-=\active%
  \global\def\code{\bgroup%
    \catcode`\-=\active \let-\codedash%
    \Rd@code}}
\def\codedash{-\discretionary{}{}{}}
\def\Rd@code#1{\texttt{\nohyphenation#1}\egroup}

\newcommand{\codefun}[1]{\code{#1()}}
\let\proglang=\textsf
\newcommand{\pkg}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\newcommand{\Rpackage}[1]{{\textsf{#1}}}

%% almost as usual
\author{Michael Friendly\\York University, Toronto}
\title{Working with categorical data with \proglang{R} and the \pkg{vcd} and \pkg{vcdExtra} packages}
\date{\footnotesize{Using \Rpackage{vcdExtra} version \Sexpr{packageDescription("vcdExtra")[["Version"]]}
  and \Rpackage{vcd} version \Sexpr{packageDescription("vcd")[["Version"]]}; Date: \Sexpr{Sys.Date()}}}


%% for pretty printing and a nice hypersummary also set:
%\Plainauthor{Michael Friendly} %% comma-separated
%\Shorttitle{vcd tutorial} %% a short title (if necessary)
%\Plaintitle{Tutorial: Working with categorical data with R and the vcd package}

%\SweaveOpts{engine=R,eps=TRUE,height=6,width=7,results=hide,fig=FALSE,echo=TRUE}
\SweaveOpts{engine=R,height=6,width=7,results=hide,fig=FALSE,echo=TRUE}
\SweaveOpts{prefix.string=fig/vcd-tut,eps=FALSE}
\SweaveOpts{keep.source=TRUE}
%\SweaveOpts{concordance=TRUE}
\setkeys{Gin}{width=0.7\textwidth}

<<preliminaries,echo=FALSE,results=hide>>=
set.seed(1071)
#library(vcd)
library(vcdExtra)
library(ggplot2)
#data(Titanic)
data(HairEyeColor)
data(PreSex)
data(Arthritis)
art <- xtabs(~Treatment + Improved, data = Arthritis)
if(!file.exists("fig")) dir.create("fig")
@



%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle

%% an abstract and keywords
\begin{abstract}
This tutorial describes the creation of frequency and contingency tables
from categorical variables, along with tests of independence, measures
of association, and methods for graphically displaying results. 
The framework is provided by the \proglang{R} package \pkg{vcd},
but other packages are used to help with various tasks.  
The \pkg{vcdExtra} package extends the graphical and statistical methods provided by
\pkg{vcd}.
\end{abstract}

%\keywords{contingency tables, mosaic plots, sieve plots, 
%categorical data, independence, conditional independence, generalized linear models,
%\proglang{R}}
%\Plainkeywords{contingency tables, mosaic plots, 
%  sieve plots, categorical data, independence, 
%  conditional independence, generalized linear models, R}

{\small
% \sloppy
% \begin{multicols}{2}
 \tableofcontents
% \end{multicols}
}


\section[Introduction]{Introduction}\label{sec:intro}
%% Note: If there is markup in \(sub)section, then it has to be escape as above.


This tutorial, part of the \pkg{vcdExtra} package,
describes how to work with categorical data in the context of
fitting statistical models in \proglang{R} and visualizing
the results using the \pkg{vcd} and \pkg{vcdExtra} packages.
It focuses first on methods and tools for creating
and manipulating \proglang{R} data objects which represent
frequency and contingency tables involving
categorical variables.

Further sections describe some simple methods for calculating 
tests of independence and measures of association amomg 
categorial variables, and also
methods for graphically displaying results.

There is much more to the analysis of categorical data than is described here,
where the emphasis is on cross-tabulated tables of frequencies (``contingency tables''),
statistical tests, associated \loglin\ models, and visualization of \emph{how}
variables are related.

A more general treatment of graphical methods for categorical data is contained
in my book, \emph{Visualizing Categorical Data} \citep{vcd:Friendly:2000}, for which
\pkg{vcd} is a partial \proglang{R} companion, covering topics not otherwise
available in \proglang{R}.  On the other hand, the implementation of graphical
methods in \pkg{vcd} is more general in many respects than what I provided in
\proglang{SAS}. Statistical models for categorical data in \proglang{R} have
been extended considerably with the \pkg{gnm} package for generalized \emph{nonlinear}
models.  The \pkg{vcdExtra} package extends \pkg{vcd} methods to models fit using
\codefun{glm} and \codefun{gnm}.

A more complete theoretical description of these statistical
methods is provided in Agresti's \citeyearpar{vcd:Agresti:2002,Agresti:2013}
\emph{Categorical Data Analysis}. For this, see
the \proglang{Splus/R} companion by Laura Thompson,
\url{https://home.comcast.net/~lthompson221/Splusdiscrete2.pdf}
and Agresti's support web page,
\url{http://www.stat.ufl.edu/~aa/cda/cda.html}.



\section[Creating frequency tables]{Creating and manipulating frequency tables}\label{sec:creating}

\proglang{R} provides many methods for creating frequency and contingency tables. Several are
described below. In the examples below, we use some real examples and some anonymous
ones, where the variables \code{A}, \code{B}, and \code{C}  represent
categorical variables, and \code{X} represents an arbitrary \proglang{R} data object.

The first thing you need to know is that categorical data can be represented in
three different forms in \proglang{R}, and it is sometimes necessary to convert
from one form to another, for carrying out statistical tests, fitting models
or visualizing the results.  Once a data object exists in \proglang{R}, 
you can examine its complete structure with
the \codefun{str} function, or view the names of its components with the
\codefun{names} function.
 
\begin{description}
  \item[case form] a data frame containing individual observations, with one or
  more factors, used as the classifying variables. In case form, there may also
  be numeric covariates. 
  The total number of observations
  is \code{nrow(X)}, and the number of variables is \code{ncol(X)}.
  	
\Example
The \data{Arthritis} data is available in case form in the \pkg{vcd} package. 
There are two explanatory factors: \code{Treatment} and \code{Sex}. \code{Age}
is a numeric covariate, and \code{Improved} is the response--- an ordered factor,
with levels 
\code{\Sexpr{paste(levels(Arthritis$Improved),collapse=' < ')}}.  
Excluding \code{Age}, we would have
a $2 \times 2 \times 3$ contingency table for \code{Treatment}, \code{Sex} and \code{Improved}.
%\code{"None" < "Some" < "Marked"}.
<<case-form,results=verbatim>>=
names(Arthritis)      # show the variables
str(Arthritis)        # show the structure
head(Arthritis,5)     # first 5 observations, same as Arthritis[1:5,] 
@
  
  \item[frequency form] a data frame containing one or more factors, and a frequency 
  variable, often called \code{Freq} or \code{count}.  The total number of observations
  is 
  \verb|sum(X$Freq)|,  
  \code{sum(X[,"Freq"])} or some equivalent form.
  The number of cells in the table is \code{nrow(X)}.

\Example 
For small frequency tables, it is often convenient to enter them in frequency form
using \codefun{expand.grid} for the factors and \codefun{c} to list the counts in a vector.
The example below, from \cite{vcd:Agresti:2002} gives results for the 1991 General Social Survey,
with respondents classified by sex and party identification.
<<frequency-form,results=verbatim>>=
# Agresti (2002), table 3.11, p. 106
GSS <- data.frame(
  expand.grid(sex=c("female", "male"), 
              party=c("dem", "indep", "rep")),
  count=c(279,165,73,47,225,191))
GSS
names(GSS)
str(GSS)
sum(GSS$count)
@

  \item[table form]  a matrix, array or table object, whose elements are the frequencies
  in an $n$-way table.  The variable names (factors) and their levels are given by
  \code{dimnames(X)}. The total number of observations
  is \code{sum(X)}.  The number of dimensions of the table is \code{length(dimnames(X))},
  and the table sizes are given by \code{sapply(dimnames(X), length)}.
  	
\Example 
The \data{HairEyeColor} is stored in table form in \pkg{vcd}.  
<<table-form1,results=verbatim>>=
str(HairEyeColor)                      # show the structure
sum(HairEyeColor)                      # number of cases
sapply(dimnames(HairEyeColor), length) # table dimension sizes
@

\Example 
Enter frequencies in a matrix, and assign \code{dimnames},
giving the variable names and category labels.  Note that, by default,
\codefun{matrix} uses the elements supplied by \emph{columns} in the
result, unless you specify \code{byrow=TRUE}.
<<table-form2,results=verbatim>>=
## A 4 x 4 table  Agresti (2002, Table 2.8, p. 57) Job Satisfaction
JobSat <- matrix(c(1,2,1,0, 3,3,6,1, 10,10,14,9, 6,7,12,11), 4, 4)
dimnames(JobSat) = list(income=c("< 15k", "15-25k", "25-40k", "> 40k"),
                satisfaction=c("VeryD", "LittleD", "ModerateS", "VeryS"))
JobSat
@
\data{JobSat} is a matrix, not an object of \code{class("table")}, and some functions
are happier with tables than matrices.
You can coerce it to a table with \codefun{as.table},
<<table-form3,results=verbatim>>=
JobSat <- as.table(JobSat)
str(JobSat)
@

\end{description}

\subsection[Ordered factors]{Ordered factors and reordered tables}\label{sec:ordered-factors}
In table form, the values of the table factors are ordered by their position in the table.
Thus in the \data{JobSat} data, both \code{income} and \code{satisfaction} represent ordered
factors, and the \emph{positions} of the values in the rows and columns reflects their
ordered nature.

Yet, for analysis, there are time when you need \emph{numeric} values for the levels
of ordered factors in a table, e.g., to treat a factor as a quantitative variable.
In such cases, you can simply re-assign the \code{dimnames} attribute of the table
variables.  For example, here, we assign numeric values to \code{income} as the middle of their
ranges, and treat \code{satisfaction} as equally spaced with integer scores.

<<relevel,results=hide,eval=FALSE>>=
dimnames(JobSat)$income<-c(7.5,20,32.5,60)
dimnames(JobSat)$satisfaction<-1:4
@

For the  \data{HairEyeColor} data, hair color and eye color are ordered arbitrarily.
For visualizing the data using mosaic plots and other methods described below, it 
turns out to be more useful to assure that both hair color and eye color are
ordered from dark to light.
Hair colors are actually ordered this way already, and it is easiest to re-order
eye colors by indexing. Again \codefun{str} is your friend.

<<reorder1,results=verbatim>>=
HairEyeColor <- HairEyeColor[, c(1,3,4,2), ]
str(HairEyeColor)
@
This is also the order for both hair color and eye color shown in 
the result of a correspondence analysis (\figref{fig:ca-haireye}) below.

With data in  case form or frequency form, when you have ordered factors
represented with character values, you must ensure that they are treated
as ordered in \proglang{R}.%
\footnote{In \proglang{SAS}, many procedures offer the option
	\code{order = data | internal | formatted} to allow character values
	to be ordered according to (a) their order in the data set, (b)
	sorted internal value, or (c) sorted formatted representation
	provided by a \proglang{SAS} format.
}

Imagine that the \data{Arthritis} data was read from a text file.  
By default the \code{Improved} will be ordered alphabetically:
\code{Marked},
\code{None},
\code{Some}--- not what we want.  In this case, the function
\codefun{ordered} (and others) can be useful. 

<<reorder2,echo=TRUE,eval=FALSE>>=
Arthritis <- read.csv("arthritis.txt",header=TRUE)
Arthritis$Improved <- ordered(Arthritis$Improved, levels=c("None", "Some", "Marked"))
@

With this order of  \code{Improved}, the response in this data,
a mosaic display of \code{Treatment} and \code{Improved} (\figref{fig:arthritis})shows a clearly
interpretable pattern.
<<Arthritis,fig=TRUE,echo=FALSE,height=6,width=7,results=hide,include=FALSE>>=
mosaic(art, gp = shading_max, split_vertical = TRUE, main="Arthritis: [Treatment] [Improved]")
@

%\setkeys{Gin}{width=0.7\textwidth}
\begin{figure}[htb]
\begin{center}
%<<Arthritis,fig=TRUE,echo=FALSE,height=6,width=7,results=hide>>=
%mosaic(art, gp = shading_max, split_vertical = TRUE, main="Arthritis: [Treatment] [Improved]")
%@
\includegraphics[width=0.7\textwidth]{fig/vcd-tut-Arthritis}
\caption{Mosaic plot for the \data{Arthritis} data, showing the marginal model of independence
	for Treatment and Improved.  Age, a covariate, and Sex are ignored here.}
\label{fig:arthritis}
\end{center}
\end{figure}

Finally, there are situations where, particularly for display purposes, you
want to re-order the \emph{dimensions} of an $n$-way table, or change the
labels for the variables or levels.
This is easy when the data are in table form: \codefun{aperm} permutes
the dimensions, and assigning to \code{names} and \code{dimnames}
changes variable names and level labels respectively.
We will use the following version of \data{UCBAdmissions} in
\secref{sec:mantel} below.%
\footnote{
Changing \code{Admit} to \code{Admit?} might be useful for display purposes, but is
dangerous--- because it is then
difficult to use that variable name in a model formula.
See \secref{sec:tips} for options \code{labeling\_args} and \code{set\_labels}
to change variable and level names 
for displays in the \code{strucplot} framework.
}
<<reorder3,results=verbatim>>=
UCB <- aperm(UCBAdmissions, c(2, 1, 3))
dimnames(UCB)[[2]] <- c("Yes", "No")
names(dimnames(UCB)) <- c("Sex", "Admit?", "Department")
ftable(UCB)
@

%There is one subtle ``gotcha'' here:  \codefun{aperm} returns an object of class \class{"array"},
%whereas \data{UCBAdmissions} is of class \class{"table"}, so methods defined for \code{table}
%objects will not work on the permuted array.
%The solution is to reassign the \code{class} of the result of  \codefun{aperm}.
%
%<<reorder4,results=verbatim, keep.source=TRUE>>=
%class(UCBAdmissions)
%class(UCB)
%str(as.data.frame(UCBAdmissions))  # OK
%str(as.data.frame(UCB))            # wrong
%
%class(UCB) <- "table"
%str(as.data.frame(UCB))            # now OK
%@
%
\subsection[structable()]{\codefun{structable}}\label{sec:structable}
For 3-way and larger tables
the \codefun{structable} function in \pkg{vcd} provides a convenient and flexible tabular display.
The variables assigned to the rows and columns of a two-way display can be specified
by a model formula.
<<structable,results=verbatim>>=
structable(HairEyeColor)                   # show the table: default
structable(Hair+Sex ~ Eye, HairEyeColor)   # specify col ~ row variables
@
It also returns an object of class \code{"structable"} which may be plotted with 
\codefun{mosaic} (not shown here).
<<structable1,eval=FALSE>>=
HSE < - structable(Hair+Sex ~ Eye, HairEyeColor)   # save structable object
mosaic(HSE)                                        # plot it
@

\subsection[table() and friends]{\codefun{table} and friends}\label{sec:table}

You  can  generate frequency  tables from factor variables  using the  \codefun{table} function,  tables  of
proportions using  the \codefun{prop.table} function,  and marginal  frequencies using
\codefun{margin.table}.

<<setup>>=
 n=500
 A <- factor(sample(c("a1","a2"), n, rep=TRUE))
 B <- factor(sample(c("b1","b2"), n, rep=TRUE))
 C <- factor(sample(c("c1","c2"), n, rep=TRUE))
 mydata <- data.frame(A,B,C)
@

<<table-ex1>>=
# 2-Way Frequency Table
attach(mydata)
mytable <- table(A,B)   # A will be rows, B will be columns
mytable                 # print table

margin.table(mytable, 1) # A frequencies (summed over B)
margin.table(mytable, 2) # B frequencies (summed over A)

prop.table(mytable)    # cell percentages
prop.table(mytable, 1) # row percentages
prop.table(mytable, 2) # column percentages
@

\codefun{table} can  also  generate  multidimensional  tables  based  on  3  or  more
categorical variables. In  this case, use  the \codefun{ftable}  or \codefun{structable}
function to print  the
results more attractively.

<<table-ex2>>=
# 3-Way Frequency Table
mytable <- table(A, B, C)
ftable(mytable)
@

\codefun{table}  ignores missing values by default. 
To include \code{NA} as a category in counts, include the
table option  \code{exclude=NULL} if  the variable  is a  vector. If  the variable is a
factor you  have to  create a  new factor  using \code{newfactor  <- factor(oldfactor,
exclude=NULL)}. 

\subsection[xtabs()]{\codefun{xtabs}}\label{sec:xtabs}

The \codefun{xtabs} function allows you to create crosstabulations of data using formula style input.
This typically works with case-form data supplied in a data frame or a matrix.
The result is a contingency table in array format, whose dimensions are determined by  
the terms on the right side of the formula.

<<xtabs-ex1>>=
# 3-Way Frequency Table
mytable <- xtabs(~A+B+C, data=mydata)
ftable(mytable)    # print table
summary(mytable)   # chi-square test of indepedence
@

If a variable is included on the left side of the formula, it is assumed to be a
vector  of  frequencies  (useful if the data have already been tabulated in frequency form).

<<xtabs-ex2,results=verbatim>>=
(GSStab <- xtabs(count ~ sex + party, data=GSS))
summary(GSStab)
@

\subsection[Collapsing over factors]{Collapsing over table factors: \codefun{aggregate}, \codefun{margin.table} and \codefun{apply}}
It sometimes happens that we have a data set with more variables or factors than 
we want to analyse, or else, having done some initial analyses, we decide that
certain factors are not important, and so should be excluded from graphic displays
by collapsing (summing) over them.  For example, mosaic plots and fourfold displays
are often simpler to construct from versions of the data collapsed over
the factors which are not shown in the plots.

The appropriate tools to use again depend on
the form in which the data are represented--- a case-form data frame, a
frequency-form data frame (\codefun{aggregate}), or a table-form array or 
table object (\codefun{margin.table} or \codefun{apply}).

When the data are in frequency form, and we want to produce another
frequency data frame, \codefun{aggregate} is a handy tool, using
the argument \code{FUN=sum} to sum the frequency variable over the 
factors \emph{not} mentioned in the formula.

\Example 
The data frame \data{DaytonSurvey} in the \pkg{vcdExtra} package represents a
$2^5$ table giving the frequencies of reported use (``ever used?'') of 
alcohol, cigarettes and marijuana in a sample of high school seniors,
also classified by sex and race.  

<<dayton1, results=verbatim>>=
str(DaytonSurvey)
head(DaytonSurvey)
@

To focus on the associations among the
substances, we want to collapse over sex and race. The right-hand side of the formula
used in the call to \codefun{aggregate} gives the factors to be retained in the
new frequency data frame, \code{Dayton.ACM.df}.

<<dayton2, results=verbatim>>=
# data in frequency form
# collapse over sex and race
Dayton.ACM.df <- aggregate(Freq ~ cigarette+alcohol+marijuana, 
                           data=DaytonSurvey, FUN=sum)
Dayton.ACM.df
@

When the data are in table form, and we want to produce another
table, \codefun{apply} with \code{FUN=sum} can be used in a similar way
to sum the table over dimensions not mentioned in the \code{MARGIN}
argument.  \codefun{margin.table} is just a wrapper for \codefun{apply}
using the \codefun{sum} function.


\Example 
To illustrate, we first convert the \data{DaytonSurvey} to a 5-way
table using \codefun{xtabs}, giving \code{Dayton.tab}.  

<<dayton3, results=verbatim>>==
# in table form
Dayton.tab <- xtabs(Freq~cigarette+alcohol+marijuana+sex+race, data=DaytonSurvey)
structable(cigarette+alcohol+marijuana ~ sex+race, data=Dayton.tab)
@
Then, use \codefun{apply} on \code{Dayton.tab} to give the
3-way table \code{Dayton.ACM.tab} summed over sex and race.
The elements in this new table are the column sums for 
\code{Dayton.tab} shown by \codefun{structable} just above.

<<dayton4, results=verbatim>>==
# collapse over sex and race
Dayton.ACM.tab <- apply(Dayton.tab, MARGIN=1:3, FUN=sum)
Dayton.ACM.tab <- margin.table(Dayton.tab, 1:3)   # same result
structable(cigarette+alcohol ~ marijuana, data=Dayton.ACM.tab)
@



Many of these operations can be performed using the \verb|**ply()| functions
in the \pkg{plyr} package.
For example, with the data in a frequency form data frame, use \codefun{ddply}
to collapse over unmentioned factors, and \codefun{plyr::summarise}%
\footnote{
Ugh. This \pkg{plyr} function clashes with a function of the same name in \pkg{vcdExtra}.
In this document I will use the explicit double-colon notation to keep them
separate.
}
as the function to be applied to each piece.
<<dayton5, eval=FALSE>>==
Dayton.ACM.df <- ddply(DaytonSurvey, .(cigarette, alcohol, marijuana), 
                       plyr::summarise, Freq=sum(Freq))
@

\subsection[Collapsing levels]{Collapsing table levels: \codefun{collapse.table}}

A related problem arises when we have a table or array and for some purpose
we want to reduce the number of levels of some factors by summing subsets
of the frequencies.  For example, we may have initially coded Age in 10-year
intervals, and decide that, either for analysis or display purposes, we
want to reduce Age to 20-year intervals.  The \codefun{collapse.table} function
in \pkg{vcdExtra} was designed for this purpose.

\Example
Create a 3-way table, and collapse Age from 10-year to 20-year intervals. 
First, we generate a $2 \times 6 \times 3$ table of random counts from a 
Poisson distribution with mean of 100.
<<collapse1,results=verbatim>>=
# create some sample data in frequency form
sex <- c("Male", "Female")
age <- c("10-19", "20-29",  "30-39", "40-49", "50-59", "60-69")
education <- c("low", 'med', 'high')
data <- expand.grid(sex=sex, age=age, education=education)
counts <- rpois(36, 100)   # random Possion cell frequencies
data <- cbind(data, counts)

# make it into a 3-way table
t1 <- xtabs(counts ~ sex + age + education, data=data)
structable(t1)
@ 
Now collapse \code{age} to 20-year intervals, and \code{education}
to 2 levels. In the arguments, levels of \code{age} and \code{education}
given the same label are summed in the resulting smaller table.
<<collapse2, results=verbatim>>=
# collapse age to 3 levels, education to 2 levels
t2 <- collapse.table(t1, 
         age=c("10-29", "10-29",  "30-49", "30-49", "50-69", "50-69"),
         education=c("<high", "<high", "high"))
structable(t2)
@


\subsection[Converting]{Converting among frequency tables and data frames}

As we've seen, a given contingency table can be represented 
equivalently in different forms,
but some \proglang{R} functions were designed for one particular representation.
\tabref{tab:convert} shows some handy tools for converting from one form to another.

\begin{table}[htb]
 \caption{Tools for converting among different forms for categorical data}\label{tab:convert}
 \begin{center}
   \begin{tabular}{llll}
	\hline
                 & \multicolumn{3}{c}{\textbf{To this}} \\
	\textbf{From this}      &     Case form         & Frequency form             &  Table form \\
	\hline
	Case form      &   noop                 & \verb|xtabs(~A+B)|        &  \verb|table(A,B)|  \\ 
	Frequency form &  \verb|expand.dft(X)|  & noop                      & \verb|xtabs(count~A+B)|\\
	Table form     &  \verb|expand.dft(X)|  & \verb|as.data.frame(X)|   &  noop \\
	\hline
   \end{tabular}
 \end{center}
\end{table}

A contingency table in table form (an object of \code{class(table)}) can be converted
to a data.frame with \codefun{as.data.frame}.%
\footnote{
Because \proglang{R} is object-oriented, this is actually a short-hand for
the function \codefun{as.data.frame.table}.
}  
The resulting
\code{data.frame} contains columns
representing the classifying factors and the table entries (as a column named by
the \code{responseName} argument, defaulting to \code{Freq}. This is the inverse of \codefun{xtabs}.

\Example 
Convert the \code{GSStab} in table form to a data.frame in frequency form.
<<convert-ex1,results=verbatim>>=
as.data.frame(GSStab)
@

\Example Convert the \code{Arthritis} data in case form to a 3-way table of
\code{Treatment} $\times$ \code{Sex} $\times$ \code{Improved}.
Note the use of \codefun{with} to avoid having to use \code{Arthritis\$Treatment} etc. within the call to \codefun{table}.%
\footnote{
\codefun{table} does not allow a \code{data} argument to provide
an environment in which the table variables are to be found.  In the 
examples in \secref{sec:table} I used \code{attach(mydata)} for this purpose,
but \codefun{attach} leaves the variables in the global environment,
while \codefun{with} just evaluates the \codefun{table} expression in a
temporary environment of the data.
}
<<convert-ex2,results=verbatim>>=
Art.tab <-with(Arthritis, table(Treatment, Sex, Improved))
str(Art.tab)
ftable(Art.tab)
@


There may also be times that you will need an equivalent case form \code{data.frame}
with factors  representing the table variables
rather than the frequency  table.
For example, the \codefun{mca} function in package \pkg{MASS}
only operates on data in this format. 
Marc Schwartz provided code for \codefun{expand.dft} on the  Rhelp
mailing list for converting a table back into a case form \code{data.frame}.
This function is included in \pkg{vcdExtra}.

\Example Convert the \data{Arthritis} data in table form (\code{Art.tab}) back to a \code{data.frame}
in case form, with factors
\code{Treatment}, \code{Sex} and \code{Improved}.
<<convert-ex3,results=verbatim>>=
Art.df <- expand.dft(Art.tab)
str(Art.df)
@

\subsection{A complex example}\label{sec:complex}

If you've followed so far, you're ready for a more complicated example.
The data file, \code{tv.dat} represents a 4-way table of size
$5 \times 11 \times 5 \times 3$ where the table variables (unnamed in the file)
are read as \code{V1} -- \code{V4}, and the cell frequency is read
as \code{V5}.  The file, stored in the \code{doc/extdata} directory
of \pkg{vcdExtra}, can be read as follows:
<<tv1,results=verbatim>>=
tv.data<-read.table(system.file("doc","extdata","tv.dat",package="vcdExtra"))
head(tv.data,5)
@
For a local file, just use \codefun{read.table} in this form:
<<tv2,eval=FALSE>>=
tv.data<-read.table("C:/R/data/tv.dat")
@
The data \code{tv.dat} came from the initial implementation of 
mosaic displays in \proglang{R} by Jay Emerson. 
In turn, they came from the initial development of mosaic displays 
\citep{vcd:Hartigan+Kleiner:1984}
that illustrated the method with data on a large sample of TV viewers
whose behavior had been recorded for the Neilson ratings.
This data set contains sample television audience data from Neilsen
Media Research for the week starting November 6, 1995.

\begin{flushleft}
The table variables are:\\
~~~\code{V1}-- values 1:5 correspond to the days Monday--Friday;\\
~~~\code{V2}-- values 1:11 correspond to the quarter hour times 8:00PM through 10:30PM;\\
~~~\code{V3}-- values 1:5 correspond to ABC, CBS, NBC, Fox, and non-network choices;\\
~~~\code{V4}-- values 1:3 correspond to transition states: turn the television Off, Switch channels, 
 or Persist in viewing the current channel.
\end{flushleft}
	 
We are interested just the cell frequencies, and rely on the facts that the
(a) the table is complete--- there are no missing cells,
so \code{nrow(tv.data)}=\Sexpr{nrow(tv.data)};
(b) the observations are ordered so that \code{V1} varies most rapidly and
\code{V4} most slowly.  From this, we can just extract the frequency column
and reshape it into an array.
<<tv2,results=hide>>=
TV <- array(tv.data[,5], dim=c(5,11,5,3))                                        
dimnames(TV) <- list(c("Monday","Tuesday","Wednesday","Thursday","Friday"), 
                c("8:00","8:15","8:30","8:45","9:00","9:15","9:30",         
                  "9:45","10:00","10:15","10:30"),                            
                c("ABC","CBS","NBC","Fox","Other"), c("Off","Switch","Persist"))
names(dimnames(TV))<-c("Day", "Time", "Network", "State")
@

More generally (even if there are missing cells), we can 
use \codefun{xtabs} (or \codefun{plyr::daply})
to do the cross-tabulation, using \code{V5} as the
frequency variable.  Here's how to do this same operation with \codefun{xtabs}:
<<tv2a,eval=FALSE>>=
TV <- xtabs(V5 ~ ., data=tv.data)
dimnames(TV) <- list(Day=c("Monday","Tuesday","Wednesday","Thursday","Friday"), 
                Time=c("8:00","8:15","8:30","8:45","9:00","9:15","9:30",         
                       "9:45","10:00","10:15","10:30"),                            
                Network=c("ABC","CBS","NBC","Fox","Other"), 
                State=c("Off","Switch","Persist"))
@

But this 4-way table is too large and awkward to work with. Among the networks,
Fox and Other occur infrequently. 
We can also cut it down to a 3-way table by considering only viewers who persist
with the current station.%
\footnote{This relies on the fact that that indexing
an array drops dimensions of length 1 by default,
using the argument \code{drop=TRUE};
the result is coerced to the lowest possible dimension.
}

<<tv3,results=verbatim>>=
TV <- TV[,,1:3,]     # keep only ABC, CBS, NBC
TV <- TV[,,,3]       # keep only Persist -- now a 3 way table
structable(TV)
@

Finally, for some purposes, we might want to collapse the 11 times into a smaller number.
Here, we use \codefun{as.data.frame.table} to convert the table back to a data frame,
 \codefun{levels} to re-assign the values of \code{Time},
 and finally, \codefun{xtabs} to give a new, collapsed frequency table.

<<tv4,results=verbatim>>=
TV.df <- as.data.frame.table(TV)
levels(TV.df$Time) <- c(rep("8:00-8:59",4),rep("9:00-9:59",4), rep("10:00-10:44",3))
TV2 <- xtabs(Freq ~ Day + Time + Network, TV.df)
structable(Day ~ Time+Network,TV2)
@
Whew! See \figref{fig:TV-mosaic} for a mosaic plot of the \code{TV2} data.


\section{Tests of Independence}

\subsection{CrossTable}

OK, now we're ready to do some analyses.  For tabular displays,
the  \codefun{CrossTable}  function in  the \pkg{gmodels}  package produces  cross-tabulations
modeled after \code{PROC FREQ} in \proglang{SAS} or \code{CROSSTABS} in \proglang{SPSS}. 
It has a wealth of options for the quantities that can be shown in each cell.

<<xtabs-ex2, results=verbatim>>=
# 2-Way Cross Tabulation
library(gmodels)
CrossTable(GSStab,prop.t=FALSE,prop.r=FALSE,prop.c=FALSE)
@
There are  options to  report percentages  (row, column,  cell), specify decimal
places, produce Chi-square,  Fisher, and McNemar  tests of independence,  report
expected  and residual  values (pearson,  standardized, adjusted  standardized),
include missing values as valid, annotate with row and column titles, and format
as \proglang{SAS} or \proglang{SPSS} style output! See \code{help(CrossTable)} for details.


\subsection{Chi-square test}

For 2-way tables you can use \codefun{chisq.test} to test independence of the row
and column variable. By default,  the $p$-value is calculated from  the asymptotic
chi-squared distribution of the test  statistic. Optionally, the $p$-value can  be
derived via Monte Carlo simulation. 

<<chisq,results=verbatim>>=
(HairEye <- margin.table(HairEyeColor, c(1, 2)))
chisq.test(HairEye)
@

\subsection{Fisher Exact Test}\label{sec:Fisher}

\code{fisher.test(X)} provides an  exact test of  independence. \code{X} must be  a two-way
contingency table in table form.  Another form, 
\code{fisher.test(X, Y)} takes two
categorical vectors of the same length.  
For tables larger than $2 \times 2$ the method can be computationally intensive (or can fail) if
the frequencies are not small.

<<fisher,results=verbatim>>=
fisher.test(GSStab)
@

But this does not work because \data{HairEye} data has $n$=592 total frequency.
An exact test is unnecessary in this case.
<<echo=TRUE,eval=FALSE>>=
fisher.test(HairEye)
@
%# <<echo=TRUE,eval=FALSE,results=verbatim>>=
%# #cat(try(fisher.test(HairEye)))
%# @
\begin{Soutput}
Error in fisher.test(HairEye) : FEXACT error 6.
LDKEY is too small for this problem.
Try increasing the size of the workspace.
\end{Soutput}

\subsection[Mantel-Haenszel test]{Mantel-Haenszel test and conditional association}\label{sec:mantel}

Use the  \code{mantelhaen.test(X)} function  to perform  a Cochran-Mantel-Haenszel 
$\chi^2$ chi
test  of   the  null  hypothesis   that  two  nominal   variables  are
\emph{conditionally independent}, $A \perp B \given C$, in each stratum,  assuming that there is no  three-way
interaction. \code{X} is  a 3 dimensional  contingency table, where  the last dimension
refers to the strata.

The \data{UCBAdmissions} serves as an example of a $2 \times 2 \times 6$ table,
with \code{Dept} as the stratifying variable.
<<mantel1,results=verbatim>>=
## UC Berkeley Student Admissions
mantelhaen.test(UCBAdmissions)
@
The results show no evidence for association between admission and gender
when adjusted for department.  However, we can easily see that the assumption
of equal association across the strata (no 3-way association) is probably
violated. For $2 \times 2 \times k$ tables, this can be examimed
from the odds ratios for each $2 \times 2$ table (\codefun{oddsratio}), and
tested 
by  using
\verb|woolf_test()| in \pkg{vcd}.

%<<mantel2,results=verbatim>>=
%oddsRatio <- function(x) (x[1,1]*x[2,2])/(x[1,2]*x[2,1])
%apply(UCBAdmissions, 3, oddsRatio)
%
%woolf_test(UCBAdmissions) 
%@ 
<<mantel2,results=verbatim>>=
oddsratio(UCBAdmissions, log=FALSE)
lor <- oddsratio(UCBAdmissions)  # capture log odds ratios
summary(lor)
woolf_test(UCBAdmissions) 
@ 
We  can visualize the  odds ratios of  Admission for
each  department  with  fourfold  displays  using  \codefun{fourfold}.  The cell
frequencies $n_{ij}$  of each  $2 \times  2$ table  are shown  as a  quarter circle whose
radius is proportional to $\sqrt{n_{ij}}$, so  that its area is proportional to  the
cell frequency.
Confidence rings for the odds ratio allow a visual test of the null of no association; 
the rings for adjacent quadrants overlap \emph{iff} the observed counts are consistent 
with the null hypothesis.  In the extended version (the default), brighter colors
are used where the odds ratio is significantly different from 1.
The following lines produce \figref{fig:fourfold1}.%
\footnote{The color values \code{col[3:4]} were modified from their default values
to show a greater contrast between significant and insignifcant associations here.}
<<fourfold1,fig=TRUE,height=6,width=9,results=hide,include=FALSE>>=
col <- c("#99CCFF", "#6699CC", "#F9AFAF", "#6666A0", "#FF0000", "#000080")
fourfold(UCB,mfrow=c(2,3), color=col)
@
 
%\setkeys{Gin}{width=0.8\textwidth}
\begin{figure}[htb]
\begin{center}
%<<fourfold1,fig=TRUE,height=6,width=9,results=hide>>=
%col <- c("#99CCFF", "#6699CC", "#F9AFAF", "#6666A0", "#FF0000", "#000080")
%fourfold(UCB,mfrow=c(2,3), color=col)
%@
\includegraphics[width=0.8\textwidth,trim=80 50 80 50]{fig/vcd-tut-fourfold1}
\caption{Fourfold display for the \data{UCBAdmissions} data. Where the odds ratio differs
	significantly from 1.0, the confidence bands do not overlap, and the circle quadrants are
	shaded more intensely.}
\label{fig:fourfold1}
\end{center}
\end{figure}

Another \pkg{vcd} function, \codefun{cotabplot}, provides a more general approach
to visualizing conditional associations in contingency tables,
similar to trellis-like plots produced by \codefun{coplot} and lattice graphics.
The \code{panel} argument supplies a function used to render each conditional 
subtable. The following gives a display (not shown) similar to \figref{fig:fourfold1}.
<<fourfold2,eval=FALSE>>=
cotabplot(UCB, panel = cotab_fourfold)
@

When we want to view the conditional
probabilities of a response variable (e.g., \code{Admit})
in relation to several factors,
an alternative visualization is a \codefun{doubledecker} plot.
This plot is a specialized version of a mosaic plot, which
highlights the levels of a response variable (plotted vertically)
in relation to the factors (shown horizontally). The following
call produces \figref{fig:doubledecker}, where we use indexing
on the first factor (\code{Admit}) to make \code{Admitted}
the highlighted level.

In this plot, the
association between \code{Admit} and \code{Gender} is shown
where the heights of the highlighted conditional probabilities
do not align. The excess of females admitted in Dept A stands out here.

<<doubledecker,fig=TRUE, height=5,width=9,include=FALSE>>=
doubledecker(Admit ~ Dept + Gender, data=UCBAdmissions[2:1,,])
@

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig/vcd-tut-doubledecker}
\caption{Doubledecker display for the \data{UCBAdmissions} data. The heights
of the highlighted bars show the conditional probabilities of \texttt{Admit},
given \texttt{Dept} and \texttt{Gender}.}
\label{fig:doubledecker}
\end{center}
\end{figure}

Finally, the there is a \codefun{plot} method for \code{oddsratio} objects.
By default, it shows the 95\% confidence interval for the log odds ratio.
\figref{fig:oddsratio} is produced by:
<<oddsratio0,eval=FALSE>>=
plot(lor, xlab="Department", ylab="Log Odds Ratio (Admit | Gender)")
@

\setkeys{Gin}{width=0.5\textwidth}
\begin{figure}[htb]
\begin{center}
<<oddsratio,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot(lor, xlab="Department", ylab="Log Odds Ratio (Admit | Gender)")
@
\caption{Log odds ratio plot for the \data{UCBAdmissions} data.}
\label{fig:oddsratio}
\end{center}
\end{figure}

\subsection[CMH tests: ordinal factors]{Cochran-Mantel-Haenszel tests for ordinal factors}\label{sec:CMH}
The standard $\chi^2$ tests for association in a two-way table
treat both table factors as nominal (unordered) categories.
When one or both factors of a two-way table are
quantitative or ordinal, more powerful tests of association
may be obtaianed by taking ordinality into account, using
row and or column scores to test for linear trends or differences
in row or column means.

More general versions of the CMH tests (Landis etal., 1978) are provided by assigning
numeric scores to the row and/or column variables. 
For example, with two ordinal factors (assumed to be equally spaced), assigning
integer scores, \code{1:R} and \code{1:C} tests the linear $\times$ linear component
of association. This is statistically equivalent to the Pearson correlation between the
integer-scored table variables, with $\chi^2 = (n-1) r^2$, with only 1 $df$
rather than $(R-1)\times(C-1)$ for the test of general association.

When only one table
variable is ordinal, these general CMH tests are analogous to an ANOVA, testing
whether the row mean scores or column mean scores are equal, again consuming
fewer $df$ than the test of general association.

The \codefun{CMHtest} function in \pkg{vcdExtra} now calculates these various
CMH tests for two possibly ordered factors, optionally stratified other factor(s).

\Example
Recall the $4 \times 4$ table, \code{JobSat} introduced in \secref{sec:creating},
<<jobsat, results=verbatim>>=
JobSat
@

Treating the \code{satisfaction} levels as equally spaced, but using
midpoints of the \code{income} categories as row scores gives the following results:
<<cmh1, results=verbatim>>=
CMHtest(JobSat, rscores=c(7.5,20,32.5,60))
@

Note that with the relatively small cell frequencies, the test for general 
give no evidence for association. However, the the \code{cor} test for linear x linear
association on 1 df is nearly significant. The \pkg{coin} contains the
functions \verb|cmh_test()| and \verb|lbl_test()|
for CMH tests of general association and linear x linear association respectively.

\subsection{Measures of Association}

There are a variety of statistical measures of \emph{strength} of association for
contingency tables--- similar in spirit to $r$ or $r^2$ for continuous variables.
With a large sample size, even a small degree of association can show a 
significant $\chi^2$, as in the example below for the \data{GSS} data.

The  \codefun{assocstats}  function in \pkg{vcd}  calculates  the   $\phi$
contingency coefficient,  and Cramer's  V for  an $r \times c$  table. 
The input must be in table form, a two-way $r \times c$ table.  
It won't work with \data{GSS} in frequency form, but by now you should know how
to convert.
<<assoc1,results=verbatim>>=
assocstats(GSStab)
@

For tables with ordinal variables, like \data{JobSat}, some people prefer the
Goodman-Kruskal $\gamma$ statistic (\citet[\S 2.4.3]{vcd:Agresti:2002})
based on a comparison of concordant
and discordant pairs of observations in the case-form equivalent of a two-way table.
<<gamma,results=verbatim>>=
GKgamma(JobSat)
@

A web article by Richard Darlington,
\url{http://www.psych.cornell.edu/Darlington/crosstab/TABLE0.HTM}
gives further description of these and other measures of association.

\subsection{Measures of Agreement}
The
\codefun{Kappa} function in the \pkg{vcd} package calculates Cohen's $\kappa$ and weighted
$\kappa$ for a square two-way table with the same row and column categories \citep{Cohen:60}.%
\footnote{ 
Don't confuse this with \codefun{kappa} in base \proglang{R} that computes something
entirely different (the condition number of a matrix).
}
Normal-theory $z$-tests are obtained by dividing $\kappa$ by its asymptotic standard
error (ASE).  A \codefun{confint} method for \code{Kappa} objects provides confidence intervals.
<<kappa,results=verbatim>>=
(K <- Kappa(SexualFun))
confint(K)
@

A visualization of agreement, both unweighted and weighted for degree of departure
from exact agreement is provided by the \codefun{agreementplot} function.
\figref{fig:agreesex} shows the agreementplot for the \data{SexualFun} data,
produced as shown below. The Bangdiwala measures represent the proportion of the
shaded areas of the diagonal rectangles, using weights $w_1$ for exact agreement,
and $w_2$ for partial agreement one step from the main diagonal.
<<agreesex,fig=TRUE,height=6,width=7,results=verbatim,include=FALSE>>=
agree <- agreementplot(SexualFun, main="Is sex fun?")
unlist(agree)
@

%\setkeys{Gin}{width=0.5\textwidth}
\begin{figure}[htb]
\begin{center}
%<<agreesex,fig=TRUE,height=6,width=7,results=verbatim>>=
%agree <- agreementplot(SexualFun, main="Is sex fun?")
%agree
%@
\includegraphics[width=0.4\textwidth,trim=50 25 50 25]{fig/vcd-tut-agreesex}
\caption{Agreement plot for the \data{SexualFun} data.}
\label{fig:agreesex}
\end{center}
\end{figure}
In other examples, the agreement plot can help to show \emph{sources}
of disagreement.  For example, when the shaded boxes are above or below the diagonal
(red) line, a lack of exact agreement can be attributed in part to
different frequency of use of categories by the two raters-- lack of
\emph{marginal homogeneity}.
	
\subsection{Correspondence analysis}
Use the \pkg{ca} package for correspondence analysis for visually exploring relationships
between rows and columns in contingency tables.  For an $r \times c$ table,
the method provides a breakdown of the Pearson $\chi^2$ for association in up to $M = \min(r-1, c-1)$
dimensions, and finds scores for the row ($x_{im}$) and column ($y_{jm}$) categories
such that the observations have the maximum possible correlations.%
\footnote{
	Related methods are the non-parametric
	CMH tests using assumed row/column scores (\secref{sec:CMH}),
	the analogous \codefun{glm} model-based methods (\secref{sec:CMH}),
	and the more general RC models which can be fit using \codefun{gnm}.
	Correspondence analysis differs in that it is a primarily
	descriptive/exploratory method (no significance tests), but is
	directly tied to informative graphic displays of the row/column categories.
}	
	


Here, we carry out a simple correspondence analysis of the \data{HairEye} data.
The printed results show that nearly 99\% of the association between hair color and eye color
can be accounted for in 2 dimensions, of which the first dimension accounts for 90\%.
<<ca1,results=verbatim>>=
library(ca)
ca(HairEye)
@

The resulting \code{ca} object can be plotted just by running the \codefun{plot}
method on the \code{ca} object, giving the result in
\figref{fig:ca-haireye}.  \codefun{plot.ca} does not allow labels for dimensions;
these can be added with \codefun{title}.
It can be seen that most of the association is accounted for by the ordering
of both hair color and eye color along Dimension 1, a dark to light dimension.
<<ca-haireye0,echo=TRUE,eval=FALSE>>=
plot(ca(HairEye), main="Hair Color and Eye Color")
title(xlab="Dim 1 (89.4%)", ylab="Dim 2 (9.5%)")
@

\setkeys{Gin}{width=0.7\textwidth}
\begin{figure}[htb]
\begin{center}
<<ca-haireye,results=hide,fig=TRUE,echo=FALSE>>=
plot(ca(HairEye), main="Hair Color and Eye Color")
title(xlab="Dim 1 (89.4%)", ylab="Dim 2 (9.5%)")
@
\caption{Correspondence analysis plot for the \data{HairEye} data.}
\label{fig:ca-haireye}
\end{center}
\end{figure}

\section{Loglinear Models}\label{sec:loglin}

You can  use the  \codefun{loglm}  function in  the \pkg{MASS}  package to fit log-linear
models.  Equivalent models can also be fit (from a different perspective) as generalized
linear models with the \codefun{glm}  function using the \code{family='poisson'} argument,
and the \pkg{gnm} package provides a wider range of generalized \emph{nonlinear} models,
particularly for testing structured associations.
The visualization methods for these models were originally developed for models fit using \codefun{loglm},
so this approach is emphasized here.  Some extensions of these methods for models
fit using \codefun{glm} and \codefun{gnm} are contained in the \pkg{vcdExtra} package
and illustrated in \secref{sec:glm}.

Assume we  have a 3-way  contingency table based  on
variables A, B, and C.  
The possible different forms of \loglin\ models for a 3-way table are shown in \tabref{tab:loglin-3way}.
The \textbf{Model formula} column shows how to express each model for \codefun{loglm} in \proglang{R}.%
\footnote{
For \codefun{glm}, or \codefun{gnm}, with the data in the form of a frequency data.frame, the same model is specified
in the form \code{glm(Freq} $\sim$ \code{..., family="poisson")}, where \texttt{Freq} is  the name
of the cell frequency variable and \texttt{...} specifies the \textbf{Model formula}.
}
In the \textbf{Interpretation} column, the symbol ``$\perp$'' is to be read as ``is independent of,''
and ``$\given$'' means ``conditional on,'' or ``adjusting for,'' or just ``given''.

\begin{table}[htb]
 \caption{Log-linear Models for Three-Way Tables}\label{tab:loglin-3way}
 \begin{center}
 \begin{tabular}{llll}
  \hline
  \textbf{Model}           & \textbf{Model formula}  & \textbf{Symbol}& \textbf{Interpretation} \\
  \hline\hline 
  Mutual independence      & \verb|~A + B + C|       & $[A][B][C]$    & $A \perp B \perp C$ \\ 
  Joint independence       & \verb|~A*B + C|         & $[AB][C]$      & $(A \: B) \perp C$ \\ 
  Conditional independence & \verb|~(A+B)*C|         & $[AC][BC]$     & $(A \perp B) \given C$ \\ 
  All two-way associations & \verb|~A*B + A*C + B*C| & $[AB][AC][BC]$ & homogeneous association  \\ 
  Saturated model          & \verb|~A*B*C|           & $[ABC]$        & 3-way association \\ 
  \hline
 \end{tabular}
 \end{center}
\end{table}



For example, the formula \verb|~A + B + C| specifies the model of \emph{mutual independence} with
no associations among the three factors.  In standard notation for the expected frequencies
$m_{ijk}$, this corresponds to
\begin{equation*}
	\log ( m_{ijk} ) = \mu + \lambda_i^A + \lambda_j^B + \lambda_k^C \equiv \texttt{A + B + C}
\end{equation*}
The parameters $\lambda_i^A , \lambda_j^B$ and  $\lambda_k^C$ pertain to the differences among the
one-way marginal frequencies for the factors A, B and C. 

Similarly, the model of \emph{joint independence},  $(A \: B) \perp C$,
allows an association between A and B, but specifies that
C is independent of both of these and their combinations,
\begin{equation*}
	\log ( m_{ijk} ) = \mu + \lambda_i^A + \lambda_j^B + \lambda_k^C + \lambda_{ij}^{AB} \equiv \texttt{A * B + C}
\end{equation*}
where the parameters $\lambda_{ij}^{AB}$ pertain to the overall association between A and B (collapsing over C).

In the literature or text books, you will often find these models expressed in shorthand symbolic notation,
using brackets, \texttt{[ ]} to enclose the \emph{high-order terms} in the model.
Thus, the joint independence model can be denoted \texttt{[AB][C]}, as shown in the \textbf{Symbol} column in 
\tabref{tab:loglin-3way}.

Models of \emph{conditional independence} allow (and fit)  two of the three possible
two-way associations.  There are three such models, depending on which variable is conditioned upon.
For a given conditional independence model, e.g., \texttt{[AB][AC]}, the given variable is the one
common to all terms, so this example has the interpretation $(B \perp C) \given A$.


\subsection[Fitting with loglm()]{Fitting with \codefun{loglm}}\label{sec:loglm}
For example, we can fit the model of mutual independence among hair color, eye color and sex 
in \data{HairEyeColor} as

<<loglm-hec1,results=verbatim>>=
library(MASS)
## Independence model of hair and eye color and sex.  
hec.1 <- loglm(~Hair+Eye+Sex, data=HairEyeColor)
hec.1
@

Similarly, the models of conditional independence and joint independence are specified as
<<loglm-hec2,results=verbatim>>=
## Conditional independence
hec.2 <- loglm(~(Hair + Eye) * Sex, data=HairEyeColor)
hec.2
@
<<loglm-hec3,results=verbatim>>=
## Joint independence model.  
hec.3 <- loglm(~Hair*Eye + Sex, data=HairEyeColor)
hec.3
@
Note that printing the model gives a brief summary of the goodness of fit.
A set of models can be compared using the \codefun{anova} function.

<<loglm-anova,results=verbatim>>=
anova(hec.1, hec.2, hec.3)
@

%Martin Theus and Stephan Lauer have written an excellent article on  Visualizing
%Loglinear Models, using  mosaic plots. There  is also great  tutorial example by
%Kevin Quinn on analyzing loglinear models via glm.

\subsection[Fitting with glm() and gnm()]{Fitting with \codefun{glm} and \codefun{gnm}}\label{sec:glm}
The \codefun{glm} approach, and extensions of this in the \pkg{gnm} package allows a
much wider class of models for frequency data to be fit than can be handled by
\codefun{loglm}.  Of particular importance are models for ordinal factors and for
square tables, where we can test more structured hypotheses about the patterns
of association than are provided in the tests of general assosiation under
\codefun{loglm}. These are similar in spirit to the 
non-parametric CMH tests described in \secref{sec:CMH}.

\Example
The data \code{Mental} in the \pkg{vcdExtra} package gives a two-way table in frequency form
classifying young people by
their mental health status and parents' socioeconomic status (SES), where
both of these variables are ordered factors.
<<mental1,results=verbatim>>=
str(Mental)
xtabs(Freq ~ mental+ses, data=Mental)   # display the frequency table
@
Simple ways of handling ordinal variables involve assigning scores to the table
categories, and the simplest cases are to use integer scores, either for the row variable (``column
effects'' model), the column variable (``row effects'' model), or both (``uniform association''
model).

<<mental2>>=
indep <- glm(Freq ~ mental + ses, family = poisson, data = Mental)  # independence model
@
To fit more parsimonious models than general association, we can define
numeric scores for the row and column categories
<<mental3>>=
# Use integer scores for rows/cols 
Cscore <- as.numeric(Mental$ses)
Rscore <- as.numeric(Mental$mental)	
@
Then, the row effects model, the column effects model, and the uniform association
model can be fit as follows:
<<mental4,results=verbatim>>=
# column effects model (ses)
coleff <- glm(Freq ~ mental + ses + Rscore:ses, family = poisson, data = Mental)

# row effects model (mental)
roweff <- glm(Freq ~ mental + ses + mental:Cscore, family = poisson, data = Mental)

# linear x linear association
linlin <- glm(Freq ~ mental + ses + Rscore:Cscore, family = poisson, data = Mental)
@
The \codefun{Summarize} in \pkg{vcdExtra} provides a nice, compact summary of
the fit statistics for a set of models, collected into a \class{glmlist} object.
Smaller is better for AIC and BIC.
<<mental4a,results=verbatim>>=
# compare models using AIC, BIC, etc
vcdExtra::LRstats(glmlist(indep, roweff, coleff, linlin))
@
For specific model comparisons, we can also carry out tests of \emph{nested} models with
\codefun{anova} when those models are listed from smallest to largest.
Here, there are two separate paths from the most restrictive (independence) model
through the model of uniform association, to those that allow only one of
row effects or column effects.
<<mental5,results=verbatim>>=
anova(indep, linlin, coleff, test="Chisq")	
anova(indep, linlin, roweff, test="Chisq")	
@
The model of linear by linear association seems best on all accounts.
For comparison, one might try the CMH tests on these data:

<<mental6,results=verbatim>>=
CMHtest(xtabs(Freq~ses+mental, data=Mental))
@

\subsection{Non-linear terms}
The strength of the \pkg{gnm} package is that it handles a wide variety of models
that handle non-linear terms, where the parameters enter the model beyond a simple
linear function.
The simplest example is the Goodman RC(1) model, which allows a multiplicative
term to account for the association of the table variables.
In the notation of generalized linear models with a log link, this can be expressed as
\begin{equation*}
\log \mu_{ij} = \alpha_i + \beta_j + \gamma_{i} \delta_{j}	
\end{equation*}
where the row-multiplicative effect parameters $\gamma_i$ and 
corresponding column parameters $\delta_j$ are estimated from the data.%
\footnote{
This is similar in spirit to a correspondence analysis with a single
dimension, but as a statistical model.
}
Similarly, the RC(2) model adds two multiplicative terms to
the independence model,
\begin{equation*}
\log \mu_{ij} = \alpha_i + \beta_j + \gamma_{i1} \delta_{j1} + \gamma_{i2} \delta_{j2}
\end{equation*}

In the \pkg{gnm} package, these models may be fit using the \codefun{Mult}
to specify the multiplicative term, and \codefun{instances} to specify several
such terms.

\Example
For the \code{Mental} data, we fit the RC(1) and RC(2) models, and compare
these with the independence model.

<<mental7,results=verbatim>>=
RC1 <- gnm(Freq ~ mental + ses + Mult(mental,ses), data=Mental, 
             family=poisson, , verbose=FALSE)
RC2 <- gnm(Freq ~ mental+ses + instances(Mult(mental,ses),2), data=Mental, 
             family=poisson, verbose=FALSE)
anova(indep, RC1, RC2, test="Chisq")
@

\section{Mosaic plots}\label{sec:mosaic}

Mosaic plots provide an ideal method both for visualizing contingency tables and for 
visualizing the fit--- or more importantly--- lack of fit of a \loglin\ model.
For a two-way table, \codefun{mosaic} fits a model of independence, $[A][B]$
or \verb|~A+B| as an \proglang{R} formula.
For $n$-way tables, \codefun{mosaic} can fit any \loglin\ model, and can also be
used to plot a model fit with \codefun{loglm}.
See \citet{vcd:Friendly:1994,vcd:Friendly:1999} for the statistical ideas behind these
uses of mosaic displays in connection with \loglin\ models.

The essential idea is to recursively sub-divide a unit square into rectangular ``tiles'' for the
cells of the table, such
that the are area of each tile is proportional to the cell frequency.
For a given \loglin\ model, the tiles can then be shaded in various ways to reflect
the residuals (lack of fit) for a given model.  The pattern of residuals can then
be used to suggest a better model or understand \emph{where} a given model fits or
does not fit.

\codefun{mosaic} provides a wide range of options for the directions of splitting,
the specification of shading, labeling, spacing, legend and many other details.
It is actually implemented as a special case of a more general
class of displays for $n$-way tables called \code{strucplot}, including
sieve diagrams, association plots, double-decker plots as well as mosaic
plots.  For details, see \code{help(strucplot)} and the ``See also'' links,
and also \citet{vcd:Meyer+Zeileis+Hornik:2006b}, which is available as
an \proglang{R} vignette via \code{vignette("strucplot", package="vcd")}.

\figref{fig:arthritis}, showing the association between
\code{Treatment} and \code{Improved} was produced with the following
call to \codefun{mosaic}.
<<Arthritis1,echo=TRUE,eval=FALSE>>=
mosaic(art, gp = shading_max, split_vertical = TRUE, 
       main="Arthritis: [Treatment] [Improved]")
@
Note that the residuals for the independence model were not large
(as shown in the legend),
yet the association between \code{Treatment} and \code{Improved}
is highly significant.
<<art1,results=verbatim>>=
summary(art)
@
In contrast, one of the other shading schemes, from \citet{vcd:Friendly:1994}
(use: \verb|gp = shading_Friendly|), 
uses fixed cutoffs of $\pm 2, \pm 4$,
to shade cells which are \emph{individually} significant
at approximately $\alpha = 0.05$ and $\alpha = 0.001$ levels, respectively.
The right panel below uses \verb|gp = shading_Friendly|.

\setkeys{Gin}{width=0.5\textwidth}
<<art21,fig=TRUE,height=6,width=7,echo=FALSE>>=
mosaic(art, gp = shading_max, split_vertical = TRUE, 
       main="Arthritis: gp = shading_max")
@
<<art22,fig=TRUE,height=6,width=7,echo=FALSE>>=
mosaic(art, gp = shading_Friendly, split_vertical = TRUE, 
       main="Arthritis: gp = shading_Friendly")
@

\subsection[Mosaics for loglinear models]{Mosaics for \loglin\ models}\label{sec:mosaic-llm}

When you have fit a \loglin\ model using \codefun{loglm}, 
and saved the result (as a \code{loglm} object) the simplest way to display the
results is to use the \codefun{plot} method for the \code{loglm} object.
Calling \code{mosaic(loglm.object)} has the same result.
In \secref{sec:loglm} above, we fit several different models to the
\data{HairEyeColor} data.  We can produce mosaic displays of each just
by plotting them:

<<hec-mosaic,eval=FALSE>>=
# mosaic plots, using plot.loglm() method
plot(hec.1, main="model: [Hair][Eye][Sex]")
plot(hec.2, main="model: [HairSex][EyeSex]")
plot(hec.3, main="model: [HairEye][Sex]")
@

\setkeys{Gin}{width=0.32\textwidth}
<<hec1,fig=TRUE,height=6,width=7,echo=FALSE>>=
plot(hec.1, main="model: [Hair][Eye][Sex]")
@
<<hec2,fig=TRUE,height=6,width=7,echo=FALSE>>=
plot(hec.2, main="model: [HairSex][EyeSex]")
@
<<hec3,fig=TRUE,height=6,width=7,echo=FALSE>>=
plot(hec.3, main="model: [HairSex][EyeSex]")
@

Alternatively, you can supply the model formula to \codefun{mosaic}
with the \code{expected} argument.  This is passed to \codefun{loglm},
which fits the model, and returns residuals used for shading in the plot. 

For example, here we examine the \data{TV2} constructed in \secref{sec:complex}
above.  The goal is to see how Network choice depends on (varies with)
Day and Time. To do this: 
\begin{itemize}
	\item We fit a model of joint independence of
\code{Network} on the combinations of \code{Day} and \code{Time},
with the model formula \verb|~Day:Time + Network|.
  \item To make the display more easily read, we place \code{Day} and \code{Time}
on the vertical axis and \code{Network} on the horizontal,
  \item The \code{Time} values overlap on the right vertical axis, so we use
  \codefun{level} to abbreviate them.  \codefun{mosaic} also supports a 
  more sophisticated set of labeling functions.  Instead of changing the data
  table, we could have used 
  \verb|labeling_args = list(abbreviate = c(Time = 2))| for a similar effect.
\end{itemize}

The following call to \codefun{mosaic} produces \figref{fig:TV-mosaic}:
<<TV-mosaic0,eval=FALSE>>=
dimnames(TV2)$Time <- c("8", "9", "10")     # re-level for mosaic display
mosaic(~ Day + Network + Time, data=TV2, expected=~Day:Time + Network, 
         legend=FALSE, gp=shading_Friendly)
@

\setkeys{Gin}{width=0.75\textwidth}
\begin{figure}[htb]
\begin{center}
<<TV-mosaic,fig=TRUE,height=6,width=6,echo=FALSE>>=
dimnames(TV2)$Time <- c("8", "9", "10")     # re-level for mosaic display
mosaic(~ Day + Network + Time, data=TV2, expected=~Day:Time + Network, 
         legend=FALSE, gp=shading_Friendly)
@
\caption{Mosaic plot for the \data{TV} data
	showing model of joint independence, \texttt{Day:Time + Network} .}
\label{fig:TV-mosaic}
\end{center}
\end{figure}

From this, it is easy to read from the display how network choice varies with day and
time. For example, CBS dominates in all time slots on Monday;
ABC and NBC dominate on Tuesday, particularly in the later time slots;
Thursday is an NBC day, while on Friday, ABC gets the greatest share.

In interpreting this mosaic and other plots, it is important to understand that
associations included in the model---here, that between day and time---are \emph{not}
shown in the shading of the cells, because they have been fitted (taken into account)
in the \loglin\ model. 

For comparison, you might want to try fitting the model of homogeneous association.
This allows all pairs of factors to be associated, but asserts that each pairwise
association is the same across the levels of the remaining factor.
The resulting plot displays the contributions to a 3-way association, but is not shown here.
<<TV-mosaic1,eval=FALSE>>=
mosaic(~ Day + Network + Time, data=TV2, 
         expected=~Day:Time + Day:Network + Time:Network, 
         legend=FALSE, gp=shading_Friendly)
@
\subsection[Mosaics for glm() and gnm() models]{Mosaics for \codefun{glm} and \codefun{gnm} models}\label{sec:mosglm}
The \pkg{vcdExtra} package provides an additional method, \codefun{mosaic.glm} 
for models fit with \codefun{glm} and \codefun{gnm}.%
\footnote{
Models fit with \codefun{gnm} are of \code{class = c("gnm", "glm", "lm")},
so all \code{*.glm} methods apply, unless overridden in the \pkg{gnm} package.
}
These are not restricted to the 
Poisson family, but only apply to cases where the response variable is non-negative. 

\Example
Here, we plot the independence and the linear-by-linear association model
for the Mental health data from \secref{sec:glm}.
These examples illustrate some of the options for labeling (variable names and
residuals printed in cells).  Note that the \code{formula} supplied to \codefun{mosaic}
for \class{glm} objects refers to the order of factors displayed in the plot, not the model.
<<mental-plots,eval=FALSE>>=
long.labels <- list(set_varnames = c(mental="Mental Health Status", ses="Parent SES"))
mosaic(indep, ~ses+mental, residuals_type="rstandard",  
     labeling_args = long.labels, labeling=labeling_residuals,
     main="Mental health data: Independence")

mosaic(linlin, ~ses+mental, residuals_type="rstandard", 
     labeling_args = long.labels, labeling=labeling_residuals, suppress=1, 
     gp=shading_Friendly, main="Mental health data: Linear x Linear")
@

\setkeys{Gin}{width=0.49\textwidth}
<<mental-plots1,fig=TRUE,echo=FALSE,height=6,width=6>>=
long.labels <- list(set_varnames = c(mental="Mental Health Status", ses="Parent SES"))
mosaic(indep, ~ses+mental, residuals_type="rstandard",  
     labeling_args = long.labels, labeling=labeling_residuals,
     main="Mental health data: Independence")
@
<<mental-plots2,fig=TRUE,echo=FALSE,height=6,width=6>>=
long.labels <- list(set_varnames = c(mental="Mental Health Status", ses="Parent SES"))
mosaic(linlin, ~ses+mental, residuals_type="rstandard", 
     labeling_args = long.labels, labeling=labeling_residuals, suppress=1, 
     gp=shading_Friendly, main="Mental health data: Linear x Linear")
@

The \pkg{gnm} package also fits a wide variety of models with nonlinear terms or terms for
structured associations of table variables.  In the following, we fit the RC(1)
model
\begin{equation*}
	\log ( m_{ij} ) = \mu + \lambda_i^A + \lambda_j^B + \phi \mu_i \nu_j 
\end{equation*}
This is similar to the linear by linear model, except that the row effect
parameters ($\mu_i$) and column parameters ($\nu_j$) are estimated from the data
rather than given assigned equally-spaced values.  The multiplicative terms
are specified by the \codefun{Mult}.

<<mental-RC1, eval=FALSE>>=
Mental$mental <- C(Mental$mental, treatment)
Mental$ses <- C(Mental$ses, treatment)
RC1model <- gnm(Freq ~ mental + ses + Mult(mental, ses),
                family = poisson, data = Mental)
mosaic(RC1model, residuals_type="rstandard", labeling_args = long.labels, 
       labeling=labeling_residuals, suppress=1, gp=shading_Friendly,
       main="Mental health data: RC(1) model")
@
Other forms of nonlinear terms are provided for the inverse of a predictor
(\codefun{Inv}) and the exponential of a predictor (\codefun{Exp}).
You should read \code{vignette("gnmOverview", package="gnm")} for further details.

\subsection{Mosaic tips and techniques}\label{sec:tips}

The \pkg{vcd} package implements an extremely general collection of graphical methods for
$n$-way frequency tables within the strucplot framework, which includes mosaic plots
(\codefun{mosaic}), as well as association plots (\codefun{assoc}), sieve diagrams
(\codefun{sieve}), as well as tabular displays (\codefun{structable}).

The graphical methods in \pkg{vcd}  support a wide of options that control almost all of the details
of the plots, but it is often difficult to determine what arguments you need to supply
to achieve a given effect from the \code{help()}.  As a first step, you should read the 
\code{vignette("strucplot")} in \pkg{vcd} to understand the overall structure of
these plot methods.  The notes below describe a few useful things that may not be obvious,
or can be done in different ways.

\subsubsection[Changing labels]{Changing the labels for variables and levels}
With data in contingency table form or as a frequency data frame, it often happens
that the variable names and/or the level values of the factors, while suitable
for analysis, are less than adequate when used in mosaic plots and other strucplot
displays.  

For example, we might prefer that a variable named \code{ses} appear
as \code{"Socioeconomic Status"}, or a factor with levels \code{c("M", "F")}
be labeled using \code{c("Male", "Female")} in a plot.  Or, sometimes we
start with a factor whose levels are fully spelled out
(e.g., \code{c("strongly disagree", "disagree", "neutral", "agree", "strongly agree")}),
only to find that the level labels overlap in graphic displays.

The structplot framework in \pkg{vcd} provides an extremely large variety of
functions and options for controlling almost all details of text labels in mosaics
and other plots.  See \code{help(labelings)} for an overview. 

For example, in \secref{sec:ordered-factors} we showed how to rearrange the dimensions
of the \code{UCBAdmissions} table, change the names of the table variables, and
relabel the levels of one of the table variables.
The code below changes the actual table for plotting purposes, but we pointed out that
these changes can create other problems in analysis.

<<reorder3a,results=hide>>=
UCB <- aperm(UCBAdmissions, c(2, 1, 3))
names(dimnames(UCB)) <- c("Sex", "Admit?", "Department")
dimnames(UCB)[[2]] <- c("Yes", "No")
@

The same effects can be achieved \emph{without} modifying the data
using the \verb|set_varnames| and
\verb|set_labels| options in \codefun{mosaic} as follows:

<<reorder3b,results=hide, keep.source=TRUE>>=
vnames <- list(set_varnames = c(Admit="Admission", Gender="Sex", Dept="Department"))

lnames <- list(Admit = c("Yes", "No"),
              Gender = c("Males", "Females"),
              Dept = LETTERS[1:6])

mosaic(UCBAdmissions, labeling_args=vnames, set_labels=lnames)
@

In some cases, it may be sufficient to abbreviate (or clip, or rotate) level names to avoid 
overlap.  For example, the statements below produce another version of
\figref{fig:TV-mosaic} with days of the week abbreviated to their first three letters.
Section 4 in the \code{vignette("strucplot")} provides many other examples.

<<TV-mosaic2,results=hide>>=
dimnames(TV2)$Time <- c("8", "9", "10")     # re-level for mosaic display
mosaic(~ Day + Network + Time, data=TV2, expected=~Day:Time + Network, 
         legend=FALSE, gp=shading_Friendly, 
         labeling_args=list(abbreviate=c(Day=3)) )
@

%\subsubsection{Fitting complex models with glm() and gnm()}



\section[Continuous predictors]{Continuous predictors}\label{sec:contin}
When continuous predictors are available---and potentially important--- in explaining a
categorical outcome, models for that outcome include:
logistic regression (binary response), 
the proportional odds model (ordered polytomous response),
multinomial (generalized) logistic regression.
Many of these are special cases of the generalized linear model using the
\code{"poisson"} or \code{"binomial"} family and their relatives.

\subsection{Spine and conditional density plots}\label{sec:spine}
I don't go into fitting such models here, but I would be remiss not to illustrate some
visualizations in \pkg{vcd} that are helpful here.
The first of these is the spine plot or spinogram \citep{vcd:Hummel:1996}
 (produced with \codefun{spine}).
These are special cases of mosaic plots with 
specific spacing and shading to show how a categorical response varies with
a continuous or categorical predictor.

They are also a generalization of stacked bar plots where not the heights  but
the \emph{widths} of the bars corresponds to the relative frequencies of \code{x}. The heights
of the  bars then  correspond to  the conditional  relative frequencies  of {y} in
every \code{x} group.

\Example
For the \data{Arthritis} data, we can see how \code{Improved} varies with \code{Age}
as follows.  \codefun{spine} takes a formula of the form 
\verb|y ~ x| with a single dependent factor and a single explanatory variable \code{x}
(a numeric variable or a factor).
The range of a numeric variable\code{x} is divided into intervals based on the
\code{breaks} argument, and stacked bars are drawn to show the distribution of
\code{y} as \code{x} varies.  As shown below, the discrete table that is visualized
is returned by the function.

<<spine1,fig=FALSE,results=verbatim>>=
(spine(Improved ~ Age, data = Arthritis, breaks = 3))
(spine(Improved ~ Age, data = Arthritis, breaks = "Scott"))
@

\setkeys{Gin}{width=0.49\textwidth}
<<spine2,fig=TRUE,echo=FALSE,height=6,width=6>>=
(spine(Improved ~ Age, data = Arthritis, breaks = 3))
@
<<spine3,fig=TRUE,echo=FALSE,height=6,width=6>>=
(spine(Improved ~ Age, data = Arthritis, breaks = "Scott"))
@


The conditional density plot \citep{vcd:Hofmann+Theus} is a further generalization.
This visualization technique is similar to spinograms, but uses a smoothing approach
rather than discretizing the explanatory variable.  As well, it uses 
the original \code{x} axis and not a distorted one.

\setkeys{Gin}{width=0.6\textwidth}
\begin{figure}[htb]
\begin{center}
<<cdplot,fig=TRUE,height=6,width=6>>=
cdplot(Improved ~ Age, data = Arthritis)
with(Arthritis, rug(jitter(Age), col="white", quiet=TRUE))
@
\caption{Conditional density plot for the \data{Arthritis} data
	showing the variation of Improved with Age.}
\label{fig:cd-plot}
\end{center}
\end{figure}

In such plots, it is useful to also see the distribution of the observations
across the horizontal axis, e.g., with a \codefun{rug} plot.
\figref{fig:cd-plot} uses \codefun{cdplot} from the \pkg{graphics} package
rather than \verb|cd_plot()| from \pkg{vcd}, and is produced with
<<cdplot1,echo=TRUE>>=
cdplot(Improved ~ Age, data = Arthritis)
with(Arthritis, rug(jitter(Age), col="white", quiet=TRUE))
@
From \figref{fig:cd-plot} it can be easily seen that the proportion
of patients reporting Some or Marked improvement increases with Age,
but there are some peculiar bumps in the distribution.
These may be real or artifactual, but they would be hard to see
with most other visualization methods.
When we switch from non-parametric data exploration to parametric
statistical models, such effects are easily missed.

\subsection[Model-based plots]{Model-based plots: effect plots and \pkg{ggplot2} plots}\label{sec:modelplots}
The nonparametric conditional density plot uses smoothing methods
to convey the distributions of the response variable,
but displays that are simpler to interpret can often be obtained by plotting
the predicted response from a parametric model.

For complex \codefun{glm} models with interaction effects, the \pkg{effects}
package  provides the most useful displays,
plotting  the predicted values for a given term, averaging over other
predictors not included in that term.  I don't illustrate this here, but
see \citet{effects:1,effects:2} and \code{help(package="effects")}.

Here I just briefly illustrate the capabilities of the \pkg{ggplot2} package
for model-smoothed plots of categorical responses in \codefun{glm} models.

\Example
The \data{Donner} data frame in \pkg{vcdExtra} gives details on the survival
of 90 members of the Donner party,
a group of people who attempted to migrate to California in 1846.
They were trapped by an early blizzard on the eastern side of the
Sierra Nevada mountains, and before they could be rescued, 
nearly half of the party had died.
What factors affected who lived and who died?

<<donner1, results=verbatim>>=
data(Donner, package="vcdExtra")
str(Donner)
@
A potential model of interest
is the logistic regression model for $Pr(survived)$, allowing separate
fits for males and females as a function of \code{age}.
The key to this is the \verb|stat_smooth()| function, using
	\code{method = "glm", family = binomial}. The \verb|formula = y ~ x|
	specifies a linear fit on the logit scale (\figref{fig:donner3}, left)

<<donner2a, fig=FALSE, eval=FALSE>>=
# separate linear fits on age for M/F
ggplot(Donner, aes(age, survived, color = sex)) +
  geom_point(position = position_jitter(height = 0.02, width = 0)) +
  stat_smooth(method = "glm", family = binomial, formula = y ~ x,
           alpha = 0.2, size=2, aes(fill = sex))
@
Alternatively, we can allow a quadratic relation with \code{age}
by specifying \verb|formula = y ~ poly(x,2)| (\figref{fig:donner3}, right).
<<donner2b, fig=FALSE, eval=FALSE>>=
# separate quadratics
ggplot(Donner, aes(age, survived, color = sex)) +
  geom_point(position = position_jitter(height = 0.02, width = 0)) +
  stat_smooth(method = "glm", family = binomial, formula = y ~ poly(x,2),
           alpha = 0.2, size=2, aes(fill = sex))
@

\setkeys{Gin}{width=0.49\textwidth}
\begin{figure}[htb]
\begin{center}
<<donner3a,fig=TRUE,echo=FALSE,height=6,width=6, results=hide>>=
ggplot(Donner, aes(age, survived, color = sex)) +
  geom_point(position = position_jitter(height = 0.02, width = 0)) +
  stat_smooth(method = "glm", family = binomial, formula = y ~ x,
           alpha = 0.2, size=2, aes(fill = sex))

@
<<donner3b,fig=TRUE,echo=FALSE,height=6,width=6, results=hide>>=
# separate quadratics
ggplot(Donner, aes(age, survived, color = sex)) +
  geom_point(position = position_jitter(height = 0.02, width = 0)) +
  stat_smooth(method = "glm", family = binomial, formula = y ~ poly(x,2),
           alpha = 0.2, size=2, aes(fill = sex))
@
\caption{Logistic regression plots for the  \data{Donner} data
	showing survival vs. age, by sex. Left: linear logistic model; right: quadratic model}
\label{fig:donner3}
\end{center}
\end{figure}
These plots very nicely show (a) the fitted $Pr(survived)$ for males and females;
(b) confidence bands around the smoothed model fits and (c) the individual
observations by jittered points at 0 and 1 for those who died and survided, respectively.

\bibliography{vcd,vcdExtra}

\end{document}
