###################################################
### chunk number 1: ICV1
###################################################
#line 278 "Lecture1_2011.Rnw"
library(apsrtable)
setwd("~/desktop/ICPSR_slides/Lecture 1/")
wvs<-read.table('Weakliem.txt', header=T) 
with(wvs, plot(secpay ~ gini, cex=1.5, xlab="GINI", ylab="Attitudes toward Inequality (mean)"))
outs <- which(rownames(wvs)%in% c("Slovakia","CzechRepublic","Chile"))
with(wvs, text(gini[outs], secpay[outs], rownames(wvs)[outs], pos=c(2,4,4)))
abline(lm(secpay ~ gini, data=wvs), lty=2, col="red", lwd=1.5)
with(wvs, lines(lowess(gini, secpay, f=.5), col="red", lwd=1.5))
legend("topright", c("Linear Model","LOWESS"), lty=c(2,1), col=c("red","red"), xjust=1, inset=.01)
 
 
###################################################
### chunk number 2: fig1
###################################################
#line 291 "Lecture1_2011.Rnw"
#line 278 "Lecture1_2011.Rnw"
library(apsrtable)
setwd("~/desktop/ICPSR_slides/Lecture 1/")
wvs<-read.table('Weakliem.txt', header=T) 
with(wvs, plot(secpay ~ gini, cex=1.5, xlab="GINI", ylab="Attitudes toward Inequality (mean)"))
outs <- which(rownames(wvs)%in% c("Slovakia","CzechRepublic","Chile"))
with(wvs, text(gini[outs], secpay[outs], rownames(wvs)[outs], pos=c(2,4,4)))
abline(lm(secpay ~ gini, data=wvs), lty=2, col="red", lwd=1.5)
with(wvs, lines(lowess(gini, secpay, f=.5), col="red", lwd=1.5))
legend("topright", c("Linear Model","LOWESS"), lty=c(2,1), col=c("red","red"), xjust=1, inset=.01)
#line 292 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 3: ICVmod1
###################################################
#line 296 "Lecture1_2011.Rnw"
mod <- lm(secpay ~ gini, data=wvs)
apsrtable(mod, model.names="", Sweave=TRUE)
 
 
###################################################
### chunk number 4: ICV2
###################################################
#line 318 "Lecture1_2011.Rnw"
with(wvs, plot(secpay ~ gini, cex=1.5, xlab="GINI", ylab="Attitudes toward Inequality (mean)", type="n"))
with(wvs, points(secpay[democrat == 0] ~ gini[democrat == 0], col="red", pch=1, cex=1.5))
with(wvs, points(secpay[democrat == 1] ~ gini[democrat == 1], col="blue", pch=2, cex=1))
with(wvs, text(gini[outs], secpay[outs], rownames(wvs)[outs], pos=c(2,4,4)))
abline(lm(secpay ~ gini, data=wvs, subset=democrat==0), col="red", lty=2, lwd=1.5)
abline(lm(secpay ~ gini, data=wvs, subset=democrat==1), col="blue", lty=2, lwd=1.5)
legend("topright", c("Non-Democracy","Democracy"), pch=c(1,2), col=c("red","blue"), xjust=1, inset=.01)
 
 
###################################################
### chunk number 5: fig2
###################################################
#line 329 "Lecture1_2011.Rnw"
#line 318 "Lecture1_2011.Rnw"
with(wvs, plot(secpay ~ gini, cex=1.5, xlab="GINI", ylab="Attitudes toward Inequality (mean)", type="n"))
with(wvs, points(secpay[democrat == 0] ~ gini[democrat == 0], col="red", pch=1, cex=1.5))
with(wvs, points(secpay[democrat == 1] ~ gini[democrat == 1], col="blue", pch=2, cex=1))
with(wvs, text(gini[outs], secpay[outs], rownames(wvs)[outs], pos=c(2,4,4)))
abline(lm(secpay ~ gini, data=wvs, subset=democrat==0), col="red", lty=2, lwd=1.5)
abline(lm(secpay ~ gini, data=wvs, subset=democrat==1), col="blue", lty=2, lwd=1.5)
legend("topright", c("Non-Democracy","Democracy"), pch=c(1,2), col=c("red","blue"), xjust=1, inset=.01)
#line 330 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 6: ICVmod2
###################################################
#line 334 "Lecture1_2011.Rnw"
mod2 <- lm(secpay ~ gini*democrat, data=wvs)
apsrtable(mod2, model.names="", Sweave=TRUE)
 
 
###################################################
### chunk number 7: ICV3
###################################################
#line 356 "Lecture1_2011.Rnw"
with(wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], plot(secpay ~ gini, cex=1.5, xlab="GINI", ylab="Attitudes toward Inequality (mean)", type="n"))
with(wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], points(secpay[democrat == 0] ~ gini[democrat == 0], col="red", pch=1, cex=1.5))
with(wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], points(secpay[democrat == 1] ~ gini[democrat == 1], col="blue", pch=2, cex=1))
abline(lm(secpay ~ gini, data=wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], subset=democrat==0), col="red", lty=2, lwd=1.5)
abline(lm(secpay ~ gini, data=wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], subset=democrat==1), col="blue", lty=2, lwd=1.5)
legend("topright", c("Non-Democracy","Democracy"), pch=c(1,2), col=c("red","blue"), xjust=1, inset=.01)
 
 
###################################################
### chunk number 8: fig3
###################################################
#line 366 "Lecture1_2011.Rnw"
#line 356 "Lecture1_2011.Rnw"
with(wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], plot(secpay ~ gini, cex=1.5, xlab="GINI", ylab="Attitudes toward Inequality (mean)", type="n"))
with(wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], points(secpay[democrat == 0] ~ gini[democrat == 0], col="red", pch=1, cex=1.5))
with(wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], points(secpay[democrat == 1] ~ gini[democrat == 1], col="blue", pch=2, cex=1))
abline(lm(secpay ~ gini, data=wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], subset=democrat==0), col="red", lty=2, lwd=1.5)
abline(lm(secpay ~ gini, data=wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),], subset=democrat==1), col="blue", lty=2, lwd=1.5)
legend("topright", c("Non-Democracy","Democracy"), pch=c(1,2), col=c("red","blue"), xjust=1, inset=.01)
#line 367 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 9: ICVmod3
###################################################
#line 371 "Lecture1_2011.Rnw"
mod3 <- lm(secpay ~ gini*democrat, data=wvs[-which(rownames(wvs) %in% c("CzechRepublic","Slovakia")),])
apsrtable(mod3, model.names="", Sweave=TRUE)
 
 
###################################################
### chunk number 10: nl1
###################################################
#line 423 "Lecture1_2011.Rnw"
set.seed(123)
y <- seq(from=1,to=100, length=100)
x <- log(y) + rnorm(100, 0, 0.4)
plot(x ~ y, axes=F, xlab="x", ylab="y", cex=1.5)
box()
lines(y, log(y))
 
 
###################################################
### chunk number 11: nonlin1
###################################################
#line 431 "Lecture1_2011.Rnw"
#line 423 "Lecture1_2011.Rnw"
set.seed(123)
y <- seq(from=1,to=100, length=100)
x <- log(y) + rnorm(100, 0, 0.4)
plot(x ~ y, axes=F, xlab="x", ylab="y", cex=1.5)
box()
lines(y, log(y))
#line 432 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 12: nl2
###################################################
#line 434 "Lecture1_2011.Rnw"
set.seed(123)
x <- seq(from=-100,to=100, length=100)
y <- x + 2*x^2 + 3*x^3 
plot(y + rnorm(100, 0, 250000) ~ x, axes=F, xlab="x", ylab="y", cex=1.5)
box()
lines(y ~ x, lwd=2)
 
 
###################################################
### chunk number 13: nonlin2
###################################################
#line 442 "Lecture1_2011.Rnw"
#line 434 "Lecture1_2011.Rnw"
set.seed(123)
x <- seq(from=-100,to=100, length=100)
y <- x + 2*x^2 + 3*x^3 
plot(y + rnorm(100, 0, 250000) ~ x, axes=F, xlab="x", ylab="y", cex=1.5)
box()
lines(y ~ x, lwd=2)
#line 443 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 14: nl3
###################################################
#line 445 "Lecture1_2011.Rnw"
set.seed(123)
x <- seq(from=-100,to=100, length=100)
y <- 20 + 50*x + x^2
plot(y+rnorm(100, 0, 1000) ~ x, axes=F, xlab="x", ylab="y", cex=1.5)
box()
lines(y ~ x, lwd=2)
 
 
###################################################
### chunk number 15: nonlin3
###################################################
#line 453 "Lecture1_2011.Rnw"
#line 445 "Lecture1_2011.Rnw"
set.seed(123)
x <- seq(from=-100,to=100, length=100)
y <- 20 + 50*x + x^2
plot(y+rnorm(100, 0, 1000) ~ x, axes=F, xlab="x", ylab="y", cex=1.5)
box()
lines(y ~ x, lwd=2)
#line 454 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 16: regspline
###################################################
#line 513 "Lecture1_2011.Rnw"
plot(c(0,10), c(0, 10), type="n", axes=F, xlab="", ylab="")
segments(0,0,0,10)
segments(0,0,10,0)
text(0.25, 9,"Election", pos=4, cex=2)
text(9, 9,"Election", pos=2, cex=2)
segments(9,8.5,6,4, lwd=2)
segments(0,8.5,6,4, lwd=2)
mtext(expression(X[0]), 1, at=6, cex=1.5)
mtext("Public Support", 2, at=5, cex=1.5)
mtext("Time X", 1, at=9.5, cex=1.5)
segments(6,1, 6,4, lty=3)
text(6,.5,"Campaign Begins", cex=1.5)
mtext("0", 1, at=0, line=0, cex=1.25)
text(6,3.8,"knot", pos=4, cex=1.25)
 
 
###################################################
### chunk number 17: splines1
###################################################
#line 548 "Lecture1_2011.Rnw"
#line 513 "Lecture1_2011.Rnw"
plot(c(0,10), c(0, 10), type="n", axes=F, xlab="", ylab="")
segments(0,0,0,10)
segments(0,0,10,0)
text(0.25, 9,"Election", pos=4, cex=2)
text(9, 9,"Election", pos=2, cex=2)
segments(9,8.5,6,4, lwd=2)
segments(0,8.5,6,4, lwd=2)
mtext(expression(X[0]), 1, at=6, cex=1.5)
mtext("Public Support", 2, at=5, cex=1.5)
mtext("Time X", 1, at=9.5, cex=1.5)
segments(6,1, 6,4, lty=3)
text(6,.5,"Campaign Begins", cex=1.5)
mtext("0", 1, at=0, line=0, cex=1.25)
text(6,3.8,"knot", pos=4, cex=1.25)
#line 549 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 18: lowfig
###################################################
#line 640 "Lecture1_2011.Rnw"
prestige <- read.table("prestige.txt", header=T)
with(prestige, plot(prestige ~ income, cex=1.5, xlab ="Pineo-Porter Prestige Score", ylab="Average Income"))
with(prestige, lines(lowess(income, prestige, f=.7), col="red"))
with(prestige, lines(lowess(income, prestige, f=.2), col="blue"))
abline(lm(prestige ~ income, data=prestige))
legend("bottomright", c("Linear","Lowess (f=.7)","Lowess (f=.2)"), lty=c(1,2,3), 
    lwd=c(1.5,1.5,1.5), col=c("black","red","blue"), inset=.01)
 
 
###################################################
### chunk number 19: lowess1
###################################################
#line 649 "Lecture1_2011.Rnw"
#line 640 "Lecture1_2011.Rnw"
prestige <- read.table("prestige.txt", header=T)
with(prestige, plot(prestige ~ income, cex=1.5, xlab ="Pineo-Porter Prestige Score", ylab="Average Income"))
with(prestige, lines(lowess(income, prestige, f=.7), col="red"))
with(prestige, lines(lowess(income, prestige, f=.2), col="blue"))
abline(lm(prestige ~ income, data=prestige))
legend("bottomright", c("Linear","Lowess (f=.7)","Lowess (f=.2)"), lty=c(1,2,3), 
    lwd=c(1.5,1.5,1.5), col=c("black","red","blue"), inset=.01)
#line 650 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 20: smsp1
###################################################
#line 696 "Lecture1_2011.Rnw"
plot(wvs$secpay ~ wvs$gini, cex=1.5, xlab="Gini", ylab="Secpay")
lmod1 <- lm(secpay ~ gini, data=wvs)
abline(lmod1)
lines(lowess(wvs$gini, wvs$secpay, f=.8), lwd=2, col="red")
mod1 <- with(wvs, smooth.spline(x=gini, y=secpay, df=5))
pred.seq <- wvs$gini[order(wvs$gini)]
mod1.fit <- predict(mod1, newdata=data.frame(gini=pred.seq))
lines(mod1.fit$y ~ mod1.fit$x, lwd=2, col="green")
legend("topright", c("Linear","Lowess","Smoothing Spline"), 
    col=c("black","red","green"), lty=c(1,1,1), xjust=1, lwd=c(2,2,2), inset=.01)
 
 
###################################################
### chunk number 21: lmsmoothlowess1
###################################################
#line 708 "Lecture1_2011.Rnw"
#line 696 "Lecture1_2011.Rnw"
plot(wvs$secpay ~ wvs$gini, cex=1.5, xlab="Gini", ylab="Secpay")
lmod1 <- lm(secpay ~ gini, data=wvs)
abline(lmod1)
lines(lowess(wvs$gini, wvs$secpay, f=.8), lwd=2, col="red")
mod1 <- with(wvs, smooth.spline(x=gini, y=secpay, df=5))
pred.seq <- wvs$gini[order(wvs$gini)]
mod1.fit <- predict(mod1, newdata=data.frame(gini=pred.seq))
lines(mod1.fit$y ~ mod1.fit$x, lwd=2, col="green")
legend("topright", c("Linear","Lowess","Smoothing Spline"), 
    col=c("black","red","green"), lty=c(1,1,1), xjust=1, lwd=c(2,2,2), inset=.01)
#line 709 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 22: multnp
###################################################
#line 734 "Lecture1_2011.Rnw"
Model.lowess<-loess(secpay~gini+gdp, data=wvs, span=.8, degree=1)
gini2<-seq(min(wvs$gini), max(wvs$gini), len=15)
gdp2<-seq(min(wvs$gdp), max(wvs$gdp), len=15)
data<-expand.grid(gini=gini2, gdp=gdp2)
secpay.fitted<-matrix(predict(Model.lowess, data), 15, 15)
persp(gini2, gdp2, secpay.fitted, theta=35, ph=20, ticktype='detailed',
  xlab='Gini Coefficient', ylab='GDP', zlab='Attitudes', shade=.5)
 
 
###################################################
### chunk number 23: multnpreg1
###################################################
#line 743 "Lecture1_2011.Rnw"
#line 734 "Lecture1_2011.Rnw"
Model.lowess<-loess(secpay~gini+gdp, data=wvs, span=.8, degree=1)
gini2<-seq(min(wvs$gini), max(wvs$gini), len=15)
gdp2<-seq(min(wvs$gdp), max(wvs$gdp), len=15)
data<-expand.grid(gini=gini2, gdp=gdp2)
secpay.fitted<-matrix(predict(Model.lowess, data), 15, 15)
persp(gini2, gdp2, secpay.fitted, theta=35, ph=20, ticktype='detailed',
  xlab='Gini Coefficient', ylab='GDP', zlab='Attitudes', shade=.5)
#line 744 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 24: gdpf
###################################################
#line 801 "Lecture1_2011.Rnw"
require(mgcv, quietly=T)
mod1 <- gam(secpay ~ s(gini) + s(gdp), data=wvs)
gdp.dat <- with(wvs, data.frame(gdp=seq(from=min(gdp), to=max(gdp), length=25), gini=mean(gini, na.rm=T)))
pred.gdp <- predict(mod1, gdp.dat, se.fit=T)
plot.gdp <- cbind(with(gdp.dat, data.frame(gdp=gdp)), with(pred.gdp, data.frame(est=fit, lower=fit-2*se.fit, upper=fit+2*se.fit)))
 
with(plot.gdp, plot(gdp,est, type="l", ylim=c(min(lower), max(upper)), xlab="GDP", ylab="Predicted Secpay"))
with(plot.gdp, lines(gdp, lower, lty=2))
with(plot.gdp, lines(gdp, upper, lty=2))
 
 
###################################################
### chunk number 25: ginif
###################################################
#line 813 "Lecture1_2011.Rnw"
gini.dat <- with(wvs, data.frame(gini=seq(from=min(gini), to=max(gini), length=25), gdp=mean(gdp, na.rm=T)))
pred.gini <- predict(mod1, gini.dat, se.fit=T)
plot.gini <- cbind(with(gini.dat, data.frame(gini=gini)), with(pred.gini, data.frame(est=fit, lower=fit-2*se.fit, upper=fit+2*se.fit)))
 
with(plot.gini, plot(gini,est, type="l", ylim=c(min(lower), max(upper)), xlab="GINI", ylab="Predicted Secpay"))
with(plot.gini, lines(gini, lower, lty=2))
with(plot.gini, lines(gini, upper, lty=2))
 
 
###################################################
### chunk number 26: gdpfig
###################################################
#line 822 "Lecture1_2011.Rnw"
#line 801 "Lecture1_2011.Rnw"
require(mgcv, quietly=T)
mod1 <- gam(secpay ~ s(gini) + s(gdp), data=wvs)
gdp.dat <- with(wvs, data.frame(gdp=seq(from=min(gdp), to=max(gdp), length=25), gini=mean(gini, na.rm=T)))
pred.gdp <- predict(mod1, gdp.dat, se.fit=T)
plot.gdp <- cbind(with(gdp.dat, data.frame(gdp=gdp)), with(pred.gdp, data.frame(est=fit, lower=fit-2*se.fit, upper=fit+2*se.fit)))
 
with(plot.gdp, plot(gdp,est, type="l", ylim=c(min(lower), max(upper)), xlab="GDP", ylab="Predicted Secpay"))
with(plot.gdp, lines(gdp, lower, lty=2))
with(plot.gdp, lines(gdp, upper, lty=2))
#line 823 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 27: ginifig
###################################################
#line 825 "Lecture1_2011.Rnw"
#line 813 "Lecture1_2011.Rnw"
gini.dat <- with(wvs, data.frame(gini=seq(from=min(gini), to=max(gini), length=25), gdp=mean(gdp, na.rm=T)))
pred.gini <- predict(mod1, gini.dat, se.fit=T)
plot.gini <- cbind(with(gini.dat, data.frame(gini=gini)), with(pred.gini, data.frame(est=fit, lower=fit-2*se.fit, upper=fit+2*se.fit)))
 
with(plot.gini, plot(gini,est, type="l", ylim=c(min(lower), max(upper)), xlab="GINI", ylab="Predicted Secpay"))
with(plot.gini, lines(gini, lower, lty=2))
with(plot.gini, lines(gini, upper, lty=2))
#line 826 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 28: outmod
###################################################
#line 848 "Lecture1_2011.Rnw"
mod1 <- lm(secpay ~ gini + gdp, data=wvs)
 
 
###################################################
### chunk number 29: inflplot1
###################################################
#line 852 "Lecture1_2011.Rnw"
library(car)
influencePlot(mod1, identify="auto")
 
 
###################################################
### chunk number 30: davisinf
###################################################
#line 865 "Lecture1_2011.Rnw"
data(Davis)
davis.mod <-lm(repwt~weight, data=Davis) 
with(Davis, plot(repwt ~ weight, xlab="Weight", ylab="Reported Weight"))
abline(davis.mod)
 
 
###################################################
### chunk number 31: davisnoinf
###################################################
#line 871 "Lecture1_2011.Rnw"
davis.mod1 <-lm(repwt~weight, data=Davis, subset = -12) 
with(Davis[-12,], plot(repwt ~ weight, xlab="Weight", ylab="Reported Weight"))
abline(davis.mod1)
 
 
###################################################
### chunk number 32: davisinffig
###################################################
#line 878 "Lecture1_2011.Rnw"
#line 865 "Lecture1_2011.Rnw"
data(Davis)
davis.mod <-lm(repwt~weight, data=Davis) 
with(Davis, plot(repwt ~ weight, xlab="Weight", ylab="Reported Weight"))
abline(davis.mod)
#line 879 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 33: Davismod1
###################################################
#line 884 "Lecture1_2011.Rnw"
apsrtable(davis.mod, model.names="", Sweave=TRUE)
 
 
###################################################
### chunk number 34: davisnoinffig
###################################################
#line 891 "Lecture1_2011.Rnw"
#line 871 "Lecture1_2011.Rnw"
davis.mod1 <-lm(repwt~weight, data=Davis, subset = -12) 
with(Davis[-12,], plot(repwt ~ weight, xlab="Weight", ylab="Reported Weight"))
abline(davis.mod1)
#line 892 "Lecture1_2011.Rnw"
 
 
###################################################
### chunk number 35: Davismod1a
###################################################
#line 897 "Lecture1_2011.Rnw"
apsrtable(davis.mod1, model.names="", Sweave=TRUE)
 
 
###################################################
### chunk number 36: huber1
###################################################
#line 944 "Lecture1_2011.Rnw"
require(MASS, quietly=T)
huber.mod <- rlm(repwt~weight, method="MM", data=Davis)
plot(Davis$repwt ~ Davis$weight, xlab="Weight", ylab="Reported Weight")
abline(davis.mod, lwd=2)
abline(huber.mod, col="red", lty=2, lwd=2)
legend("topleft", 
c("OLS","Huber"), col=c("black","red"), lty=c(1,2), inset=.01)
 
 
###################################################
### chunk number 37: huber_w
###################################################
#line 970 "Lecture1_2011.Rnw"
w <- huber.mod$w
names(w) <- rownames(Davis)[
	-which(is.na(Davis$repwt))]
w <- w[order(w)]
pdf("huber_weights.pdf", 
	height=6, width=6)
plot(1:length(w), w, 
	xlab="index", 
	ylab="Weight")
text((1:length(w))[w<.3], 
	jitter(w[w < .3],1), 
	names(w)[w<.3], pos=4)
invisible(dev.off())
 
 
###################################################
### chunk number 38: florida1
###################################################
#line 995 "Lecture1_2011.Rnw"
data(Florida)
mod1 <- lm(BUCHANAN ~ GORE, data=Florida)
mod2 <- lm(BUCHANAN ~ GORE, data=Florida, subset=-which(rownames(Florida) =="PALM.BEACH"))
mod3 <- lm(BUCHANAN ~ GORE, data=Florida, subset=-which(rownames(Florida) %in% c("DADE", "BROWARD")))
mod4 <- lm(BUCHANAN ~ GORE, data=Florida, subset=-which(rownames(Florida) %in% c("PALM.BEACH","DADE", "BROWARD")))
gpb <- which(rownames(Florida) =="PALM.BEACH")
gd <- which(rownames(Florida) %in% c("DADE"))
gb <-  which(rownames(Florida) %in% c("BROWARD"))
plot(BUCHANAN ~ GORE, data=Florida)
with(Florida, text(GORE[gpb], BUCHANAN[gpb], rownames(Florida)[gpb], pos=1))
with(Florida, text(GORE[gb], BUCHANAN[gb], rownames(Florida)[gb], pos=2))
with(Florida, text(GORE[gd], BUCHANAN[gd], rownames(Florida)[gd], pos=1))
 
# Add the linear model results with different colors and line-types
abline(mod1)
abline(mod2, col="red", lty=2)
abline(mod3, col="blue", lty=3)
abline(mod4, col="green", lty=4)
 
# Add a legend in the upper-left corner
legend("topleft", c("OLS", "No Palm Beach", "No Dade or Broward", "No PB, Dade or Broward"), 
    col=c("black", "red", "blue", "green"), lty=c(1,2,3,4), inset=.01)
 
 
###################################################
### chunk number 39: rrji1
###################################################
#line 1049 "Lecture1_2011.Rnw"
mod1a <- rlm(BUCHANAN ~ GORE, data=Florida, method="MM")
names(mod1a$w) <- rownames(Florida)
plot(mod1a$w[mod1a$w < .5], xlim=c(0.5,8.5), ylim=c(0,.3), pch=16)
text((1:sum(mod1a$w < .5)), mod1a$w[mod1a$w<.5], names(mod1a$w)[mod1a$w < .5], pos=3)	
 
 
###################################################
### chunk number 40: bugo
###################################################
#line 1056 "Lecture1_2011.Rnw"
lowwt <- names(mod1a$w)[mod1a$w < .5]
pchs <- rep(1, nrow(Florida))
cols <- rep("black", nrow(Florida))
pchs[which(rownames(Florida) %in% lowwt)] <- 16
cols[which(rownames(Florida) %in% lowwt)] <- "red"
plot(BUCHANAN ~ GORE, data=Florida, pch = pchs, col=cols)
abline(lm(BUCHANAN ~ GORE, data=Florida), lty=1, col="black", lwd=1.5)
abline(lm(BUCHANAN ~ GORE, data=Florida, subset = !(rownames(Florida) %in% lowwt)), lty=2, col="red", lwd=1.5)
legend("topleft", c("All Obs", "Red Points Removed"), col=c("black", "red"), lty=c(1,2), inset=.01)
 
 
###################################################
### chunk number 41: meds
###################################################
#line 1113 "Lecture1_2011.Rnw"
library(boot)
sds <- aggregate(wvs$gdp, list(wvs$democrat), sd)
ns <- aggregate(wvs$gdp, list(wvs$democrat), length)
meds <- aggregate(wvs$gdp, list(wvs$democrat), median)
means <- aggregate(wvs$gdp, list(wvs$democrat), mean)
theor.ci <- diff(meds[,2]) + 2*sqrt(sum((1.25^2 * sds[,2]^2)/ns[,2]))*c(-1,1)
theor.ci2 <- diff(means[,2]) + 2*sqrt(sum((sds[,2]^2)/ns[,2]))*c(-1,1)
 
diff.medians <- function(data, i){
	tmp <- data[i,]
	m1 <- quantile(tmp[which(tmp$democrat == 0), "gdp"], .5)
	m2 <- quantile(tmp[which(tmp$democrat == 1), "gdp"], .5)
	diff <- m2-m1
	diff
}
diff.means <- function(data, i){
	tmp <- data[i,]
	m1 <- mean(tmp[which(tmp$democrat == 0), "gdp"])
	m2 <- mean(tmp[which(tmp$democrat == 1), "gdp"])
	diff <- m2-m1
	diff
}
 
myboot <- boot(wvs, diff.medians, R=2500)
myboot2 <- boot(wvs, diff.means, R=2500)
med.cis <- rbind(theor.ci, boot.ci(myboot, type="bca")$bca[4:5])
mean.cis <- rbind(theor.ci2, boot.ci(myboot2, type="bca")$bca[4:5])
plot.dat <- data.frame(
	points = c(diff(means[,2]), diff(means[,2]), 
		diff(meds[,2]), diff(meds[,2])),
	lower = c(mean.cis[,1], med.cis[,1]),
	upper = c(mean.cis[,2], med.cis[,2]),
	stat = factor(c(1,1,2,2), levels=1:2, labels=c("Mean", "Median")), 
	kind = factor(c(1,2,1,2), levels=1:2, labels=c("Theoretical", "Bootstrap")))
library(lattice)
yrg <- range(c(plot.dat[,c("lower", "upper")]))
 
trellis.par.set(strip.background=list(col="white"))
print(
xyplot(points ~ kind | stat, data=plot.dat, xlab = "", ylab = "",
	ylim = yrg + c(-1,1)*.1*diff(yrg), 
	panel = function(x,y,subscripts){
		panel.points(x,y,pch=16, col="black")
		panel.segments(x,plot.dat$lower[subscripts], x, plot.dat$upper[subscripts], lty=1, col="black")
	}))
 
 
###################################################
### chunk number 42: fmod
###################################################
#line 1206 "Lecture1_2011.Rnw"
library(flexmix)
dmod <- lm(prestige ~ type*income*education, data=Duncan)
fmod <- flexmix(prestige ~ type, data=Duncan, 
	model = FLXMRglmfix(nested=list(k=c(1,1), formula=c(~ .*education, ~.*income))))
 
b1 <- dmod$coef
b2a <- fmod@components$Comp.1[[1]]
b2a <- b2a@parameters$coef
b2b <- fmod@components$Comp.2[[1]]
b2b <- b2b@parameters$coef
 
mat <- round(cbind(b1, b2a[match(names(b1), names(b2a))], b2b[match(names(b1), names(b2b))]), 2)
colnames(mat) <- c("OLS", "FMM: C1", "FMM: C2")
 
 
###################################################
### chunk number 43: fmod_tab
###################################################
#line 1223 "Lecture1_2011.Rnw"
library(xtable)
xtable(mat)
 
 
###################################################
### chunk number 44: fmod_preds
###################################################
#line 1236 "Lecture1_2011.Rnw"
pred.dmod <- fitted(dmod)
preds.fmod <- fitted(fmod)
pred.fmod1 <- preds.fmod[cbind(1:nrow(preds.fmod), fmod@cluster)] 
pred.fmod2 <- diag(preds.fmod %*% t(fmod@posterior$scaled))
cormat <- cor(cbind(pred.dmod, pred.fmod1, pred.fmod2), Duncan$prestige)
colnames(cormat) <- "Correlation"
rownames(cormat) <- c("OLS", "FMM1", "FMM2")
 
plot.dat2 <- data.frame(
	pred = c(pred.dmod, pred.fmod1, pred.fmod2), 
	orig = rep(Duncan$prestige, 3), 
	gp = factor(rep(1:3, each=length(pred.dmod)), 
		levels=1:3, labels=c("OLS", "FMM1", "FMM2")))
pdf("fmpred.pdf", height=6, width=9)
trellis.par.set(strip.background = list(col="white"))
xyplot(pred ~ orig | gp, data=plot.dat2, 
	xlab = "Original", ylab = "Predicted", layout=c(3,1))
invisible(dev.off())
 
 
###################################################
### chunk number 45: post_dens
###################################################
#line 1264 "Lecture1_2011.Rnw"
plot(density(fmod@posterior$scaled), xlab="Pr(Member of Group 1)", ylab="Density", main="")