#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

library(MASS)

head(cement)

#postscript("../../WriteUp/Graphics/Chapter6/prob_2_pairs_plot.eps", onefile=FALSE, horizontal=FALSE)
pairs(cement)
#dev.off()

# Form an initial model 
m0 = lm( y ~ ., data=cement )

# Look at what the AIC recommends removing: 
stepAIC( m0 )
m1 = lm( y ~ x1 + x2 + x4, data=cement )
summary(m1)

m2 = lm( y ~ x1 + x2 + x3, data=cement )
summary(m2)

# How well does this linear model fit?
par(mfrow=c(1,4))
plot(m2)
par(mfrow=c(1,1))

summary( lm( y ~ x1 + x2, data=cement ) ) # remove x3 

# Lets try the suggested transformation:
#
log_cement = cement
log_cement$x1 = log( cement$x1 / (100 - cement$x1 ) )
log_cement$x2 = log( cement$x2 / (100 - cement$x2 ) )
log_cement$x3 = log( cement$x3 / (100 - cement$x3 ) )
log_cement$x4 = log( cement$x4 / (100 - cement$x4 ) )

pairs( log_cement ) 

# Again x2 and x4 are highly correlated
log_cement.lm = lm( y ~ x1 + x2 + x3, data=log_cement )

# Summary plots indicate that we have more outliers with the log transformed model