# R code for model building
# We analyze the surgical unit data
# using multiple linear regression


# Save the data file into a directory and 
# use the full path name:

surg.data <- read.table(file = "z:/My Documents/teaching/stat_704/surgicalunitdata.txt", 
header=FALSE, col.names = c('x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'y', 'lny'))

# attaching the data frame:

attach(surg.data)

# scatter plot matrix:

pairs(surg.data)

# pairwise correlations:

cor(surg.data)

# To install the "leaps" package:
# Go to "Packages" menu and choose
# "Install package(s) from CRAN"
# scroll to the "leaps" package and choose it
# An Internet connection is needed to install it.

# Then you can proceed with the code...

# Load the "leaps" package:

library(leaps)

# There is no easy equivalent to SAS's "STEPWISE" in R (I don't think, anyway)
# The step function does something similar...

step(lm(lny ~ 1), data=surg.data, scope = ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)

# But the leaps function provides the all-subsets analysis:

# Printing the 2 best models of each size, using the Cp criterion:

leaps(x = cbind(x1, x2, x3, x4, x5, x6, x7, x8), y = lny, nbest=2, method="Cp")

# Printing the 2 best models of each size, using the adjusted R^2 criterion:

leaps(x = cbind(x1, x2, x3, x4, x5, x6, x7, x8), y = lny, nbest=2, method="adjr2")

# Calculating the AIC of any particular model:

# For example, here is the AIC of the 8-predictor model:

surg.fit.full <- lm(lny ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)

AIC(surg.fit.full)

## a version of BIC or Schwarz' BC :
AIC(surg.fit.full, k = log(length(lny)))

##############################################################################