chi2 should support categorical data other than binary or document count

[altanova]

Describe the bug

When looking for correlation between features (for feature selection), I found that sklearn implementation of Chi2 test of independence produce significantly different results from scipy.stats implementation.

My sample data contains 300 records, with 6 anonymized categorical features and the label. My focus is on the feature A. This data is available in this folder in github . The file sample300.csv has the file, while the file chi2_showcase.ipynb has the code demonstrating the mismatch.

For the feature A, sklearn's SelectKBest() returned the lowest ranking, suggesting there is no correlation between A and the target. But scipy.stats.chi2_contingency() returned very different result, suggesting the correlation is very high.

Because of mismatch between the two, I went a long way performing a number of different tests described in detail in this article The results suggest that the scipy implementation is correct, while sklearn implementation is incorrect.

Steps/Code to Reproduce

Please see the two links given above, where I provided the full source code and the results.
The piece of code is quite standard:

fs = SelectKBest(score_func=skfs.chi2, k = 'all')
X, y = df[cat_feature_cols], df[label]
selector = fs.fit(X, y)
kbest = pd.DataFrame({'feature': X.columns, 'score': fs.scores_})
kbest.sort_values(by = 'score', ascending = False).reset_index()

Expected Results

I would expect that sklearn.feature_selection.SelectKBest(score_func=skfs.chi2) returns same, or at least similar results (p-value and chi2 statistics) as scipy.stats.chi2_contingency() . In the particular case of feature A from my set, these expected results are:

chi2 = 127.497517
p-value = 1.445816e-29

Actual Results

For the feature A in my set, sklearn.feature_selection.chi2() (encapsulated inside SelectKBest(score_func=skfs.chi2)) returned lowest rank of all features, suggesting no correlation. Feature A has score 1.412797, while other features score between 1647 and 24.

In contrast, scipy.stats.chi2_contingency() gave highest rank to feature A, suggesting high correlation. The other tests described in the article suggest that the latter is correct.

Versions

System:
    python: 3.7.7 (default, May  6 2020, 11:45:54) [MSC v.1916 64 bit (AMD64)]
executable: C:\Users\pplaszczak\AppData\Local\Continuum\anaconda3\python.exe
   machine: Windows-10-10.0.18362-SP0

Python dependencies:
          pip: 21.0.1
   setuptools: 52.0.0.post20210125
      sklearn: 0.24.1
        numpy: 1.19.2
        scipy: 1.6.2
       Cython: 0.29.22
       pandas: 1.2.3
   matplotlib: 3.3.4
       joblib: 1.0.1
threadpoolctl: 2.1.0

Built with OpenMP: True

[glemaitre]

So our implementation of sklearn.feature_selection.chi2 and the scipy.stats.chisquare provides the same results:

# %%
import pandas as pd

df = pd.read_csv('sample300.csv')
df = df.rename(columns = {'A': 'education', 
                          'E': 'expertise', 
                          'label': 'success'})
label = 'success'
X, y = df.drop(columns=label), df[label]

# %%
from sklearn.feature_selection import chi2

scores, pvalues = chi2(X, y)
pd.Series(scores, index=X.columns)

# %%
import numpy as np
from scipy.sparse import issparse
from sklearn.utils import check_array
from sklearn.preprocessing import LabelBinarizer
from sklearn.utils.extmath import safe_sparse_dot

def preprocess_X_y(X, y):
    X = check_array(X, accept_sparse="csr")
    if np.any((X.data if issparse(X) else X) < 0):
        raise ValueError("Input X must be non-negative.")

    Y = LabelBinarizer().fit_transform(y)
    if Y.shape[1] == 1:
        Y = np.append(1 - Y, Y, axis=1)

    observed = safe_sparse_dot(Y.T, X)  # n_classes * n_features

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

f_obs, f_exp = preprocess_X_y(X, y)

# %%
from scipy.stats import chisquare

pd.Series(chisquare(f_obs, f_exp).statistic, index=X.columns)

chi2

education       1.412797
B              69.178778
C             232.577148
D             116.422861
expertise    1647.696011
F              24.250652
dtype: float64

chisquare

education       1.412797
B              69.178778
C             232.577148
D             116.422861
expertise    1647.696011
F              24.250652
dtype: float64

So I have to look at how the f_obs and f_exp are constructed because this is the only things that is different here.

[glemaitre]

OK so you are building the contingency matrix of the categories for each matrix while here we are computing the sum of the event, so equivalently, we have:

df.groupby("success").sum()

So we don't compute at all the same statistics. The reason is this feature selection is supposed only to work with binary indicators or words counts as mentioned in the documentation:

which must contain only non-negative features such as booleans or frequencies
(e.g., term counts in document classification), relative to the classes.

[altanova]

Thank you, first of all I appreciate that you read very thoroughly my case and even run your code on my own data. I think this was a lot of work.

In the same time, I note that your response seems to base on a logical error. I don't think it is correct to compute the statistics using this code: df.groupby("success").sum() . Let me give a simplistic example, based on my case. Suppose in Education column, we have 500 records with value 1 (the coded value 1 may indicate Oxford University), all of them Success = 0. And we also have 100 records with value 5 (indicating Cambridge), with Success = 0. Your contingency value will return the following result:

0: 500
1: 500

...indicating that there is no difference between the two, while in fact the difference is dramatic. This is because in my view it is a logical error to add (sum) numbers that in fact merely represent numeric codes of categorical variables.

But also if the data contains factual numbers, this logic will not work: as another example, let's consider the Iris dataset: imagine in one category we have 100 specimens with sepal length 5, and in the other 500 specimens with sepal lenght 1. In both case the sum is 500. So again, the contingency table defined as above will return the same results for both cases, suggesting there is no difference.

You also say this code supposed to only work on binary indicators or word counts. Indeed, I agree if all features were binary, it would perhaps work. But this is not what the documentation says. It says the data should be positive, not binary. Instead of "only binary", it says "such as binary":
https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.chi2.html

I believe it is a common understanding in the industry that selectKBest (based on sklearn's chi implementation) can be used to identify the best categorical feature (typically not binary). But if we adopt the logic of building the RC table as presented above, then this means this is not reliable if features are not binary. And this makes almost all known usages of SelectKBest () invalid. For example:

https://www.datatechnotes.com/2021/02/seleckbest-feature-selection-example-in-python.html

https://www.kaggle.com/jepsds/feature-selection-using-selectkbest?utm_campaign=News&utm_medium=Community&utm_source=DataCamp.com

https://programtalk.com/python-examples/sklearn.feature_selection.SelectKBest/

https://medium.com/mlearning-ai/how-to-select-features-using-selectkbest-in-python-c5a5239969f0

So if my understanding is correct, the sklearn's selectKBest() correctly performs the chi2 computation but before doing so it creates a RC table in an incorrect way, and hence the results are wrong. I would propose that building the RC table on df.groupby("success").sum() is a logical fault. But it is beyond my capability to look inside the sklearn source code and verify whether this problem is really responsible for the discrepancy in the results I am seing. However, from your description I conclude that we have hit an important problem that possibly impacts numerous users. Perhaps you could confirm or reject this idea?

[glemaitre]

I checked a bit more the (past) documentation, and the chi2 was really intended to work with binary features our the results of a CountVectorizer. Now, this is not super explicit and I agree that this is indeed an issue.

I think that we can reopen this issue and:

1. Improve the documentation in the meanwhile to mention that it is only for the above use case.
2. Implement an improvement that would build the proper RC table with a potential new keyword to switch how the matrix should be built. I assume that the current way is really efficient but restricted to the above cases.
3. Improve the user guide where we should explicitly show when to use this statistic.

[altanova]

Thank you, this is appreciated. I agree that until this is implemented, meanwhile the mention in the docs will be a good solution and I am happy to help propagate the information too. When it comes to the implementation, you may review my code for building the RC table, and if it is any useful feel free to use, improve it or incorporate it. Alternatively if you want me to improve my stuff, let me know. I know there are some deficiencies: https://github.com/altanova/stuff/tree/main/projs/2021-10-chi2

[glemaitre]

When it comes to the implementation, you may review my code for building the RC table, and if it is any useful feel free to use, improve it or incorporate it.

We will need a pure NumPy implementation :)

[altanova]

Agreed, my implementation uses Pandas so it is not usable here.

However I looked carefully and spotted one more issue with the existing code. One issue is what we already discussed (the statistics are implemented as if the features were binary). But there is another issue too. To correctly compute the chi square, one needs to implement the following condition: at least 80% of the RC contingency table cells need the expected count of 5 or more. Authors differ here, some say 80%, some say 100%, however there is general consensus that you cannot compute the chi square if too many cells have low expected frequencies:

https://libguides.library.kent.edu/spss/chisquare
https://ismayc.github.io/teaching/sample_problems/chi-square.html

Then what can be done? I believe that the RC rows with too low expected frequencies should be aggregated into one row, before calculating the chi2. However, the current implementation of SelectKBest does not meet this condition - this sort of calculation is not present. I checked experimentally that this makes quite a big difference in the calculation results.

Here is my implementation, now improved so that it now contains the correct (I hope!) calculation of this condition (majority of the cells >= 5 expected count), and it also has a parameter switch to disable this condition:
https://github.com/altanova/stuff/blob/main/projs/2021-10-chi2/chi2_util.py

And here is the notebook that demonstrates how to use it, and it also shows what happens if you disable the condition that majority of the cells must be >= 5. It shows that the results are quite different.
https://github.com/altanova/stuff/blob/main/projs/2021-10-chi2/chi2_usage.ipynb

Edit: having written the above, I realized that the problem I describe above might not apply to binary features. This problem will appear only if we work with categorical features with more possible values (thus more cells in the RC table, and then more possibility of cells with count <5). So in the end I think this bug still falls into category "chi2 should support categorical features".

[altanova]

also, to make this clear, from your code snippet I understood that the bug lies deeper than in SelectKBest.
The problem concerns sklearn.feature_selection.chi2, and that's where the fix needs to be applied.
The fact that it also exhibits in SelectKBest is secondary, as the latter method is a wrapper to the former.
This implies that until the fix is applied, the correct bypass is to use scipy.
This wasn't clear to me at first.