require(leaps)
swiss <- read.delim("swiss.txt",row.names=1)

## Model selection criteria
fit <- lm(Fertility~Agriculture+Education,data=swiss)
fit.full <- lm(Fertility~.,data=swiss)
n <- nrow(swiss)
sig2 <- summary(fit.full)$sigma^2
extractAIC(fit)
extractAIC(fit,k=log(n))
extractAIC(fit,scale=sig2)
sum((fit$resid/(1-hatvalues(fit)))^2)

## Best subsets
regs <- regsubsets(Fertility~.,data=swiss)
plot(regs)

regs <- regsubsets(Fertility~.,data=swiss,nbest=10)
plot(regs)
plot(regs,scale="Cp")

summary(lm(scale(Fertility)~0+scale(Examination),data=swiss))

## Stepwise
regs <- regsubsets(Fertility~.,data=swiss)
regs.f <- regsubsets(Fertility~.,data=swiss,method="forward")
par(mfrow=c(1,2))
plot(regs)
plot(regs.f)

x1 <- rnorm(100)
x2 <- rnorm(100)
x <- (x1+x2)/2+rnorm(100,sd=0.5)
y <- x1+x2+rnorm(100)
df <- data.frame(y,x,x1,x2)
regs.bs <- regsubsets(y~.,data=df)
regs.f <- regsubsets(y~.,data=df,method="forward")
par(mfrow=c(1,2))
plot(regs.bs)
mtext("Best subsets",line=1)
plot(regs.f)
mtext("Forward stepwise",line=1)

regs.bs <- regsubsets(Fertility~.*.,data=swiss)
regs.f <- regsubsets(Fertility~.*.,data=swiss,method="forward")
par(mfrow=c(1,2))
plot(regs.bs)
plot(regs.f)