HOW to add the leaps package to R:  (UNIX systems)

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) # in isds or under windows if package has been added
library(leaps, lib.loc="~/Rpackages")  
help(leaps)
help(regsubsets)

source("bayesreg.all.R")
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.zs <- bayesreg(ozoneX,ozone[,4],prior="ZSnullLAP", nbest=50,nvmax=35,regmeth="seqrep")

#out.oz65.bic <- bayesreg(ozoneX,ozone[,4],prior="BIC", nbest=50,nvmax=20,regmeth="seqrep")

names(out.oz65.zs)

 [1] "post.order" "postprob"   "label"      "r2"         "bic"       
 [6] "size"       "which"      "which.int"  "margs"      "probne0"   
[11] "postmean"   "postsd"     "ols"        "se"   


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