\documentstyle[11pt]{article}
\renewcommand{\baselinestretch}{2.1}
\textwidth 6.75in
\textheight 8.50in
\oddsidemargin 0.0in
\evensidemargin 0.0in
\topmargin -0.25in


\begin{document}

\begin{center}
\bf \large 
A Simple Hazards Model for Ordered Categorical Outcomes 
\end{center}


\begin{center}
Philip J. Smith,\footnote{Address for correspondence:
Centers for Disease Control and Prevention,
Division of Diabetes Translation,
4770 Buford Highway N.E., Mail Stop K-10,
Atlanta, GA 30341 U.S.A.}
Alula Hadgu, \& Susan J. Kenny
\end{center}

\begin{center}
Centers for Disease Control and Prevention, U.S.A.
%Division of Diabetes Translation
%\end{center}
%\begin{center}
%4770 Buford Highway N.E., Mail Stop K-10,
%Atlanta, GA 30341
\end{center}

\begin{center}
\today
\end{center}


%\vskip .5in
\section*{Summary}

A simple hazards model belonging to the class
of generalized linear models
is proposed for the analysis of 
ordered categorical data. 
%Our finding contrasts with the notion conveyed by
%Aitkin et al.~(1989, p.240--242)
%that models for ordered categorical data
%do not belong to the class of generalized linear models.
In comparison to
McCullagh's~(1980) multivariate methods
and William's and Grizzle's~(1972) weighted least squares
approach
for analyzing ordered categorical data,
our method can use continuous covariates as predictors
and can be easily fitted to data
using existing software 
implementing binary regression with
logistic or complementary log--log links.
Informal diagnostic methods for assessing goodness of fit are
sketched
and a test for assessing the proportional hazards assumption is
given
when predictors are categorical.
Data from the Wisconsin Epidemiologic Study of Diabetic
Retinopathy (Klein, Klein, Moss, et al. 1989)
and the 
Survey of Physician Practice Behaviors conducted by the
National Institute of Diabetes and Digestive and Kidney Disease
are used to illustrate the methods.

\vskip .5in

\noindent
{\sc Key words:}
Binomial kernel,
complementary log--log link,
deviance residuals,
%disease progression,
logistic link,
proportional hazards model,
generalized linear models,
goodness of link test,
hazard function,
iteratively reweighted least squares.



\section{Introduction}


An ordinal categorical outcome 
is defined by an underlying continuous
or unobservable latent random variable $A$
that has been discretized into adjacent 
intervals.
Let $a \sb 0 < a \sb 1 < \ldots < a \sb K $ denote defined ``cut
points''
of the underlying continuous variable
and let $I \sb j = ( a \sb {j-1} , a \sb j ]$
denote the $j$th ordered category of the discretized variable,
$j=1, \ldots , K.$
For generality, let $a \sb 0 = -\infty$
and $a \sb K = +\infty.$
%
Finally, let 
$\{ d \sb j \}$ denote counts observed within 
ordered categories $\{ I \sb j \}$ 
from a random sample of size $n = \sum d \sb j .$

The objective of ordered categorical data analysis
is to summarize how the  distribution of the counts $\{ d \sb j
\}$ 
observed across the ordered categories
depend upon 
cut point ``effects'' $\{ \theta \sb j  \}$ 
that are functions
of the underlying 
continuous  variable's unknown distribution 
at their associated cut points $\{ a \sb j \}.$ 
Further, when the distribution of the counts 
across the ordered categories 
is believed to depend upon predictor variables
the objective of ordered categorical data analysis 
is to model how the relationship between the counts
and the cut point effects
depends upon the predictor variables.

As a first example,
Table~1
lists sampled distributions of U.S. physicians' agreement
on whether educational programs 
should be provided for persons
with insulin dependent diabetes melitus.
Each responding physician registered agreement
on an ordered scale that ranged
from strongly agree (SA), agree (A), neutral (N), disagree (D),
to strongly
disagree (SD).
Because education of diabetic persons is an objective
of the Year 2000 Goals of the U.S. Public Health Service
(1990),
an analysis of physician's current attitudes 
toward the objective
provides information
on persons who provide primary care to persons with diabetes
and insight on problems the Service
needs to anticipate in pursuing its goal.
Our analysis sought to assess the extent to which 
the distribution across the attitude scale
depends upon medical specialty and physician age.




As a second example, 
Table~2 gives data from the Wisconsin 
Epidemiological Study of Diabetic Retinopathy (WESDR).
This table lists the empirical distribution 
across an ordered scale measuring the 
degree of diabetic retinopathy 
of persons at the first 4 year follow up of the scale.
Counts are given for persons who at the inception of the 
longitudinal study had high or low glycemic control
as measured by glycosylated hemoglobin.
The purpose our analysis was to model 
progression toward proliferative diabetic retinopathy (PDR)
and
death 
on the continuous variate,
duration of diabetes, 
after
adjusting for baseline glycemic control.

The notion of accounting for ordering in analyzing
a categorical response variable
goes back at least a century
and includes investigations by
Pearson~(1901), Yule~(1912),
Yates~(1948),
Stuart~(1953),
Goodman and Kruskal~(1954),
and Williams~(1952).
Perhaps the most widely used method for analyzing
data from ordered categorical responses is
based on Armitage's~(1955) method of
assigning increasing scores to the ordered categories
and comparing the average scores for each 
categorical predictor variable.
Agresti~(1988)  provided other methods 
of using scores based on Poisson regression.
However, Graubard and Korn~(1987)
showed that inference may be sensitive to
the choice of scores.
Williams and Grizzle~(1972)  formulated a 
proportional odds model for accounting for an ordered categorical
outcome within the weighted least squares (WLS) approach. 
However, Zeger~(1988)  pointed out that when
the frequencies within categories are small 
estimates of the covariance weighting matrix are unstable
and WLS methods are inefficient. Further, the WLS approach
cannot account for continuous predictors.
%
More recently, McCullagh~(1980) 
applied a
{\sl multivariate} generalized linear model 
to obtain both proportional hazards and proportional odds models.
These methods seem to be specially tailored 
for applications where predictor variables are 
restricted to be categorical.
Anderson~(1984) 
provided further insight for these models.
Other methods for modeling ordered categorical data include those
described by
Thompson and Baker~(1981)
who used composite link functions
and
Hastie, Botha, and Schnitzler~(1986) who used
nonparametric smoothing methods.

This paper presents a simple method for analyzing
ordered categorical data
when covariate predictors are either categorical or continuous
within the class of generalized linear models.
Pregibon~(1980) and Tosteson and Begg~(1988)
have noted the absence of methods for ordinal categorical data
using 
continuous predictors.
The  approach we describe is
related to methods for
analyzing grouped survival data
discussed by
Bartlett~(1978),
Prentice and Gloeckler~(1978),
and Aitkin, Anderson, Francis, and Hinde~(1989).
In contrast to McCullagh's methods implemented 
in the computer package PLUM~(1988)
the statistical methods we propose 
can be 
modeled easily using many statistical software program
packages including GLIM~(1986) and S~(1992).

%Prentice and Gloeckler \cite{Pre_Glo}
%nd Aitkin, Anderson, Francis, and Hinde \cite{AAFH}
%have provide discussions of the analysis of
%grouped survival data.

Although the proportional odds model for 
analyzing ordered categorical data has grown in popularity,
McCullagh has pointed out 
that there is little reason
{\sl a priori} to believe that the proportional odds model should
fit ordered categorical data better than the proportional hazards
model does.
His example using family income data 
in which the proportional hazards provides a better fit to the
data reinforces this claim well.
In comparison to the proportional odds model,
the proportional hazards model has the added advantage that  
statistical methods including tests 
and refined residual analyses are available for
assessing the validity of the proportional hazards assumption.
Although it is easy to conceive of informal graphical methods
for assessing goodness of fit for the proportional odds model,
we know of no formal test for assessing the validity of the
link function in this model. 
Our paper presents goodness of link tests
as a starting point for evaluating
both the complementary log-log and logit links
for analyzing ordered categorical data.

In Section~2, the discrete proportional hazards model for
analyzing ordered categorical outcomes is described
and computational methods are sketched.
In Section~3, we review informal graphical diagnostic checks for
evaluating the fit of the model and give a goodness of link test 
for assessing the validity of the proportional hazards
assumption.
Sections~4 and 5 illustrate the methods 
using data
from the Survey of Physician Practice Behaviors
and the WESRD.


\section{The discrete proportional hazards model}

Let
$f(\cdot),$ 
$S(\cdot),$
and
$h(\cdot)$
denote 
the density function,
the complement of the cumulative distribution function,
and
the hazard function of the underlying continuous random variable
$A.$
Given the ordered categories $\{ I \sb j \}$ defined in
Section~1,
we define the probability mass function on the ordered
categories by
\begin{eqnarray*}
f \sb j =
P(A \in I \sb j ) =
s \sb j - s \sb {j+1}
\end{eqnarray*}
\noindent
where
\begin{eqnarray*}
s \sb j = P(A > a \sb {j-1}) = \sum \sb {k=j} \sp K f \sb k
\end{eqnarray*}
\noindent
is the probability of exceeding the $j-1$st cut point, and
\begin{eqnarray*}
h \sb j = P (A \in I \sb j | A > a \sb {j-1} ) = f \sb j / s \sb
j 
\end{eqnarray*}
\noindent
is the hazard for interval $I \sb j .$


Assuming that the $i$th sampled unit's underlying continuous
outcome $A$ 
belongs to interval $I \sb {i \sb j},$ 
the contribution to the likelihood of the $i$th sampled unit is
\begin{eqnarray*}
L \sb i
=
f \sb {i \sb j}
=
h \sb {i \sb j} s \sb {i \sb j}
=
h \sb {i \sb j}  \Pi \sb {j=1} \sp {{i \sb j}-1} (1 - h \sb j).
\end{eqnarray*}

Letting
$
w \sb {ij} =
     1 $ if $A \in I \sb j$ for unit $i$
     and
     $ = 0, $ otherwise,
the likelihood over all sampled units is
\begin{eqnarray}
L 
& = &
\Pi \sb {i=1} \sp n L \sb i
\nonumber
\\
& = &
\Pi \sb {i=1} \sp n 
     \Pi \sb {j=1} \sp {i \sb j}
     \left[
     h \sb j  \sp {w \sb {ij} }
      (1 - h \sb j) \sp {1- w \sb {ij}}
     \right].
\label{eq:like0}
\end{eqnarray}
\noindent
More generally, letting hazard for interval $I \sb j$
depend upon characteristics of th $i$th unit (\ref{eq:like0})
is
\begin{eqnarray}
L 
=
\Pi \sb {i=1} \sp n 
     \Pi \sb {j=1} \sp {i \sb j}
     \left[
     h \sb {ij}  \sp {w \sb {ij} }
      (1 - h \sb {ij}) \sp {1- w \sb {ij}}
     \right].
\label{eq:like1}
\end{eqnarray}

Letting $x \sb i = ( x \sb {i1} , \ldots , x \sb {ip} )$
denote a vector of $p$ continuous and/or categorical
covariates for the $i$th unit,
for the Cox's~(1972) proportional hazards model,
\begin{eqnarray*}
h(a, x\sb i )
=
\lambda (a) e \sp { \beta \sp T x\sb i}
\end{eqnarray*}
\noindent
and
\begin{eqnarray*}
S(a, x\sb i )
=
\exp ( - \Lambda (a)  e \sp { \beta \sp T x\sb i} )
\end{eqnarray*}
\noindent
where
$  \Lambda (a) = \int \sb {-\infty} \sp a \lambda ( u ) du $
is the integrated hazard function.

Then,
\begin{eqnarray*}
s \sb {ij} 
& = &
\exp ( - \Lambda ( a \sb {j-1} ) \: e \sp { \beta \sp T x\sb i}
),
\\
h \sb {ij} 
& = &
1 -
\exp \left[ -  e \sp { \beta \sp T x\sb i} \:
     \{
     \Lambda ( a \sb {j} ) 
     -
     \Lambda ( a \sb {j-1} ) 
     \}
      \right]
\\
& = &
1 -
\exp \left[  -  e \sp { \theta \sb j + \beta \sp T x\sb i}
\right]
\end{eqnarray*}
where
$
\theta \sb j
= 
\log 
\left\{ 
     \Lambda ( a \sb {j} ) 
     -
     \Lambda ( a \sb {j-1} ) 
\right\} .
$

To obtain maximum likelihood estimates of $\beta$ 
and the cut point effects $\{ a \sb j \}$
note that (\ref{eq:like1})
is the product of $\sum i \sb j$ Bernoulli densities.
In the parlance of generalized linear models,
the random component of
each density is $ w \sb {ij}$
having a Bernoulli distribution
with
``success'' probability $h \sb {ij}. $
The 
link function $g ( \cdot ) $ connecting
the linear predictor
$ \theta \sb j + \beta \sp T x\sb i $
to the mean of $E ( w \sb {ij} ) = h \sb {ij}$
is the complementary log--log link
\begin{eqnarray}
\label{eta}
\log \{ - \log ( 1 - h \sb {ij} ) \}
= 
g( h \sb {ij} )
= 
\theta \sb j + \beta \sp T x\sb i .
\end{eqnarray}
\noindent
To obtain inferences on the parameters $\beta$
the $ w \sb {ij}$ may be treated as 
independent random variables. 
The iteratively reweighted least squares algorithm
implemented in many standard packages
such as GLIM and S may be used to obtain
maximum likelihood estimates, standard errors,
and deviance tests.

As McCullagh~(1980) has pointed out, 
subject matter considerations may dictate
may dictate the choice of an alternative link,
such as the logit, for the analysis. 
%However, there is often little reason to 
%prefer the complementary log--log link to the logistic link,
%{\sl a priori.}
Subsequent sections discuss the issue of link selection.


For the important case when the objective of the analysis 
is to compare the distribution of
counts across each of $s$ subpopulations
(\ref{eq:like1}) can be reexpressed as
\begin{eqnarray}
\label{eq:like3}
L
=
\Pi \sb {i=1} \sp s
\Pi \sb {j=1} \sp K
\left[
h \sb {ij} \sp {d \sb {ij}}
(1 - h \sb {ij} ) \sp { n \sb {ij} - d \sb {ij} }
\right]
\end{eqnarray}
\noindent
where
$ d \sb {ij} $
is the number of sampled units from the $i$th subpopulation 
belonging to intervals $I \sb j;$
$n \sb {ij}$
is the number of  sampled units from the $i$th subpopulation
belonging to interval 
$I \sb j, $
$I \sb {j+1}, \ldots, $
and
$I \sb n ;$
and
$x\sb i$
is a indicator vector of population membership
and interaction terms. Section~4 gives an example illustrating
this case.


\section{Informal diagnostics and a  test for proportional
hazards}

\subsection{Informal diagnostics}

Perhaps the most fundamental diagnostic check 
for assessing systematic departures from any model
involves inspecting residuals. 
For the binomial kernels in (\ref{eq:like3})
letting $\hat h \sb {ij}$ denote the fitted mean of $d \sb {ij}$
deviance residuals are defined by
\begin{eqnarray}
\label{dres}
e \sb {ij} 
=
\mbox{sign}
(d \sb {ij} - \hat  h  \sb {ij} )
\left[
2
\left\{
d \sb {ij} \log ( d \sb {ij} / \hat  h  \sb {ij} )
+
( n \sb {ij} - d \sb {ij} )
\log
\left[
( n \sb {ij} - d \sb {ij} )
/
(n \sb {ij} - \hat  h  \sb {ij} )
\right]
\right\}
\right]
\sp {1/2}
\end{eqnarray}
\noindent
$j = 1 , \ldots , i \sb j,$
$i = 1 , \ldots ,K.$
Provided the model is correctly specified
the deviance residuals are normally distributed and
a normal probability plot of the deviance residuals
will provide a preliminary assessment of goodness of fit,
as Pierce and Schafer~(1986) showed.
A further check for systematic departures from the model
can be obtained by plotting
the fitted linear predictors (\ref{eta}) 
versus the deviance residuals (\ref{dres}).
A null pattern in this plot
provides evidence supporting goodness of fit.
However, 
systematic trend in this plot may indicate
the wrong choice of link function
or  scale for an explanatory variable,
or point to a missing quadratic term in the linear predictor.
Finally,
a plot of the absolute value of the deviance residuals
against fitted values $\hat h \sb {ij}$
will provide a check on the adequacy of the assumed
variance function:
the null pattern shows no trend
but a poorly chosen variance function will.
In particular,
a positive trend provides evidence of underdispersion
and
a negative trend provides evidence of overdispersion.
To correct for either effect,
a more appropriate variance function may be selected
and quasi-likelihood methods employed.
Alternatively, compensation for departures from nominal
dispersion
can be made using the methods described by
Smith and Heitjan~(1993).
Methods for checking for other isolated departures 
in generalized linear models are described by
Pregibon~(1981),
Williams~(1987),
and
McCullagh and Nelder~(1989).

Kalbfleisch and Prentice~(1980)
described informal diagnostic methods 
designed specifically
for assessing 
the proportional hazards assumption.
In particular, 
if the proportional hazards model holds
then
the effects of 
predictive covariates on
$
\log 
\left[
-\log 
\left\{ 
S(a, x \sb i)
\right\}
\right]
$
are additive
and plots of 
$
\log 
\left[
-\log 
\left\{ 
d \sb {ij} / n \sb {ij}
\right\}
\right]
$
versus the cut points $\{ a \sb j \}$
should appear to be roughly parallel
for each subpopulation
$i=1 , \ldots , s.$
Aitkin and Clayton~(1980)
describe other informal graphical methods
for assessing the proportional hazards assumption.


\subsection{A test for proportional hazards}

In survival analysis 
Schoenfeld~(1980),
Harrell~(1986),
and Gill and Schumacher~(1986)
describe formal tests for assessing 
the validity of the proportional hazards assumption.
For ordered categorical data as would be obtained from
(\ref{eq:like3})
we present a new and simple test for assessing
the assumption.
This test examines the appropriateness of the
complementary log--log link 
that yields the proportional hazards model 
with binomial data.
A convenient advantage to our method is that
it can be obtained as a part of the usual fitting procedure 
in statistical software packages such as GLIM and S that
implement methods for generalized linear models.
%
In particular,
consider the one parameter link function family 
discussed by Pregibon~(1980):
\begin{eqnarray}
\label{eq:link}
g( h \sb {ij} , \alpha ) 
=
\log
\left[
     \left.
\left\{
     \left(
     \frac{1}{1-h \sb {ij}}
     \right)
     \sp \alpha
     -
     1
\right\}
     \right/
     \alpha
\right]  .
\end{eqnarray}
\noindent
When $\alpha \rightarrow 0$ (\ref{eq:link}) 
corresponds to the 
complementary log--log link 
and 
when $\alpha = 1 $ (\ref{eq:link}) 
corresponds to the logistic link.
Letting 
$\alpha \sb *$ 
and
$\alpha \sb 0 = 0$
denote the ``true'' 
and hypothesized values for $\alpha ,$
Pregibon shows that 
the correct link can be approximated by
\begin{eqnarray}
\label{eq:true}
\left.
g ( h \sb {ij} , \alpha \sb * ) 
\doteq
g ( h \sb {ij} , \alpha \sb 0 )
+
(\alpha \sb * - \alpha \sb 0 )
\frac{\partial}{\partial \alpha}
g( h \sb {ij} , \alpha )
\right| \sb {\alpha = \alpha \sb 0} .
\end{eqnarray}
\noindent
Using (\ref{eq:true}) the ``true'' linear predictor
can be approximated by
\begin{eqnarray*}
g( h \sb {ij} , \alpha \sb 0 )
=
\theta \sb j + \beta \sp T x \sb i + \gamma z \sb {ij}
\end{eqnarray*}
\noindent
where 
$\gamma = \alpha \sb * - \alpha \sb 0$
is a parameter to be estimated 
\begin{eqnarray}
\label{eq:test}
z \sb {ij}
& = &
\lim \sb {\alpha \rightarrow 0}
\frac{\partial}{\partial \alpha}
g( \hat h \sb {ij} , \alpha )
\nonumber
\\
& = &
-
\frac 1 2
\log ( 1 - \hat h \sb {ij} )
\end{eqnarray}
\noindent
where $\{ \hat h \sb {ij} \}$ are the fitted hazard probabilities
obtained from the binary regression of
$\{ d \sb {ij} \}$ on $\{ x \sb i \}$
using the complementary log--log link.
Thus, 
the problem of assessing the proportional hazards assumption
can be formulated as a goodness of link test.
Within this framework
the proportional hazards assumption 
can be tested simply
by assessing the significance of the additional covariate
$z \sb {ij}$ 
using a deviance test. 

Similarly, to assess the adequacy of the logit link 
a deviance test may be used to assess the significance of  
$z \sb {ij} = \log ( \hat h \sb {ij} / ( 1 - \hat h \sb {ij} ))$
in a model using the logit link
where $\{ \hat h \sb {ij} \}$ are the fitted hazard probabilities
obtained from the binary regression of
$\{ d \sb {ij} \}$ on $\{ x \sb i \}$
using the same link.



\section{Analysis of physician attitudes}

Due to the paucity of data 
among persons who strongly disagreed 
to the notion of providing educational programs to insulin
dependent 
persons with diabetes, 
the adjacent D and SD categories in Table~1  were combined.
The purpose of the analysis was 
to assess how the distribution 
of sampled counts observed across 
the ordered agreement scale (CAT with 4 levels) 
 depended upon medical specialty (SPEC with 4 levels) 
and physician age (AGE with 3 levels).
Table~3 shows the analysis 
of deviance table for our analysis. 

Models 1-3 in Table~3
indicate that  medical specialty 
is a significant predictor of the distribution across the ordered
agreement scale: 
removal of the interaction between specialty and category
in model 2
results in the deviance increasing by 19.04
on 9 degrees of freedom corresponding to a p-value of 0.025.
Subsequent addition of AGE yields model 4  
indicating that age is a significant predictor.
However, the interaction of CAT and AGE in model  5 is not
significant.


Figure~1 shows the fitted probabilities obtained from model 6.
This figure shows that pediatricians are more  likely to strongly
agree to the notion of education for insulin dependent diabetics.
Also, this figure shows that younger physicians are more prone to
strongly agree than older physicians, regardless of specialty
type.
Figure 2 shows a normal probability plot of the deviance
residuals from model 6.
This plot 
provides informal evidence that the residuals behave as one would
expect from a sample of normally distributed variates: 
the tails of the empirical distribution are somewhat shorter than

normal
as would be expected from a small to moderate sample size.
Figures 3 and 4 provide further graphical diagnostics using
Cleveland's~(1979) lowess smoother to
aid in identifying trends.
Lack of a  trend in the plot of the deviance residuals versus the
fitted linear predictor in Figure 3 indicates an adequate
selection of terms in the linear predictor.
Also, lack of a  trend in the plot of the absolute values of the
deviance residuals in Figure 4 indicates that the data is
nominally dispersed according to the 
binomial distribution assumed in (\ref{eq:like3}).

To assess the proportional hazards assumption, Figure 5 shows 
a plot of the empirical values of
$\log ( -\log ( h \sb {ij}  ) )$  versus the ordered attitude
scale.
Because differences between the empirical transformed hazards are
roughly
constant for each medical specialty, 
the plot provides evidence that effects 
are additive on the transformed scale 
and that the proportional hazards assumption is not unreasonable.

Model~7 in Table~3 provides 
a more formal assessment of the goodness of the choice 
of the complementary log--log link and the proportional hazards
model. 
In particular, the 
change in deviance resulting from adding (\ref{eq:test}) 
to the model is -1 on 1 degree of freedom corresponding 
to a p-value of 0.31. 
This provides 
information that the proportional hazards 
model is adequate for the physician attitude data.




\section{The diabetic retinopathy data}

The purpose of the analysis of the diabetic retinopathy 
data was to build
statistical models
for the purpose of predicting the distribution across the
diabetic retinopathy
scale given in Table~2 for persons who have high glycemic
control and for persons who manage their glycemic levels poorly.
Since the distribution across the retinopathy scale (RET) is 
known be strongly associated with duration of diabetes (DUR), 
a continuous variate, 
and glycemic control (GLY), these factors were used in the model.


Table~4 gives the analysis of deviance table for the WESDR
data. For numerous analyses that we have conducted, 
the deviance statistics have been {\sl very} close for the logit
and complementary log--log link. For the analysis of the WESDR
data
deviance statistics were identical to {\sl at least} one digit to
the
right of the decimal point. This empirical evidence
suggests that in there may be very little information from which
to
choose between the logistic and complementary log--log link
in many applied situations. 
In related problems,
Tosteson and Begg~(1988) indicate that an appropriate selection
between these links is not the most crucial issue in conducting
a statistical analysis of ordered categorical data.
McCullagh~(1980) recommends that the choice of a link function
to be guided by issues specific to the nature of the application.

For the WESDR data our analysis indicates that the first order
interactions
of RET with GLY and DUR are highly significant, indicating that
the 
distribution across the retinopathy scale depends upon both
glycemic control and duration of diabetes. In particular, 
fitted values from this model indicate that for a given duration,
persons with high glycemic control are less likely to progress to
PDR and death. Also, for either glycemic group, longer durations
were associated with a greater likelihood of progression toward
PDR and death.  The second order interaction of DUR.GLY.RET was
not significant ($p=0.20$) indicating that these factors do not
necessarily
act synergistically on either the logit or complementary log--log
scale.


\section*{References}

\medskip \noindent
Agresti, A. (1988).
{\sl Analysis of Ordered Categorical Data.}
New York: John Wiley and Sons.

\medskip \noindent
Aitkin, M., Anderson, D., Francis, B., and Hinde, J. (1989)
{\sl Statistical Modelling in GLIM.}

\medskip \noindent
Aitkin, M. and Clayton, D. (1980).
Fitting the exponential, Weibull and extreme value distributions
to complex
survival data using GLIM.
{\sl Journal of the Royal Statistical Society, Series C,}
{\bf 29}, 156-163.

\medskip \noindent
Anderson, J.A. (1984).
Regression and ordered categorical variables (with discussion).
{\sl Journal of the Royal Statistical Society, Series B,}
{\bf 46}, 1-30.

\medskip \noindent
Armitage, P. (1955).
Tests for linear trends in proportions and frequencies.
{\sl Biometrics,}
{\bf 11}, 375-386.

\medskip \noindent
Bartlett, N.R. (1978).
A survival model for a wood preservation trial.
{\sl Biometrics,}
{\bf 34}, 673-679.

\medskip \noindent
Chambers, J.M. and Hastie, T.J. (1992).
{\sl Statistical Models in S,}
Wadsworth and Brooks: Pacific Grove, California.

\medskip \noindent
Cox, D.R. (1972).
Regression models and life tables (with discussion).
{\sl Journal of the Royal Statistical Society, Series B,}
{\bf 34}, 187-220.

\medskip \noindent
Cleveland, W. (1979).
Robust locally weighted regression and smoothing scatterplots.
{\sl Journal of the American Statistical Association,}
{\bf 74}, 829-836.

\medskip \noindent
Department of Health and Human Services. (1990).
{\sl Healthy People 2000:
National Health Promotion and Disease Prevention Objectives.}
U.S. Government Printing Office: Washington, D.C.
Oxford: Clarendon Press.

\medskip \noindent
Gill, R. and Schumacher, M.
A simple test of the proportional hazards assumption.
{\sl Biometrika,}
{\bf 74}, 289-300.

\medskip \noindent
Goodman, L.A. and Kruskal, W.H. (1954).
Measures of association for cross classifications, Part I.
{\sl Journal of the American Statistical Association,}
{\bf 49}, 732.

\medskip \noindent
Graubard, B.I. and Korn, E.L. (1987).
Choice of column scores for testing independence in ordered 
$2 \times k$ contingency tables.
{\sl Biometrics,}
{\sl 43}, 471-476.

\medskip \noindent
Harrell, F. (1986).
The PHGLM procedure,
In
{\sl SUGI Supplemental Library User's Guide, Version 5 Edition,}
437-466.

\medskip \noindent
Hastie, T.J., Botha, J.L., and Schnitzler, C.M. (1986).
Regression with an ordered categorical response.
{\sl Statistics in Medicine,}
{\bf 8}, 785-794.

\medskip \noindent
Kalbfleisch, J.D. and Prentice, R.L. (1980).
{\sl The Analysis of Failure Time Data.}
John Wiley and Sons: New York.

\medskip \noindent
Klein, R., Klein, E.K., Moss. S.E., Davis, D.D., DeMets, D.L.
(1989).
The Wisconsin Epidemiologic Study of Diabetic Retinopathy,
{\sl Archives of Ophthamology,}
{\bf 107}, 224-249.

\medskip \noindent
McCullagh, P. (1988).
PLUM: an interactive computer package for analysis of ordinal
data.
Technical report:
Department of Statistics, University of Chicago.


\medskip \noindent
McCullagh, P. (1980).
Regression models for ordinal data (with discussion).
{\sl Journal of the Royal Statistical Society, Series B,}
{\bf 42}, 109-142.

\medskip \noindent
McCullagh, P. and Nelder, J.A. (1989).
{\sl Generalized Linear Models.}
Chapman and Hall: New York.

%\medskip \noindent
%Liang, K.Y. and Zeger, S.L. (1986).
%Longitudinal data analysis using generalized linear models.
% {\sl Biometrika,}
%{\bf 73}, 13-22.


\medskip \noindent
Payne, C.D. (ed.) (1986).
{\sl The GLIM Manual, Release 3.77.}
Numerical Algorithms Group: Oxford.

\medskip \noindent
Pearson, K. (1901). 
Mathematical contributions to the theory of evolution.
{\sl Philosophical Transactions of the Royal Society,}
{\bf 195}, 79-150.

\medskip \noindent
Pierce, D.A. and Schafer, D.W. (1986).
Residuals in generalized linear models.
{\sl Journal of the American Statistical Association,}
{\bf 81}, 997-983.

\medskip \noindent
Pregibon, D. (1980). Discussion commentary of 
``Regression models for ordinal data'' by Peter McCullagh,
{\sl Journal of the Royal Statistical Society, Series B,}
{\bf 42}, 138--139.

\medskip \noindent
Pregibon, D. (1980).
Goodness of link tests for generalized linear models.
{\sl Journal of the Royal Statistical Society, Series C,}
{\bf 29}, 15-24.

\medskip \noindent
Pregibon, D. (1981).
Logistic regression diagnostics.
{\sl Annals of Statistics,}
{\bf 9}, 705-724.


\medskip \noindent
Prentice, R.L. and Gloeckler, L.A. (1978).
Regression analysis of grouped survival data with application
to breast cancer data.
{\sl Biometrics,}
{\bf 34}, 57-67.


\medskip \noindent
Schoenfeld, D. (1980).
Chi-squared goodness-of-fit tests for the proportional hazards
regression model.
{\sl Biometrika,}
{\bf 67}, 145-153.

\medskip \noindent
Smith, P.J. and Heitjan, D. (1993).
Testing and adjusting for overdispersion in 
generalized linear models.
{\sl Journal of the Royal Statistical Society, Series C,}
{\bf 42}, 31-42.

\medskip \noindent
Stuart, A. (1953).
The estimation and comparison of strengths of association in
contingency
tables.
{\sl Biometrika,}
{\bf 40},
105.

\medskip \noindent
Tosteson, A.N. and Begg, C.B. (1988).
A general regression methodology for ROC curve estimation.
{\sl Medical Decision Making,}
{\bf 8}(3), 204--215.

\medskip \noindent
Thompson, R. and Baker, R.J. (1981).
Composite link functions in generalized linear models.
{\sl Journal of the Royal Statistical Society, Series C,}
{\bf 30}, 125-131.

\medskip \noindent
Yates, F. (1948). 
The analysis of contingency tables with groupings based on
quantitative characters.
{\sl Biometrika,}
{\bf 35}, 176.

\medskip \noindent
Yule, G.U. (1912).
On methods of measuring association between two attributes
(with discussion).
{\sl Journal of the Royal Statistical Society,}
{\bf 75}, 579-652.

\medskip \noindent
Williams, D.A. (1987).
Generalized linear model diagnostics using the deviance and
single-case deletions.
{\sl Journal of the Royal Statistical Society, Series C,}
{\bf 36}, 181-191.

\medskip \noindent
Williams, E.J. (1952).
Use of scores for the analysis of association in contingency
tables.
{\sl Biometrika,}
{\bf 39}, 274.

\medskip \noindent
Williams, O.D. and Grizzle, J.E. (1972).
Analysis of contingency tables having ordered response
categories.
{\sl Journal of the American Statistical Association,}
{\bf 67}, 55-63.

\medskip \noindent
Zeger, S.L. (1988).
Commentary.
{\sl Statistics in Medicine,}
{\bf 7}, 161-168.

\eject

\begin{table}[h]
%\label{tab:data}
\caption
{Distributions of physicians attitudes}
\begin{center}
\begin{tabular}{ l c  r r r r r  }
            &     & \multicolumn{5}{c}{agreement scale} \\
\cline{3-7} 
medical specialty & age & SA & A & N & D & SD            \\ 
\hline\hline
Pediatricians & $\le$ 39 & 47   &     45   &      8   &      4  
&      0   \\
           & 40-54   &75   &     73   &     14   &      4   &    

1   \\
           & 55+     &27   &     36   &      9   &      3   &    

0   \\
\hline
Internists    & $\le$ 39 &31   &     67   &     26   &      7   &

    0   \\
           & 40-54   &45   &     67   &     46   &     20   &    

1   \\
           & 55+     &26   &     37   &     32   &     13   &    

2   \\
\hline
Family Practitioners & $\le$ 39 &45   &     71   &     31   &    
23   &      0   \\
             & 40-54   &28   &     62   &     22   &     21   &  

  1   \\
             & 55+     &18   &     66   &     33   &     22   &  

  2   \\
\hline
General Practitioners & $\le$ 39 & 9   &     19   &     12   &   

 5   &      0   \\
                & 40-54   &13   &     43   &      5   &     13  
&      0   \\
                & 55+     &28   &     68   &     31   &     26  
&      1   \\
\hline \hline
\end{tabular}
\end{center}
\end{table}

\eject

\begin{table}[h]
%\label{tab:WESDR}
\caption
{Distributions of diabetic retinopathy scores}
\begin{center}
\begin{tabular}{ c c c c c c  }
            & \multicolumn{5}{c}{retinopathy scale}  \\
\cline{2-6} 
glycemic control & none & mild & moderate & PDR & dead  \\  
\hline\hline
high             &   82 &   72 &        9 &  30 &   10  \\
low              &   17 &   64 &       28 &  60 &   16  \\
\hline \hline
\end{tabular}
\end{center}
\end{table}

\eject
\vskip 4in
\begin{table}[h]
%\label{tab:anodev}
\caption
{Analysis of deviance table for the physicians attitude analysis}
\begin{center}
\begin{tabular}{ l c  r r r r r  }
\multicolumn{1}{c}{model} & 
added/deleted terms & 
deviance & 
d.f. & 
$\Delta(\mbox{deviance})$ & 
$\Delta(\mbox{d.f.})$ & 
p-value            \\  \hline\hline
1: CAT + SPEC + CAT.SPEC &  CAT*SPEC    & 38.12 & 32 &        &  

&       \\
2: CAT + SPEC            & -CAT.SPEC    & 57.16 & 41 & +19.04 &
+9 & 0.025 \\ 
3: CAT + SPEC + CAT.SPEC &  +CAT.SPEC   & 38.12 & 32 &        &  

&       \\
4: CAT*SPEC + AGE        &    +AGE      & 27.74 & 30 & -10.38 &
-2 & $<0.01$ \\ 
5: CAT*(SPEC+AGE)        &    +CAT.AGE  & 23.02 & 24 &  -4.72 &
-6 &   0.58    \\ 
6: CAT*SPEC + AGE        &    -AGE      & 27.74 & 30 &        &  

&     \\ 
7: CAT*SPEC + AGE + z    &  +z          & 26.74 & 29 & -1.00  &
-1 & 0.31 \\
\hline \hline 
\end{tabular}
\end{center}
\end{table}
\begin{table}[h]
%\label{tab:andevW}
\caption
{Analysis of deviance table for the WESDR data analysis}
\begin{center}
\begin{tabular}{ l c  r r r r r  }
\multicolumn{1}{c}{model} & 
added/deleted terms & 
deviance & 
d.f. & 
$\Delta(\mbox{deviance})$ & 
$\Delta(\mbox{d.f.})$ & 
p-value            \\  \hline\hline
1: RET                   &                 & 1133.1 & 967 &      

&    &           \\
2: RET * GLY             & + GLY + GLY.RET & 1065.2 & 962 & -67.9

& -9 & $\ll 0.01$\\ 
3: RET * (GLY + DUR )    & + DUR + DUR.RET &  864.2 & 957 &-201.0

& -5 & $\ll 0.01$\\
\hline \hline 
\end{tabular}
\end{center}
\end{table}



\end{document}
