###################################################
### chunk number 1: vif1
###################################################
#line 166 "ModelSelection_2011.Rnw"
options(useFancyQuotes=F)
library(car)
data(Prestige)
mod <- lm(prestige ~ income + education, data=Prestige)
summary(mod)


###################################################
### chunk number 2: vif2
###################################################
#line 182 "ModelSelection_2011.Rnw"
vif(mod)
sqrt(vif(mod))


###################################################
### chunk number 3: gvif1
###################################################
#line 229 "ModelSelection_2011.Rnw"
mod <- lm(prestige ~ income + education + type, data=Prestige)
summary(mod)


###################################################
### chunk number 4: gvif2
###################################################
#line 241 "ModelSelection_2011.Rnw"
vif(mod)


###################################################
### chunk number 5: nntest
###################################################
#line 680 "ModelSelection_2011.Rnw"
library(car)
data(Prestige)
mod1 <- lm(prestige ~ income + women, data=na.omit(Prestige), y=T)
mod2 <- lm(prestige ~ education + type + women, data=na.omit(Prestige), y=T)
smod1 <- summary(mod1)
smod2 <- summary(mod2)
library(lmtest)
jtest(mod1, mod2, Prestige)
encomptest(mod1, mod2, Prestige)


###################################################
### chunk number 6: nntest2
###################################################
#line 699 "ModelSelection_2011.Rnw"
coxtest(mod1, mod2, Prestige)
vuong <- (logLik(mod2) - logLik(mod1))-((mod2$rank - mod1$rank)/2)*log(nrow(na.omit(Prestige)))
vuong
# distribution-free test
il1 <- log(dnorm(mod1$y, mod1$fitted, smod1$sigma))
il2 <- log(dnorm(mod2$y, mod2$fitted, smod2$sigma))
cuts <- qbinom(c(.05,.95), size=length(mod1$y), prob=.5)
d <- sum(il1 > il2)
pbinom(d, length(il1), .5)
AIC(mod1)
BIC(mod1)


###################################################
### chunk number 7: step1
###################################################
#line 751 "ModelSelection_2011.Rnw"
data(Ericksen)
mod <- lm(undercount ~ minority + crime + poverty + language +
    highschool + housing, data=Ericksen)
mod.step.aic <- step(mod)
mod.step.bic <- step(mod, k=log(length(Ericksen$undercount)))


###################################################
### chunk number 8: modsel
###################################################
#line 812 "ModelSelection_2011.Rnw"
library(leaps)
X <- Ericksen[,c("minority", "crime", "poverty", "language", "highschool", "housing")]
y <- Ericksen$undercount
rmods <- regsubsets(x=X, y=y, method="exhaustive", all.best=TRUE, nbest=10)


###################################################
### chunk number 9: modsel2
###################################################
#line 823 "ModelSelection_2011.Rnw"
sqrt(vif(mod))


###################################################
### chunk number 10: subsetfig
###################################################
#line 853 "ModelSelection_2011.Rnw"
library(car)
pdf("subsetfig.pdf", height=6, width=6)
subsets(rmods, statistic="cp", legend=F)
dev.off()


###################################################
### chunk number 11: runmods
###################################################
#line 950 "ModelSelection_2011.Rnw"
c1 <- t(combn(4,1))
c2 <- t(combn(4,2))
c3 <- t(combn(4,3))
c4 <- t(combn(4,4))
l <- list()
k <- 1
for(i in 1:nrow(c1)){
	l[[k]] <- c1[i,]
k <- k+1
}
for(i in 1:nrow(c2)){
	l[[k]] <- c2[i,]
k <- k+1
}
for(i in 1:nrow(c3)){
	l[[k]] <- c3[i,]
k <- k+1
}
for(i in 1:nrow(c4)){
	l[[k]] <- c4[i,]
k <- k+1
}
P <- Prestige[,c("income", "education", "women", "type", "prestige")]
dats <- lapply(l, function(x)P[,c(x,5)])

mods <- lapply(dats, function(x)lm(
	as.formula(paste("prestige ~ ", ifelse(length(names(x)) > 2, paste(names(x)[1:(ncol(x)-1)], collapse=" + "), colnames(x)[1]), sep="")), data=x))
AICs <- sapply(mods, AIC)
delta.aic <- AICs - min(AICs)
w <- exp(-delta.aic/2)/sum(exp(-delta.aic/2))


###################################################
### chunk number 12: mmsv
###################################################
#line 989 "ModelSelection_2011.Rnw"
coefs <- vars <- matrix(0, ncol=6, nrow=length(mods))
colnames(coefs) <- names(coef(mods[[15]]))
for(i in 1:length(mods)){
	coefs[i,match(names(coef(mods[[i]])),colnames(coefs))] <- coef(mods[[i]])
	vars[i,match(names(coef(mods[[i]])),colnames(coefs))] <- diag(vcov(mods[[i]]))
}

theta.bar.hat <- w %*% coefs
part.2 <- sweep(coefs, 2, theta.bar.hat)^2
var.theta.bar.hat <- (w %*% sqrt(vars + part.2) )^2
## t-statistics
theta.bar.hat/var.theta.bar.hat


###################################################
### chunk number 13: impvars
###################################################
#line 1025 "ModelSelection_2011.Rnw"
impvars <- (w %*% (coefs != 0))[,2:5]
names(impvars) <- c("Income", "Education", "Women", "Type")
impvars <- impvars[order(impvars)]
pdf("impdot.pdf", height=6, width=6)
dotchart(impvars, pch=16)
invisible(dev.off())