Is there a way to optimize regression according to a specific criterion? I am curious, is there a way to optimize a regression according to a specific statistic? 
Let's say I am interested in a model with the best possible AIC statistic (or MSE or whatever measurement I am interested in) - could I somehow direct the regression to give me the top X models that would do this? (Of course I would not ignore the other measures, but would it be possible to ask for this?) What software supports this, or would you write your own code (in say R)?
Also in general, when multiple regression results are displayed (let's say an all-possible or best-subset regression is done), is there a ranking and if so, according to what measure/criteria?
I'm not saying that this is the best way to evaluate models, but perhaps it would be a way to explore candidate models? (This is really not my main question though.)
Thanks .. learning a lot reading this site.
 A: Statistical appropriateness aside, R provides some nice functions that allow for this type of analysis. You can take a look at the leaps() function within the leaps package. The leaps() function returns the top n models and some fit statistics for a given set of parameters. stepAIC() within MASS is another handy function for this type of analysis.
There's a decent tutorial on the statmethods site describing these and some other techniques. 
A: Try wle.cp from package wle:
http://cran.r-project.org/web/packages/wle/index.html
It's based on Mallow's Cp:
Mallow's Cp
Here's the example given in the reference manual.
library(wle)
x.data <- c(runif(60,20,80),runif(5,73,78))
e.data <- rnorm(65,0,0.6)
y.data <- 8*log(x.data+1)+e.data
y.data[61:65] <- y.data[61:65]-4
z.data <- c(rep(0,60),rep(1,5))
plot(x.data,y.data,xlab="X",ylab="Y")

xx.data <- cbind(x.data,x.data^2,x.data^3,log(x.data+1))
colnames(xx.data) <- c("X","X^2","X^3","log(X+1)")
result <- wle.cp(y.data~xx.data,boot=10,group=10,num.sol=2)
summary(result)
plot(result,num.max=15)

result <- wle.cp(y.data~xx.data+z.data,boot=10,group=10,num.sol=2)
summary(result)
plot(result,num.max=15)

The output from the last summary(result) statement is:
Call:
wle.cp(formula = y.data ~ xx.data + z.data, boot = 10, group = 10, 
    num.sol = 2)


Weighted Mallows Cp:
      (Intercept) xx.dataX xx.dataX^2 xx.dataX^3 xx.datalog(X+1) z.data   wcp
 [1,]           0        0          0          0               1      1 1.570
 [2,]           1        0          0          0               1      1 2.372
 [3,]           0        1          0          0               1      1 2.510
 [4,]           0        0          1          0               1      1 2.564
 [5,]           0        0          0          1               1      1 2.570
 [6,]           1        1          1          1               0      1 4.088
 [7,]           0        1          1          1               1      1 4.289
 [8,]           1        0          1          1               1      1 4.530
 [9,]           1        1          0          1               1      1 4.710
[10,]           1        1          1          0               1      1 4.888 
[11,]           1        1          1          1               1      1 6.000

Printed the first  11  best models 

The top model (row [1,] above), which uses log(X+1) from xx.data (see the above 1 under xx.datalog(X+1)) and z.data (see the above 1 under z.data) has the lowest Mallow's Cp value (wcp = 1.57).
The final plot(result,num.max=15) statement provides the following graph where any green model under the black line follows Mallow's criteria.  The blue model in the lower left area is the "best Mallow's Cp model" (see the above list).

A: At the simple end of the spectrum: Minitab will do "best subsets regression", which will find the best one predictor model, the best two predictor model, the best 3 predictor model, the best 4 predictor model, etc. The criterion is r-squared.
