General principles

• The larger n is, then the larger a value of p that can be used before variance increase overwhelms bias reduction.

• However, that bias-variance tradeoff transition point for p depends on the application. There is no magic formula for best p as a function of n.

• Nevertheless, a commonly-used rule of thumb is to limit p to at most n^0.5. Some theoretical analyses even suggest the more conservative rule p < log(n). Again, these are oversimplifications, but you may find them useful as guides to intuition.

• It should be kept in mind that categorical features can greatly add to the value of p. See below.

Feature selection methods should produce an ordered sequence of candidate models

A nethod should produce some kind of ordered list, for instance, best single predictor; best pair of predictors; best triplet etc. We can evaluate each of the p feature sets via cross-validation–i.e. noting the testAcc component in the output of qeML functions–and then choosing the one that peforms best. (A method may simply rank features according to importance, but we still would assess the performance of the first feature, the first two features together, the first three features together and so on.)

In the all possible subsets (APS) approach, by contrast, one looks at all the 2^p possible feature sets, a prohibitively large number computationally.

Moreover, remember, we are working with randomness. We choose holdout sets randomly, and usually the dataset itself is considered to be a random sample from some (real or conceptual) population. Due to this randomness, it’s possible that some combination of features will appear to perform well, just by accident. The more feature sets we look at, the greater the chance of this occurring.

Note that some of the methods discussed below do their own internal cross-validation, and ultimately give us the “best” feature set according to that criterion. Nevertheless, they still provide an ordered feature set.

This last point is important because we may wish to use the selected feature set sequence as input to a different ML method. There is no theoretical reason to believe the best feature set for one ML method is also best for another, but as noted, feature selection is an ad hoc business.

Methods based on p-values

A classic approach for linear and generalized linear models is to first fit a full model, then retain only those features that are “statistically significant”. In the nyctaxi taxi data shown above, we would retain, for instance PULocationID33 but not PULocationID35.

This would violate our Desired Property A above. Use of p-values in geneeral is problematic (see the “NoPVals” vignette), and this is no exception. Just because β_i is nonzero does not mean that Feature i has strong predictive power; what if β_i is nonzero but tiny?

One could generate an ordered sequence of feature sets in the above manner by first setting the threshold for a p-value very low, then progressively higher This is similar to stepwise regression: In the forward version, one starts with no features, but adds them one at a time. At each step, one tests for the β coefficient being 0, and then chooses the feature that is most “significant.” When no new feature is found to be “significant,” the process stops, and we use whatever features managed to make it in to the model so far. The backward version starts with a full model, and removes features one at a time. Though this one may seem to be an improvement, in that it does take into account that a feature may be useful individually but not if similar features are already in the equation, it still suffers from the fundamental problems of p-values.

The LASSO

Here we fit a linear model, subject to the constraint that no estimated β_i is larger (in absolute value) than λ, a hyperparameter. Starting at 0, we increase λ to ever-larger values, having the effect that one estimated β_i becomes nonzero at a time. This gives us an ordered sequence of candidate features, as desired.

This is similar to a forward stepwise method, but NOT based on p-values; instead, the decision of which feature to use next is based on cross-validated mean squared prediction error.

Here is the LASSO approach on the nyctaxi data:

> lassout <- qeLASSO(nyctaxi,'tripTime') 
> lassout$lambda.whichmin
[1] 51
> lassout$whenEntered
  trip_distance PULocationID.132 PULocationID.186      DayOfWeek.1
               2               29               31               31
DOLocationID.132 DOLocationID.124      DayOfWeek.5  DOLocationID.33
              33               34               34               35
DOLocationID.189 DOLocationID.225 PULocationID.213  PULocationID.76
              35               35               37               38
PULocationID.263   DOLocationID.1 DOLocationID.177      DayOfWeek.2
              38               38               38               38
PULocationID.124  DOLocationID.35  DOLocationID.36  DOLocationID.89
              39               39               39               39
DOLocationID.231 DOLocationID.263 PULocationID.239 DOLocationID.155
              39               39               40               40
 PULocationID.40  PULocationID.43 PULocationID.146 PULocationID.161
              41               41               41               41
PULocationID.230  DOLocationID.22  DOLocationID.37 DOLocationID.161
              41               41               41               41
DOLocationID.162 DOLocationID.251 PULocationID.163 PULocationID.211
              41               41               42               42
PULocationID.257  DOLocationID.49  DOLocationID.81 DOLocationID.179
              42               42               42               42
...

For each value of λ run by the code, one or more features joined the model. The cross-valued mean squared prediction error reached its minimum at the 51st λ. Along the way, trip_distance was the first feature added, at the 2nd λ, followed by PULocationID132 (29th), PULocationID186 (31st), DayOfWeek.1 (31st) and so o.

Methods based on measures of feature importance

Consider the Random Forests (RF) ML method. It has various hyperparameters, such as maximum number of levels in a tree, and we wish to find a “good” combination of those settings. As a byproduct, though, the fact that at each level RF is entering a new feature, many RF implemntations compute some kind of variable importance ranking–thus giving us our desired ordered sequence.

We can then run RF on each feature set in the sequence, or as mentioned, try the sequence on some other ML method. We would then use whichever feature set yielded the best cross-validated performance.

Here is an example with the nyctaxi data:

> z <- qeRFranger(nyctaxi,'tripTime')
holdout set has  1000 rows
Loading required namespace: ranger
Warning message:
In eval(tmp, parent.frame()) :
  9 rows removed from test set, due to new factor levels
> z$variable.importance
trip_distance  PULocationID  DOLocationID     DayOfWeek
   2447411690     212569502     201650058      85733194

So, our feature set sequence might be:

trip_distance
trip_distance and PULocationID
trip_distance and PULocationID and DOLocationID
trip_distance and PULocationID and DOLocationID
trip_distance and PULocationID and DOLocationID and DayOfWeek

Note that the algorithm chose the categorical variables as a whole here, not breaking down into their dummy variable components.

Another example is Prnciple Components Analysis (PCA). Here our p features are replaced by p linear combinations (the “principle components,” PCs) of the original features. The PCs are uncorrelated.

The PCs form our new features, and are ordered by variance, largest to smallest. Since a variable with small variance is essentially constant and thus has little predictive value, we take only the PCs with larger variance.

For the first m PCs, say, we compare the sum of their variances (which is also the variance of their sum) to the total variance of the response variable Y (i.e. the variable we are trying to predict). By varying the proportion, we produce an ordered sequence of feature sets, as desired.

The qeML function qePCA first finds the PCs, then applies whatever ML method is specified by the user. The function qeUMAP does the same for UMAP, a kind of nonlinear analog of PCA.

One more measure of variable importance is Shapley values. This concept describes an attempt to use game theory to apportion credit for a win among several players in a team; the “players” here are the features. While it is a clever idea, it has been the subject of much criticism in terms of practical value.

Feature Ordering by Conditional Independence (FOCI)

This is my personal “go to” method for feature selection. Its theoretical foundations are complex, but it boils down to measuring the reduction in variance in predicting Y from X, versus predicting Y simply from its mean.

FOCI is implemented in a CRAN package of the same name, and qeML’s qeFOCI function provides a convenient interface. Example:

> w <- qeFOCI(nyctaxi,'tripTime')
> w$selectedVar
   index            names
1:     1    trip_distance
2:   202  DOLocationID.75
3:   348      DayOfWeek.1
4:    78 PULocationID.151

Trip distance, dropoff location 72, Mondays and pickup location 151 were chosen, in that order, after which the algorithm decided that adding any further features was counterproductive. Again, we could simply take that 4-feature set for whatever ML method we wish, or we could use the above to set up an ordered sequence of feature sets,

trip_distance
trip_distance and DOLocationID.75
trip_distance and DOLocationID.75 and DayOfWeek.1
trip_distance and DOLocationID.75 and DayOfWeek.1 and PULocationID.151

running our desired ML method on each feature set, then choosing the one with best cross-validation performance.

Direct dimension reduction for categorical data

Rather than doing feature selection per se, we might transform the data to summary variables, or embeddings. We could group the pickup/dropoff locations by ZIP Code, for instance.

Or, we could consider only the more frequently used locations. The qeML function factorToTopLevels would help here. We will use it to reduce pickup location to just a few places that appear a lot in our data; all others will be lumped together as ‘other’.

We invoke the function as follows:

> factorToTopLevels(nyctaxi$PULocationID,lowCountThresh=0)

alt text

Ah, yes, lots of locations have small counts, less than 50 and probably many in the single-digits range. So, PULocationID is a prime candidate for simplication. Let’s go for a cutoff of 75.

> newPUlocs <- factorToTopLevels(nyctaxi$PULocationID,lowCountThresh=75)

So, now we have a new, simpler version of the pickup location data. Let’s take a look at old versus new.

> head(newPUlocs,50)
 [1] 236   238   100   161   236   161   238   107   170   other 170   137  
[13] 239   143   164   164   151   239   161   100   142   107   238   43   
[25] 237   170   238   238   236   161   236   other 230   138   other 79   
[37] 186   132   other 100   234   142   151   263   236   142   142   144  
[49] 236   230  
43 Levels: 100 107 113 114 13 132 137 138 140 141 142 143 144 148 151 ... other
> head(nyctaxi$PULocationID,50)
 [1] 236 238 100 161 236 161 238 107 170 211 170 137 239 143 164 164 151 239 161
[20] 100 142 107 238 43  237 170 238 238 236 161 236 41  230 138 261 79  186 132
[39] 152 100 234 142 151 263 236 142 142 144 236 230
224 Levels: 1 3 4 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 ... 265

So for example pickup location in row 10 of our data was location 211, but it was one of the “rare” locations, so it was recoded to “other”.

How much dimension reduction did we accomplish?

> length(levels(newPUlocs))
[1] 43

That’s much less than the original 265. Of course, we could reduce even further by using a larger threshold.

So, let’s try out the new data:

> newDOlocs <- factorToTopLevels(nyctaxi$DOLocationID,lowCountThresh=75)
> nyctaxi$PULocationID <- newPUlocs
> nyctaxi$DOLocationID <- newDOlocs
> z <- qeLin(nyctaxi,'tripTime')

Dimension reduction at least got us a stable solution. Does it predict well?

> z$testAcc
[1] 240.5642
> z$baseAcc
[1] 430.4696

Mean Absolute Prediction Error is much less than what one would have by simply guessing all cases to be the overall mean trip time.

By varying lowCountThresh, we could generate an ordered sequence of feature sets, as before.