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

DF = load_exercise_3_7_data()

# Part (a):
#
m = lm( P ~ AGE, data=DF )

# Lets plot the data and models:
#
#postscript("../../WriteUp/Graphics/Chapter3/ex_3_7_part_a.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,3))
plot( m$fitted.values, m$residuals, type='p', pch=19, cex=1.5, xlab='y_hat', ylab='residual' )
grid()
plot(m, which=2)
grid()
plot(m$fitted.values, DF$P, xlab='y_hat', ylab='y')
grid()
par(mfrow=c(1,1))
#dev.off()


# Part (b):
#
library(MASS)

#postscript("../../WriteUp/Graphics/Chapter3/ex_3_7_part_b.eps", onefile=FALSE, horizontal=FALSE)
bc = boxcox(m)
#dev.off()

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

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

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


# Part (c):
#
DF$LOGP = log( DF$P )
m_log = lm( LOGP ~ AGE, data=DF )

DF$one_over_sqrt_P = 1/sqrt( DF$P )
m_one_over = lm( one_over_sqrt_P ~ AGE, data=DF )

# Make some predictions with different models:
p1 = predict( m, newdata=data.frame( AGE=2 ), level=0.95, interval="prediction" )
p = predict( m_log, newdata=data.frame( AGE=2 ), level=0.95, interval="prediction" )
p2 = exp( p )
p = predict( m_one_over, newdata=data.frame( AGE=2 ), level=0.95, interval="prediction" )
p3 = 1/(p^2)
t = p3[2]; p3[2] = p3[3]; p3[3] = t; # swap upper and lower bounds

# Put everything together:
#
R = rbind( p1, p2, p3 )
rownames(R) = 1:dim(R)[1]
R = data.frame( R )
R$pi_range = R$upr - R$lwr
print( R )


# Plot the data and the three models:
#
#postscript("../../WriteUp/Graphics/Chapter3/ex_3_7_part_c.eps", onefile=FALSE, horizontal=FALSE)
plot( DF$AGE, DF$P, type='p', pch=19, cex=1.5, xlab='AGE', ylab='P' )
all_ages = unique( DF$AGE )
points( all_ages, predict( m, newdata=data.frame( AGE=all_ages ) ), type='l', col='blue' )
points( all_ages, exp( predict( m_log, newdata=data.frame( AGE=all_ages ) ) ), type='l', col='purple' )
p = predict( m_one_over, newdata=data.frame( AGE=all_ages ) )
points( all_ages, 1/(p^2), type='l', col='green' )
legend( 'topright',
        c('P = beta_0 + beta_1 AGE', 'P = exp(beta_0 + beta_1 AGE)', 'P=1/sqrt(beta_0 + beta_1 AGE)'),
        col=c('blue', 'purple', 'green'),
        lty=c(1, 1, 1))
grid()
#dev.off()