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

save_plots = F

library(AER)
data(CPS1988)
attach(CPS1988)

# A general plotting function: 
#
plot_diagnostics = function(fitLm){
  
  par(mfrow=c(3,2))
  resid = rstudent(fitLm)
  plot(fitLm$fit,resid,ylim=c(-3,+3),main="(a)")
  lines(lowess(fitLm$fit,resid,f=0.2),lwd=5,col="red")
  abline(h=0,col="blue",lwd=5)
  
  plot(fitLm$fit,abs(resid),ylim=c(0,2.0),main="(b)")
  lines(lowess(fitLm$fit,abs(resid),f=0.2),lwd=5,col="red")
  abline(h=mean(abs(resid)),col="blue",lwd=5)

  qqnorm(resid,datax=F,main="(c)")
  qqline(resid,datax=F,lwd=5,col="blue")

  plot(education,resid,ylim=c(-5,5),main="(d)")
  lines(lowess(education,resid),lwd=5,col="red")
  abline(h=0,col="blue",lwd=5)
  
  plot(experience,resid,ylim=c(-3.0,+3.0),main="(e)")
  lines(lowess(experience,resid),lwd=5,col="red")
  abline(h=0,col="blue",lwd=5)

}

# The initial model:
# 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter13/cps1988_model_plots.eps", onefile=FALSE, horizontal=FALSE) }

fitLm1 = lm( wage ~ education + experience + ethnicity ) 
plot_diagnostics(fitLm1) 

if( save_plots ){ dev.off() }


# The log model:
# 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter13/cps1988_log_model_plots.eps", onefile=FALSE, horizontal=FALSE) }

fitLm2 = lm( log(wage) ~ education + experience )
plot_diagnostics(fitLm2) 

if( save_plots ){ dev.off() }


# The experience^2 model:
# 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter13/cps1988_log_model_w_experience2_plots.eps", onefile=FALSE, horizontal=FALSE) }

fitLm3 = lm( log(wage) ~ education + experience + I(experience^2) ) 
plot_diagnostics(fitLm3) 

if( save_plots ){ dev.off() }


# The model thus far when we add in ethnicity:
# 
fitLm4 = lm( log(wage) ~ education + experience + I(experience^2) + ethnicity ) 
plot_diagnostics(fitLm4) 


library(faraway)
par(mfrow=c(1,3))
plot(hatvalues(fitLm4))
plot(sqrt(cooks.distance(fitLm4)))
halfnorm(sqrt(cooks.distance(fitLm4)))