###################################################
### chunk number 1: plot1
###################################################
#line 446 "Lecture9_2010.Rnw"
library(boot)
weakliem <- read.table("Weakliem2.txt")
plot(secpay ~ gini, data=weakliem)
outs <- which(rownames(weakliem) %in% 
	c("Slovakia","CzechRepublic"))
with(weakliem, text(gini[outs], secpay[outs], 
	rownames(weakliem)[outs], pos=4))
 
 
###################################################
### chunk number 2: modsum1
###################################################
#line 460 "Lecture9_2010.Rnw"
mod <- lm(secpay ~ gini, data=weakliem)
library(apsrtable)
apsrtable(mod, model.names = "", Sweave=T, digits=4)
 
 
###################################################
### chunk number 3: infl
###################################################
#line 471 "Lecture9_2010.Rnw"
library(car)
influencePlot(mod, identify="auto")
 
 
###################################################
### chunk number 4: modsum2
###################################################
#line 483 "Lecture9_2010.Rnw"
mod2 <- lm(secpay ~ gini, data=weakliem, subset=-c(7,26))
apsrtable(mod2, model.names="", Sweave=T, digits=4)
 
 
###################################################
### chunk number 5: rob1
###################################################
#line 493 "Lecture9_2010.Rnw"
library(MASS)
mod3 <- rlm(secpay ~ gini, data=weakliem) 
summary(mod3)
 
 
###################################################
### chunk number 6: plotres
###################################################
#line 513 "Lecture9_2010.Rnw"
plot(residuals(mod3))
text(outs, residuals(mod3)[outs], 
	rownames(weakliem)[outs], pos=2)
 
 
###################################################
### chunk number 7: loww1
###################################################
#line 533 "Lecture9_2010.Rnw"
loww <- with(mod3, which(w < 1))
with(mod3, plot(w))
text(loww, mod3[["w"]][loww], 
	rownames(weakliem)[loww], pos=2)
 
 
###################################################
### chunk number 8: rob2
###################################################
#line 543 "Lecture9_2010.Rnw"
mod4 <- rlm(secpay ~ gini, data=weakliem, 
	psi= psi.bisquare)
summary(mod4)
 
 
###################################################
### chunk number 9: plotres4
###################################################
#line 562 "Lecture9_2010.Rnw"
with(mod4, plot(residuals))
text(outs, mod4[["residuals"]][outs], 
	rownames(weakliem)[outs], pos=2)
 
 
###################################################
### chunk number 10: wcompare
###################################################
#line 580 "Lecture9_2010.Rnw"
plot(mod4[["w"]], mod3[["w"]], xlim=c(0,1),
    ylim=c(0,1), xlab="bisquare weights",
    ylab="Huber Weights")
abline(0,1)
text(mod4[["w"]][loww], mod3[["w"]][loww], 
	rownames(weakliem)[loww], pos=2)
 
 
###################################################
### chunk number 11: robmm
###################################################
#line 714 "Lecture9_2010.Rnw"
mod5 <- rlm(secpay ~ gini, data=weakliem, 
	method="MM")
summary(mod5)
 
 
###################################################
### chunk number 12: compmod
###################################################
#line 725 "Lecture9_2010.Rnw"
mat <- matrix(NA, ncol=3, nrow=5)
mat[1,1:2] <- c(coef(mod)[2], sqrt(diag(vcov(mod))[2]))
mat[2,1:2] <- c(coef(mod2)[2], sqrt(diag(vcov(mod2))[2]))
mat[3,1:2] <- c(coef(mod3)[2], sqrt(diag(vcov(mod3))[2]))
mat[4,1:2] <- c(coef(mod4)[2], sqrt(diag(vcov(mod4))[2]))
mat[5,1:2] <- c(coef(mod5)[2], sqrt(diag(vcov(mod5))[2]))
mat[,3] <- mat[,1]/mat[,2]
colnames(mat) <- c("Estimate", "SE", "t")
rownames(mat) <- c("OLS (with outliers)", "OLS (no outliers)", 
	"M (Huber)", "M (bisquare)", "MM (Huber)")
library(xtable)
xtable(mat, digits=4)
 
 
###################################################
### chunk number 13: compfig
###################################################
#line 752 "Lecture9_2010.Rnw"
plot(secpay ~ gini, data=weakliem)
abline(mod, lty=1)
abline(mod2, lty=2)
abline(mod3, lty=3)
abline(mod4, lty=4)
abline(mod5, lty=5)
legend("topright", rownames(mat), lty=1:5, inset=.01)
 
 
###################################################
### chunk number 14: rr1
###################################################
#line 782 "Lecture9_2010.Rnw"
plot(residuals(mod) ~ residuals(mod3))
abline(lm(residuals(mod) ~ residuals(mod3)), lty=3)
abline(0,1, col="gray60")
legend("bottomright", c("LM Line", "45-degree"), lty=c(3,1), 
	col=c("black", "gray60"), inset=.01)
 
 
###################################################
### chunk number 15: rr2
###################################################
#line 789 "Lecture9_2010.Rnw"
plot(residuals(mod) ~ residuals(mod4))
abline(lm(residuals(mod) ~ residuals(mod4)), lty=3)
abline(0,1, col="gray60")
legend("bottomright", c("LM Line", "45-degree"), lty=c(3,1), 
	col=c("black", "gray60"), inset=.01)
 
 
###################################################
### chunk number 16: rr3
###################################################
#line 796 "Lecture9_2010.Rnw"
plot(residuals(mod) ~ residuals(mod5))
abline(lm(residuals(mod) ~ residuals(mod5)), lty=3)
abline(0,1, col="gray60")
legend("bottomright", c("LM Line", "45-degree"), lty=c(3,1), 
	col=c("black", "gray60"), inset=.01)
 
 
###################################################
### chunk number 17: assocmod1
###################################################
#line 847 "Lecture9_2010.Rnw"
library(foreign)
assoc <- read.spss("assoc.sav", 
	use.value.labels=T, to.data.frame=T)[1:1000, ]
assoc.mod <- glm(ASSOC ~ SEX + AGE + SES, 
	data=assoc, family=poisson(link=log))
apsrtable(assoc.mod, model.names="", Sweave=T, digits=4)
 
 
###################################################
### chunk number 18: glmcook
###################################################
#line 864 "Lecture9_2010.Rnw"
plot(cookd(assoc.mod))
bigd <- which(cookd(assoc.mod) > .05)
text(bigd, cookd(assoc.mod)[bigd], 
	as.character(bigd), pos=1)
 
 
###################################################
### chunk number 19: glmdfb
###################################################
#line 870 "Lecture9_2010.Rnw"
dfb <- dfbetas(assoc.mod)
plot(dfb[,6], ylab = "DFBetas for Unskilled")
text(bigd, dfb[bigd, 6], 
	as.character(bigd), pos=1)
 
 
###################################################
### chunk number 20: bigdobs
###################################################
#line 877 "Lecture9_2010.Rnw"
assoc[bigd, ]
 
 
###################################################
### chunk number 21: glmrob
###################################################
#line 898 "Lecture9_2010.Rnw"
library(robustbase)
glm2 <- glmrob(ASSOC ~ SEX + AGE + SES, data=assoc, family=poisson(link=log))
summary(glm2)
 
 
###################################################
### chunk number 22: exptab
###################################################
#line 916 "Lecture9_2010.Rnw"
mat <- cbind(exp(coef(assoc.mod)), exp(coef(glm2)))
colnames(mat) <- c("GLM", "Robust GLM")
xtable(mat)
 
 
###################################################
### chunk number 23: rreg_b
###################################################
#line 936 "Lecture9_2010.Rnw"
mod2 <- rlm(secpay ~ gini, data=weakliem, method="MM")
summary(mod2)
 
 
###################################################
### chunk number 24: booth1
###################################################
#line 964 "Lecture9_2010.Rnw"
boot.huber <- function(data, inds, maxit=20){
	assign(".inds", inds, envir=.GlobalEnv)
  	mod <- rlm(secpay ~ gini, data=data[.inds,],
	 	maxit=maxit)
	remove(".inds", envir=.GlobalEnv)
  	coefficients(mod)
}
 
 
###################################################
### chunk number 25: boot1
###################################################
#line 986 "Lecture9_2010.Rnw"
boot.ranx <- boot(data=weakliem, statistic=boot.huber,
	R=1500, maxit=100)
boot.ranx
 
 
###################################################
### chunk number 26: boot.ci1
###################################################
#line 1002 "Lecture9_2010.Rnw"
boot.ci(boot.ranx, type=c("normal","perc", "bca"), index=2)
 
 
###################################################
### chunk number 27: boot2
###################################################
#line 1012 "Lecture9_2010.Rnw"
fit <- fitted(mod2)
e <- residuals(mod2)
boot.huber.fixed <- function(data, inds, maxit=20){
	assign(".inds", inds, envir=.GlobalEnv)
  	y.star <<- fit + e[.inds]
	mod <- update(mod2, y.star ~ .)
  	remove(".inds", envir=.GlobalEnv)
  	coefficients(mod)
}
boot.fixx <- boot(weakliem, boot.huber.fixed, 
	R=1500, maxit=100)
boot.fixx
 
 
###################################################
### chunk number 28: bootdist
###################################################
#line 1033 "Lecture9_2010.Rnw"
d1 <- density(boot.ranx$t[,2])
d2 <- density(boot.fixx$t[,2])
xl <- range(c(d1$x, d2$x))
yl <- range(c(d1$y, d2$y))
pdf("boot_rlm_dist.pdf", height=6, width=6)
plot(d1, type="l", xlim=xl, ylim=yl, 
	xlab = "Bootstrap Coefficient")
lines(d2, lty=2)
legend("topleft", c("Random X", "Fixed X"), 
	lty=c(1,2), inset=.01)
invisible(dev.off())
 
 
###################################################
### chunk number 29: cis
###################################################
#line 1055 "Lecture9_2010.Rnw"
ranx.ci <- boot.ci(boot.ranx, index=2, 
	type=c("perc", "bca"))
fixx.ci <- boot.ci(boot.fixx, index=2, 
	type=c("perc", "bca"))
analytical.ci <- coef(mod2)[2] + qt(.975, 
	nrow(weakliem)-mod2$rank)*
	sqrt(diag(vcov(mod2)))[2] * c(-1,1)
cis <- rbind(ranx.ci$bca[4:5], 
	  fixx.ci$bca[4:5], 
	  analytical.ci)
rownames(cis) <- c("randomX", "fixedX", "analytical")
colnames(cis) <- c("lower", "upper")
round(cis, 3)