#  This is code to read the file containing
#  the BPD data as a data frame in Splus.

bpd <- read.table("c:/temp/bpd.dat",header=T)

#  Code for making a scatter plot matrix
#  The par( ) function is used to specify
#  features of the plot such as the shape 
#  and size of plotting symbols and the size
#  of the plot.  The following options are
#  specified

#  fin=c(w,h)  requests a plot that us w inches
#              wide and h inches high
#  pch=18      requests a filled rectangle as a 
#              plotting symbol
#  mkh=b       requests plotting symbols that
#              that are b inched high
#  cex=1.2     increases size of printed letters
#  lwd=3       controls widths of lines and curves 


par(din=c(7,7), pch=18, mkh=.15, cex=1.2, lwd=3) 
pairs(bpd, panel=function(x,y){ points(x,y)})


# Add the natural logagrithms of the exposure
# variables and squares and crossproducts to 
# the data frame

bpd$lnL <- log(bpd$Low+1)
bpd$lnM <- log(bpd$Medium+1)
bpd$lnH <- log(bpd$High+1)

# Use the  glm  function to fit a
# logistic regression model  

  bpd2 <- glm(BPD ~ lnL + lnM + lnH, 
              family=binomial, data=bpd)

# Print a summary of the results

summary(bpd2)
anova(bpd2, test="Chisq")


# Search the simplest model that 
# adequately describes the trends
# in the data

# Forward selection: try entering the
#   explanatory variables in different
#   orders

 bpd2 <- glm(BPD ~ lnM + lnH + lnL, 
              family=binomial, data=bpd)
summary(bpd2)
anova(bpd2, test="Chisq")


 bpd2 <- glm(BPD ~ lnH + lnL + lnM, 
              family=binomial, data=bpd)
summary(bpd2)
anova(bpd2, test="Chisq")