s you can see in the commands, a logistic regression is fit to
y (the last column in the data file) and the first four columns (the
predictor variables).
The command is 'glm').

Say the output is stored in 'logitmodel'

>logitmodel = glm(y ~ x, family="binomial")

>summary(logitmodel)       #gives the ANOVA table.

also
>logitmodel$aic    #spits out the AIC value.


You can repeat the above steps for different combinations of predictor variables and compute the AIC in each case and select the model that gives you the 'minimum' AIC.

---
To test the model following the steps in Chapter 15 of
Helsel and Hirsch book

>logitmodel$null.deviance    gives Null Deviance = G02
>logitmodel$deviance         gives the Residual Deviance = G2


The above two are output from the summary
LR0 = logitmodel$null.deviance - logitmodel$deviance   #G02 - G2
k = length(x[1,])   # the number of predictor variables.
if
>1-pchisq(LR0, k)
is less than 'alpha'
then reject the Null hypothesis - which is that all the coefficients
corresponding to the variables are = 0 (i.e., model is no good)
in favour of alternate hypothesis that the coeffecients are different
from 0 (i.e., model is good)

-- 
ANOVA table is performed on the fitted model. This is similar to the
ANOVA table from the linear regression.
Then see which variables are insignificant and drop them out and fit the
logistic regression again with the reduced set of
variables. This is your 'best' model.

---------------------------------------------------------------------------

Extra commands.....

-----------------------------------------------

Below is also commands to perform the logistic regression using the
equations..

Here is what you do.
>zz=glm(y ~ x, family=binomial(link=logit))
(if you don't give link=logit wihin binomial it should still be fine -
this just being explicit)

>n=length(y)
>X=cbind(rep(1,n),x)       #create the X matrix with the first column as 1

> ypred=predict.glm(zz)       #this will predict the 'log of odds ratio'

>Now get this into the original space..
> ypredo=exp(ypred)/(1+exp(ypred))      #this will have values between
0 and 1

> vmat=diag(ypredo*(1-ypredo))      (this is a diagonal matrix of size
n x n - i.e. off diagonal elements are 0)

>varbetas=solve(t(X) %*% vmat %*% X)      #recall this is similar to
solve(t(X) %*% X)*errorvar in the case of multiple

                #regression

>sqrt(diag(varbetas))          #will spit out the standard errors of
the coefficients.

(http://www.roguewave.com/support/docs/sourcepro/analyticsug/3-3.html)
This link has the equations..

>summary(zz)                     #this will print out the summary table
a la 'ls.print'
The standard error column from this should match up with that from the
above command (sqrt(diag(varbetas))