Should Sklearn's `chi2` function use the Chi-Square statistic or p-value for feature selection? TLDR: From my understanding, as Chi-Square is a statistic dependent on the degrees of freedom, it cannot be used to compare reliably across different features, but it seems to be used in Sklearn for example. In my view, p-values should be used instead.
Hi,
I am juggling my mind over the following:


*

*The p-value for a given Chi-Square statistic depends on the degrees of freedom (namely (columns - 1) * (rows - 1)).

*This means that, for a given value of the Chi-Square statistic (say 2), it might or might not translate into statistical significance (say at alpha =  0.5), depending on the degrees of freedom.

*When applied to different features against the target, this could result in different degrees of freedom.


*

*To better illustrate, let's say the target is binary, and feature 1 is categorical with 3 distinct values, while feature 2 is also categorical with 10 potential values. The degrees of freedom for feature 1 would be 2, and for feature 2 would be 9.


*When I look into Sklearn's chi2 code and documentation, I conclude that the Chi-Square statistic is in fact used to sort the features for subsequent selection.

*This could result in situations in which we are leaving in some features that are not significantly related to the target (as per the statistical test) while dropping others that are significantly related to the target.

*Using p-values instead would solve this issue, as these are already adjusted for the degrees of freedom.


Could you help me by identifying where my reasoning is failing (if it is indeed)?
Thanks!
Refs: Sklearn docs, Related question
 A: I was actually thinking it will not be an issue with sklearn because you need to do onehot encoding for categorical variables anyway, but as you rightly noted in your comments, it is not a contingency table. For binary cases, it tests whether there's an excess of feature values in  you are trying to predict.
We can use an example dataset:
import sklearn
from sklearn.datasets import load_iris
X, y = load_iris(return_X_y=True)
# use only 2 label
y = (y==1).astype(int)

And run the relevant part from the source code, this part does a dot product between binarized label and features, basically to get the sum per feature, for rows with labels 1:
Y = sklearn.preprocessing.LabelBinarizer ().fit_transform(y)
observed = sklearn.utils.extmath.safe_sparse_dot(Y.T, X)          

So observed looks like this:
array([[296.8, 138.5, 213. ,  66.3]])

In this example 1/3 of the labels are 1, we expect the sum of each feature to be 1/3 of the total, and this is used as expected: 
feature_count = X.sum(axis=0).reshape(1, -1)

array([[876.5, 458.6, 563.7, 179.9]])

class_prob = Y.mean(axis=0).reshape(1, -1)
expected = np.dot(class_prob.T, feature_count)

array([[292.16666667, 152.86666667, 187.9       ,  59.96666667]])

And then the chisquare is done using a function defined in sklearn, to test observed and predicted.
When you have a k-class prediction (k>2), the observed and predicted will be a kxn matrix, and the chi-square will need to be done on k-1 degree of freedom.
