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