#
# Epage 
#
# 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.
# 
#-----

# version$language == "R" for R version$language == NULL for SPlus

if(is.null(version$language) == FALSE){ 
  require(alr3)
}else{
  library(alr3)
}

data(wool)
attach(wool)

# 7.6.1:
# 
postscript("../../WriteUp/Graphics/Chapter7/prob_6_orig_scatter_plot.eps", onefile=FALSE, horizontal=FALSE)
pairs( ~Len+Amp+Load+Cycles ) 
dev.off()

# 7.6.2 create a model that has two parts: a main effects part (linear factors/terms) and the interaction terms
# 
# see epage 144
# 
lenF <- factor( Len, ordered=FALSE ) 
ampF <- factor( Amp, ordered=FALSE )
loadF <- factor( Load, ordered=FALSE )

m2 <- lm( Cycles ~ lenF + ampF + loadF + lenF:ampF + lenF:loadF + ampF:loadF ) 
summary(m2)


# 7.6.3: create a main effects model to predict Cycles no interaction terms and
#        use an inverse fitted plot to determine a tranformation of Cycles 
# 
m1 <- lm( Cycles ~ lenF + ampF + loadF )
summary(m1)

anova(m2,m1) # <- shows that there is a statistically significant reduction in variance in using the more complicated model.

ans <- bctrans( Cycles ~ Len + Amp + Load )
summary(ans)
#plot(ans)

woolT <- data.frame( logLen=log(Len), logAmp=log(Amp), logLoad=log(Load), Cycles=Cycles )

# look for a tranformation of Cycles: 
mt1 <- lm( Cycles ~ logLen + logAmp + logLoad, data=woolT ) 

cycPower <- inv.tran.plot( Cycles, predict(mt1) )

unlist( inv.tran.estimate( Cycles, predict(mt1) ) # <- gives the standard error 

woolTT <- data.frame( logLen=log(Len), logAmp=log(Amp), logLoad=log(Load), logCycles=log(Cycles) )

mt2 <- lm( logCycles ~ logLen + logAmp + logLoad, data=woolTT )    
       
# 7.6.4: refit the interactions model from 7.6.2 for the responce log(Cycles) and
#        now show that none of the interaction terms are needed ...
#
# I'm a bit confused as to how to do this with factors ... 
#
logLenF <- factor(woolT$logLen,ordered=FALSE)
logAmpF <- factor( woolT$logAmp, ordered=FALSE )
logLoadF <- factor( woolT$logLoad, ordered=FALSE )

mp1 <- lm( log(Cycles) ~ logLenF + logAmpF + logLoadF + logLenF:logAmpF + logLenF:logLoadF + logAmpF:logLoadF ) 
summary(mp1)
       
summary( lm( log(Cycles) ~ logLenF + logAmpF + logLoadF )  )