十一、非参数检验：使用python进行卡方两样本独立性检验

在单样本问题中, 人们想要检验的是总体的中心是否等于一个已知的值. 但在实际问题中, 更受注意的往往是比较两个总体的位置参数; 比如, 两种训练方法中哪一种更出成绩, 两种汽油中哪一种污染更少, 两种市场营销策略中哪种更有效等等.

1. ${\chi }^2$独立性检验的原理

若随机变量$X,Y$的分布函数分别为$F_1(x)和F_2(y)$, 且联合分布为$F(x,y)$, 则X与Y的独立性归结为假设检验问题：
$$H_0:F(x,y)=F_1(x)F_2(y)H_1:F(x,y)\ne F_1(x)F_2(y)$$
若X与Y为分类变量，其中X的取值为$X_1,X_2,...,X_r$, Y的取值为$Y_1,Y_2,...,Y_s$, 将X与Y的各种情况组合用一张$r\times s$列联表表示，称为$r\times s$二维列联表，如下图所示：

表中$n_{ij}$表示n个随机试验的结果中X取$X_i$及Y取$Y_j$的频数，${\sum }_{i=1}^r{\sum }_{j=1}^s{n}_{ij}=n$.
$$n_{i.}=\sum _{j=1}^s{n}_{ij},i=1,2,...,r,表示各行之和$$
$$n_{.j}=\sum _{i=1}^r{n}_{ij},j=1,2,...s,表示各列之和$$
令$p_{ij}=P(X=X_i,Y=Y_j),p_{i.}=P(X=X_i),p_{.j}=P(Y=Y_j),i=1,2,...,r;j=1,2,...,s,$ 则X与Y的独立性检验等价于下述检验：
$$H_0:p_{ij}=p_{i.}{p}_{.j},\forall 1\le i\le r,1\le j\le s{H}_1:\exists (i,j),p_{ij}\ne p_{i.}{p}_{.j}$$
注: 若X与Y 为连续型随机变量, 这时将它们的取值范围分成r个及s个互不相交的小区间, 用$n_{ij}$ 表示n个随机试验的结果中“X属于第i个小区间, Y 属于第k个小区间”的频数 $(i=1,2,...,r;k=1,2,...,s)$ 这时可将X与Y 的独立性转化为列联表的独立性检验问题.
由于$p_{i.}$的极大似然估计为${\hat{p}}_{i.}=n_{i.}/n$, $p_{.j}$的极大似然估计为${\hat{p}}_{.j}/n$, 因此若$H_0$成立，则$p_{ij}$的极大似然估计为${\hat{p}}_{i.}{\hat{p}}_{.j}=n_{i.}{n}_{.j}/n^2$. 从而X取代$X_i$, Y取代$Y_j$(试验数据落入第$(i,j)$个类)的理论频数为$n\times n_{i.}{n}_{.j}/n^2=n_{i.}{n}_{.j}/n$. 由此构造检验统计量为：
$${\chi }^2=\sum _{i=1}^r\sum _{k=1}^s[n_{ij}-\frac{n_{i.}{n}_{.j}}{n}{]}^2/\frac{n_{i.}{n}_{.j}}{n}$$
可以证明在原假设成立时， ${\chi }^2$近似服从${\chi }^2((r-1)(s-1))$

2. 使用Python展示和计算

在python中的pandas库中有一个DataFrame的数据结构，可以用来保存二维数据，并进行计算。
例1： 对表1所示频数分布表，以59%显著水平，检验色觉与性别是否相关。

2. 1. 初始化数据

import numpy as np
import pandas as pd
from scipy import stats
d=np.array([[442, 514],[38,6]])
r,s = d.shape
df = pd.DataFrame(d, columns=['male', 'female'], index=['normal', 'blindness'])

上述代码已将初始化数据存入df，显示如下：

2.2. 计算边际频数

首先在df基础上建立一个数据框（数据一致）。并在新建的数据框基础继续计算。

df1 = df[:]
df1= pd.DataFrame(df1)
df1['r_tot']=np.sum(df1,axis=1)
df1.loc['c_tot']=np.sum(df1, axis=0)

经过上面的处理，df1就是带有编辑频数的数据框。显示如下：

2.3. 计算理论频数分布表

r_tot = df1['r_tot'][:-1]
c_tot=df1.loc['c_tot'][:-1]
total = np.sum(r_tot)
data = np.zeros((r,s))
for i in range(len(r_tot)):
    for j in range(len(c_tot)):
        data[i,j]=r_tot[i]*c_tot[j]/total
df2=pd.DataFrame(data, index=['normal', 'blindness'], columns=['male', 'female'])

至此，我们将计算得到的理论频数分布表保存在了df2中，而原始数据在df中，如下图所示：

2.4. 统计量计算

为了清晰整个计算过程，这里添加一步（对计算结果来说是多余的）产生统计量数据框的代码：

(df-df2)**2/df2

其产生的统计量数据框如下图：

其实，下面一行代码即可获得卡方统计量的值并获取pval：

chi_square=np.sum((df-df2)**2/df2).sum()
stats.chi2.sf(chi_squre, df=(r-1)*(s-1))

得到${\chi }^2=27.13874340443378,pval=1.8936459120177876e-07$
因为$pval<\alpha =0.05$, 所以拒绝原假设，可以认为色盲和性别有关系