select.cp.aic.bic.keep.fun_function(x,y,keepit)
{
  ## this program works for a search of models with main effects and no interactions
  ## For the sex discrim data, the command to search all models within the SSOM *and*
  ## keep the 4 main effects in is:
  ## select.cp.aic.bic.keep.fun(y=log(bsal),x=cbind(senior,age,educ,exper,senior2,age2,educ2,exper2,senior.i.age,senior.i.educ,senior.i.exper,age.i.educ,age.i.exper,educ.i.exper),keep=1:4)
  ## the output is ordered according to BIC as in display 12.12 of the book.
  n_length(y)
  fullmod <- lm(y~x)
  s2full <- (summary(fullmod)$sigma)^2
  sstot_crossprod(y-mean(y))
  nullmod <- lm(y~1)
  null.p <- 1
  null.cp <- 1
  null.s2 <- summary(nullmod)$sigma^2
  null.bic <- n*log(null.s2)+null.p*log(n)
  null.aic <- n*log(null.s2)+2
  null.r2 <- summary(nullmod)$r.squared
  null.r2adj <- 1-null.s2/(sstot/(n-1))
  r2info_leaps(x,y,method="r2",keep=keepit)
  mstot <- sstot/n
  p_r2info$size
  r2 <- r2info$r2/100
  r2adj <- (r2-(1-r2)*(p-1)/(n-p))
  ssres_(1-r2info$r2/100)*sstot
  sigma2hat_ssres/(n-p)
  Cp <- n-p*sigma2hat/s2full+(2*p-n)
  AIC_n*log(sigma2hat)+2*p
  BIC_n*log(sigma2hat)+p*log(n)
  ord <- order(BIC)
  output_data.frame(model=r2info$label[ord],p=r2info$size[ord],Cp=round(Cp,3)[ord],AIC=round(AIC,3)[ord],BIC=round(BIC,3)[ord],R2=round(r2info$r2/100,4)[ord],R2adj=round(r2adj,4)[ord])
  print("Subset selection criteria results:")
  nullmod <- data.frame("NULL",round(null.p,0),round(null.cp,3),round(null.aic,3),round(null.bic,3),round(null.r2,4),round(null.r2adj,4))
  names(nullmod) <- names(output)
return(rbind(output,nullmod))
}