teach.df <- data.frame(read.table("http://jgill.wustl.edu/data/spector.mazzeo.data",
        col.names=c("GPA","TUCE","PSI","GRADE")))

 teach.df$GRADE <- factor(teach.df$GRADE)
 levels(teach.df$GRADE) <- c("FAIL","PASS")

 teach.df$PSI <- factor(teach.df$PSI)
 levels(teach.df$PSI) <- c("OLD","NEW")
 attach(teach.df)

summary(teach.df)
       GPA             TUCE        PSI      GRADE   
 Min.   :2.060   Min.   :12.00   OLD:18   FAIL:21  
 1st Qu.:2.812   1st Qu.:19.75   NEW:14   PASS:11  
 Mdian :3.065   Median :22.50                     
 Mean   :3.117   Mean   :21.94                     
 3rd Qu.:3.515   3rd Qu.:25.00                     
 Max.   :4.000   Max.   :29.00                     

postscript("teach.ps")
par(mfrow=c(1,3),mar=c(2,2,2,2),oma=c(4,3,2,2))
boxplot(GPA~GRADE,col="green",main="GPA")
boxplot(TUCE~GRADE,col="green",main="TUCE")
plot(PSI[GRADE=="FAIL"],col="gray",main="gray for FAIL",ylim=c(0,15))
par(new=T)
dev.off()

teach.logit.fit <- glm(GRADE ~ GPA + TUCE + PSI, 
	family=binomial(link=logit), data = teach.df)
summary(teach.logit.fit)

Call:
glm(formula = GRADE ~ GPA + TUCE + PSI, family = binomial(link = logit), data = teach.df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9551  -0.6453  -0.2570   0.5888   2.0966  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept) -13.02119    4.91635  -2.649  0.00808 
GPA           2.82609    1.26051   2.242  0.02496  
TUCE          0.09516    0.14135   0.673  0.50083   
PSINEW        2.37866    1.06242   2.239  0.02516  
---

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 41.183  on 31  degrees of freedom
Residual deviance: 25.779  on 28  degrees of freedom
AIC: 33.779

Number of Fisher Scoring iterations: 4

glm.summary(teach.logit.fit)

            Coefficient Std. Error 0.95 CI Lower 0.95 CI Upper
(Intercept)   -13.02119    4.91635     -22.65706      -3.38533
GPA             2.82609    1.26051       0.35554       5.29663
TUCE            0.09516    0.14135      -0.18189       0.37220
PSINEW          2.37866    1.06242       0.29636       4.46097

# FIRST DIFFERENCE PLOT
par(mfrow=c(1,2),cex=1.3)
ruler <- seq(min(GPA),max(GPA),length=500)
xbeta <-  1/(1+exp(-teach.logit.fit$coef[1] 
        - teach.logit.fit$coef[2]*ruler
        - teach.logit.fit$coef[3]*mean(teach.df$TUCE)))
plot(ruler,xbeta,type="b",xlab="Levels of GPA",
        ylab="Probability of Passing",ylim=c(0,1))
text(3.4,0.2,"OLD Method")
xbeta2 <-  1/(1+exp(-teach.logit.fit$coef[1] 
        - teach.logit.fit$coef[2]*ruler 
        - teach.logit.fit$coef[3]*mean(teach.df$TUCE) 
        - teach.logit.fit$coef[4]))
lines(ruler,xbeta2)
text(3.1,0.8,"NEW Method")

ruler <- seq(min(TUCE),max(TUCE),length=500)
xbeta <-  1/(1+exp(-teach.logit.fit$coef[1] 
        - teach.logit.fit$coef[3]*ruler
        - teach.logit.fit$coef[2]*mean(teach.df$GPA)))
plot(ruler,xbeta,type="b",xlab="Levels of TUCE",
        ylab="Probability of Passing",ylim=c(0,1))
text(22,0.18,"OLD Method")
xbeta2 <-  1/(1+exp(-teach.logit.fit$coef[1] 
        - teach.logit.fit$coef[3]*ruler 
        - teach.logit.fit$coef[2]*mean(teach.df$GPA) 
        - teach.logit.fit$coef[4]))
lines(ruler,xbeta2)
text(19,0.6,"NEW Method")

# RESIDUALS PLOT
par(mfrow=c(2,2))
qqnorm(resid(teach.logit.fit,type="response"),ylab="Response")
mtext(side=3,"Normal Quantiles vs. Response Residuals")
abline(0,0.33)
qqnorm(resid(teach.logit.fit,type="pearson"),ylab="Pearson")
mtext(side=3,"Normal Quantiles vs. Pearson Residuals")
abline(0,0.75)
qqnorm(resid(teach.logit.fit,type="deviance"),ylab="Deviance")
mtext(side=3,"Normal Quantiles vs. Deviance Residuals")
abline(0,1)
qqnorm(resid(teach.logit.fit,type="working"),ylab="Working")
mtext(side=3,"Normal Quantiles vs. Working Residuals")
abline(0,2)