\documentstyle[11pt,fullpage,psfig,times]{article}

\parindent=0pt
\parskip=0.5em

\long\def\comment#1{}

\newcommand{\combos}[2]{\left(\begin{array}{c}#1\\#2\end{array}\right)}

\let\bar\overline

\def\eq#1{(\ref{#1})}

\def\logit{\mbox{logit}\,}

\def\tp{^T}

\def\Var{\mbox{Var}\,}
\let\var\Var
\def\sigsq{{\hat\sigma^2_\epsilon}}

\def\Varsigma{{\mathbf{\not{\hspace*{-0.33em}\Sigma}}}}
\def\tr{\mbox{tr\,}}

\def\diag{\mbox{\it diag\,}}

\let\vc\Varsigma
\let\phi\varphi

\def\Corr{\mbox{Corr\,}}

\let\ul\mb
\let\by\times

\def\ul#1{\underline{#1}}
\def\ulem#1{\ul{\em #1}}

% \def\mb#1{{\mathbf #1}}
\let\mb\ul

\def\xtx{X\tp X}

\def\Re{{\cal R}}

\newenvironment{boxed}{\begin{center}\begin{tabular}{|p{5.5in}|}
\hline}{\\ \hline\end{tabular}\end{center}}

\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{lemma}{Lemma}

%\slides

\begin{document}

\begin{center}\Large\bf
A Little Generalized Linear Models and a Lot of Logistic Regression
\end{center}

\section{Some References}

Around 1972, the analysis of categorical data was revolutionized by
Nelder and Wedderburn who recognized that many disparate models for
discrete data could be unified as special cases of the {\em
  generalized linear model (GLM)}.  A big advantage of GLM's is that
they make the analysis of discrete categorical data (which had been a
bag of mostly seemingly-unrelated tricks) look much more like linear
regression and ANOVA.  For example, we can use our intutions about
fitting models with {\tt lm()} and {\tt aov()} to inform our analysis
using {\tt glm()}\footnote{Not to be confused with PROC GLM and the
  ``general linear model'' in SAS.  PROC GLM simply fits linear
  regression and ANOVA models.  The SAS procedure that corresponds to
  {\tt glm()} in SPLUS is called PROC GENMOD.} in SPLUS.

You can learn about log-linear models and logistic regression models
from the following sources (among others!).

\begin{description}

\item Agresti, A. (1990). {\em Categorical data analysis.}  New York:
  Wiley.
  
  A fully-detailed textbook on categorical data, generalized linear
  models, etc.

\item Chambers, J. M. and Hastie, T. J. (1992).  {\em Statistical
    models in S.}  Pacific Grove, CA: Wadsworth \& Brooks/Cole.  {\sl
    [Especially Chapter 6.]}
  
  A nice summary of what {\tt glm()} can do.  

% \item Fienberg, S. E. (1987).  {\em The analysis of cross-classified
%     categorical data, 2$^{nd}$ Edition.}  Cambridge MA: MIT Press.
%   
%   Can be thought of as an abbreviated version of Agresti (1990).  Much
%   more emphasis on tables of counts, much less on logistic regression
%   and generalized linear models. Many of the examples and intuitions
%   for this part of the course will come from this book.

\item Fox, J. (1997). {\em Applied regression analysis, linear models,
    and related methods.}  SAGE publications, Thousand Oaks CA.

  Chapter 15 contains a nice introduction to logit models for
  dichotomous and polytomous data, and GLM's.

\item Hosmer, D. W., Jr. and Lemeshow, S. (1989).  {\em Applied
logistic regression.} New York: Wiley.

  From a customer review at {\tt amazon.com}: ``This book is widely
  considered the `bible' of logistic regression analysis. It
  provides an accessible introduction to the theory of logistic
  modeling, and gives in-depth coverage of the proper use of the
  method, including interpretation, diagnostics, and practical
  considerations.'' {\em An updated edition appeared in 2000.}

\item McCullagh, P. and Nelder, J. A. (1989).  {\em Generalized
    linear models, 2$^{nd}$ edition.} London: Chapman and Hall.
  
  This book really codified the revolution in discrete data analysis,
  by two of the main developers of the generalized linear model
  approach.

\item Rodr\'{i}gues, G. (Fall 2001).  {\em WWS509 Generalized Linear
  Models.}  {\tt http://data.princeton.edu/}

  Online notes for a course in GLM's from a research demographer at
  Princeton.  Very nice.

\item Venables, W. N. and Ripley, B. D. (1994).  {\em Modern
    applied statistics with Splus.}  New York: Springer-Verlag.
    {\sl [Especially Chapter 7.]}

  As usual, Venables and Ripley is better after you sort of know what
  is going on.  But then it is a great reference and reminder of what
  is possible.  {\em Current edition is 3${^rd}$, some page number and
  chapter number references have changed\ldots}

\end{description}

\section{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, for a total of six O-rings per
flight.  There was data about the number of O-rings out of six that
failed on each of 23 previous space shuttle flights.  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\footnote{Some engineers did look at more complete data, or at
least saw that the likelihood of failure at $31^\circ$F was high; and
implored management to delay the flight.  The reasons that they were
not listened to, and instead the plot of the failure-only data was
used to justify liftoff, makes a fascinating study in (the failure of)
scientific writing.}.  They produced a plot that looked something like
Figure~\ref{fail-only}.

\begin{verbatim}
> chal <- read.table("chal.dat", head=T)
> attach(chal)
> names(chal)
[1] "Temp" "Fail"

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

\begin{figure}
\centerline{\psfig{figure=chal1.eps,height=3.5in}}
\caption{Plot of the number of O-ring failures including only space
shuttle flights prior to Challenger in which failures
occured.}\label{fail-only}
\end{figure}

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
focused on a plot of all the  data, as in Figure~\ref{all-data}, they
might have behaved differently.

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

\begin{figure}
\centerline{\psfig{figure=chal2.eps,height=3.5in}}
\caption{Plot of the number of O-ring failures including all space
shuttle flights prior to Challenger, whether or not failures
occured.}\label{all-data}
\end{figure}

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.




% \subsection{How to look at these data}

\subsection{The logistic model}

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).  

\comment{

In the past when $y$ (or for that matter $x$) has been count,
fraction or proportion data we have looked at models of the form
\begin{eqnarray*}
        \log(\mb y) & = & X\mb\beta + \mb\epsilon ~, \\
        \logit(\mb y) & = & X\mb\beta + \mb\epsilon ~,\\
        {\rm probit}\,(\mb y) & = & X\mb\beta + \mb\epsilon ~,~{\rm etc.}
\end{eqnarray*}
With the Challenger data this is not likely to work well {\em
directly}, because of all the 0's in the data (and in other settings,
the high frequency of 0\% and 100\% proportions that we might see).

}

Instead of modelling $\mb y$ directly, we will model $p(t)$, 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.

\comment{

\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 (since it is bounded above by zero).  

Let's 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)& =& \combos{n}{r} p(t)^r (1-p(t))^{n-r}  \\
          & = & \combos{n}{r} {\left ( \frac{p(t)}{1-p(t)} \right )}^r (1-p(t))^{n} \\
          & = & \mbox{(Fcn of $r$)}\times\mbox{(Fcn of $p$ and $r$)} \times
                 \mbox{(Fcn 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{\combos{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 
\[\mbox{logit}\,(p) = \log{\left ( \frac{p}{1-p}\right )} 
\]
is called the
logistic transformation.  If 
\[ x= \mbox{logit}\,(p) = \log{\left ( \frac{p}{1-p}\right )}
\] 
then the inverse function is
\[ p = \mbox{expit}\,(x) = \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.


\subsection{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.


\subsection{The model}

The S-Plus command to fit logistic regressions is {\tt glm}, which
also fits members of the more general class of models known as {\em
generalized linear models}.  To fit a logistic regression, {\tt glm}
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, we 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.   

The {\tt family} option tells Splus what distribution the individual
$y$'s should have, and we have chosen {\tt binomial} since that is the
way our model was developed above.  There are other {\tt family}
values for other kinds of data (poisson, Gaussian, inverse Gaussian,
Gamma, etc.).  By default the {\tt binomial} family uses the {\tt
logit link function}, i.e. it builds a linear model for $\logit
p(y=1|x)$.  If you wanted to fit a {\em probit} model, for example,
you would say {\tt family=binomial(link = probit)} instead of just
{\tt family=binomial}.

Now let's 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}

So what does all that mean?  And what's a ``deviance''??

\subsection{Maximum likelihood, standard errors, IRLS}

\paragraph{Likelihood equations.}

{\tt glm} is not just doing a least-squares fit, but actually
maximizing the likelihood.    Let's re-cap the model, and see how the
fit should work.  Let
\[
        X = \left[\begin{array}{ll}
        1 & t_1 \\
        1 & t_2 \\
        \vdots & \vdots \\
        1 & t_{23} 
        \end{array}\right]
        = \left[\begin{array}{l}
        \mb x_1\tp \\
        \mb x_2\tp \\
        \vdots \\
        \mb x_{23}\tp \\
        \end{array}\right]
\]
and recall that our model has the form
\[
        p_i \equiv P[{\em failure}|\mb x_i\tp] 
        =  p(\mb x_i\tp\mb\beta) = \frac{\exp[\mb x_i\tp\mb\beta]}
                {1 + \exp[\mb x_i\tp\mb\beta]}
\]
and therefore if we experience $y_i$ failures in $n_i$ attempts
($n_i\equiv 6$ for us), $i=1,\ldots,N=23$, we have
\begin{eqnarray*}
        P(R=y_i | \mb x_i\tp)& =& \combos{n_i}{y_i} p_i^{y_i }
                        (1-p_i))^{n_i-y_i}  \\
          & = & \combos{n_i}{y_i} \left ( 
        \frac{p_i}{1-p_i)} \right )^{y_i} (1-p_i)^{n_i} \\
          & = & \combos{n_i}{y_i} \exp[y_i\cdot \mb x_i\tp\mb\beta] 
         (1 + \exp[ \mb x_i\tp\mb\beta])^{-n_i}
\end{eqnarray*}
where in the last line I have applied the useful identity $1/(1-p_i) = 1
+ p_i/(1-p_i)$.  It follows that the log-likelihood for all $N=23$ shuttle
flights is
\begin{eqnarray*}
        \log L(\ul y | X) & = & \log \prod_{i=1}^N P(R=y_i | \mb x_i\tp)\\
        & = & \sum_i \log \combos{n_i}{y_i} + \sum_i y_i\cdot \mb x_i\tp\mb\beta
                - \sum_i n_i\log(1 + \exp[  \mb x_i\tp\mb\beta])
\end{eqnarray*}
Differentiating with respect to $\mb\beta$ we obtain
\[
        \frac{\partial\log L}{\partial\mb\beta} ~ = ~
        \sum_i y_i \,\mb x_i ~ - ~ \sum_i n_i \, 
        \frac{\exp[ \mb x_i\tp\mb\beta]}{1 + \exp[ \mb x_i\tp\mb\beta]} 
        \, \mb x_i
\]
Setting this equal to zero to obtain the likelihood equations, and
noting the definition of $ p_i = p(\mb x_i\tp\mb\beta)$ above, we get
\[
        \sum_i n_i \, p_i \, \mb x_i
        ~ = ~ 
        \sum_i y_i \,\mb x_i 
\]
Since $n_i p_i  = E[y_i|\mb x_i\tp]$, we have the usual
maximum likelihood intepretation that this equates the observed
sufficient statistics $\sum_i y_i \mb x_i$ with their expected
values.  

It is also worth noting that the likelihood equations above can be
written as
\begin{equation}
        X\tp \mb f(X\mb\beta) = X\tp y,
        \label{l-e}
\end{equation}
where 
\[
        \mb f(X\mb\beta) = \left[\begin{array}{c}
        n_1 p(\mb x_1\tp\mb\beta) \\
        n_2 p(\mb x_2\tp\mb\beta) \\
        \vdots \\
        n_N p(\mb x_N\tp\mb\beta) 
        \end{array}\right]
\]
is the vector of ``fitted values'' ($\hat y_i = f_i = p(\mb
x_i\tp\mb\beta)$), and this is essentially like the normal equations
of OLS regression,
\[
        X\tp X\mb\beta = X\tp y,
\]
where $\hat{\mb y} = X\mb\beta$ is the vector of fitted values for the
linear regresion.

\paragraph{Asymptotic standard errors.}

Once we have found the MLE's $\hat{\mb\beta}$ by solving \eq{l-e}, their
asymptotic information matrix can be calculated as usual as
\[
        I(\mb\beta) = - E\left[
        \frac{\partial^2 \log L}{\partial\mb\beta\partial\mb\beta\tp}
        \right]
        = \mbox{ [many steps omitted] } =
        X\tp V X
\]
where 
\[
        V ~ = \left[\begin{array}{cccc}
        n_1 p_1(1-p_1) & 0 & \cdots & 0 \\
        0 & n_2 p_2(1-p_2) & \cdots & 0 \\
        \vdots & \vdots & \ddots & \vdots \\
        0 & 0 & \cdots & n_N p_N(1-p_N) 
        \end{array}\right]
\]
and hence the asymptotic variance-covariance matrix of
$\hat{\mb\beta}$ is
\[
        \Sigma_\beta = \Var(\hat{\mb\beta}) ~ = ~ (X\tp V X)^{-1}
\]
This looks exactly like the variance-covariance matrix for
$\hat{\mb\beta}$ in a WLS regression.  It leads immediately to
asymptotic confidence ellipsoids for linear functions $L\mb\beta$ of
the parameters, based on the so-called {\em Wald statistics} [see
Appendix D.6.2 in Fox, or section 15.3 of Rawlings et al., for example]
\[
        (L\mb\beta - L\mb\beta_0)\tp 
        [L\Sigma_\beta L\tp]^{-1}
        (L\mb\beta - L\mb\beta_0) ,
\]
and similarly for hypothesis tests, confidence intervals for
individual parameters, approximate Scheffe' corrections for multiple
inferences, etc., etc.

\paragraph{Finding the MLE's.}

Solving the likelihood equations \eq{l-e} requires some sort of
iterative method.  ``Newton's method'' from elementary calculus for
finding a root $\hat\beta$ of the equation $\ell(\beta) = 0$ by
iterating
\[
        \beta_{k+1} = \beta_k - \ell(\beta_k)/\ell'(\beta_k)
\]
can be generalized to the multivariate-$\beta$ case, and the case in
which $\ell(\mb\beta) = \partial\log L(\mb\beta)/\partial \mb\beta$ is
the derivative of the log-likelihood.  In the case of logistic
regression, this boils down (after much fooling around) to iterating
\[
        \mb\beta_{k+1} 
        = \mb\beta_k + (X\tp V_k X)^{-1} X\tp (\mb y - f(X\mb\beta_k))
\]
where $V_k$ is defined as $V$ above, but using $\mb\beta_k$.  Note
that this is a plausible thing to do, since when $\beta_k$ has
converged to a constant value $\mb\beta_k\equiv\hat{\mb\beta}$, we get
\begin{equation}
        \hat{\mb\beta} 
        = \hat{\mb\beta} + (X\tp VX)^{-1} X\tp (\mb y - \mb f(X\hat{\mb\beta}))
        \label{irls-start}
\end{equation}
so that \eq{l-e}, $X\tp f(X\hat{\mb\beta}) = X\tp \mb y$, is in fact
satisfied.

Note that \eq{irls-start} may be rewritten as
\[
        \hat{\mb\beta} = (X\tp V X)^{-1} X\tp V \mb y^*
\]
where
\[
        \mb y^* = X\hat{\mb\beta} + V^{-1}(\mb y - \mb f(X\hat{\mb\beta}))
\]
which suggests that an alternate way to iterate a solution of the
likelihood equations is
\[
\begin{array}{rrcl}
        \mbox{Let} & \mb y_k^* 
                & = & X\mb\beta_k + V_k^{-1}(\mb y - \mb f(X\mb\beta_k))
                        \\[0.5em]
        \mbox{and then} & \mb\beta_{k+1} 
                & = & (X\tp V_k X)^{-1} X\tp V_k \mb y_k^*
\end{array}
\]
Note that at each stage, $\mb\beta_{k+1}$ are the estimated regression
coefficients from a WLS problem with weight matrix $V_k^{-1/2}$.  This
method is called {\em iteratively reweighted least squares} (IRLS) and
is a standard way to calculate MLE's for generalized linear models; it
is in fact exactly what {\tt glm} is doing in Splus.




\subsection{Evaluating the model}

\paragraph{Model deviance.}

Since we are out of the realm of least-squares (except as a
computational tool, as with IRLS above), the $F$-testing approach
using $c\cdot(RSS_{small}-RSS_{big})/RSS_{big}$ to model fit is of
limited use in logistic regression (and glm's more generally).  A
useful generalization of RSS (it is essentially equivalent for OLS
regression with normal errors) is the {\em model deviance}
\[
        G^2 = -2\log L
\]
where $\log L$ is the log-likelihood we were working with above.
Exploiting the analogy between $G^2$ and RSS, 
\[
        R^2 = 1 - \frac{G^2}{G^2_0}
\]
is a reasonable analogy for $R^2$, where $G^2$ is the deviance of the
model we are interested in, and $G^2_0$ is the deviance of the
intercept-only model.

Using classical MLE theory we also know that for two nested models,
``big'' and ``small'',
\[
        G^2_{small} - G^2_{big} = -2\log\frac{L_{small}}{L_{big}}
        \sim \chi^2_{df}
\] 
where $df$ is the difference in the number of parameters in the two
models (really, the difference in ranks of the corresponding $X$
matrices). 

\paragraph{Residuals.}

A basic reference for evaluating the fit of logistic regression models
is Pregibon (1981, {\em Annals of Statistics, 9,} 705--724).

The raw residuals 
\[
        E_i  = y_i - n_i p_i
\]
are not very useful in logistic regression, since their variance
$n_ip_i(1-p_i)$ varies greatly with $p_i=p(\mb x_i\tp\hat{\mb\beta})$
[and $n_i$].  The {\em Pearson residuals}
\[
        Z_i = \frac{E_i}{\sqrt{n_ip_i(1-p_i)}}
\]
have more stable variances; and the sum of their squares
\[
        Z^2 = \sum_i Z_i^2
\]
provides a measure of goodness-of-fit analogous to $G^2$.  Recall that
\[
        G^2 = -2\log L = -2 \sum_i [ y_i \log p_i + (n_i - y_i)\log(1-p_i) ] ;
\]
this is a measure of the goodness-of-fit of the model also, since it
implicitly contrasts the fit of the model with the fit of a model that
fits the data perfectly (so that $L=1$ and $\log L=0$).  The
contribution of individual observations to $G^2$ is sometimes written
as
\[
        G_i 
        = {\tt sign}(E_i) \sqrt{-2[ y_i \log p_i + (n_i - y_i)\log(1-p_i) ]} ;
\]
these $G_i$ are called {\em deviance residuals}.

Pearson residuals and deviance residuals are the most commonly used
for diagnosing the fit of a logistic regression.  The default
diagnostic plots for logistic regression (and any glm) in Splus
involve plots of both deviance residuals and Pearson residuals.
\begin{verbatim}
> par(mfrow=c(2,2))
> plot(mymod)
\end{verbatim}
The result is plotted in Figure~\ref{std-diag-chal}.  

\begin{figure}
\centerline{\psfig{figure=chal3.eps,angle=-90,height=3.5in}}
\caption{Splus's standard logistic regression diagnostic plots for the
Challenger data.}\label{std-diag-chal}
\end{figure}

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 the evidence is less than clear.  The places where
there are big residuals are for temperatures where there were both no
failures, and failures. 

Another kind of residual that is sometimes used is the {\em working
residuals}; these are the terms $V^{-1}(\mb y - \mb
f(X\hat{\mb\beta}))$ in the working $y$ vector 
\[
        \mb y^* = X\hat{\mb\beta} + V^{-1}(\mb y - \mb f(X\hat{\mb\beta}))
\]
from IRLS.  All three kinds of residuals can be obtained in Splus
using the {\tt residuals(glm.object, \\ type=[your choice])} command,
where ``{\tt[your choice]}'' should be one of {\tt response} (for raw
residuals), {\tt pearson}, {\tt deviance} or {\tt working}.  Fox
(p. 458) calls working residuals ``logit residuals''.

\paragraph{Leverage and ``delete-one'' diagnostics.}

An analogue of the hat matrix,
\[
        H = V^{1/2} X(X\tp V X)^{-1}V^{1/2}
\]
is available (for example, $(I-H)\mb Z = \mb Z$), and the diagonal
elements of it can be interpreted the same way as $h_{ii}$ in linear
regression, as measures of leverage of each observation.

Recall that in OLS we could obtain studentized residuals by fitting
the usual model, but with an additional indicator variable for
observation $i$.  The same thing can be done in logistic regression:
we fit the model
\[
        logit(\mb p) = X\mb\beta + \gamma D_{{\rm obs~} i}
\]
and then the same model without $\gamma$.  The log-likelihood ratio
statistic (or deviance statistic)
\[
        G^2_{(-i)} = -2\log(L_{\gamma}/L_{\gamma=0})
\]
is asymptotically $\chi^2_1$ and its square root $G_{(-i)}$ (with sign
chosen to agree with $E_i$ is analogous to a studentized residual.


It can be expensive to refit the model $N$ times to obtain each
studentized residual in this way, but Pregibon (1981) shows that
\[
        G^2_{(-i)} \approx G_i^2 + \frac{Z^2_i h_{ii}}{1-h_{ii}} ;
\]
substituting $G_i$ for $Z_i$ and taking the square root we obtain
approximate studentized residuals
\[
        G_i^* = \frac{G_i}{\sqrt{1-h_{ii}}} .
\]
The $G_i^*$ are approximately normally distributed when the model
holds (hence {\tt qqnorm} plots etc. are useful for checking the
model, checking for outliers, erc; though discreteness of the response
can still make interpretation difficult).

Similarly $DFBETA$ plots and Cook's distances can be obtained exactly
by refitting the model $N$ times and looking at $\hat{\mb\beta} -
\hat{\mb\beta}_{(-i)}$.  Approximations are again provided by
Pregibon (1981): For $DFBETA$,
\[
        \hat{\mb\beta} - \hat{\mb\beta}_{(-i)} 
        \approx (X\tp V X)^{-1} \mb x_i \frac{E_i}{1-h_{ii}}
\]
where, as usual, $\mb x_i$ is the $i^{th}$ row of the $X$ matrix,
written as a column vector; and for Cook's $D$,
\[
        D_i \approx \frac{E_i^2}{n_ip_i(1-p_i)} \cdot
                \frac{h_{ii}}{(1-h_{ii})^2}
        = Z^2_i \cdot \frac{h_{ii}}{(1-h_{ii})^2}
\]
These are all very much like the analogous formulas for OLS regression
(cf. e.g. p. 290 of Fox); some more details on their derivation are
given on pp. 457--462 of Fox; see also Pregibon (1981) and Landwehr et
al. (1980, {\em Proceedings of the ASA Business and Economics
Section,} 15--20).

Figure~\ref{delete-ones} illustrates some of these diagnostics.  
\begin{verbatim}
> X _ model.matrix(mymod)
> V _ diag(mymod$weights)
> W _ diag(sqrt(mymod$weights))
> H _ W%*%X%*%solve(t(X)%*%V%*%X)%*%t(X)%*%W
> hats _ diag(H)
> studres _ residuals(mymod,type="deviance")/sqrt(1-hats)
> CooksD _ residuals(mymod,type="pearson")^2 * hats/(1-hats)^2
> source("../week10/qqenv.splus")
> par(mfrow=c(2,2))
> plot(predict(mymod),studres)
> abline(h=0) ; identify(predict(mymod),studres)
> qqenv(studres) 
> plot(hats) 
> plot(CooksD)
\end{verbatim}


\begin{figure}
\centerline{\psfig{figure=chal-del1.eps,height=3.5in}}
\caption{Some additional diagnostic plots for the Challenger
data.}\label{delete-ones}
\end{figure}

\subsection{Fitted values and extrapolations}

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 of 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.




\subsection{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}
> newmod <- glm(Bn ~ Temp, family=binomial)
> summary(newmod)

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(newmod)
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

> plot(newmod)

> X _ model.matrix(newmod)
> V _ diag(newmod$weights)
> W _ diag(sqrt(newmod$weights))
> H _ W%*%X%*%solve(t(X)%*%V%*%X)%*%t(X)%*%W
> hats _ diag(H)
> studres _ residuals(newmod,type="deviance")/sqrt(1-hats)
> CooksD _ residuals(newmod,type="pearson")^2 * hats/(1-hats)^2
> par(mfrow=c(2,2))
> plot(predict(newmod),studres)
> abline(h=0) ; identify(predict(newmod),studres)
> qqenv(studres) 
> plot(hats) 
> plot(CooksD)
\end{verbatim}
Figure~\ref{new-diag} gives the diagnostic plots.

\begin{figure}

\centerline{\psfig{figure=std-diag.eps,height=3.5in}}
\centerline{\psfig{figure=new-del1.eps,height=3.5in}}


\caption{Standard and additional diagnostic plots for the 0/1 model
for the Challenger data.}\label{new-diag}
\end{figure}



\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



