重采样方法与 scipy.stats.chi2_contigency 的卡方检验 P 值

Question

本题参考《O'Relly Practical Statistics for Data Scientist 2nd Edition》一书第 3 章，卡方检验部分。

本书提供了一个卡方测试用例的示例，其中假设一个网站具有三个不同的标题，由 1000 名访问者运行。结果显示每个标题的点击次数。

观察到的数据如下：

Headline   A    B    C
Click      14   8    12
No-click   986  992  988

期望值的计算公式如下：

Headline   A        B        C
Click      11.13    11.13    11.13
No-click   988.67   988.67   988.67

皮尔逊残差定义为：

表现在的位置：

Headline   A        B        C
Click      0.792    -0.990   0.198
No-click   -0.085   0.106   -0.021

卡方统计量是 Pearson 残差平方和：。这是 1.666

到目前为止，一切都很好。现在是重采样部分：

1. Assuming a box of 34 ones and 2966 zeros
2. Shuffle, and take three samples of 1000 and count how many ones(Clicks)
3. Find the squared differences between the shuffled counts and expected counts then sum them.
4. Repeat steps 2 to 3, a few thousand times.
5. The P-value is how often does the resampled sum of squared deviations exceed the observed.

书中提供的重采样python测试代码如下：（可以从https://github.com/gedeck/practical-statistics-for-data-scientists/tree/master/python/code下载）

## Practical Statistics for Data Scientists (Python)
## Chapter 3. Statistial Experiments and Significance Testing
# > (c) 2019 Peter C. Bruce, Andrew Bruce, Peter Gedeck

# Import required Python packages.

from pathlib import Path
import random

import pandas as pd
import numpy as np

from scipy import stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats import power

import matplotlib.pylab as plt

DATA = Path('.').resolve().parents[1] / 'data'

# Define paths to data sets. If you don't keep your data in the same directory as the code, adapt the path names.

CLICK_RATE_CSV = DATA / 'click_rates.csv'

...

## Chi-Square Test
### Chi-Square Test: A Resampling Approach

# Table 3-4
click_rate = pd.read_csv(CLICK_RATE_CSV)
clicks = click_rate.pivot(index='Click', columns='Headline', values='Rate')
print(clicks)

# Table 3-5
row_average = clicks.mean(axis=1)
pd.DataFrame({
    'Headline A': row_average,
    'Headline B': row_average,
    'Headline C': row_average,
})

# Resampling approach
box = [1] * 34
box.extend([0] * 2966)
random.shuffle(box)

def chi2(observed, expected):
    pearson_residuals = []
    for row, expect in zip(observed, expected):
        pearson_residuals.append([(observe - expect) ** 2 / expect
                                  for observe in row])
    # return sum of squares
    return np.sum(pearson_residuals)

expected_clicks = 34 / 3
expected_noclicks = 1000 - expected_clicks
expected = [34 / 3, 1000 - 34 / 3]
chi2observed = chi2(clicks.values, expected)

def perm_fun(box):
    sample_clicks = [sum(random.sample(box, 1000)),
                     sum(random.sample(box, 1000)),
                     sum(random.sample(box, 1000))]
    sample_noclicks = [1000 - n for n in sample_clicks]
    return chi2([sample_clicks, sample_noclicks], expected)

perm_chi2 = [perm_fun(box) for _ in range(2000)]

resampled_p_value = sum(perm_chi2 > chi2observed) / len(perm_chi2)

print(f'Observed chi2: {chi2observed:.4f}')
print(f'Resampled p-value: {resampled_p_value:.4f}')

chisq, pvalue, df, expected = stats.chi2_contingency(clicks)
print(f'Observed chi2: {chi2observed:.4f}')
print(f'p-value: {pvalue:.4f}')

现在，我运行了 perm_fun(box) 2,000 次，并获得了 0.4775 的重采样 P 值。然而，如果我运行 perm_fun(box) 10,000 次和 100,000 次，我两次都能获得 0.84 的重采样 P 值。在我看来，P 值应该在 0.84 左右。为什么 stats.chi2_contigency 显示如此小的数字？

我运行2000次得到的结果是：

Observed chi2: 1.6659
Resampled p-value: 0.8300
Observed chi2: 1.6659
p-value: 0.4348

如果我运行它 10,000 次，结果是：

Observed chi2: 1.6659
Resampled p-value: 0.8386
Observed chi2: 1.6659
p-value: 0.4348

软件版本：

pandas.__version__:         0.25.1
numpy.__version__:          1.16.5
scipy.__version__:          1.3.1
statsmodels.__version__:    0.10.1
sys.version_info:           3.7.4

Answer 1

I ran your code trying 2000, 10000, and 100000 loops, and all three times I got close to .47. I did, however, get an error at this line that I had to fix:

resampled_p_value = sum(perm_chi2 > chi2observed) / len(perm_chi2)

这里perm_chi2 is a list and chi2observed is a float, so I wonder how this code ever ran for you (perhaps whatever you did to fix it was the source of the error). In any case, changing it to the intended

resampled_p_value = sum([1*(x > chi2observed) for x in perm_chi2]) / len(perm_chi2)

allowed me to run it and get close to .47.

确保在更改迭代次数时，仅更改 2000，而不更改其他数字。