source('../../Data/data_loaders.R')

DF = load_appendix_polymer_data()

# Part (a):
#
m = lm( IV ~ MW, data=DF )

# Lets plot the data and models:
#
#postscript("../../WriteUp/Graphics/Chapter3/ex_3_10_part_a.eps", onefile= FALSE, horizontal=FALSE)
par(mfrow=c(1, 3))
plot( DF$MW, DF$IV, type='p', pch=19, cex=1.5, xlab='MW', ylab='IV', main='data and linear model' )
abline(m, col='red')
grid()
plot(m, which=1)
grid()
plot(m, which=2)
grid()
par(mfrow=c(1,1))
#dev.off()

# Part (b):
#
library(MASS)
bc = boxcox(m)

# What lambda value gives the largest log-likelihood:
#
print(sprintf("Boxcox lambda estimate: %10.6f", bc$x[ which.max(bc$y)]))

# Lets try Atkinsons method:
#
source('utils.R')

atk_res = Atkinsons( DF, 'IV', c( 'MW' ) )
print(sprintf("Atkinsons lambda estimate: %10.6f (t-value %10.6f)", atk_res$lambda_hat, atk_res$gamma_t))


# Part (c):
#
DF$LOGIV = log(DF$IV)
m_log = lm( LOGIV ~ MW, data=DF )

#postscript("../../WriteUp/Graphics/Chapter3/ex_3_10_part_c.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,3))
plot(DF$MW, DF$LOGIV, type='p', pch=19, cex=1.5, xlab='MW', ylab='log(IV)', main='data and linear model' )
abline(m_log, col='red')
grid()
plot(m_log, which=1 )
grid()
plot(m_log, which=2 )
grid()
par(mfrow=c(1,1))
#dev.off()

# Part (d): Implement the Box-Tidwell method to study transformation of MW
#
bt_res = Box_Tidwell( DF, 'LOGIV', c( 'MW' ) ) 
print(sprintf("Box-Tidwell %s alpha estimate: %.3f (t-value %.3f)", bt_res$variable_name[1], bt_res$alpha[1], bt_res$t_value[1]))

# Build this model (and look at the regression diagnostics for this model):
#
m_bt = lm( LOGIV ~ I(MW^(-0.1)), data=DF )
#plot(m_bt)


# Part (e):
#
DF$LOGMW = log(DF$MW)
m_log_log = lm( LOGIV ~ LOGMW, data=DF )

#postscript("../../WriteUp/Graphics/Chapter3/ex_3_10_part_e.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,3))
plot(DF$LOGMW, DF$LOGIV, type='p', pch=19, cex=1.5, xlab='log(MW)', ylab='log(IV)', main='data and linear model' )
abline(m_log_log, col='red')
grid()
plot(m_log_log, which=1 )
grid()
plot(m_log_log, which=2 )
grid()
par(mfrow=c(1,1))
#dev.off()


#
# As an aside, this is the model suggested by the box-cox transformation above:
#
DF$IV_to_LAMBDA = DF$IV^1.4
m_alt = lm( IV_to_LAMBDA ~ MW, data=DF )

#postscript("../../WriteUp/Graphics/Chapter3/ex_3_10_part_e_aux.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,3))
plot(DF$MW, DF$IV_to_LAMBDA, type='p', pch=19, cex=1.5, xlab='MW', ylab='IVA(1ambda)', main='data and linear model' )
abline(m_alt, col='red')
grid()
plot(m_alt, which=1)
grid()
plot(m_alt, which=2)
grid()
par(mfrow=c(1,1))
#dev.off()


# Lets see if we then want to apply a power transformation to MW
#
bt_res = Box_Tidwell( DF, 'IV_to_LAMBDA', c( 'MW' ) )
print(sprintf("Box-Tidwell %s alpha estimate: %.3f (t-value %.3f)", bt_res$variable_name[1], bt_res$alpha[1], bt_res$t_value[1]))