Pearson chi2 tests of independence: differences between Scipy and Scikit-learn

Question

Scipy and Scikit-learn both implement Pearson's $\chi^2$ test of independence, but they give different results. The former matches what you might expect when computing this test "by hand", the latter does not. What nonstandard assumption or extension is used in the Scikit-learn version that causes the discrepancy?

---

This issue was raised in a mailing list thread, but the resolution was a bit hazy to me:

Lars Buitinck wrote:

The difference seems (thinking out loud) to stem from assumptions about the input. feature_selection.chi2 (implicitly) assumes a multinomial event model, so each X[i, j] is the frequency with which event j was observed when drawing X[i].sum() times from a multinomial. A zero input value is interpreted as the absence of an event, rather than a separate 0 event.

Christian Jauvin wrote (paraphrased for clarity):

If I understand you correctly, one way to reconcile the difference between the two interpretations (multinomial vs binomial) would be to x = np.append(1 - x, x, axis=1). Summing the chi-squared value of each feature (1.5 + 0.5) then yields the same result as obtained with scipy.stats.chi2_contingency. Does that make sense?

Lars never replied to this post, and Christian never followed up, so I assume that Christian was satisfied and/or they both forgot about the thread.

---

Here are both implementations, for reference:

• Scipy v1.7.1
• Scikit-learn v1.0

---

This is how I always learned to compute the "expected" number of data points in a sample with category $i$ in the first variable and category $j$ in the second:

$$ E_ij = N p_{i \cdot} p_{\cdot j} $$

Where:

• $N$ is the number of observations in the sample
• $p_{i \cdot}$ is the marginal proportion of category $i$ in the first variable $X$
• $p_{\cdot j}$ is the marginal proportion of category $j$ in the second variable $Y$

We use $p_{i \cdot}$ to estimate $\Pr \left( X = i \right)$, and we use $p_{\cdot j}$ to estimate $\Pr \left(Y = j \right)$, and we use $N p_{i \cdot} p_{\cdot j}$ as our estimate of $\Pr \left( X = i \right) \Pr \left( Y = j \right)$, which under the independence assumption will be $\Pr \left( X = i \land Y = j \right)$.

---

Here is a cleaned-up version of the code in the mailing list thread, showing the different results:

import numpy as np
import pandas as pd
from sklearn.feature_selection import chi2
from sklearn.preprocessing import LabelBinarizer
from scipy.stats import chi2_contingency


## Construct a sample of 2 binary variables

data = pd.DataFrame(np.vstack((
    [[0, 0]] * 18,
    [[0, 1]] * 7,
    [[1, 0]] * 42,
    [[1, 1]] * 33
)), columns=['x', 'y'])
x = data['x']
y = data['y']


## Compute a cross-tab

xtab_xy = pd.crosstab(x, y)

# Scipy chi2 test
sp_chi2_val, sp_chi2_p, sp_chi2_dof, sp_chi2_exp = chi2_contingency(xtab_xy)
print((sp_chi2_val, sp_chi2_p, sp_chi2_dof, sp_chi2_exp))
# Output:
#  (
#    1.3888888888888888,
#    0.2385928293164321,
#    1,
#    array([[15., 10.], [45., 30.]])
#  )



## Scikit-learn chi2 test

sk_chi2_val, sk_chi2_p = chi2(x.to_frame(), y)
print((sk_chi2_val, sk_chi2_p))
# Output:
# (
#   array([0.5]),
#   array([0.47950012])
# )
```

Answer 1

This is not entirely a statistics question, and mainly a programming one. To answer the statistics aspect, this is very simple once you understand what sklearn is doing: the chi2 function performs a goodness-of-fit test on your data, not a chi-squared test of independence. This is why you end up with different results.

But since you probably need some additional explanation not to take my answer at face value, and since we're at it anyway, here is a short explanation on the programming aspect (i.e. how sklearn performs this goodness-of-fit test): the observed values are the observed counts of 0/1 in $y$ filtered by $x=1$, and the expected values are the proportions of 0/1 in $y$ multiplied by the sum of 1s in $x$.

Below is a piece of code to reconcile scipy and sklearn, which will probably make the explanation above much clearer. The variable data comes from the code you provided in your question, and the scipy chisquare function performs a goodness-of-fit test:

> from scipy import stats
> observed = data[data["x"] == 1]["y"].value_counts()
> print(observed)
0    42
1    33
Name: y, dtype: int64

>expected = data["y"].value_counts(normalize=True)*observed.sum()
>print(expected)
0    45.0
1    30.0
Name: y, dtype: float64

>result = stats.chisquare(observed, f_exp=expected)
>print(result)
Power_divergenceResult(statistic=0.5, pvalue=0.47950012218695337)

which is in line with the sklearn output that you mention in your question.

Answer 2

χ² (chi-squared) statistic of scipy.stats.chi2_contingency vs sklearn.feature_selection.chi2

It appears from reading Scikit-learn χ² (chi-squared) statistic and corresponding contingency table that sklearn does not perform a standard contingency table analysis when calculating the χ² statistic between two categorical variables.

For example, given the data below

# libraries
import scipy.stats as sps
import pandas as pd

# data
data = pd.DataFrame({'gender': ['female']*60+['female']*54+['female']*46+['female']*41+['male']*40+['male']*44+['male']*53+['male']*57,
                     'education level': ['high school']*60+['bachelors']*54+['masters']*46+['phd']*41+['high school']*40+['bachelors']*44+['masters']*53+['phd']*57})

# contingency table
contingency_table = pd.crosstab(index=data['gender'],
                                columns=data['education level'])

the χ² statistic

print sps.chi2_contingency(contingency_table, correction=False)

is 8.0060 using scipy.

If we were to instead use sklearn

# libraries
from sklearn.preprocessing import LabelEncoder
from sklearn.feature_selection import chi2

# label encoder for categorical features
le = LabelEncoder()

# transforming the categorical features
data['education level le'] = pd.DataFrame(le.fit_transform(data['education level']))
data['gender le'] = pd.DataFrame(le.fit_transform(data['gender']))

we would obtain

chi_2, p_value = chi2(data['education level le'].values.reshape(-1, 1),
                   data['gender le'].values)

print chi_2

a value of 4.6557.

I am not entirely sure which of these two methods is the appropriate one to use in order to determine whether the two categorical features are independent or not. Related question 549371.

Update (3/5/23): It looks like the contingency matrix is a 2x2 matrix and you get a single p-val for the set of all the features in producing the target, whereas the feature_selection chi2 shows the p-vals for each feature independently of the others in producing the target. I believe that is why the two functions do not produce the same results.