卡方检验及其Python实现

• 分类数据的 拟合优度检验

• 独立性检验

分类数据的 拟合优度检验

前面我已经写了关于几种常见的假设检验内容，而 检验主要是测试样本分类数据的分布是否符合预期分布。相信大家如果学过高中生物，都知道孟德尔——遗传学之父，当时他根据颜色和形状把豌豆分为四类：黄圆、绿圆、黄皱和绿皱.孟德尔根据遗传学原理判断这四类的比例应为9:3:3:1.为做验证,孟德尔分别统计了这四类豌豆的个数，正是利用 检验证明了这令人激动的结论

在处理分类数据时，这些类别值本身对统计检验没有多大用处，比如像“男性”、“女性”和“其他”这样的类别数据没有任何数学意义。所以处理分类变量的检验是基于变量计数，而不是变量本身的实际值。

下面通过生成一些虚假的人口统计数据，并通过 检验来检验它们是否不同：

import numpy as np
import pandas as pd
import scipy.stats as stats

national = pd.DataFrame(["white"]*100000 + ["hispanic"]*60000 +\
                        ["black"]*50000 + ["asian"]*15000 + ["other"]*35000)
           

minnesota = pd.DataFrame(["white"]*600 + ["hispanic"]*300 + \
                         ["black"]*250 +["asian"]*75 + ["other"]*150)

national_table = pd.crosstab(index=national[0], columns="count")
minnesota_table = pd.crosstab(index=minnesota[0], columns="count")

print( "National")
print(national_table)
print(" ")
print( "Minnesota")
print(minnesota_table)

National
col_0      count
0               
asian      15000
black      50000
hispanic   60000
other      35000
white     100000
 
Minnesota
col_0     count
0              
asian        75
black       250
hispanic    300
other       150
white       600

检验是基于 检验统计量。使用以下公式计算检验统计量的值：

observed = minnesota_table

national_ratios = national_table/len(national)  # 实际值

expected = national_ratios * len(minnesota)   # 理论值

chi_squared_stat = (((observed-expected)**2)/expected).sum()

print(chi_squared_stat)

col_0
count    18.194805
dtype: float64

检验假设所有预期计数均不小于5,如果某一类别的个数小于5，就将相邻的某些类别合成为一类。

拒绝域：W={ },其实r为类别数，a为显著性水平

crit = stats.chi2.ppf(q = 0.95, # 找到95%置信度的临界值
                      df = 4)   # 自由度个数

print("Critical value")
print(crit)

p_value = 1 - stats.chi2.cdf(x=chi_squared_stat,  # P值
                             df=4)
print("P value")
print(p_value)

Critical value
9.487729036781154
P value
[0.00113047]

由于检验统计量大于P值，所以得出结论，有95%的把握认为上述两个总体的分布不是相同的。

当然也可以使用scipy.stats.chisquare()函数，十分快捷！

stats.chisquare(f_obs= observed,   # 观察值
                f_exp= expected)   # 理论值

Power_divergenceResult(statistic=array([18.19480519]), pvalue=array([0.00113047]))

独立性检验

独立性检验是统计学的另一种检验方式，它是根据次数判断两类变量彼此相关或相互独立的假设检验。下面生成一些虚假的选民投票数据并进行独立性测试，用于确定教育、政治观点和其他偏好等变量是否因性别、种族和宗教等人口因素而有所不同：

np.random.seed(10)

voter_race = np.random.choice(a= ["asian","black","hispanic","other","white"],
                              p = [0.05, 0.15 ,0.25, 0.05, 0.5],
                              size=1000)

voter_party = np.random.choice(a= ["democrat","independent","republican"],
                              p = [0.4, 0.2, 0.4],
                              size=1000)

voters = pd.DataFrame({"race":voter_race, 
                       "party":voter_party})

voter_tab = pd.crosstab(voters.race, voters.party, margins = True)

voter_tab.columns = ["democrat","independent","republican","row_totals"]

voter_tab.index = ["asian","black","hispanic","other","white","col_totals"]

observed = voter_tab.iloc[0:5,0:3]   
print(voter_tab)

            democrat  independent  republican  row_totals
asian             21            7          32          60
black             65           25          64         154
hispanic         107           50          94         251
other             15            8          15          38
white            189           96         212         497
col_totals       397          186         417        1000

对于独立性测试，使用与拟合优度检验相同的检验统计量。主要区别在于，独立性检验必须在二维表格中计算每个单元格的预期计数，而不是一维表格。要获得单元格的预期计数，需要将该单元格的行总计乘以该单元格的列总计，然后除以观察的总数。可以通过np.outer()除以总的观察数快速获得表中所有单元格的理论值

expected =  np.outer(voter_tab["row_totals"][0:5],
                     voter_tab.loc["col_totals"][0:3]) / 1000

expected = pd.DataFrame(expected)

expected.columns = ["democrat","independent","republican"]
expected.index = ["asian","black","hispanic","other","white"]

print(expected)

          democrat  independent  republican
asian       23.820       11.160      25.020
black       61.138       28.644      64.218
hispanic    99.647       46.686     104.667
other       15.086        7.068      15.846
white      197.309       92.442     207.249

现在可以按照之前相同的步骤来计算检验统计量，临界值和p值：

chi_squared_stat = (((observed-expected)**2)/expected).sum().sum()

print(chi_squared_stat)

7.169321280162059

注意：调用此处使用sum()方法两次：第一次是获取列和，第二次是将列和相加，返回整个二维表的总和。

crit = stats.chi2.ppf(q = 0.95, #找到95%置信度的临界值
                      df = 8)   

print("Critical value")
print(crit)

p_value = 1 - stats.chi2.cdf(x=chi_squared_stat,  # P值
                             df=8)
print("P value")
print(p_value)

Critical value
15.50731305586545
P value
0.518479392948842

独立性测试的自由度等于每个变量中类别数减去1的乘积。在本例中，有一个5x3表，因此df=4x2=8。

同样可以使用scipy快速进行独立性测试

stats.chi2_contingency(observed= observed)

(7.169321280162059, 0.518479392948842, 8, array([[ 23.82 ,  11.16 ,  25.02 ],
        [ 61.138,  28.644,  64.218],
        [ 99.647,  46.686, 104.667],
        [ 15.086,   7.068,  15.846],
        [197.309,  92.442, 207.249]]))

输出检验统计量的值、p值和自由度以及理论值矩阵。

7.169321280162059<15.50731305586545,落入接受域，故认为上述两变量之间无显著关系。