\section*{End of Stroke example}

\subsection*{Influence of a point}

There {\bf are} good measures of influence for glm models, but they
don't seem to be implemented (i.e., I could not find them) in the
standard packages.  Rather than just give up, I quickly hacked
together a ``Poor Man's'' influence function.   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.  This plot is surprisingly effective and quickly shows us
where the pathology in the data lies.

\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 following plot.

\centerline{\psfig{figure=strokeinf.eps,width=6in}}

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.  


\section*{More Examples}

Here I will pull together a number of examples, drawn directly from
the Book by Venables and Ripley.

First, the budworm data.  The following 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

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

\centerline{\psfig{figure=budplot2.eps,angle=-90,width=5in}}

Given that we have more than 0-1 data, it makes a lot of sense to at
least try to plot the (observed) logits,

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

\centerline{\psfig{figure=budplot3.eps,angle=-90,width=5in}}

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{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, but unless they show something very odd, we are done.

\begin{verbatim}
> plot(mymod) 
\end{verbatim}
\centerline{\psfig{figure=budresid.eps,angle=-90,width=5in}}

That 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 book).   The book also goes on to talk about LD50,
and how to estimate it


%\section*{Low Birth Weight example}

%This example consists of a dataset of 189 births at a US hospital,
%with the main interest in low birth weight.  The available variables,
%are

%\begin{tabular}{l l}
%{\tt low} & Birth Weight less than 2.5kg (0/1) \\
%{\tt age} & age of mother in years \\
%{\tt lwt} & weight of mother (lbs) at last menstrual period \\
%{\tt race} & white/black/other \\
%{\tt smoke} & smoking status during pregnancy (0/1) \\
%{\tt ptl} & number of previous premature labours \\
%{\tt ht} & history of hypertension (0/1) \\
%{\tt ui} & has uterine irritability (0/1) \\
%{\tt ftv} & number of physician visits in first trimester \\
%{\tt bwt} & actual birth weight (grams) 
%\end{tabular}

%and the first few lines of the data are
%\begin{verbatim}
%> birthwt
%low age lwt race smoke ptl ht ui ftv  bwt 
%  0  19 182    2     0   0  0  1   0 2523
%  0  33 155    3     0   0  0  0   3 2551
%  0  20 105    1     1   0  0  0   1 2557
%  0  21 108    1     1   0  0  1   2 2594
%  0  18 107    1     1   0  0  1   0 2600
%  0  21 124    3     0   0  0  0   0 2622
%  0  22 118    1     0   0  0  0   1 2637
%  0  17 103    3     0   0  0  0   1 2637
%  0  29 123    1     1   0  0  0   1 2663
%  0  26 113    1     1   0  0  0   0 2665
%\end{verbatim}

%Although the actual birth weights are available, we shall concentrate
%on estimating the 0/1 variable of low birth weight or not.  
%Both {\tt ptl} and {\tt ftv} are (in some sense) ordinal, but their
%exact values are unlikely to be important.  We reduce {\tt ptl} to a
%0/1 variable and {\tt ftv} to a variable with just three levels
%(0/1/2+).
%\begin{verbatim}
%> ptd <- factor(ptl)
%> ftv <- factor(ftv)
%> ftv
%  [1] 0 3 1 2 0 0 1 1 1 0 0 1 0 2 0 0 0 3 0 1 2 3 1 0 2 0 0 2 0 1 1 1 1 1 0 2 2
% [38] 0 2 1 2 2 1 0 0 0 4 0 2 0 1 0 0 2 0 0 0 0 0 0 0 2 0 0 0 1 2 6 1 2 0 2 1 0
% [75] 0 0 1 4 0 0 1 0 0 0 0 0 0 0 0 1 0 2 0 0 0 1 1 0 0 1 1 0 0 1 0 0 1 0 2 4 2
%[112] 1 2 1 0 1 0 0 2 1 1 0 1 0 2 2 1 0 1 1 0 2 0 0 0 0 1 1 0 1 0 0 0 1 0 2 2 0
%[149] 0 0 1 2 0 0 0 0 3 1 0 0 0 1 0 0 0 0 4 0 1 0 1 0 0 0 0 0 1 3 0 2 1 3 0 0 2
%[186] 2 0 0 3
%> levels(ftv)[ -(1:2) ] <- "2+"
%> ftv
% [1] 0  2+ 1  2+ 0  0  1  1  1  0  0  1  0  2+ 0  0  0  2+ 0  1  2+ 2+ 1  0  2+
%[26] 0  0  2+ 0  1  1  1  1  1  0  2+ 2+ 0  2+ 1  2+ 2+ 1  0  0  0  2+ 0  2+ 0 
%[51] 1  0  0  2+ 0  0  0  0  0  0  0  2+ 0  0  0  1  2+ 2+ 1  2+ 0  2+ 1  0  0 
%[76] 0  1  2+ 0  0  1  0  0  0  0  0  0  0  0  1  0  2+ 0  0  0  1  1  0  0  1 
%[101]1  0  0  1  0  0  1  0  2+ 2+ 2+ 1  2+ 1  0  1  0  0  2+ 1  1  0  1  0  2+
%[126]2+ 1  0  1  1  0  2+ 0  0  0  0  1  1  0  1  0  0  0  1  0  2+ 2+ 0  0  0 
%[151]1  2+ 0  0  0  0  2+ 1  0  0  0  1  0  0  0  0  2+ 0  1  0  1  0  0  0  0 
%[176]0  1  2+ 0  2+ 1  2+ 0  0  2+ 2+ 0  0  2+
%\end{verbatim}

%In the interim, I have created a data frame ({\tt bwt}), with these
%modified variables.


% 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: 


