\section*{Logistic Regression---Again}


\subsection*{The Challenger Data}

You are probably all familiar with the Challenger Space Shuttle
accident of January 20, 1986.  I will leave out all the interesting
substantive details about the accident and focus instead on the data
about O-rings that was available at the time of the accident.

Well, I do need a little background, so here goes.  The space shuttle
uses two booster rockets and each rocket is made of several pieces
whose joints are sealed by rubber O-rings.   The O-rings are supposed
to stop the hot gases that are inside the booster from escaping.  Each
booster contains 3 primary O-rings.   There was data about the O-rings
on 23 previous flights of space shuttles.  Here is that data.
\begin{verbatim}
Temp Fail          Temp Fail
53   2             70   1
57   1             70   1
58   1             72   0
63   1             73   0
66   0             75   0
67   0             75   2
67   0             76   0
67   0             76   0
68   0             78   0
69   0             79   0
70   0             81   0
70   0
\end{verbatim}

On the day of the fated flight, the forecast launch temperature was
$31^\circ$F.  The O-rings were known to be sensitive to temperature.
The previous coldest launch temperature was $53^\circ$.   There was
discussion about the cold temperature before launch, and the engineers
did look at the data, but {\em only\/} the data for the failures.
They produced a plot that looked something like this.

\begin{verbatim}
> chal <- read.table("chal.dat", head=T)
> attach(chal, pos=1)
> ls()
[1] ".Last.value" "Fail"        "Temp"       
> n <- rep(6, length(Fail))

> plot(Temp[Fail >0], Fail[Fail>0])
\end{verbatim}

\centerline{\psfig{figure=chal1.eps,angle=-90,width=3in}}

Based on the erroneous assumption that the flights with no failures
did not add any information, they noted that there did not appear to
be any relationship between failures and temperature.

Had they produced a plot of the full data, which follows, they might
have behaved differently.

\begin{verbatim}
> plot(Temp, Fail)
\end{verbatim}

\centerline{\psfig{figure=chal2.eps,angle=-90,width=3in}}

In addition before each flight the O-rings were given a pressure
test.  The pressure at which that test was conducted increased over
time, and it is possible that the higher pressure test could have
caused more failures.  That data, and a lot more discussion, is
available in the paper {\em Risk Analysis of the Space Shuttle:
  Pre-{\em Challenger\/} Prediction of Failure}, by S.R~Dalal,
E.B.~Fowlkes, and B.~Hoadley, in V.~84, pp.945--957, of the Journal of
the American Statistical Association (JASA).   {\bf You} should be
able to read most of that article.  It is well written and easy to
follow.




\subsubsection*{How to look at these data}

Obviously, we are most interested how probability of failure relates
to temperature, and other quantities known before a flight.  Given our
data, it is not immediately obvious what response variable we should
be looking at.  For example, should we look at the number of failures
of O-rings in a flight, or just the fact that an O-ring failed?  These
might lead to different results (but fortunately for us they can be
looked at using the same analyses).  

We have in the past looked at models for count data that were of the
form $$
\log{m} = \mbox{some combination of factor effects} $$

Although these loglinear models did not require discrete (categorical)
predictors, the only examples we saw were of that form.  The current
data can also be modelled using the standard loglinear model, but an
equivalent and more natural approach is to set up the model as a
logistic regression.

In this case, allow  $p(t)$ to be the probability of failure given the
temperature $t$.  If one wanted to be really careful, one might write
this as
$$
  p(F|t) = \mbox{probability of failure given temperature was } t
$$
but the shorthand notation will probably suffice.

A first approach is to consider that the failures are a binomial
process with the probability of failure dependent on $t$, but with the
events themselves independent.  Thus using the first form of the data
(where we explicitly model the number of O-rings that failed) we are
assuming that
\begin{itemize}
\item The 6 O-rings on any flight are fail independently (but with a
  common probability), and
\item The 23 flights are independent.
\end{itemize}

The second assumption is more tenable than the first, but even the
first seems quite reasonable.  In order to fret about the first
assumption, we will also analyze a version of the data that just
records whether there was a failure or not.


\subsection*{The Logistic model}

Clearly, we can not get away with modelling
$$
  p(t) = \mbox{some arbitrary function of } t $$
unless we are very careful about the arbitrary function.  One would
quickly find oneself in the position where the models produced
negative probabilities, or probabilities greater than 1.  

One might think of the model (based on our loglinear experience)
$$
\log{p(t)} = \mbox{some arbitrary function of } t $$
but it suffers from the same problems.  

Lets go back to first principles and think about the data we are
trying to model.

We believe (at least I am suggesting) that the data can be modelled as
follows.  Let $n$ be the number of trials that are done for each run
of the experiment (in the Challenger case this is either 6 if each
O-ring is considered or 1 if just overall success/failure is
considered) and let $r$ be the number of failures in those $n$
trials.  Let $R$ be the corresponding random variable.  Then we
believe that
\begin{eqnarray*}
P(R=r | t)& =& {\binom n r} p(t)^r (1-p(t))^{n-r}  \\
          & = & {\binom n r} {\left ( \frac{p(t)}{1-p(t)} \right )}^r (1-p(t))^{n} \\
          & = & \mbox{Fn of r }\times\mbox{Fn of $p$ and $r$ } \times
                 \mbox{Fn of only $p$}
\end{eqnarray*}

Likelihood theory (and common sense) suggests that the only thing of
any importance here is the part that depends on both $r$ and $p$.  In
other words, if you want to maximize the likelihood only the middle
bit counts.

If we go one step further, and consider the log-likelihood 
$$
\log{{\binom n r}} + r \times \log{\left ( \frac{p(t)}{1-p(t)}
  \right )} + n \times \log{(1-p(t))} $$

then we see that modelling might be easy and natural (and the
asymptotic statistics might work) if we form models that are linear in
$$
\log{\left ( \frac{p(t)}{1-p(t)}\right )} $$

That is exactly what we will do.

The function 
$$\log{\left ( \frac{p}{1-p}\right )} $$
is called the
logistic transformation.  If 
$$ x= \log{\left ( \frac{p}{1-p}\right )}$$ 
then the inverse function is
$$ p = \frac{\exp(x)}{1 + \exp(x)} $$

We will (obviously, from the name) model the logit of the
probabilities as a function of the predictors.  


\subsection*{Back to the Challenger Data}

How do we want to model the challenger data?   We (with the current
data) only have temperature as a predictor, so we are stuck with
models of the form,
$$
\mbox{logit}(p(t)) = \mbox{some function of} t $$

The obvious first step is to plot the data (how can we do that?), then
fit some simple models, then think about transformations.   In this
example we are particularly interested in extrapolations (dangerous
though that is!) to temperatures that are less than we have seen.


\subsubsection*{Plotting the data}

This is harder than one might imagine.  One (mostly) can not just plot
the logits, because much of the data is at the extremes where $p=0$ or
$(1-p) = 0$.  Often the best we can do is the plots that we did
earlier, perhaps adding some sort of smooth line to the data.  (Next
example, we will try that).  We will be able to do better once we have
a model and can consider residuals---but be careful which residuals we
think about.


\subsubsection*{The model}

The S-Plus logistic regression command needs to know both the number
of successes and the number of failures, so we need to start out by
creating such a variable,
\begin{verbatim}
> FN <- cbind(Fail, n-Fail)
\end{verbatim}

Now, just use the {\tt glm} command, with the {\tt family=binomial}
option, to fit the model, viz.
\begin{verbatim}
> mymod <- glm( FN ~ Temp, family=binomial)
> mymod <- glm( cbind(Fail, n-Fail) ~ Temp, family=binomial)
\end{verbatim}

Either of the above forms will work, but the first is easier in the
long run.   Now look at the results of the fit,
\begin{verbatim}
> summary(mymod)

Call: glm(formula = FN ~ Temp, family = binomial)
Deviance Residuals:
        Min         1Q     Median          3Q      Max 
 -0.9522773 -0.7829968 -0.5411832 -0.04378972 2.651492

Coefficients:
                 Value Std. Error   t value 
(Intercept)  5.0848620  3.0477855  1.668379
       Temp -0.1155992  0.0469317 -2.463136

(Dispersion Parameter for Binomial family taken to be 1 )

    Null Deviance: 24.23036 on 22 degrees of freedom

Residual Deviance: 18.08633 on 21 degrees of freedom

Number of Fisher Scoring Iterations: 4 

Correlation of Coefficients:
     (Intercept) 
Temp -0.9932182 
> anova(mymod)
Analysis of Deviance Table

Binomial model

Response: FN

Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev 
NULL                    22   24.23036
Temp  1 6.144035        21   18.08633
\end{verbatim}

The model fits well (although, just like $R^2$ is regression, one
should be cautious about this, and the temperature coefficient is
certainly significant.

At this point we need to
\begin{itemize}
\item Get predictions for the low-temperature extrapolations, and
\item Evaluate the fit of the model
\end{itemize}

\paragraph{Extrapolating}


Lets work on that, in that order (which is really back-to-front, but
useful for teaching purposes).

The fitted values (on the probability scale) are obtained via,
\begin{verbatim}
> fitted(mymod)
         1         2         3          4          5          6          7 
 0.2607849 0.1817866 0.1652202 0.09994062 0.07278457 0.06535816 0.06535816
          8          9         10         11         12         13         14 
 0.06535816 0.05864156 0.05257638 0.04710711 0.04710711 0.04710711 0.04710711
         15         16         17         18         19         20        21 
 0.03775043 0.03376847 0.02698616 0.02698616 0.02411117 0.02411117 0.0192299
         22         23 
 0.01716665 0.01367159
\end{verbatim}

The predicted values (on the logit scale) are,
\begin{verbatim}
> predict(mymod)
         1        2         3         4         5         6         7 
 -1.041893 -1.50429 -1.619889 -2.197885 -2.544682 -2.660281 -2.660281
         8        9        10        11        12        13        14 
 -2.660281 -2.77588 -2.891479 -3.007079 -3.007079 -3.007079 -3.007079
        15        16        17        18        19        20        21 
 -3.238277 -3.353876 -3.585074 -3.585074 -3.700673 -3.700673 -3.931872
        22        23 
 -4.047471 -4.278669
> fitted(mymod) - unlogit(predict(mymod))
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 
 0 0 0 0 0 0 0 0 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
\end{verbatim}

(unlogit is the following simple function),

\begin{verbatim}
> unlogit
function(x)
exp(x)/(1 + exp(x))
\end{verbatim}

For communication purposes, it makes sense to work on the probability
scale (and hence use {\tt fitted}), but for all other purposes the
predicted values are what we want.   In particular, the {\tt predict}
function has an option that will give the Standard Errors of the
predicted values (on the logit scale), viz.,
\begin{verbatim}
> predict(mymod,se=T)
$fit:
         1        2         3         4         5         6         7 
 -1.041893 -1.50429 -1.619889 -2.197885 -2.544682 -2.660281 -2.660281
         8        9        10        11        12        13        14 
 -2.660281 -2.77588 -2.891479 -3.007079 -3.007079 -3.007079 -3.007079
        15        16        17        18        19        20        21 
 -3.238277 -3.353876 -3.585074 -3.585074 -3.700673 -3.700673 -3.931872
        22        23 
 -4.047471 -4.278669

$se.fit:
         1         2        3         4         5         6         7 
 0.7682786 0.5943294 0.556385 0.4298909 0.4298809 0.4441498 0.4441498
         8         9        10        11        12        13        14 
 0.4441498 0.4647341 0.4908397 0.5216385 0.5216385 0.5216385 0.5216385
        15        16        17        18        19        20        21 
 0.5942932 0.6348851 0.7222049 0.7222049 0.7682357 0.7682357 0.8637766
        22       23 
 0.9129226 1.013299

$residual.scale:
[1] 1.189907

$df:
[1] 21
\end{verbatim}

So, a 95\% confidence interval for the first predicted value on the
logit scale would be
$$
-1.041893  \pm 2 \times 0.7682786 $$
Transferred to the probability scale, this is an interval
$(0.070,0.621)$, with a center at $0.260$.  

But what about predictions nearer the take-off temperature?

\begin{verbatim}
 predict(mymod, data.frame(Temp=c(30,31,32,33)), se=T)
$fit:
        1        2        3       4 
 1.616887 1.501288 1.385689 1.27009

$se.fit:
        1        2        3       4 
 1.972255 1.917739 1.863302 1.80895
\end{verbatim}

(That was pretty sneaky, but the {\tt predict} function requires a
data frame predictor values).   

The bottom line is that at $31^\circ$F our model extrapolates to
$p = 0.83$, with a 95\% confidence interval of
$(0.089, 0.996)$, which is VERY wide.   Even with this wide interval,
the point estimate is something to worry about.

\paragraph{Evaluating the Model}

\begin{verbatim}
> plot(mymod)
\end{verbatim}

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

Notice that the residual plot has at least two recognizable
components, but it couldn't help but look like this.  There is perhaps
some suggestion that a function of temperature might be needed, like
perhaps $t^2$, but don't be fooled.  The places where there are big
residuals are for temperatures where there were both no failures, and
failures.   I tried a number of additional functions of temperature,
and nothing much worked.  

\subsubsection*{The Zero-One version of the model}

Very quickly, lets look at the model when we just look at a failure or
not.

\begin{verbatim}
> Bn <- cbind(sign(Fail), 1 - sign(Fail))
> Bn
   [,1] [,2] 
 1    1    0
 2    1    0
 3    1    0
 4    1    0
 5    0    1
 6    0    1
 7    0    1
 8    0    1
 9    0    1
10    0    1
11    0    1
12    0    1
13    1    0
14    1    0
15    0    1
16    0    1
17    0    1
18    1    0
19    0    1
20    0    1
21    0    1
22    0    1
23    0    1
\end{verbatim}
then
\begin{verbatim}
> mymod <- glm(Bn ~ Temp, family=binomial)
> summary(mymod)

Call: glm(formula = Bn ~ Temp, family = binomial)
Deviance Residuals:
       Min       1Q     Median        3Q      Max 
 -1.061112 -0.76128 -0.3782822 0.4523828 2.217506

Coefficients:
                  Value Std. Error   t value 
(Intercept)  15.0422911  7.3366528  2.050293
       Temp  -0.2321537  0.1076141 -2.157279

(Dispersion Parameter for Binomial family taken to be 1 )

    Null Deviance: 28.26715 on 22 degrees of freedom

Residual Deviance: 20.31519 on 21 degrees of freedom

Number of Fisher Scoring Iterations: 4 

Correlation of Coefficients:
     (Intercept) 
Temp -0.99716   
> anova(mymod)
Analysis of Deviance Table

Binomial model

Response: Bn

Terms added sequentially (first to last)
     Df Deviance Resid. Df Resid. Dev 
NULL                    22   28.26715
Temp  1  7.95196        21   20.31519
\end{verbatim}

So far things look very similar, including the residual plot (believe
me).


\subsection*{Another Example, Stroke}


\subsubsection*{Background}

The data are from a prospective study whose purpose was to
characterize the spontaneous recovery of language functioning in
patients during the first 3 to 4 months following a stroke.  Details
of the design and results of this study may be found in Holland et
al. (1985) and Holland, Greenhouse, Fromm, and Swindell (1989).
Briefly, patients admitted to a hospital with the diagnosis of stroke
and who met the study inclusion criteria were seen daily (6 days/week)
for 15 minutes beginning 48 to 72 hours poststroke and continuing
until hospital discharge.  The daily visits consisted of
semistructured conversational interactions and systemic scoring of the
patient's communication behavior.  At hospital discharge, patients
were formally tested using the Western Aphasia Battery (WAB).  If the
WAB indicated aphasic dysfunction (AQ $<$ 93.8), the patient was seen
and assessed at 1 month postdischarge. If this second WAB demonstrated
continued difficulty the patient was assessed again at 2 months
postdischarge.  The data set to be discussed here consists of 61
patients who had unilateral (right- or left-hemisphere) stroke,
ischemic or hemorrhagic etiology, and impaired language functioning at
time of hospital discharge.  We will use these data to investigate the
correlates of language recovery during the two month period
post-hospital discharge.  Definitions and the codes for the variables
of interest are presented in Table~\ref{tab:defn}.

The multiple logistic regression model will allow us to investigate
simultaneously the effect of a number of different variables on the
likelihood of spontaneous language recovery.  A feature of the
logistic regression model is that it allows us to use explanatory
variables that are either quantitative or discrete.  Using the
variables defined and coded in Table~\ref{tab:defn} for the stroke
recovery study, 
the specification of the deterministic component of the multiple
linear logistic regression model on the logit scale is
\begin{center}
\hfill{}\small{
$ln(\frac{\pi}{1-\pi})  =  \beta_{0} + \beta_{1}\cdot (SIDE) + \beta_{2}\cdot
(GENDER) + \beta_{3}\cdot (AGE) + \beta_{4}\cdot (LOSD) 
 + \beta_{5}\cdot (WAB1)$}. \hfill{}(1)
\end{center}
Note that if we were to dichotomize length of stay, age and WAB1 then
a multiple cross-classification table of these data would consist of
$2^{6}=64$ cells or subgroups and that with 61 patients such a table
would have numerous sparse or empty cells which would be a problem for
a chi-square analysis.  An advantage of a statistical model such as
the one given in equation (1), however, is that the model provides a
structure that smooths the fit of the probability of a response to all
the cells in such a table.

\begin{table}[htbp]
\begin{verbatim}
Side  Gender  Age  LOSD  WAB1  Status

  0     1      74   17   92.5    1
  1     1      75   30   87.8    0
  0     1      72   17   65.9    0
  0     1      60   21   86.4    1
  0     0      69   11   90.0    1
  0     0      62   13   6.2     0
\end{verbatim}
  \caption{Example of what the data file 'stroke.dat' looks like}
  \label{tab:data}
\end{table}



\begin{table}[htb]
\begin{center}
    \leavevmode
\begin{tabular}{lll} \hline \hline \\
{\em Variable} & {\em Description} & {\em Values/Codes}\\ 
\hline\\
STATUS & Status at last & 0 = Not Recovered \\
& observation & 1 = Recovered\\
\\
SIDE & Side of lesion & 0 = Left\\
& & 1 = Right\\
\\
GENDER & & 0 = Female\\
& & 1 = Male\\
\\
AGE & Age at stroke & In years\\
& onset &\\
\\
LOSD & Length of hospital & In days\\
& stay &\\
\\
WAB1 & WAB at  & Range: 0 - 100\\
& hospital discharge\\
\hline \hline
\end{tabular}
\end{center}
  \caption{Variables and Codes from the Stroke Recovery Study }
  \label{tab:defn}
\end{table}


\subsubsection*{Interpretation of Regression Coefficients}

The output below presents part of the S-Plus output listing the
estimates of the regressions coefficients and their standard errors.
We see that the 
coefficients for gender, length-of-stay, and the WAB score at
discharge are all significantly related to the probability of
spontaneous language recovery at the 0.05 level.  The coefficient for
length-of-stay is negative indicating that the longer patients are in
the hospital the smaller the chance of spontaneous language recovery
relative to patients who have shorter lengths-of-stay even after
controlling for the effect of the WAB1 score.  The coefficient for
WAB1 is positive as we would expect.  Males, after controlling for the
other variables in the model, are $exp(2.88)=17.8$ times as likely to
recover spontaneous language function as females.  A 95\% confidence
interval for $\beta_{2}$, the adjusted log odds ratio of recovery for
males relative to females, however, is large and ranges from 0.47 to
5.29.  We also note that side of stroke and subject's age do not have
a statistically significant association with spontaneous language
recovery after controlling for the effects of the other variables.  
Details about these tests can be found in Hosmer and
Lemeshow (1989, Chapter 5).

\begin{verbatim}
> stroke <- read.table("stroke.dat", head=T) 
> stroke[1:10,]
   Side Gender Age LOSD WAB1 Status 
 1    0      0  79   27 18.2      0
 2    0      0  82   17 86.0      0
 3    1      0  72   42 46.3      0
 4    0      0  59   23 35.6      0
 5    0      0  60   28  0.0      0
 6    1      1  75   12 81.9      1
 7    0      0  75   53 56.0      0
 8    0      0  85   46 17.3      0
 9    0      0  79   58 74.9      0
10    0      1  82   34 18.2      0

> mymod <- glm(Status ~ Side + Gender + Age + LOSD +WAB1, family=binomial,
+ data=stroke)  

> summary(mymod)

Call: glm(formula = Status ~ Side + Gender + Age + LOSD + WAB1, family = binomial, 
        data = stroke)
Deviance Residuals:
       Min         1Q       Median            3Q      Max 
 -1.317524 -0.2999725 -0.003778562 -0.0002212266 2.862154

Coefficients:
                    Value Std. Error    t value 
(Intercept) -11.072759221 6.28633143 -1.7614024
       Side  -1.369026127 1.10318267 -1.2409786
     Gender   2.876569185 1.25325078  2.2952862
        Age  -0.006938071 0.03062231 -0.2265691
       LOSD  -0.196137069 0.09595028 -2.0441532
       WAB1   0.164467235 0.07172681  2.2929674

(Dispersion Parameter for Binomial family taken to be 1 )

    Null Deviance: 65.71913 on 60 degrees of freedom

Residual Deviance: 28.9878 on 55 degrees of freedom

Number of Fisher Scoring Iterations: 7 

Correlation of Coefficients:
       (Intercept)       Side     Gender        Age       LOSD 
  Side  0.2553975                                             
Gender -0.5007016  -0.2793548                                 
   Age -0.3373404  -0.1265137  0.2006563                      
  LOSD  0.1199385   0.2046374 -0.5309336 -0.1894954           
  WAB1 -0.9292116  -0.3055970  0.5278826  0.0616042 -0.2969773
> anova(mymod)
Analysis of Deviance Table

Binomial model

Response: Status

Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev 
  NULL                    60   65.71913
  Side  1  0.10011        59   65.61902
Gender  1  3.36010        58   62.25892
   Age  1  1.45449        57   60.80444
  LOSD  1 12.25425        56   48.55018
  WAB1  1 19.56238        55   28.98780
> 
\end{verbatim}


\subsubsection*{Analysis of the Fit of the Model}


Figure~\ref{fig:resid} presents a residual analysis for this logistic
regression. 

One relatively easy method for the analysis of the fit of the logistic
regression model is based on calculating the predicted probability of
spontaneous language recovery based on the model in (1) compared to
the observed probability of language recovery 
Specifically, substitute the estimated regression coefficients given
in above into the multiple logistic regression equation in (1).
For each patient in the study calculate his or her individual
predicted probability of language recovery by calculating the linear
combination of each patient's individual set of explanatory variables
with the estimated regression coefficients and use the inverse of the
logit transformation to solve for the probability of language
recovery, $\hat{\pi}$.
For example, consider a patient whose stroke was on the left side
(SIDE=0), was a 61 year old male (GENDER=1, AGE=61), stayed in the
hospital 19 days (LOSD=19), and whose WAB at discharge was 81.6
(WAB1=81.6).  The estimated probability of spontaneous language
recovery based on the model for such a patient is
\begin{displaymath}
  \hat{\pi}  =  \frac{1}{1 +
e^{-(-11.1-1.37\cdot(0)+2.88\cdot(1)-0.0069\cdot(61)-0.196\cdot(19)+0.165\cdot(81.6))}} = 0.77 .
\end{displaymath}
This value can be interpreted to mean that 77 patients out of 100 with
the above characteristics can be expected to recover normal language
functioning spontaneously within the first 2 months post-hospital
discharge.  To study the goodness of fit of the model, we obtain the
frequency distribution of the values of $\hat{\pi}$, and determine, say,
the quintiles of the distribution.  The observed number of patients
who had normal language functioning are tallied within each quintile,
and the expected number of patients who had normal language
functioning (predicted by the model) is determined by summing the
$\hat{\pi}$'s for all the patients in the quintile.
Table~\ref{tab:fit} presents 
the expected and observed number of patients who had normal language
functioning for each quintile of risk of language recovery.
Comparisons of the expected number predicted by the model to the
observed number suggests that this model fits the data reasonably
well.

A test for whether there is a significant difference between the
observed and expected numbers in Table~ref{tab:fit} can be constructed
based on a 
Pearson chi-square test.  This chi-square test would have
degrees-of-freedom equal to the number of groups (e.g., 5 for the
number of quintiles) minus 2 (Hosmer and Lemeshow, 1989, chapter 5).
In our example, this chi-square statistics is equal to 3.60 with 3
degrees-of-freedom and has a non-significant $p$-value=0.69, indicating
that the model adequately fits the data. A small p-value would
indicate a significant lack of fit of the model.


\begin{table}[htb]
  \begin{center}
    \leavevmode
\begin{tabular}{rccc}
 &                 & Expected  & Observed  \\
 & Quintile of risk & Number    & Number \\ \hline
Highest & 5    & 10.7   & 11 \\
        & 4    &  3.1   &  2  \\     
        & 3    &  0.2   &  1  \\
        & 2    &  0.0   &  0  \\ 
Lowest  & 1    &  0.0   &  0  \\ \hline
Total   &      & 14.0   & 14  \\ \hline
\end{tabular}
    \caption{Analysis of the fit of the logit model in Equation (1)}
    \label{tab:fit}
  \end{center}
\end{table}


\begin{figure}[htb]
  \begin{center}
    \leavevmode
    \psfig{figure=stroke.eps,angle=-90,width=6in}
    \caption{Residual plots}
    \label{fig:resid}
  \end{center}
\end{figure}





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


