###################################################
### chunk number 1: ss1
###################################################
#line 173 "Lecture7_2010.Rnw"
library(pspline)
library(foreign)
jacob <- read.dta("jacob.dta")
perot <- with(jacob, perotvote - min(perotvote))
perot <- perot/max(perot)
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, spar=0), lwd=3, col="red"))
with(jacob, sm.spline(perot, chal_vote, spar=0))
with(jacob, sm.spline(perot, chal_vote, spar=0))
 
 
###################################################
### chunk number 2: ss2
###################################################
#line 183 "Lecture7_2010.Rnw"
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, spar=0.1), lwd=3, col="red"))
 
 
###################################################
### chunk number 3: ss3
###################################################
#line 188 "Lecture7_2010.Rnw"
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, spar=1), lwd=3, col="red"))
 
 
###################################################
### chunk number 4: ss4
###################################################
#line 192 "Lecture7_2010.Rnw"
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, spar=10), lwd=3, col="red"))
 
 
###################################################
### chunk number 5: ss1a
###################################################
#line 304 "Lecture7_2010.Rnw"
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, df=2), lwd=3, col="red"))
 
 
###################################################
### chunk number 6: ss2a
###################################################
#line 309 "Lecture7_2010.Rnw"
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, df=5), lwd=3, col="red"))
 
 
###################################################
### chunk number 7: ss3a
###################################################
#line 314 "Lecture7_2010.Rnw"
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, df=10), lwd=3, col="red"))
 
 
###################################################
### chunk number 8: ss4a
###################################################
#line 318 "Lecture7_2010.Rnw"
with(jacob, plot(perot, chal_vote, cex=.5))
with(jacob, lines(sm.spline(perot, chal_vote, df=20), lwd=3, col="red"))
 
 
###################################################
### chunk number 9: prestgam
###################################################
#line 506 "Lecture7_2010.Rnw"
library(mgcv)
library(car)
data(Prestige)
prestige.gam <- gam(prestige ~ s(income) + s(education), data=Prestige)
summary(prestige.gam)
 
prestige.gam1 <- gam(prestige ~ log(income) + s(education, k=4, fx=T), data=Prestige)
prestige.gam1a <- gam(prestige ~ s(income) + s(education, k=4, fx=T), data=Prestige)
prestige.gam1a <- gam(prestige ~ income + education, data=Prestige)
 
 
 
 
###################################################
### chunk number 10: visgam1
###################################################
#line 522 "Lecture7_2010.Rnw"
vis.gam(prestige.gam, theta=-40)
 
 
###################################################
### chunk number 11: visgam2
###################################################
#line 541 "Lecture7_2010.Rnw"
vis.gam(prestige.gam, se=2, theta=-20)
 
## Added in class - 3d confidence bounds (rotatable)
eg <- expand.grid(
	income = seq.range(Prestige$income), 
	education = seq.range(Prestige$education))
 
preds <- predict(prestige.gam, eg)
preds2 <- predict(prestige.gam, eg, se.fit=T)
lower <- preds2$fit - 2*preds2$se.fit
upper <- preds2$fit + 2*preds2$se.fit
 
newdat <- cbind(rbind(eg, eg, eg), c(preds2$fit, lower, upper))
newdat$fac <- factor(rep(1:3, each=625), levels=1:3, labels=
	c("fit", "lower", "upper")
 
scatter3d(eg[,1], preds, eg[,2], fit="smooth", df.smooth=100)
scatter3d(newdat[,1], newdat[,3], newdat[,2], fit="smooth", df.smooth=100,
groups = newdat[,4], parallel=F)
 
 
prestige.gam2 <- gam(prestige ~ s(income) + s(education) + s(women), data=Prestige)
vis.gam(prestige.gam2, view=c("income", "women"))
 
 
###################################################
### chunk number 12: twofigs1
###################################################
#line 563 "Lecture7_2010.Rnw"
ed.seq <- with(Prestige, seq(min(education), max(education), 
	length=25))
ed.pred <- predict(prestige.gam, newdata=data.frame(
	education=ed.seq, income=mean(Prestige$income)), 
	se.fit=T)
lower <- with(ed.pred, fit - 2*se.fit)
upper <- with(ed.pred, fit + 2*se.fit)
plot(ed.seq, ed.pred$fit, ylim=range(c(lower,upper)), 
	xlab="Education", ylab="Predicted Prestige", type="l")
lines(ed.seq, lower, lty=2)
lines(ed.seq, upper, lty=2)
with(Prestige, rug(education))
 
 
###################################################
### chunk number 13: twofigs2
###################################################
#line 585 "Lecture7_2010.Rnw"
inc.seq <- with(Prestige, seq(min(income), max(income), 
	length=25))
inc.pred <- predict(prestige.gam, 
	newdata=data.frame(education=mean(Prestige$education), 
	income= inc.seq), se.fit=T)
lower <- with(inc.pred, fit - 2*se.fit)
upper <- with(inc.pred, fit + 2*se.fit)
plot(inc.seq, inc.pred$fit, ylim=range(c(lower,upper)), 
	xlab="Income", ylab="Predicted Prestige", type="l")
lines(inc.seq, lower, lty=2)
lines(inc.seq, upper, lty=2)
with(Prestige, rug(income))
 
 
###################################################
### chunk number 14: lintest
###################################################
#line 695 "Lecture7_2010.Rnw"
library(xtable)
prestige.ols<-gam(prestige~income+education, data=Prestige)
xtable(anova(prestige.ols, prestige.gam, test="F"))
 
 
###################################################
### chunk number 15: gamcheck
###################################################
#line 716 "Lecture7_2010.Rnw"
gam.check(prestige.gam)
 
 
###################################################
### chunk number 16: gaminter
###################################################
#line 781 "Lecture7_2010.Rnw"
type.bc<-with(Prestige, as.numeric(type=="bc"))
type.prof<-with(Prestige, as.numeric(type=="prof"))
type.wc<-with(Prestige, as.numeric(type=="wc"))
inter.gam<-gam(prestige~type+s(income,by=type.bc)
        +s(income,by=type.prof)+s(income,by=type.wc), 
		data=Prestige)
summary(inter.gam)
 
 
###################################################
### chunk number 17: bcinc
###################################################
#line 798 "Lecture7_2010.Rnw"
plot.gam(inter.gam, select=1)
title("Blue Collar")
 
 
###################################################
### chunk number 18: profinc
###################################################
#line 814 "Lecture7_2010.Rnw"
plot.gam(inter.gam, select=2)
title("Professional")
 
 
###################################################
### chunk number 19: wcinc
###################################################
#line 830 "Lecture7_2010.Rnw"
plot.gam(inter.gam, select=3)
title("White Collar")
 
 
###################################################
### chunk number 20: demogam
###################################################
#line 862 "Lecture7_2010.Rnw"
weakliem <- read.table("~/Desktop/ICPSR_Slides/Lecture 1/weakliem.txt")
nondemo <- with(weakliem, as.numeric(democrat == 0))
demo.gam<-gam(secpay~democrat + s(gini,by=democrat)
        +s(gini,by=nondemo), data=weakliem)
 
 
###################################################
### chunk number 21: sumdemogam
###################################################
#line 873 "Lecture7_2010.Rnw"
summary(demo.gam)
 
 
###################################################
### chunk number 22: dgfig1
###################################################
#line 881 "Lecture7_2010.Rnw"
plot(demo.gam, select=1)
title("Inequality in Democracies")
 
 
###################################################
### chunk number 23: dgfig2
###################################################
#line 893 "Lecture7_2010.Rnw"
plot(demo.gam, select=2)
title("Inequality in Non-democracies")
 
 
###################################################
### chunk number 24: demolin
###################################################
#line 907 "Lecture7_2010.Rnw"
demo.lm<-gam(secpay~democrat*gini, data=weakliem)
xtable(anova.gam(demo.lm, demo.gam))
 
 
###################################################
### chunk number 25: anes
###################################################
#line 1043 "Lecture7_2010.Rnw"
nes <- read.dta("anes1992.dta")
mod <- gam(votedem ~ s(pid, k=7) + s(age) + female + educ + 
	black + retnat + south, family=binomial, data=nes)
 
 
###################################################
### chunk number 26: modsum
###################################################
#line 1053 "Lecture7_2010.Rnw"
summary(mod)
 
 
###################################################
### chunk number 27: pid
###################################################
#line 1064 "Lecture7_2010.Rnw"
plot(mod, select=1)
title("Effect of Party ID")
 
 
###################################################
### chunk number 28: age
###################################################
#line 1077 "Lecture7_2010.Rnw"
plot(mod, select=2)
title("Effect of Age")
 
 
###################################################
### chunk number 29: agetest
###################################################
#line 1091 "Lecture7_2010.Rnw"
moda <- gam(votedem ~ s(pid, k=7, fx=T) + s(age) + female + educ + 
	black + retnat + south, family=binomial, data=nes)
agelin <- gam(votedem ~ s(pid, k=7, fx=T) + age + female + educ + 
	black + retnat + south, family=binomial, data=nes)
xtable(anova(moda, agelin, test='Chisq'))
 
 
###################################################
### chunk number 30: pidlintest
###################################################
#line 1128 "Lecture7_2010.Rnw"
mod2 <- gam(votedem ~ s(pid, k=7, fx=F) + age + female + educ + 
	black + retnat + south, family=binomial, data=nes)
piddum <- 	gam(votedem ~ as.factor(pid) + age + female + educ + 
		black + retnat + south, family=binomial, data=nes)
xtable(anova(mod2, piddum, test='Chisq'))
 
 
###################################################
### chunk number 31: pidpreds
###################################################
#line 1148 "Lecture7_2010.Rnw"
newdata <- with(nes, data.frame(pid = 1:7, age = mean(age, na.rm=T), 
	female = 0, educ = mean(educ, na.rm=T), black = 0, 
	retnat = factor(2, levels=1:3, labels=levels(retnat)), south = 0))
preds <- predict(mod2, newdata=newdata, se.fit=T, type="link")
lower <- plogis(with(preds, fit - 2*se.fit))
upper <- plogis(with(preds, fit + 2*se.fit))
plot(plogis(preds$fit) ~ newdata$pid, type="l", xlab="Party ID", 
	ylab="Pr(Vote Democrat)", ylim=c(min(lower), max(upper)))
lines(newdata$pid, lower, lty=2)
lines(newdata$pid, upper, lty=2)