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

% \input info.tex

\newcommand{\Var}{\mbox{Var}\,}

\newcommand{\eq}[1]{(\ref{#1})}

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

%\slides

\parindent=0pt
\parskip=0.5em

\begin{document}

\begin{center}\Large\bf
Two Examples of Logistic Regression
\end{center}


\section{Budworms}

This data set is analyzed by Venables and Ripley also, pp. 230 ff.
{\em Check page ref's in current edition.}

The data is the result of an experiment on the toxicity of the tobacco
budworm to doses of particular drug, that the moths are beginning to
show resistance to.  Batches of 20 moths of each sex were exposed for
3 days to the chemical and the number of dead (or disabled) moths in
each group were recorded.  The data is,
\begin{verbatim}
Dose Response Gender
1       1       M
2       4       M
4       9       M
8       13      M
16      18      M
32      20      M
1       0       F
2       2       F
4       6       F
8       10      F
16      12      F
32      16      F
\end{verbatim}

The doses were in micro-grams, and the dose level when up by powers of
two.  It is a judgement call as to whether you will try to model with
dose as the response or $\log_2(\mbox{dose})$.  We will try both.

\begin{verbatim}
> bud <- read.table("budworm.dat", head=T)
> bud
   Dose Response Gender 
 1    1        1      M
 2    2        4      M
 3    4        9      M
 4    8       13      M
 5   16       18      M
 6   32       20      M
 7    1        0      F
 8    2        2      F
 9    4        6      F
10    8       10      F
11   16       12      F
12   32       16      F
> attach(bud,pos=1)
> N <- 20
\end{verbatim}

First, lets plot the data, and see if anything obvious shows up.

\begin{verbatim}
> plot( Dose, Response/N)         
> points( Dose[Gender=="F"], (Response/N)[Gender=="F"], pch="F")             
> points( Dose[Gender=="M"], (Response/N)[Gender=="M"], pch="M")
\end{verbatim}

The same plot on the log(Dose) scale, is shown in Figure~\ref{freq-vs-logdose}.

\begin{verbatim}
> plot( Dose, Response/N, log="x")                              
> points( Dose[Gender=="F"], (Response/N)[Gender=="F"], pch="F")
> points( Dose[Gender=="M"], (Response/N)[Gender=="M"], pch="M")
\end{verbatim}

\begin{figure}
\centerline{\psfig{figure=budplot2.eps,angle=-90,height=3.5in}}
\caption{Frequency of budworm death (mortality) vs. log-dose of
pesticide.}\label{freq-vs-logdose} 
\end{figure}

Given that we have more than 0-1 data, it makes a lot of sense to at
least try to plot the (observed) logits; Figure~\ref{dose-vs-logdose}
does this for both dose and log(dose) of the pesticide.

\begin{verbatim}
> par(mfrow=c(2,1))
> plot( Dose, logit(Response/N), log="x")
> title("Logit vs Log(Dose)")            
> plot( Dose, logit(Response/N))         
> title("Logit vs Dose")
\end{verbatim}

\begin{figure}
\centerline{\psfig{figure=budplot3.eps,angle=-90,width=5in}}
\caption{Logit-mortality vs Dose and Log-dose of
pesticide.}\label{dose-vs-logdose}
\end{figure}

This plot makes it quite clear, that the log scale has a much better
chance of working.

So, proceed to fit the model.  We surely want a model with different
slopes for the two sexes
\begin{verbatim}
> ldose <- log(Dose)
> SF <- cbind(Response, N - Response)
> mymod <- glm(SF ~ Gender * ldose, family=binomial)   
> summary(mymod)

 ..........
Coefficients:
                  Value Std. Error    t value 
 (Intercept) -2.9935414  0.5525295 -5.4178852
      Gender  0.1749865  0.7781556  0.2248733
       ldose  1.3071341  0.2410133  5.4234939
Gender:ldose  0.5091459  0.3894474  1.3073547


    Null Deviance: 124.8756 on 11 degrees of freedom

Residual Deviance: 4.993727 on 8 degrees of freedom

> anova(mymod,test="Chisq")
Analysis of Deviance Table

Binomial model

Response: SF

Terms added sequentially (first to last)
             Df Deviance Resid. Df Resid. Dev   Pr(Chi) 
        NULL                    11   124.8756          
      Gender  1   6.0770        10   118.7986 0.0136955
       ldose  1 112.0415         9     6.7571 0.0000000
Gender:ldose  1   1.7633         8     4.9937 0.1842088
> 
\end{verbatim}
Now, what is going on here?  {\tt Gender} seems to have no effect
(based on the t-value), but the deviance is quite large?  Even more,
our plots seem to imply what looks like a gender effect.  The best way
to find out is to distrust the parameter stuff (it can get easily
confused), and try fitting a model without the interaction.
This makes sense anyway, because there is little/no evidence of an
interaction.


\begin{figure}
\centerline{\psfig{figure=budresid.eps,angle=-90,height=3.5in}}
\caption{Standard Splus residual diagnostic plots for {\tt glm(SF \~{}
Gender + ldose, family=binomial)}.}\label{diag-mf+ldose}
\end{figure}

\begin{verbatim}
> anova(mymod <- glm(SF ~ Gender + ldose, family=binomial))        

Terms added sequentially (first to last)
       Df Deviance Resid. Df Resid. Dev 
  NULL                    11   124.8756
Gender  1   6.0770        10   118.7986
 ldose  1 112.0415         9     6.7571
> summary(mymod)
Coefficients:
                Value Std. Error   t value 
(Intercept) -3.473154  0.4682939 -7.416612
     Gender  1.100743  0.3557226  3.094385
      ldose  1.535336  0.1890119  8.122959

    Null Deviance: 124.8756 on 11 degrees of freedom

Residual Deviance: 6.757064 on 9 degrees of freedom
\end{verbatim}

That is a little more like it, and seems to be a VERY good fitting
model.  We should at least look at some diagnostics, and some other
plots (see Figure~\ref{diag-mf+ldose}), but unless they show something
very odd, we are done.

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

Figure~\ref{diag-mf+ldose} looks very good.  We should finish by
plotting the data and fitted curves on the same graph, and that is
left as an exercise for you (or look it up in the Venables and
Ripley).  Venables and Ripley also go on to talk about LD50, and how
to estimate it.

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


\section{An Extended Example: Stroke}


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


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


\subsection{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}
\centerline{\psfig{figure=stroke.eps,angle=-90,width=6in}}
\caption{Splus's standard residual plots for {\tt glm(Status \~{} Side
+ Gender + Age + LOSD +WAB1, family=binomial, +
data=stroke)}.}\label{fig:resid}
\end{figure}


\subsection{How to plot these data}

Look at the stroke data.  This is explicitly 0--1 data, and it is very
hard to think about plotting it.  A simple plot of, say, {\tt Status}
against {\tt WAB1}, is shown in Figure~\ref{raw-stroke}, and it not
all that informative.  Even if I jitter the points a little, it is
still not very easy to see any trends or patterns.
\begin{verbatim}
> plot(WAB1, Status)
> plot(WAB1, Status+0.1*runif(length(Status)))
\end{verbatim}

\begin{figure}
\centerline{
\psfig{figure=stroke1.eps,angle=-90,width=3in}
\psfig{figure=stroke2.eps,angle=-90,width=3in}}
\caption{Stroke status (0/1) vs WAB at intake (WAB1) (unjittered on
the left; jittered on the right).}\label{raw-stroke}
\end{figure}




A somewhat more useful approach is to plot the data and add some sort
of smooth curve that will allow you to at least look at some of the
trends in the data.  The particular approach I took with these data
was to superimpose a (black-box) smooth curve generated by loess.  The
various pictures in Figure~\ref{loess-wab1} result from different
smoothing parameters.
\begin{verbatim}
> scatter.smooth(WAB1, Status,span=0.2) 
> title("Span = 0.2")
> scatter.smooth(WAB1, Status,span=0.35)
> title("Span = 0.35")                  
> scatter.smooth(WAB1, Status,span=0.5) 
> title("Span = 0.5")                  
> scatter.smooth(WAB1, Status,span=0.65)
> title("Span = 0.65")                  
\end{verbatim}

\begin{figure}
\centerline{\psfig{figure=stroke3.eps,angle=-90,height=3.5in}}
\caption{Various loess smooths of stroke status vs.
WAB1.}\label{loess-wab1}
\end{figure}

This is still not perfect, but you certainly get the sense that the
probability goes up, probably more or less monotonely with WAB1.  

Compare this to the similar picture(s) for Age (Figure~\ref{loess-age}).

\begin{verbatim}
> for (i in c(0.2, 0.35, .5, 0.65)){  
+ scatter.smooth(Age, Status, span=i) 
+ title(paste("Span = ", i))
+ }
\end{verbatim}

\begin{figure}
\centerline{\psfig{figure=stroke4.eps,angle=-90,height=3.5in}}
\caption{Various loess smooths of stroke status vs. age.}\label{loess-age}
\end{figure}

If anything, there is a bit of evidence for some type of quadratic
effect of Age, but the probable story is that Age is not that
important, {\bf on its own}.  We should always remember that
variables that do not seem to be important (from these marginal
plots), may be quite important when paired with other variables (that
is viewed conditionally).

With that hint about the Age variable, I re-ran the logistic
regression model for these data and obtained,
\begin{verbatim}
> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age + Age^2, family=bin>
> summary(mymod)

Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2, family = 
        binomial)
Deviance Residuals:
       Min         1Q        Median            3Q      Max 
 -1.317304 -0.2637303 -0.0002506098 -6.014239e-07 2.770847

Coefficients:
                   Value   Std. Error   t value 
(Intercept) -42.13652444 17.240544854 -2.444037
       Side  -1.70612965  1.320176507 -1.292350
     Gender   4.03384776  1.847166535  2.183803
       LOSD  -0.28388626  0.130491084 -2.175522
       WAB1   0.24891576  0.104667163  2.378165
        Age   0.92453559  0.421765743  2.192059
   I(Age^2)  -0.00791642  0.003578441 -2.212254

(Dispersion Parameter for Binomial family taken to be 1 )

    Null Deviance: 65.71913 on 60 degrees of freedom

Residual Deviance: 21.7621 on 54 degrees of freedom

> 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
    LOSD  1 13.46451        57   48.79441
    WAB1  1 19.75552        56   29.03889
     Age  1  0.05110        55   28.98780
I(Age^2)  1  7.22570        54   21.76210

\end{verbatim}

This is somewhat surprising (because I didn't think there was an Age
effect).  Could we have seen the same trend from a residual analysis
of these data?

\begin{figure}
\centerline{\psfig{figure=stroke5.eps,angle=-90,height=3.5in}}
\caption{Residual diagnostic plots for {\tt glm(Status \~{} Side
+ Gender + Age + LOSD +WAB1, family=binomial, +
data=stroke)}.}\label{fig:resid-again}
\end{figure}

Let us go back to the original model and look at residuals
(Figure~\ref{fig:resid-again}).  There is at least one large (Pearson)
residual, and it is the same point as the large (deviance) residual.
What about a plot of (deviance) residual against {\tt Age} (when {\tt
Age} is in the model).

\begin{verbatim}
> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age, family=binomial)
> scatter.smooth(Age, (residuals(mymod)))                                 
> identify(Age, (residuals(mymod)))                                       
[1] 59
> scatter.smooth(Age[-59], (residuals(mymod))[-59])                       
\end{verbatim}

We see that observation {\tt 59} is the one with the large residual.
If I plot residuals of my model against age (the lower plot in
Figure~\ref{59-inf}), the quadratic trend seems to go away, at least a
bit.

\begin{figure}
\centerline{\psfig{figure=stroke6.eps,angle=-90,height=3.5in}}
\caption{Loess smooths of residuals of {\tt glm(Status \~{} Side +
Gender + LOSD +WAB1 + Age, family=binomial)}, with and without
observation 59.}\label{59-inf}
\end{figure}

Note that this shows that observation \#59 is highly influentialon the
coefficient of {\tt Age\^{}2} in the model.  This should show up in
Cook's D or DFBETA plots as well.  Observation {\tt 59} is now
suspect, so let us remove that observation, and see what the resulting
model looks like.

\begin{verbatim}
> WW <- rep(1, length(Side))
> WW[59] <- 0
> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age,
+ weight=WW, family=binomial)                            
Warning messages:
  linear convergence not obtained in  10  iterations. in: glm.fitter(x = X, y =
        Y, w = w, start = start, offset = offset, family = family, ....
\end{verbatim}

Oh-oh, what does that mean?   (Basically, the optimization algorithm
is converging slowly, which is an indication that there are some
fitted probabilities that are very close to 0 or 1).  We can try a few
more iterations, via.

\begin{verbatim}
> mymod <- glm(Status ~ Side + Gender + LOSD +WAB1 + Age,
+ weight=WW, maxit=30,family=binomial)                   
\end{verbatim}

Now, what does the summary and anova look like?

\begin{verbatim}
> summary(mymod)

Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age, family = binomial, 
        weights = WW, maxit = 30)
Deviance Residuals:
      Min           1Q        Median            3Q      Max 
 -1.21685 -0.006138687 -2.107342e-08 -2.107342e-08 2.120409

Coefficients:
                    Value  Std. Error     t value 
(Intercept) -46.034665461 23.60416082 -1.95027757
       Side  -7.479606389  4.79366143 -1.56031178
     Gender   7.889726328  4.28891089  1.83956406
       LOSD  -0.677382165  0.43418335 -1.56012928
       WAB1   0.662724777  0.34770876  1.90597667
        Age  -0.001289562  0.04096666 -0.03147832

(Dispersion Parameter for Binomial family taken to be 1 )

    Null Deviance: 62.71879 on 60 degrees of freedom

Residual Deviance: 13.32727 on 54 degrees of freedom

Number of Fisher Scoring Iterations: 10 

Correlation of Coefficients:
       (Intercept)       Side     Gender       LOSD       WAB1 
  Side  0.8469550                                             
Gender -0.9187997  -0.8860396                                 
  LOSD  0.7944690   0.9174434 -0.8555061                      
  WAB1 -0.9831416  -0.9061715  0.9317554 -0.8758313           
   Age -0.0387868   0.0615410  0.0243350 -0.0022418 -0.0596445
Warning messages:
  1 rows with zero weights not counted in: summary(mymod)
\end{verbatim}

If I now start adding the quadratic term in {\tt Age}, I get myself
into all sorts of trouble.  

\begin{verbatim}
> formula(mymod)
Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2
> anova(mymod)
Analysis of Deviance Table

Binomial model

Response: Status

Terms added sequentially (first to last)
         Df Deviance Resid. Df Resid. Dev 
    NULL                    60   62.71879
    Side  2  0.62858        58   62.09021
  Gender  1  2.34693        57   59.74328
    LOSD  1 13.67051        56   46.07277
    WAB1  1 32.74451        55   13.32826
     Age  1  0.00099        54   13.32727
I(Age^2)  1 13.32727        53    0.00000
Warning messages:
1: linear convergence not obtained in  10  iterations. in: glm.fit(x, y, w,
        offset = offset, family = family, maxit = control$maxit,  ....
2: linear convergence not obtained in  10  iterations. in: glm.fit(x, y, w,
        offset = offset, family = family, maxit = control$maxit,  ....
\end{verbatim}

Whoa.  What is going on here?  We suddenly have a model that fits
perfectly?

The summary looks even weirder?
\begin{verbatim}
> summary(mymod)

Call: glm(formula = Status ~ Side + Gender + LOSD + WAB1 + Age + Age^2, family = 
        binomial, weights = WW, maxit = 50)
Deviance Residuals:
         Min            1Q        Median            3Q          Max 
 -3.9065e-05 -2.107342e-08 -2.107342e-08 -2.107342e-08 3.281089e-05

Coefficients:
                    Value  Std. Error       t value 
(Intercept) -2858.6085408 2227320.001 -0.0012834297
       Side  -164.0618144  165491.082 -0.0009913635
     Gender   273.5752941  237195.052  0.0011533769
       LOSD   -16.2464405   13633.473 -0.0011916582
       WAB1    24.1122280   17153.825  0.0014056473
        Age    39.8782980   35880.527  0.0011114190
   I(Age^2)    -0.3549147     302.566 -0.0011730157

(Dispersion Parameter for Binomial family taken to be 1 )

    Null Deviance: 62.71879 on 60 degrees of freedom

Residual Deviance: 0 on 53 degrees of freedom

Number of Fisher Scoring Iterations: 29 

Correlation of Coefficients:
         (Intercept)       Side     Gender       LOSD       WAB1        Age 
    Side  0.3661039                                                        
  Gender -0.8397332  -0.5132934                                            
    LOSD  0.7676113   0.5867586 -0.7405370                                 
    WAB1 -0.9678428  -0.5023038  0.8384746 -0.7822973                      
     Age -0.9537356  -0.2802197  0.7986261 -0.8085521  0.8643352           
I(Age^2)  0.9585630   0.3372172 -0.8129406  0.8428222 -0.8863173 -0.9963403
Warning messages:
1: fitted values close to 0 or 1 in: family$deviance(mu, y, w, residuals = T)
2: 1 rows with zero weights not counted in: summary(mymod)
> 
\end{verbatim}

On reflection, it is clear what is going on.  Look at the first plot
of Status against WAB1 (Figure \ref{raw-stroke}).  Here we see, that
absent the one observation in the middle (which is 59), we have an
almost perfect discrimination between the two levels of Status.  That
is the clue to what is ``broken''.  With the 7 parameters (anything
like WAB1 squared would have done), we can find a function that
perfectly discriminates between the two levels of status, and then
maximum likelihood modelling breaks down. 


\subsection{Influence of a point on the logistic regression
coefficients: DFBETA}

Pregibon's (1981) approximation to DFBETA is not pre-programmed into
Splus, and for a small data set like this it is just easy to compute
DFBETA's exactly.  The basic idea of the following code is to fit the
{\tt glm} model with each data point in turn removed.  The summary is
then a plot of each of the coefficients, with as many points on the
plot as there are data points.  As in OLS regression, the plot is
surprisingly effective and quickly shows us where the pathology in the
data lies.

{ % \small
\begin{verbatim}
# Heavily annotated, for class

# Generate a vector of all ones, of length the number of observations.
W <-  rep(1,length(Side))
# Generate current weights.
WW <- W
# Fit the base model
mymod <- glm(Status ~ Gender + LOSD + Side + Age + WAB1, weight=WW,
   family=binomial)

# Save the coefficients, just in case
base.coef <- coefficients(mymod)

# Set up a matrix to hold the coefficients for each model with one
# point deleted.
infs <- matrix(0, nrow=length(Side), ncol=length(coefficients(mymod)))

# Loop through the data, setting one weight equal to 0.
# Save the coefficients in the infs matrix
for (i in 1:length(Side)){               
        WW<-W; WW[i] <- 0;                       
        infs[i,] <- coefficients(update(mymod, . ~ .))
}
## Can't think of an easy, smart summary, so just plot the
## coefficients.

#Set up a 2 x 3 plotting frame
par(mfrow=c(2,3))

#Do the 6 plots.  After each plot, allow the user to
# identify the points which look interesting.
for(i in 1:6){
        plot(1:length(Side), infs[,i])
        title(names(coefficients(mymod))[i])
        identify(1:length(Side), infs[,i])
}
\end{verbatim}
}

After running that code, and interacting with the various plots, you
end up with the plot in Figure~\ref{dfbeta}.
There are other things we could have plotted (like following the
deviance as the models changed), but this plot alone tells us much of
the story.  

\begin{figure}
\centerline{\psfig{figure=strokeinf.eps,width=6in}}
\caption{Exact DFBETA plots, obtained by refitting the model with each
observation deleted in turn.}\label{dfbeta}
\end{figure}





\end{document}

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





\end{document}

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



