#########################
#
# File: 20111004.R
#
#       Prediction Graphs
#
#########################


# Preamble

library(RFS)


# Let's return to the GDP data set

gdp <- read.csv("http://courses.kvasaheim.com/pols6123/data/gdpcap.csv")
attach(gdp)
names(gdp)
summary(gdp)

# Let's create our model:

model <- lm(gdpcap ~ hig * govt)
summary(model)


# Change the base category for the govt variable

govt <- set.base(govt,"Democracy",gdp)
model <- lm(gdpcap ~ hig * govt)
summary(model)

govt <- set.base(govt,"Autocracy",gdp)
model <- lm(gdpcap ~ hig * govt)
summary(model)


# Plot the data

plot(hig,gdpcap)

# Create a slew of predictions:

newHIG <- seq(1,10, length=1000)
prDemo <- predict( model, newdata=data.frame(hig=newHIG, govt="Democracy" ) )
prAno  <- predict( model, newdata=data.frame(hig=newHIG, govt="Anocracy" ) )
prAuto <- predict( model, newdata=data.frame(hig=newHIG, govt="Autocracy" ) )

# Graph the predictions

lines(newHIG,prDemo, col=3, lwd=2)
lines(newHIG,prAno, col=2, lwd=2)
lines(newHIG,prAuto, col=4, lwd=2)

# Provide a legend

legend("topleft", bty="n", col=c(3,2,4), lwd=2, c("Democracy","Anocracy","Autocracy")  )



#### New example

# Use region in lieu of government type
# There are six regions, so the graph of all predictions may
# be a bit crowded. You have to decide that on your own


model <- lm(gdpcap~hig*region)
summary(model)

region <- set.base(region,"Other",gdp)
model <- lm(gdpcap~hig*region)
summary(model)


newHIG <- seq(1,10,length=1000)
prOth  <- predict( model, newdata=data.frame(hig=newHIG, region="Other") )
prAfr  <- predict( model, newdata=data.frame(hig=newHIG, region="Africa") )
prEast <- predict( model, newdata=data.frame(hig=newHIG, region="Eastern") )
prIsl  <- predict( model, newdata=data.frame(hig=newHIG, region="Islamic") )
prWest <- predict( model, newdata=data.frame(hig=newHIG, region="Western") )
prLA   <- predict( model, newdata=data.frame(hig=newHIG, region="Latin America") )

plot(hig,gdpcap, las=1, pch=16, col=8)
lines(newHIG,prOth, col=1, lwd=2)
lines(newHIG,prAfr, col=2, lwd=2)
lines(newHIG,prEast, col=3,lwd=2)
lines(newHIG,prIsl, col=4, lwd=2)
lines(newHIG,prWest, col=8, lwd=2)
lines(newHIG,prLA, col=6, lwd=2)

legend("topleft", col=c(8,3,6,4,2,1), lwd=2, 
	c("Western","Eastern","Latin American","Islamic","African","Other"), 
	bty="n")