我有以下数据:
array([[33, 250, 196, 136, 32],
[55, 293, 190, 71, 13]])
我可以从 stats.chi2_contingency(data) 中获取 p 值。
有没有类似于此 R 对象的东西 - data.chisq$residuals 来获取 Pearson 残差和标准化残差？
请您参考如下方法:
如果你不介意依赖，statsmodels有一个模块 contingency table calculations .例如，
In [2]: import numpy as np
In [3]: import statsmodels.api as sm
In [4]: F = np.array([[33, 250, 196, 136, 32], [55, 293, 190, 71, 13]])
In [5]: table = sm.stats.Table(F)
In [6]: table.resid_pearson # Pearson's residuals
Out[6]:
array([[-1.77162519, -1.61362277, -0.05718356, 2.96508777, 1.89079393],
[ 1.80687785, 1.64573143, 0.05832142, -3.02408853, -1.92841787]])
In [7]: table.standardized_resids # Standardized residuals
Out[7]:
array([[-2.62309082, -3.0471942 , -0.09791681, 4.6295814 , 2.74991911],
[ 2.62309082, 3.0471942 , 0.09791681, -4.6295814 , -2.74991911]])
如果您不想依赖 statsmodels，可以使用 scipy.stats.chi2_contingency 的结果在几行中实现这些计算。这是一个为这些残差定义函数的简短模块。它们采用观察到的频率和预期频率(由 chi2_contingency 返回)。请注意，虽然 chi2_contingency 和以下 residuals 函数适用于 n 维数组，但此处实现的 stdres 仅适用于二维数组。
from __future__ import division
import numpy as np
from scipy.stats.contingency import margins
def residuals(observed, expected):
return (observed - expected) / np.sqrt(expected)
def stdres(observed, expected):
n = observed.sum()
rsum, csum = margins(observed)
# With integers, the calculation
# csum * rsum * (n - rsum) * (n - csum)
# might overflow, so convert rsum and csum to floating point.
rsum = rsum.astype(np.float64)
csum = csum.astype(np.float64)
v = csum * rsum * (n - rsum) * (n - csum) / n**3
return (observed - expected) / np.sqrt(v)
根据您的数据，我们得到:
>>> F = np.array([[33, 250, 196, 136, 32], [55, 293, 190, 71, 13]])
>>> chi2, p, dof, expected = chi2_contingency(F)
>>> residuals(F, expected)
array([[-1.77162519, -1.61362277, -0.05718356, 2.96508777, 1.89079393],
[ 1.80687785, 1.64573143, 0.05832142, -3.02408853, -1.92841787]])
>>> stdres(F, expected)
array([[-2.62309082, -3.0471942 , -0.09791681, 4.6295814 , 2.74991911],
[ 2.62309082, 3.0471942 , 0.09791681, -4.6295814 , -2.74991911]])
这里是 R 中用于比较的计算:
> F <- as.table(rbind(c(33, 250, 196, 136, 32), c(55, 293, 190, 71, 13)))
> result <- chisq.test(F)
> result$residuals
A B C D E
A -1.77162519 -1.61362277 -0.05718356 2.96508777 1.89079393
B 1.80687785 1.64573143 0.05832142 -3.02408853 -1.92841787
> result$stdres
A B C D E
A -2.62309082 -3.04719420 -0.09791681 4.62958140 2.74991911
B 2.62309082 3.04719420 0.09791681 -4.62958140 -2.74991911
