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

if(is.null(version$language) == FALSE){ 
  require(alr3)
  require(car)
}else{
  library(alr3)
  library(car)
}

data(drugcost)
attach(drugcost)

postscript("../../WriteUp/Graphics/Chapter9/prob_12_scatter_plot_matrix.eps", onefile=FALSE, horizontal=FALSE)
pairs( ~GS+RI+F+AGE+RXPM+COPAY+MM+COST ) 
dev.off()

cor( data.frame( COST=COST, RXPM=RXPM, GS=GS, RI=RI, COPAY=COPAY, AGE=AGE, F=F, MM=MM ) )

m0 <- lm( COST ~ RXPM+GS+RI+COPAY+AGE+F+MM )
summary(m0) 

postscript("../../WriteUp/Graphics/Chapter9/prob_12_avp.eps", onefile=FALSE, horizontal=FALSE)
avp( m0, ask=FALSE, one.page=TRUE, identify.points=FALSE)
dev.off()

# lets see if we can remove RI:
m1 <- lm( COST ~ RXPM+GS+COPAY+AGE+F+MM )
summary(m1)

anova(m0,m1) # yes we can remove it ...

# can we remove MM:
m2 <- lm( COST ~ RXPM+GS+COPAY+AGE+F )
summary(m2)

anova(m1,m2) # yes we can remove it ...

# can we remove F:
m3 <- lm( COST ~ RXPM+GS+COPAY+AGE )
summary(m3)

anova(m2,m3) # lets say yes ...

# can we remove COPAY ...
m4 <- lm( COST ~ RXPM+GS+AGE )
summary(m4)

anova(m3,m4) # lets say yes ...

# can we remove AGE:
m5 <- lm( COST ~ RXPM+GS+COPAY )
summary(m5)

anova(m4,m5) # not really ... 

# what is the total loss from removing all variables ... still looks like m4 is a valid model i.e. we can remove all of these terms 
anova(m0,m4)

m4p <- update( m4, ~ . + RI ) # lets add RI back in and see if we are "sure" of the sign of its coefficient ...
summary(m4p)