HOW to add the leaps package to R:

1)  Download leaps package from http://cran.r-project.org/src/contrib/PACKAGES.html

2) tar -xzf leaps_2.6.tar.gz  

- or on some systems -

2) gunzip leaps_2.6.tar.gz

   tar -xf leaps_2.6.tar   


3) R CMD INSTALL leaps -l Rpackages

4) run R as usual

##########################################
R 

library(leaps, lib.loc="~/Rpackages")
help(leaps)

help(regsubsets)


ozone <- read.table("ozonefit.txt", header=F)
ozone.pred <- read.table("ozonepred.txt", header=F)
nfit <- dim(ozone)[1]
npred <- dim(ozone.pred)[1]
X <- poly(rbind(as.matrix(ozone[,c(1,2,3,5,6,7,8,10,11,13)]),
                as.matrix(ozone.pred[,c(1,2,3,5,6,7,8,10,11,13)])),
          degree=2)

xnames <- paste("x", as.character(1:10), sep="")
intnamesmat <- outer(xnames, xnames, paste, sep=".")
intnames <- intnamesmat[1,2]
for (i in 2:9) {
  intnames <- c(intnames, intnamesmat[1:i, i+1])}

ozoneX <- X[1:nfit,c(cumsum(1:10),cumsum(2:11), (1:65)[-c(cumsum(1:10),cumsum(2:11))])]
colnames(ozoneX) <- c(xnames, paste(xnames,".2", sep=""), intnames)
X.pred <- X[-(1:nfit),c(cumsum(1:10),cumsum(2:11), (1:65)[-c(cumsum(1:10),cumsum(2:11))])]
colnames(X.pred) <- c(xnames, paste(xnames,".2", sep=""), intnames)

ozone.lm <- lm(ozone[,4] ~ ozoneX, x=T)
out.oz65.s <- bayesreg(ozoneX,ozone[,4],prior="ZSnullLAP",nbest=50,nvmax=35,regmeth="seqrep")
names(out.oz65.f)

postscript("ozone-comp.ps")
par(mfrow=c(2,2))
plot(out.oz65.s$r2, out.oz65.s$margs, xlab="R2", ylab="Log Marginal Zellner-Siouw")
plot(out.oz65.s$bic, out.oz65.s$margs, xlab="BIC", ylab="Log Marginal Zellner-Siouw")
plot(out.oz65.s$size, out.oz65.s$margs, xlab="model dimension",
ylab="Log Marginal")
plot(out.oz65.s$size, out.oz65.s$bic, xlab="model dimension", ylab="BIC")
dev.off()