了解排列测试

Question

我正在尝试优化 scipy.stats 中实现的排列测试的性能。我的数据集包含 500,000 个观察值，每个观察值与 2,000 个二元协变量相关。我对这些观测值应用了大约 10,000 个不同的过滤器，因此需要在 500,000 个观测值中计算 20,000,000 次排列测试。

我的方法涉及使用 CuPy (500,000 x 2,000) 生成 2D 布尔数组。然后，我将此二维“协变量矩阵”按元素乘以测试统计量（浮点值）的一维向量，并沿适当的轴求和，以获得每个协变量的“参考统计量”。为了执行排列测试，我将此参考统计量与从 2D 协变量矩阵的 n 排列获得的统计量分布进行比较。

为了提高计算效率，我实现了一种生成置换协变量矩阵的替代方法。我没有使用 cp.random.permutation(covariates)，而是使用 cp.random.rand(*covariates.shape) < proportion_of_ones, which is 2-3 times faster. I hypothesized that this approach would be valid because, while the number of 'ones' in each permuted matrix might vary, the distribution should be symmetric around the original number of ones (assuming a sufficiently large n)。

我根据经验验证了这一假设，并发现在许多情况下，我的假设似乎成立，两种方法之间的估计 p 值通常相差不到 0.01。然而，当我对检验统计量的分布进行轻微修改时，我的假设不再成立，两种方法的结果出现了显着差异。我无法确定造成这种差异的原因，并发现自己陷入了僵局。

这是我的参考代码（将 cp 替换为 np，它也应该可以工作）：

from scipy.stats import ttest_ind

import pandas as pd
import numpy as np
import cupy as cp

from tqdm import tqdm


def test_permu_1(
    score: cp.ndarray, state: cp.ndarray, permu: int = 1000000
) -> cp.ndarray:

    value = cp.abs(cp.sum(score * state, axis=1))

    fraction_affected = cp.sum(state, axis=1) / state.shape[1]

    stats = []
    for _ in tqdm(range(permu), total=permu, desc="Permuting 1", ncols=120):
        mask = cp.random.rand(*state.shape, dtype=cp.float32)
        mask = mask < fraction_affected[:, None]
        stats.append(cp.sum(score * mask, axis=1))

    stats = cp.array(stats)

    return cp.sum(value < cp.abs(stats), axis=0) / permu


def test_permu_2(
    score: cp.ndarray, state: cp.ndarray, permu: int = 1000000
) -> cp.ndarray:

    value = cp.abs(cp.sum(score * state, axis=1))

    # NOTE, cp.random.permutation doesn't support axis = 1
    state = state.transpose()
    score = score[:, None]

    stats = []
    for _ in tqdm(range(permu), total=permu, desc="Permuting 2", ncols=120):
        mask = cp.random.permutation(state)
        stats.append((score * mask).sum(axis=0))

    stats = cp.array(stats)

    return cp.sum(value < cp.abs(stats), axis=0) / permu


def ttest_ind_0(score: cp.ndarray, state: cp.ndarray) -> cp.ndarray:
    np_score = score.get()
    np_state = state.get()

    ttests = []
    for i in range(np_state.shape[0]):
        tp = np_score[np_state[i]]
        tn = np_score[~np_state[i]]

        pv = ttest_ind(tp, tn, equal_var=False)
        ttests.append(pv.pvalue)  # type: ignore

    return np.array(ttests)


def dumpy(p0: np.ndarray, p1: np.ndarray, p2: np.ndarray) -> pd.DataFrame:
    df = pd.DataFrame(
        {
            "p0": p0s,
            "p1": p1s,
            "p2": p2s,
            "-": ["."] * ysize,
            "d10": p0s - p1s,
            "d12": p1s - p2s,
        }
    ).T

    with pd.option_context("float_format", "{:.5f}".format):
        print(df)

    return df


xsize = 100000
ysize = 18

# Pull scores from a normal distribution
score = cp.random.randn(xsize) + 10
state = cp.random.rand(ysize, xsize, dtype=cp.float32) < 0.001

# Perform the tests
permutations = 10000
p0s = ttest_ind_0(score, state)
p1s = test_permu_1(score, state, permutations).get()
p2s = test_permu_2(score, state, permutations).get()

dumpy(p0s, p1s, p2s)


print("Done!")

注意写着

score = cp.random.randn(xsize) + 10

的行...如果我删除

+ 10

，一切都会“有效”，正如我的假设似乎成立一样。有了这条线，两种不同的方法就“分歧”了......

这是一个不使用

+ 10

运行的示例：

Permuting 1: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [00:30<00:00, 329.59it/s]
Permuting 2: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [01:08<00:00, 146.35it/s]
          0        1       2       3        4        5        6        7        8        9        10       11      12       13       14       15       16      17
p0   0.58190  0.27232 0.34934 0.21264  0.85781  0.09202  0.20862  0.24606  0.01390  0.59029  0.28032  0.71163 0.96128  0.86587  0.26078  0.77700  0.70671 0.62766
p1   0.61460  0.22130 0.32700 0.20900  0.82670  0.05520  0.20360  0.25450  0.02340  0.59310  0.30230  0.71890 0.93790  0.89820  0.27080  0.79350  0.71250 0.60850
p2   0.61070  0.22810 0.32330 0.19780  0.83530  0.05580  0.20810  0.25300  0.02360  0.58910  0.29630  0.72530 0.93780  0.89760  0.26720  0.79580  0.72110 0.60660
-          .        .       .       .        .        .        .        .        .        .        .        .       .        .        .        .        .       .
d10 -0.03270  0.05102 0.02234 0.00364  0.03111  0.03682  0.00502 -0.00844 -0.00950 -0.00281 -0.02198 -0.00727 0.02338 -0.03233 -0.01002 -0.01650 -0.00579 0.01916
d12  0.00390 -0.00680 0.00370 0.01120 -0.00860 -0.00060 -0.00450  0.00150 -0.00020  0.00400  0.00600 -0.00640 0.00010  0.00060  0.00360 -0.00230 -0.00860 0.00190

这是使用

+ 10

进行的两次运行。请注意，这两次运行是在开始时生成的相同虚拟数据的复制（我刚刚在调试模式下运行了两次排列步骤）。我想看看这是否只是由于将 10 添加到 -1 到 1 范围内的值而导致较小数据范围内的较高变异性，但排列结果在每种方法中都是一致的。我讨厌这样无知！

Permuting 1: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [00:30<00:00, 329.59it/s]
Permuting 2: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [01:08<00:00, 146.27it/s]
          0        1        2       3        4        5       6        7       8       9        10       11       12      13       14       15       16      17
p0   0.29352  0.25426  0.33900 0.47577  0.30116  0.97216 0.81252  0.67735 0.57792 0.98254  0.38642  0.01012  0.56308 0.60787  0.46944  0.10419  0.57390 0.65886
p1   0.44620  0.45500  0.52370 0.46790  0.53630  0.50050 0.48880  0.51110 0.47350 0.49980  0.45840  0.59830  0.52310 0.46390  0.52090  0.55280  0.51700 0.48080
p2   0.13040  0.15300  0.80300 0.26190  0.87330  0.51890 0.40320  0.68480 0.28030 0.49180  0.17170  0.99550  0.69950 0.31450  0.75990  0.93470  0.73290 0.31700
-          .        .        .       .        .        .       .        .       .       .        .        .        .       .        .        .        .       .
d10 -0.15268 -0.20074 -0.18470 0.00787 -0.23514  0.47166 0.32372  0.16625 0.10442 0.48274 -0.07198 -0.58818  0.03998 0.14397 -0.05146 -0.44861  0.05690 0.17806
d12  0.31580  0.30200 -0.27930 0.20600 -0.33700 -0.01840 0.08560 -0.17370 0.19320 0.00800  0.28670 -0.39720 -0.17640 0.14940 -0.23900 -0.38190 -0.21590 0.16380
Permuting 1: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [00:30<00:00, 329.57it/s]
Permuting 2: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [01:08<00:00, 145.83it/s]
          0        1        2       3        4        5       6        7       8       9        10       11       12      13       14       15       16      17
p0   0.29352  0.25426  0.33900 0.47577  0.30116  0.97216 0.81252  0.67735 0.57792 0.98254  0.38642  0.01012  0.56308 0.60787  0.46944  0.10419  0.57390 0.65886
p1   0.45250  0.45430  0.53450 0.47370  0.53260  0.49320 0.48450  0.51640 0.47310 0.49860  0.45260  0.59260  0.51230 0.46290  0.51770  0.55190  0.51670 0.47340
p2   0.12270  0.14980  0.79630 0.25370  0.87600  0.51530 0.40710  0.67740 0.28180 0.48010  0.17010  0.99530  0.70230 0.31050  0.77380  0.93760  0.73010 0.31800
-          .        .        .       .        .        .       .        .       .       .        .        .        .       .        .        .        .       .
d10 -0.15898 -0.20004 -0.19550 0.00207 -0.23144  0.47896 0.32802  0.16095 0.10482 0.48394 -0.06618 -0.58248  0.05078 0.14497 -0.04826 -0.44771  0.05720 0.18546
d12  0.32980  0.30450 -0.26180 0.22000 -0.34340 -0.02210 0.07740 -0.16100 0.19130 0.01850  0.28250 -0.40270 -0.19000 0.15240 -0.25610 -0.38570 -0.21340 0.15540

我刚刚将 

n

的值增加了 10 倍并重新运行（只有 3 个速度条件）：

Permuting 1: 100%|███████████████████████████████████████████████████████████████| 10000/10000 [00:54<00:00, 185.17it/s]
Permuting 2: 100%|████████████████████████████████████████████████████████████████| 10000/10000 [03:20<00:00, 49.85it/s]
          0        1        2
p0  0.84096  0.06036  0.28752
p1  0.48310  0.43160  0.54530
p2  0.43590  0.03200  0.85530
-         .        .        .
d10 0.35786 -0.37124 -0.25778
d12 0.04720  0.39960 -0.31000

Answer 1

排列检验的零假设本质上是观察样本的顺序没有什么特殊之处（因此零分布是通过随机对观察值重新排序而生成的）。

蒙特卡罗检验的零假设是您的样本是伯努利分布的一系列观察结果（根据观察到的数据估计成功概率）。

由于您正在测试不同的原假设，因此可能会得到不同的结果。 （顺便说一句，您可能对含义提出的论点也不同，因为排列测试遵循“随机化”推理模型，而蒙特卡洛测试遵循“总体”推理模型 - 请参阅

“为什么排列测试是这样的”的 3.3优于生物医学研究中的 t 和 F 检验”

。） 假设代码中没有错误，为什么两者在您的情况下会产生不同的 p 值（尽管看起来相似）更多的是一个统计问题。考虑描述数学并将该部分发布在

交叉验证

。