bic.sim<-function(x,y){
# chooses from FAS using minimum BIC
require(leaps)
ok<-complete.cases(x,y)
x<-x[ok,]                            # get rid of na's
y<-y[ok]                             # since regsubsets can't handle na's
m<-ncol(x)
n<-nrow(x)
as.matrix(x)->x                      # in case x is a data frame
colnames(x)<- as.character(1:m)
regsubsets(x,y,method="forward")->out.x
# gic gives model size including intercept of min BIC model
# we use 2:(m+2) because sigma is estimated
gic<-function(rss,m,n,z){which.min(z*(2:(m+2))+n*log(rss/n)+n+2*n*log(sqrt(2*pi)))}
if(m > n) out.x$rss<-out.x$rss[1:61]
if(m > n) m<-60
bic<-gic(rss=out.x$rss,m=m,n=n,z=log(n))   # model size including intercept
x<-x[,colnames(x)[(out.x$vorder-1)[2:bic]]]
if(bic>1) x.ind<-(out.x$vorder-1)[2:bic]  else x.ind<-0
if (bic==1) {mod = lm(y~1)} else {mod = lm(y~x)}
return(list(mod=mod,size=bic-1,x.ind=x.ind))
}