fas.leaps<-function(x,y,ret=F){
# gets forward sequence using leaps package
# x is matrix or data frame of independent variables
# y is vector of responses
# ret=T returns printed stuff for plotting
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)                           # number of ind. var.
n<-nrow(x)                           # sample size
vm<-1:m
as.matrix(x)->x                      # in case x is a data frame
pvmax<-rep(0,m)
regsubsets(x,y,method="forward")->out.x
pv.orig<-1-pf((out.x$rss[vm]-out.x$rss[vm+1])*(n-vm-1)/out.x$rss[vm+1],1,n-vm-1)
for (i in 1:m){pvmax[i]<-max(pv.orig[1:i])}  # sequential max of pvalues
Cp<-out.x$rss*(n-m-1)/out.x$rss[m+1]-n+2*(c(vm,m+1))
Cp.adj<-Cp-2*(m-c(0,vm))/(n-m-3)          # Gilmour's adjusted Cp
BIC<-n*log(2*pi)+n+n*log(out.x$rss/n)+log(n)*(c(vm,m+1)+1)   # true BIC
cat("Forward addition sequence",fill=T)
cat("",fill=T)
print(out<-data.frame(p=c(vm,m+1),s=c(vm-1,m),add.var=out.x$vorder-1,
         pval=round(c(NA,pv.orig),4),pvmax=round(c(NA,pvmax),4),
         Cp=round(Cp,3),Cp.adj=round(Cp.adj,3),
         BIC=round(BIC,3),MSE=round(out.x$rss/(n-c(vm,m+1)),3)))
cat("",fill=T)
cat("p=number of effects in model, p=1 is just the intercept",fill=T)
cat("s=number of independent variables in model",fill=T)
cat("add.var=next variable to add, current model includes",fill=T)
cat("        this variable and all variables above it",fill=T)
cat("pval=p-value to add variable",fill=T)
cat("pvmax=max p-value from above, FS(alpha)=move down until pvmax > alpha",fill=T)
cat("Cp=usual Cp, note p=s+1 where s=number of x variables in model",fill=T)
cat("Cp.adj=adjusted Cp, choose smallest model with Cp.adj < p=s+1",fill=T)
cat("BIC=official BIC, number of parameters=s+2, intercept & sigma^2",fill=T)
if(ret) return(out)
}