Pearson chi2 tests of independence: differences between Scipy and Scikit-learn 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])
# )
```


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.
