# 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")