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

DF = load_appendix_deforestation_data()
DF$DEF[ DF$DEF == min(DF$DEF) ] = min(DF$DEF[ DF$DEF>0 ])/2 # boxcox cannot handle a zero response

# Part (a):
#
m = lm( DEF ~ POP + GNP, data=DF )
print( summary(m) )

# Part (b):
#
#postscript("../../WriteUp/Graphics/Chapter4/ex_4_7_residuals_vs_fitted.eps", onefile=FALSE, horizontal=FALSE)
plot( m, which=1 )
#dev.off()

# Part (c):
#
library(MASS)

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

source('../Chapter3/utils.R')

atk_res = Atkinsons( DF, 'DEF', c('POP', 'GNP') )
print(sprintf("Atkinsons lambda estimate: %.3f (t-value %.3f)", atk_res$lambda_hat, atk_res$gamma_t))

bt_res = Box_Tidwell( DF, 'DEF', c('POP', 'GNP') )
print(bt_res)

# Part (d):
#
DF$DEFT = 1/sqrt(DF$DEF) # transformed DEF
DF$LGNP = log (DF$GNP)

m_tran = lm( DEFT ~ POP + LGNP, data=DF ) # the suggested transformations
print( summary( m_tran ) )


# Part (e):
#
cm = data.frame( t(colMeans(DF)) )
cm = cm[ c('POP', 'GNP') ] # just select the two inputs
print(cm)
cm$LGNP = log(cm$GNP)

p = predict( m, newdata=cm, interval='prediction' )
p_tran = predict( m_tran, newdata=data.frame( POP=cm['POP'], LGNP=cm['LGNP'] ), interval='prediction' )
p_tran = (1/p_tran)^2

lwr = p_tran[3]
upr = p_tran[2]
p_tran[2] = lwr
p_tran[3] = upr
R = rbind( p, p_tran )
rownames(R) = c('orig_model', 'tran_model')
R = cbind( R, range=(R[,'upr'] - R[,'lwr']) )
print(R)

# Part (f):
#
#postscript("../../WriteUp/Graphics/Chapter4/ex_4_7_residuals_vs_fitted_trans_model.eps", onefile=FALSE, horizontal=FALSE)
plot( m_tran, which=1 )
#dev.off()