4

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.

whuber
  • 281,159
  • 54
  • 637
  • 1,101
Levon
  • 433
  • 7
  • 16
  • It is not clear what is your question. If you have set of regression models ranking them according to some criteria is trivial. The problem is producing the set of regressions. To do that you need some kind of data. No method can create regressions out of thin air. – mpiktas Mar 29 '11 at 13:59
  • 2
    @mpiktas - I don't think your comment is on-point. levon9 never said he wants to create regressions without data. He's looking for macros or code that constitute shortcuts to doing the work of ranking according to some criterion. This may be "trivial" in the sense of straightforward, but as I see it he's looking for a way to shorten the process. – rolando2 Mar 29 '11 at 14:22
  • @rolando2, ok my comment about thin air is not appropriate. Still I think the question is too vague. – mpiktas Mar 29 '11 at 14:31
  • @mpiktas - Fair enough. – rolando2 Mar 29 '11 at 18:15
  • @rolando2 - thanks for your comment to mpiktas, you got it exactly right. I'll be the first to admit I am a newbie when it comes to statistics, but I didn't think you could do regression without data so I didn't specifically mention it. If one can, then I guess mpiktas has a point and I should have mentioned it. – Levon Mar 29 '11 at 23:35

3 Answers3

5

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.

Chase
  • 3,055
  • 2
  • 19
  • 28
2

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).

enter image description here

bill_080
  • 3,458
  • 1
  • 20
  • 21
  • Thank you, that looks very much along the lines of what I was trying to do. I'll have to dig into this some more. – Levon Mar 29 '11 at 23:38
1

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.

zbicyclist
  • 3,363
  • 1
  • 29
  • 34
  • Thanks. Just curious, and only if you know of the top of your head, do you know if it would be possible to specify a different statistic as criteria? I.e., rank according to Adjusted $R^2$ or AIC? – Levon Mar 30 '11 at 21:49