8
$\begingroup$

If I have a dependent variable and $N$ predictor variables and wanted my stats software to examine all the possible models, there would be $2^N$ possible resulting equations.

I am curious to find out what the limitations are with regard to $N$ for major/popular statistic software since as $N$ gets large there is a combinatorial explosion.

I've poked around the various web pages for packages but not been able to find this information. I would suspect a value of 10 - 20 for $N$?

If anyone knows (and has links) I would be grateful for this information.

Aside from R, Minitab, I can think of these packages SAS, SPPS, Stata, Matlab, Excel(?), any other packages I should consider?

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ @levon9 - I was able to generate all possible subsets using 50 variables in SAS. I do not believe there is any hard limitation other than memory and CPU speed. -Ralph Winters $\endgroup$
    – Ralph Winters
    Mar 1, 2011 at 14:56
  • $\begingroup$ For what size of dataset? $\endgroup$
    – M. Tibbits
    Mar 1, 2011 at 16:17
  • $\begingroup$ Thanks very useful information. Just curious, did this this take long? $\endgroup$
    – Levon
    Mar 1, 2011 at 21:34
  • 2
    $\begingroup$ @levon9 This question generated a lot of sound answers and comments so that I've +1. But, please, forget about Excel for doing serious job in model selection... $\endgroup$
    – chl
    Mar 1, 2011 at 22:21
  • $\begingroup$ @chl .. is that because Excel is slow, or just incapable (ie would give inaccurate results?). $\endgroup$
    – Levon
    Mar 1, 2011 at 23:56

4 Answers 4

Reset to default
12
$\begingroup$

I suspect 30--60 is about the best you'll get. The standard approach is the leaps-and-bounds algorithm which doesn't require fitting every possible model. In $R$, the leaps package is one implementation.

The documentation for the regsubsets function in the leaps package states that it will handle up to 50 variables without complaining. It can be "forced" to do more than 50 by setting the appropriate boolean flag.

You might do a bit better with some parallelization technique, but the number of total models you can consider will (almost undoubtedly) only scale linearly with the number of CPU cores available to you. So, if 50 variables is the upper limit for a single core, and you have 1000 cores at your disposal, you could bump that to about 60 variables.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ leaps is great, I like the plots from it, +1. In real applications some averaging techniques work faster (and better) than pretested estimators even found from all regression models. So I would suggest to go for Bayesian model averaging (BMA package) or the algorithm I like the most - Weighted average least squares (WALS[1]) developed by J.R.Magnus et al. The Matlab code is easily transformable to R code. The good thing for WALS is $N$ computational difficulty instead of $2^N$. [1]: tilburguniversity.edu/research/institutes-and-research-groups/… $\endgroup$
    – Dmitrij Celov
    Mar 1, 2011 at 11:14
  • 1
    $\begingroup$ @Dmitrij, thanks for your comments. I've tried to remain fairly agnostic in my response regarding the utility of all-subsets regression. It seems to me that there is nearly always a better solution, but I felt that might seem like too trite of a response to the OP's question. $\endgroup$
    – cardinal
    Mar 1, 2011 at 12:34
  • $\begingroup$ @Dmitrij, BMA over main-effects models would still have the same computational complexity as all-subsets regression. No? The main advantage in BMA seems to me to be in trying to figure out which covariates are likely to be influencing the response. BMA does this by essentially averaging the log likelihoods over $2^{n-1}$ submodels. $\endgroup$
    – cardinal
    Mar 1, 2011 at 12:38
  • $\begingroup$ Thanks for the pointer to the leaps R package! I didn't know about it and it might come in handy in the future. If I could get some information on specific limitations for N for other popular packages it would be very helpful. $\endgroup$
    – Levon
    Mar 1, 2011 at 12:44
  • 1
    $\begingroup$ @levon9, I doubt it will vary much at all by package. The algorithm that leaps uses has been state of the art for at least 20 years or so. Even if you found an implementation that was twice as fast, that would mean you got to increment the number of variables you can handle by one. For every doubling of the speed, you get one more variable. Hardware, not algorithmic, limitations are your bottleneck in this case. $\endgroup$
    – cardinal
    Mar 1, 2011 at 13:00
10
$\begingroup$

Just a caveat, but feature selection is a risky business, and the more features you have, the more degrees of freedom you have with which to optimise the feature selection criterion, and hence the greater the risk of over-fitting the feature selection criterion and in doing so obtain a model with poor generalisation ability. It is possible that with an efficient algorithm and careful coding you can perform all subsets selection with a large number of features, that doesn't mean that it is a good idea to do it, especially if you have relatively few observations. If you do use all subsets selection, it is vital to properly cross-validate the whole model fitting procedure (so that all-subset selection is performed independently in each fold of the cross-validation). In practice, ridge regression with no feature selection often out-performs linear regression with feature selection (that advice is given in Millar's monograph on feature selection).

Improve this answer
$\endgroup$
6
  • 2
    $\begingroup$ @Dikran, (+1) good comments. I was trying to avoid going there since it didn't directly address the OP's question. But, I agree. All-subsets is rarely the way to go. And, if you do go that way, you need to understand all the implications. $\endgroup$
    – cardinal
    Mar 1, 2011 at 13:19
  • $\begingroup$ @Dirkan thank you for the comments, I am a real stats newbie. I realize the danger of overfitting the model when too many variables are in play, so I am just looking at various automated ways (ie without much benefit of insight) such as the stepwise approach (which might get caught in a local maxima) and the exhaustive all subsets model - and the computational limits it faces (and the external limitations imposed by packages) $\endgroup$
    – Levon
    Mar 1, 2011 at 13:48
  • 3
    $\begingroup$ @levon9, you can get over-fitting that is just as serious when you choose the features, so feature selection does not guard against over-fitting. Consider a logistic regression model used to predict the outcome of flipping a fair coin. The potential inputs are the outcome of flipping a large number of other fair coins. Some of these inputs will be positively correlated with the target, so the best all-subsets model will select inputs (even though they are useless) and you will get a model that appears to have skill, but in reality is no better than guessing. $\endgroup$
    – Dikran Marsupial
    Mar 1, 2011 at 14:39
  • $\begingroup$ @Dikran (+1) the same as @cardinal, I first wrote a similar text, but then decided it is not what @levon9 asked, because he simply was curious about the complexity :) $\endgroup$
    – Dmitrij Celov
    Mar 1, 2011 at 16:06
  • $\begingroup$ @Dikran +1 because I like such advice. $\endgroup$
    – chl
    Mar 1, 2011 at 22:22
5
$\begingroup$

I was able to generate all possible subsets using 50 variables in SAS. I do not believe there is any hard limitation other than memory and CPU speed.

Edit

I generated the 2 best models for N=1 to 50 variables for 5000 observations.

@levon9 - No, this ran in under 10 seconds. I generated 50 random variables from (0,1)

-Ralph Winters

Improve this answer
$\endgroup$
3
  • $\begingroup$ For what size of dataset? $\endgroup$
    – M. Tibbits
    Mar 1, 2011 at 16:17
  • $\begingroup$ Thanks very useful information. Just curious, did this this take long? $\endgroup$
    – Levon
    Mar 1, 2011 at 21:34
  • 1
    $\begingroup$ I have undeleted this post (and merged another of your comments in an edit) because the OP found it useful and others might too. Thank you for your contribution; please keep it up! (If you really think it should be deleted, go ahead and do so; I won't contravene your wishes again.) $\endgroup$
    – whuber ♦
    Mar 2, 2011 at 5:53
3
$\begingroup$

As $N$ gets big, your ability to use maths becomes absolutely crucial. "inefficient" mathematics will cost you at the PC. The upper limit depends on what equation you are solving. Avoiding matrix inverse or determinant calculations is a big advantage.

One way to help with increasing the limit is to use theorems for decomposing a large matrix inverse from into smaller matrix inverses. This can often means the difference between feasible and not feasible. But this involves some hard work, and often quite complicated mathematical manipulations! But it is usually worth the time. Do the maths or do the time!

Bayesian methods might be able to give an alternative way to get your result - might be quicker, which means your "upper limit" will increase (if only because it gives you two alternative ways of calculating the same answer - the smaller of two, will always be smaller than one of them!).

If you can calculate a regression coefficient without inverting a matrix, then you will probably save a lot of time. This may be particularly useful in the Bayesian case, because "inside" a normal marginalisation integral, the $X^{T}X$ matrix does not need to be inverted, you just calculate a sum of squares. Further, the determinant matrix will form part of the normalising constant. So "in theory" you could use sampling techniques to numerically evaluate the integral (even though it has an analytic expression) which will be eons faster than trying to evaluate the "combinatorical explosion" of matrix inverses and determinants. (it will still be a "combinatorical explosion" of numerical integrations, but this may be quicker).

This suggestion above is a bit of a "thought bubble" of mine. I want to actually test it out, see if it's any good. I think it would be (5,000 simulations + calculate exp(sum of squares) + calculate weighted average beta should be faster than matrix inversion for a big enough matrix.)

The cost is approximate rather than exact estimates. There is nothing to stop you from using the same set of pseudo random numbers to numerically evaluate the integral, which will again, save you a great deal of time.

There is also nothing stopping you from using a combination of either technique. Use exact when the matrices are small, use simulation when they are big. This is because in this part of the analysis. It is just different numerical techniques - just pick the technique which is quickest!

Of course this is all just a bit of "hand wavy" arguments, I don't exactly know the best software packages to use - and worse, trying to figure out which algorithms they actually use.

Improve this answer
$\endgroup$
14
  • 2
    $\begingroup$ @probabilityislogic, while your response is interesting, maybe it could be refocused to better address the OP's question. Also, no software for computing least-squares solutions performs matrix inversion, much less a determinant. Ever. Unless it's inverting a $1 \times 1$ matrix. $\endgroup$
    – cardinal
    Mar 1, 2011 at 12:28
  • 1
    $\begingroup$ @probabilityislogic, handling the $2^n$ cases efficiently quickly far outstrips the $O(n^3)$ issues of an efficient least-squares solution. That's where the leaps-and-bounds algorithm comes in. $\endgroup$
    – cardinal
    Mar 1, 2011 at 12:30
  • $\begingroup$ Thanks for the post. "Do the maths or do the time!" :-) .. I'm actually not even trying to figure out the underlying algorithms used by the packages (thought that is interesting to know), at this point really looking for specific information regarding the limitations of N by the major packages. $\endgroup$
    – Levon
    Mar 1, 2011 at 12:47
  • $\begingroup$ @cardinal The updating and downdating algorithms also exist for various matrix decomposition procedures, I suspect that is what was meant by "matrix inverse" etc. $\endgroup$
    – Dikran Marsupial
    Mar 1, 2011 at 13:04
  • $\begingroup$ @Dikran, several efficient and numerically stable approaches to least-squares exist, including methods for augmenting or reducing a design matrix by one column at a time. Sometimes it's good to understand what is happening under the surface, even if on most days you don't need to care. $\endgroup$
    – cardinal
    Mar 1, 2011 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.