##############################
#
# Script: assignment08a.R
#   Answers to the eighth-
#   week homework assignment
#
##############################


# Read in the data

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


# Perform the four modelings

m.ls <- lm( gdpcap ~ democracy * region)
summary(m.ls)

m.ll <- lm( log(gdpcap) ~ democracy * region)
summary(m.ll)

m.gs <- glm( gdpcap ~ democracy * region)
summary(m.gs)

m.gl <- glm( gdpcap ~ democracy * region, family=gaussian(link=log))
summary(m.gl)



# Determine which is best

# Well, we know:

# m.gl better than m.gs (AIC)
# m.gs = m.ls
# m.ll better than m.ll (aR2)

# This means we know m.ll and m.gl are better than the 
# other two, but cannot be compared. Thus, we conclude 
# that the two log models are appropriate and optimal.

# At least from what we now know.




# Pediction graph for m.ll

png("clm.png",height=4,width=6,units="in",res=900)
par(mar=c(4,4,1,1)+0.3)
par(cex=0.8, cex.lab=0.8, cex.axis=0.8)
plot(democracy,gdpcap, ylab="GDP per Capita ($)", xlab="Level of Democracy", yaxt="n",ylim=c(0,50000))
axis(2, c(0,"10,000","50,000"), at=c(0,10000,50000), las=1 )
newDem <- seq(-10,10,length=1000)

prW <- exp( predict(m.ll, newdata=data.frame(democracy=newDem,region="Western")) )
prE <- exp( predict(m.ll, newdata=data.frame(democracy=newDem,region="Eastern")) )
prI <- exp( predict(m.ll, newdata=data.frame(democracy=newDem,region="Islamic")) )
prL <- exp( predict(m.ll, newdata=data.frame(democracy=newDem,region="Latin America")) )
prA <- exp( predict(m.ll, newdata=data.frame(democracy=newDem,region="Africa")) )
prO <- exp( predict(m.ll, newdata=data.frame(democracy=newDem,region="Other")) )

lines(newDem,prW, col=4, lwd=2)
lines(newDem,prE, col="orange", lwd=2)
lines(newDem,prI, col=3, lwd=2)
lines(newDem,prL, col=2, lwd=2)
lines(newDem,prA, col="purple", lwd=2)
lines(newDem,prO, col=1, lwd=2)

legend("topright",c("Latin America","Western","Islamic","Eastern","Africa","Other"),col=c(2,4,3,"orange","purple",1), lwd=2, bty="n",cex=0.8)
dev.off()




# Pediction graph for m.gl

png("glm.png",height=4,width=6,units="in",res=900)
par(mar=c(4,4,1,1)+0.3)
par(cex=0.8, cex.lab=0.8, cex.axis=0.8)
plot(democracy,gdpcap, ylab="GDP per Capita ($)", xlab="Level of Democracy", yaxt="n",ylim=c(0,50000))
axis(2, c(0,"10,000","50,000"), at=c(0,10000,50000), las=1 )
newDem <- seq(-10,10,length=1000)

prW <- exp( predict(m.gl, newdata=data.frame(democracy=newDem,region="Western")) )
prE <- exp( predict(m.gl, newdata=data.frame(democracy=newDem,region="Eastern")) )
prI <- exp( predict(m.gl, newdata=data.frame(democracy=newDem,region="Islamic")) )
prL <- exp( predict(m.gl, newdata=data.frame(democracy=newDem,region="Latin America")) )
prA <- exp( predict(m.gl, newdata=data.frame(democracy=newDem,region="Africa")) )
prO <- exp( predict(m.gl, newdata=data.frame(democracy=newDem,region="Other")) )

lines(newDem,prW, col=4, lwd=2)
lines(newDem,prE, col="orange", lwd=2)
lines(newDem,prI, col=3, lwd=2)
lines(newDem,prL, col=2, lwd=2)
lines(newDem,prA, col="purple", lwd=2)
lines(newDem,prO, col=1, lwd=2)

legend("topright",c("Latin America","Western","Islamic","Eastern","Africa","Other"),col=c(2,4,3,"orange","purple",1), lwd=2, bty="n",cex=0.8)
dev.off()