bes.leaps<-function(x,y,ret=F){
# gets backward 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
pvmin<-rep(0,m)
regsubsets(x,y,method="backward")->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)
# sequential min of p-values down from full model
for (i in m:1){pvmin[i]<-min(pv.orig[m:i])}
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("Backward elimination sequence in reverse order",fill=T)
cat("Remove variables from bottom of list to top",fill=T)
cat("",fill=T)
print(out<-data.frame(p=c(vm,m+1),s=c(vm-1,m),rem.var=out.x$vorder-1,
         pval=round(c(NA,pv.orig),4),pvmin=round(c(NA,pvmin),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("rem.var=next variable to remove, current model includes",fill=T)
cat("        this variable and all variables above it",fill=T)
cat("pval=p-value to remove variable",fill=T)
cat("pvmin=min p-value from below, BS(alpha)=move down until pvmin > 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)
}