Find correlation between continuous predictive features and a continuous target feature

Question

In my data, I have about 10K predictive features (genes), and one target feature (age). I want to predict the ages according to the genes. The rows in the data are the patients. To do so I plan to use Regression Random Forest.

I don't want to use this many predictive features, so I want to do some feature selection first.

There is no linear correlation between the predictive features and the target feature (at least I didn't find any relationship for the features that I checked).

For binary features, when I want to predict gender for example, I can just use Wilcoxon test to find the most significant features that separate the two classes. Here, I can't use such a test.

How can I find the most important features for age prediction? Can I just run a random forest algorithm and then just check the most important feature? would that work without creating noise?

Here is a subset of the training set, including the age feature:

dput(train_scaled[1:20,c(1,2,3,4,5,dim(train_scaled)[2])])
structure(list(A1BG = c(1.81619824260442, 1.9986779809134, 1.91171736562985, 
1.87425799530611, 1.95720931978885, 1.68534041055052, 1.89237252718096, 
1.67216783026329, 1.94555622783709, 2.05581255682001, 1.89803035420513, 
1.7563466972377, 1.85448031100116, 1.90469081497093, 1.82958626152702, 
1.80639351405546, 1.94904037078298, 1.88121448353727, 1.90265126862802, 
1.89685997844076), ADA = c(1.25649644487907, 1.22729343759834, 
1.27344838192825, 1.25955928072103, 1.26370991138808, 1.20355435132166, 
1.23956642505305, 1.25589256673664, 1.15163992141014, 1.20146398841983, 
1.09375020345131, 1.19284479092003, 1.18821270400345, 1.15707902340534, 
1.29848225592125, 1.2563306911831, 1.29923301554395, 1.22251152311355, 
1.22795303612616, 1.48761789143517), CDH2 = c(0.53688090267567, 
0.493919738045297, 0.560208693940622, 0.588029409349587, 0.559643640625794, 
0.570599153392745, 0.562110779919758, 0.54921119370662, 0.507086211915313, 
0.496614809627379, 0.581539495325737, 0.597444486905757, 0.560166965896316, 
0.579972731871132, 0.583039148505923, 0.581924465154048, 0.566420208700464, 
0.576395012253254, 0.575907185558433, 0.453946904680819), AKT3 = c(0.917211707537678, 
0.892003590486357, 0.969818024729793, 0.978292068213014, 0.913032018184228, 
0.948312269441081, 0.947709935054217, 0.83611701240751, 0.912172816373717, 
0.98719118237761, 1.02711099335984, 0.922819275258826, 0.933697725060485, 
0.996194969362905, 0.971300509819334, 0.851048415219854, 0.9156277536571, 
0.982369058418409, 0.832254764434006, 0.905941809264712), MED6 = c(2.02291559929045, 
2.08170269351807, 2.04355176601994, 2.05526765226102, 1.93189920401206, 
2.03859461894252, 1.97348257053102, 1.9229558498545, 1.95605272086482, 
2.06298256427372, 2.11184798077237, 1.99810309844712, 2.01005618200693, 
2.06589538426559, 2.1372244020894, 1.967894127866, 2.01416144921981, 
2.02184221220218, 1.90343367987094, 1.9634446015096), age = c(69, 
30, 64, 65, 61, 70, 48, 73, 40, 58, 62, 53, 75, 68, 52, 67, 50, 
70, 78, 53)), row.names = c("Patient12", "P10", "P11", 
"PX123", "PX77", "P1", "ER45", "ER30", "Patient8", 
"Patient9", "Patient10", "EA6327611", "EA6329802", "EA6839018", "EA6389069", 
"EA6359107", "EA6359120", "EA6391391", "EA6399146", "EA6391153"), class = "data.frame")

Answer 1

This is a problematic approach.

If you just want to know about general relationships between two variables, not restricting to linear (Pearson correlation) or even monotonic (Spearman correlation) relationships, you could use a value like mutual information. The R function JMI::JMI is one way to calculate mutual information between two variables. (My experience with this function is that it is slow.)

That is probably the answer to the posted question: use mutual information to calculate/estimate the overall relationships between pairs of variables so you are not restricted to the particular relationships detected by, for instance, Pearson or Spearman correlation. Then, a flexible model like a random forest will figure out the nonlinear, nonmonotonic relationships in the regression.

However, a flexible model like a random forest also looks at interactions between variables and their nonlinear transformations. By only considering one feature at a time in the mutual information calculations, you miss all of those. In fact, any kind of feature-by-feature screening is going to miss these interactions. In terms of information theory, two variables can be independent (zero mutual information) yet be conditionally dependent, conditional on the value of a third variable. Considering just two variables at a time will always miss that conditional dependence, yet a random forest model will be able to discover such relationships and use them to make accurate predictions (subject to the usual concerns about overfitting).

Answer 2

The bootstrap can help in these types of problems and is a valuable procedure in exposing the true difficulty of the task by taking a lot of uncertainties into account if you do one-at-a-time feature selection or a joint model (e.g. elastic net or better ridge regression). The idea is to get bootstrap confidence intervals for the importance ranking of all candidate features simultaneously. Importance can be measured in the one-at-a-time case by Wilcoxon statistics (or its equivalent Somers' $D_{xy}$ rank correlation or concordance probability $c$), Spearman $\rho$, $\chi^2$, ordinary correlation, or anything you want.

Simulating data like yours is another approach, where you'll see that for non-large samples there is not even a correlation between the true and estimated gene-specific effects across genes and show a scatterplot.

These methods are described and exemplified here and here.

We are mainly fooling ourselves about the ability to answer gene-level questions in the $N << p$ case.

Aside: not sure that age in general is a meaningful target feature, depending on what even causes age to be assessed.