Python中的F-检验和P值

Question

R允许我们计算两个总体之间的F检验：

> d1 = c(2.5579227634, 1.7774243136, 2.0025207896, 1.9518876366, 0.0, 4.1984191803, 5.6170403364, 0.0)
> d2 = c(16.93800333, 23.2837045311, 1.2674791828, 1.0889208427, 1.0447584137, 0.8971380534, 0.0, 0.0)
> var.test(d1,d2)

    F test to compare two variances

data:  d1 and d2
F = 0.0439, num df = 7, denom df = 7, p-value = 0.000523
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.008789447 0.219288957
sample estimates:
ratio of variances 
        0.04390249 

请注意，这里也报告了P值。

另一个例子，R给出了以下内容：

> x1 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 68.7169110318)
> x2 = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.1863361211)
> var.test(x1,x2)
#p-value = 1.223e-09

在Python中的等价代码是什么？我查看了这个文档，但似乎没有给出我想要的。

这段代码给出了不同的P值（特别是第二个例子）：

import statistics as stats
import scipy.stats as ss
def Ftest_pvalue(d1,d2):
    """docstring for Ftest_pvalue"""
    df1 = len(d1) - 1
    df2 = len(d2) - 1
    F = stats.variance(d1) / stats.variance(d2)
    single_tailed_pval = ss.f.cdf(F,df1,df2)
    double_tailed_pval = single_tailed_pval * 2
    return double_tailed_pval

Python给出了以下结果：

In [45]: d1 = [2.5579227634, 1.7774243136, 2.0025207896, 1.9518876366, 0.0, 4.1984191803, 5.6170403364, 0.0]
In [20]: d2 = [16.93800333, 23.2837045311, 1.2674791828, 1.0889208427, 1.0447584137, 0.8971380534, 0.0, 0.0]
In [64]: x1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 68.7169110318]
In [65]: x2 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.1863361211]

In [69]: Ftest_pvalue(d1,d2)
Out[69]: 0.00052297887612346176

In [70]: Ftest_pvalue(x1,x2)
Out[70]: 1.9999999987772916

- pdubois

---

scipy 的 p 值似乎大约是 R 的一半（考虑到浮点表示和四舍五入问题）。这似乎意味着单尾检验与双尾检验的区别。 - lvc

@lvc的评论完全正确。如果您查看var.test的文档， 默认备择假设是双侧检验，而当您像Python一样计算cdf时，这本质上是一个单侧检验。 - Amit Kumar Gupta

@lvc：在我尝试了另一个示例之后，似乎并不是这样（请参见更新）。 - pdubois

在这种情况下，你不应该乘以2，而是应该先从1中减去一部分，然后再乘以2。在四舍五入后，2-1.9999999987772916等于1.223e-09。统计量F的单侧p值取决于F在平均值的哪一侧，它可以是cdf(F)或者1-cdf(F)。你试图衡量统计量比观察到的更"极端"的概率，如果F在平均值的左侧，"更极端"意味着"向左更远"，因此使用cdf(F)。如果F大于平均值，则"更极端"意味着"向右更远"，所以使用 1-cdf(F)。 - Amit Kumar Gupta

我不知道为什么scipy中不存在这个测试，但一个可能的原因（也是一个很好的替代方案）是人们更喜欢进行Levene检验，因为它被认为更加健壮。 - PlasmaBinturong

Answer 1

一个基于 rpy2 开发的实现：

import rpy2.robjects as robjects
def Ftest_pvalue_rpy2(d1,d2):
    """docstring for Ftest_pvalue_rpy2"""
    rd1 = (robjects.FloatVector(d1))
    rd2 = (robjects.FloatVector(d2))
    rvtest = robjects.r['var.test']
    return rvtest(rd1,rd2)[2][0]

有了这个结果：

In [4]: x1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 68.7169110318]
In [5]: x2 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.1863361211]
In [6]: Ftest_pvalue_rpy2(x1,x2)
Out[6]: 1.2227086010341282e-09

- neversaint

Answer 2

我应该提到xalglib是一个充满统计方法的包，可以实现以下功能： http://www.alglib.net/ http://www.alglib.net/hypothesistesting/variancetests.php 虽然它比基于scipy的原始方法不够灵活。

我应该提到，在variancetests.c中可以找到正确的双尾计算过程：

stat = ae_minreal(xvar/yvar, yvar/xvar, _state); *bothtails = 1-(fdistribution(df1, df2, 1/stat, _state)-fdistribution(df1, df2, stat, _state))

而@Amit Kumar Gupta在他的评论中所描述的是错误的（如果你仅仅将单侧p值与1之间的差异加倍，你可能会得到大于1的值）。

- Ando Jurai