Best subset selection uses the branch and bound algorithm to speed up

-
Answering my own question in the last post: best-subset selection is not stepwise. The speed was achieved due to the specific algorithm applied. Same results from R stepwise and SAS best subset were obtained due to coincidence. After I remove some variables the outputs are no longer the same. The output from R best subset is still not exactly the same as SAS best subset. The differences are very small, only on several predictors with trivial effects.

I discussed this question on SAS community, got very good comments that remind me obviously R did not search for all the combinations either. Some speed-up algorithm must have been applied.

Then I sent an email to SAS technical help (follow this page), and got a reply in a day. The algorithm is called "Furnival-Wilson leaps-and-bounds" (1974), which seems quite basic in computer science after I figured out what it is. This article is fairly easy to read: Branch and bound in statistical data analysis D.J.Hand.
It scans for combinations given each allowed numbers of variables to use (1, 2, 3, ..., p) and find an optimal solution for each given level (thus Cp is bounded). This is why the program always reports the optimal model based on the number of predictors chosen: e.g., the best 1-predictor model, the best 2-predictor model, the best 10-predictor model... etc.


I hereby attach the reply from SAS:

--------------------------

 Yang,
With selection-RSQUARE, ADJRSQ, and CP and n=number of regressors >=11, by default REG will only DISPLAY the best n subset models for each number of regressors. The best n one-variable models, best n two-variable models, etc. These can be computed (using the Furnival and Wilson algorithm) without examining every possible model of every possible size and so this is typically much faster than if all models of each size need to be displayed.

By default, when you run PROC REG with the SELECTION=CP and the STOP=10 option, and you have 20 regressors in the model, PROC REG will display at most 20 models for each of the 1-variable models, 2-variable models, 3-variable models, ...through 10-variable models.  In other words, the maximum number of models displayed will be equal to the number of predictor variables in the MODEL statement (if the number of predictors listed in the MODEL statement is greater than 11).   

To obtain more models than the number displayed by default, you will need to add the BEST= option to the MODEL statement.  For example, if you have 20 predictors in your MODEL statement, but you want to see up to 35 models in each of the possible subsets, then your PROC REG step would need to look something like:

---------------------
proc reg data=test;
  model y = x1-x20 / selection=rsquare stop=10 best=35;
run;
quit;
----------------------

I hope the above information is helpful.

Kathleen Kiernan
Senior Principal Technical Support Statistician

Comments

Post a Comment

Popular posts from this blog

Power-law distribution (Pareto)& Zipf's Law: connection and how to fit the distribution of global city population

-
Image
- This post counts for 3 days. - In short,  We use the histogram to outline the distribution, use rank-frequency plot (log-log plot) to identify the power-law distribution, and use maximum likelihood estimation to obtain the parameters. -  I will show 1 st , how the power-law function CDF and PDF can be derived from Zipf CDF; and 2 nd , this least square estimation by excel is wrong and how to get the right one. Historical background on the power-law distribution. - Power-law distribution is as common as normal distribution in the real world. In the 19 th century, Italian economist Pareto realized that 20% of the population owns 80% the wealth, which is the famous 80/20 rule. He later named the power-law distribution as Pareto distribution, which is what we studied in the probability course. - In the 20 th century, another main contributor, Zipf made a similar discovery in linguistics about the frequency of words used. - Besides social wealth, many things follow a hig
Read more

How to Draw Heatmap with Colorful Dendrogram

-
Image
Heatmap with Dendrogram: Color Branches and Labels This data visualization example include: * Hierarchical clustering, dendrogram, and heat map based on normalized odds ratios * The dendrogram was built separately to give color to dendrogram’s branches/labels based on clusters using  dendextend * Heatmap is made by  heatmap.2  from  gplots  using the built dendrogram * The rows are sorted by means from highest to lowest, it can be done in either the dendrogram or the  heatmap.2 * Use color palettes from  colorspace Data YRBSS (Youth Risk Behavior Surveillance System) survey data from CDC. This example plots the odds ratios of 63 behavioral questions on one question regarding disordered eating. The odds ratios are median centered by column and log2 transformed Version 1: Color both the branches and labels library (curl) # read file from google drive library (gplots) # heatmap.2 library (dendextend) # make and color dendrogram library (colorspace) #
Read more

The Weighted Standard Deviation

-
Image
I have revised the figures in the Conflicts in Colombia paper to add weighted s.e. Together with confidence intervals. Turns out to be not that complicated. Weighted s.d.: Add { geom_errorbar(aes(ymin=..., ymax=...), width = 0.2 } in ggplot2 to add the confidence interval into plots: The paper discussed how conflicts in Colombia makes cities develop more sluggish and denser. I will work more on it in the later posts. The complete R code for the figure # weighted means for global 200 sample ggplot(data=m4_mean_g_g200, aes(x=year, y=wm, group=Group, color=Group, shape=Group))+ geom_point(size= 4 )+ geom_line(size= 1 )+ geom_errorbar(aes(ymin=CI_low, ymax=CI_high), width= .2 )+ ggtitle( "Urban Density Change Over Time in 200 Cities" ) + labs(x= "Year" , y= "Average City Density" , color = "" , shape= "" ) + scale_y_continuous(limits = c( 30 , 150 ))
Read more