为什么scipy.stats.mstats.pearsonr的结果与scipy.stats.pearsonr不一致？

Question

我原本以为，使用带有掩码的数组输入时，scipy.stats.mstats.pearsonr的结果应该和scipy.stats.pearsonr对未掩码数据的结果一样，但实际上并不是这样的：

from pylab import randn,rand
from numpy import ma
import scipy.stats

# Normally distributed data with noise
x=ma.masked_array(randn(10000),mask=False)
y=x+randn(10000)*0.6

# Randomly mask one tenth of each of x and y
x[rand(10000)<0.1]=ma.masked
y[rand(10000)<0.1]=ma.masked

# Identify indices for which both data are unmasked
bothok=((~x.mask)*(~y.mask))

# print results of both functions, passing only the data where 
# both x and y are good to scipy.stats
print "scipy.stats.mstats.pearsonr:", scipy.stats.mstats.pearsonr(x,y)[0]
print "scipy.stats.pearsonr:", scipy.stats.pearsonr(x[bothok].data,y[bothok].data)[0]

每次运行这个代码，结果会有一点不同，但对我来说，差异大约是0.1。而且，掩码的比例越大，结果的差异就越明显。

我发现，如果对x和y使用相同的掩码，那么这两个函数的结果是一样的，也就是说：

mask=rand(10000)<0.1
x[mask]=ma.masked
y[mask]=ma.masked
...

这算是一个bug吗？还是说我应该提前处理输入数据，确保x和y的掩码是完全一样的（这肯定不应该是这样吧）？

我使用的是numpy版本'1.8.0'和scipy版本'0.11.0b1'

Answer 1

我们有两种处理缺失值的方法，分别是完全删除缺失的情况和成对删除。

在你使用 scipy.stats.pearsonr 时，如果任何变量有缺失值，就会完全删除这些情况。

numpy.ma.corrcoef 也会给出相同的结果。

查看 scipy.stats.mstats.pearsonr 的源代码，它在计算方差或均值时并不会完全删除缺失的情况。

>>> xm = x - x.mean(0)
>>> ym = y - y.mean(0)
>>> np.ma.dot(xm, ym) / np.sqrt(np.ma.dot(xm, xm) * np.ma.dot(ym, ym))
0.7731167378113557

>>> scipy.stats.mstats.pearsonr(x,y)[0]
0.77311673781135637

不过，完全删除和成对删除在均值和标准差上的差别很小。

主要的差异似乎是因为没有对非缺失元素的数量进行修正。如果不考虑自由度的修正，我得到的结果是

>>> np.ma.dot(xm, ym) / bothok.sum() / \
    np.sqrt(np.ma.dot(xm, xm) / (~xm.mask).sum() * np.ma.dot(ym, ym) / (~ym.mask).sum())
0.85855728319303393

这个结果接近于完全删除缺失情况的结果。

Answer 2

这看起来像是 scipy.stats.mstats.pearsonr 里的一个错误。它似乎是说，x 和 y 这两个值应该按索引成对出现，也就是说，如果其中一个被屏蔽了，另一个也应该被忽略。也就是说，如果 x 和 y 看起来像这样（用 -- 表示被屏蔽的值）：

x = [1, --,  3,  4,  5]
y = [9,  8, --,  6,  5]

那么 (--, 8) 和 (3, --) 都应该被忽略，结果应该和 scipy.stats.pearsonr([1, 4, 5], [9, 6, 5]) 一样。

在 这个 mstats 版本 中的错误是，计算 x 和 y 的平均值的代码没有使用共同的屏蔽。

我在 scipy 的 GitHub 网站上为此创建了一个问题：https://github.com/scipy/scipy/issues/3645