source('../../Data/data_loaders.R') DF = load_exercise_3_8_data() # Part (a): # m = lm( Y ~ X, data=DF ) # Lets plot the data and models: # #postscript("../../WriteUp/Graphics/Chapter3/ex_3_8_part_a.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1,4)) plot( DF\$X, DF\$Y, type='p', pch=19, cex=1.5, xlab='X', ylab='Y', main='data and linear model' ) abline(m, col='red') grid() plot( m\$fitted.values, m\$residuals, type='p', pch=19, cex=1.5, xlab='y_hat', ylab='residual', main='residual vs. fitted' ) grid() plot( m, which=2 ) grid() plot( m\$fitted.values, DF\$Y, xlab='y_hat', ylab='y', main='y vs. y_hat' ) 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, 'Y', c( 'X' ) ) print(sprintf("Atkinsons lambda estimate: %10.6f (t-value %10.6f)", atk_res\$lambda_hat, atk_res\$gamma_t)) # Part (c): # DF\$one_over_Y = 1/DF\$Y m_one_over = lm( one_over_Y ~ X, data=DF ) # Make some predictions with different models: p1 = predict( m, newdata=data.frame( X=3 ), level=0.95, interval="prediction" ) p = predict( m_one_over, newdata=data.frame( X=3 ), level=0.95, interval="prediction" ) p2 = 1/p t = p2[2]; p2[2] = p2[3]; p2[3] = t; # swap upper and lower bounds # Put everything together: # R = rbind( p1, p2 ) 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_8_part_c.eps", onefile=FALSE, horizontal=FALSE) plot( DF\$X, DF\$Y, type='p', pch=19, cex=1.5, xlab='X', ylab='Y' ) all_Xs = unique( DF\$X ) points( all_Xs, predict( m, newdata=data.frame( X=all_Xs ) ), type='l', col='blue' ) p = predict( m_one_over, newdata=data.frame( X=all_Xs ) ) points( all_Xs, 1/p, type='l', col='green' ) legend( 'topleft', c('Y = beta_0 + beta_1 X', 'Y=1/(beta_0 + beta_1 X)'), col=c('blue', 'green'), lty=c(1, 1) ) grid() #dev.off()