基于排列的 scipy.stats.ttest_1samp 替代方案

Question

我想使用基于排列的替代方案来测试我的观察平均值是否显着大于零。 我偶然发现了

scipy.stats.ttest_1samp

，但我不确定这是否也可以用于我的情况？我还偶然发现了

scipy.stats.permutation_test

，它似乎可以满足我的要求，但如果可以的话，我想坚持使用 mne.stats.permutation_t_test

。

示例：scipy

可以使用 

import numpy as np
from scipy import stats

# create data
np.random.seed(42)
rvs = np.random.normal(loc=5,scale=5,size=100)

# compute one-sample t-test 
t,p = stats.ttest_1samp(rvs,popmean=0,alternative='greater')

Answer 1

permutation_test

，它会“排列”观察结果的符号。假设数据已按上述方式生成，测试可以按如下方式进行

permutation_type='samples'如果您只关心 p 值，则使用

from scipy import stats
def t_statistic(x, axis=-1):
    return stats.ttest_1samp(x, popmean=0, axis=axis).statistic

res = stats.permutation_test((rvs,), t_statistic, permutation_type='samples')
print(res.pvalue)

 作为统计量而不是 

np.mean

可以获得相同的结果。

不可否认，只有一个样本的

t_statistic

的这种行为有点隐藏在文档中。

因此，如果数据仅包含一个样本，则通过独立更改每个观测值的符号来形成零分布。

但是产生相同 p 值的测试也可以作为双样本测试来执行，其中第二个样本是数据的负值。为了避免特殊情况，这实际上就是

permutation_type='samples'

在幕后所做的事情。

在这种情况下，上面的示例代码比现在的

permutation_test

快很多。不过，我会尝试在 SciPy 1.10 中改进这一点。

根据当前的

docs

Answer 2

permutation_test

函数似乎无法实现相当于单样本 t 检验的效果。但可以使用 permutation_test 来实现它，如下所示。这是基于 R 实现（在

here

找到）和 Cross Validated 上的

this

线程，并提供了进行单方面测试和针对添加的特定均值进行测试的选项。 numpy 示例 1（针对均值 0 的原假设的双边检验）：

import numpy as np

def permutation_ttest_1samp(
    data, popmean, n_resamples, alternative='two-sided', random_state=None
):

    assert alternative in ('two-sided', 'less', 'greater'), (
        "Unrecognized alternative hypothesis"
    )

    n = len(data)

    data = np.asarray(data) - popmean
    dbar = np.mean(data)
    
    absx = np.abs(data)
    z = []

    rng = np.random.RandomState(random_state)

    for _ in range(n_resamples):
        mn = rng.choice((-1,1), n, replace=True)
        xbardash = np.mean(mn * absx)
        z.append(xbardash)
    z = np.array(z)

    if alternative == 'greater':
        return 1 - (np.sum(z <= -np.abs(dbar)) / n_resamples)
    elif alternative == 'less':
        return np.sum(z <= -np.abs(dbar)) / n_resamples
    return (
        (np.sum(z >= np.abs(dbar)) + np.sum(z <= -np.abs(dbar))) / n_resamples
    )

与参数化 t 检验相比：

rng = np.random.RandomState(42)
rvs = rng.normal(loc=0, scale=0.01, size=1000)

pval = permutation_ttest_1samp(rvs, 0, 100_000, alternative='two-sided', random_state=42)
print(pval)
# 0.53206

示例 2（针对非 0 均值零假设的单边检验）

from scipy.stats import ttest_1samp

stat, pval = ttest_1samp(rvs, popmean=0, alternative='two-sided')
print(pval)
# 0.5325672436623021

与参数化 t 检验相比：

rng = np.random.RandomState(42)
rvs = rng.normal(loc=0, scale=3, size=1000)

pval = permutation_ttest_1samp(rvs, 0.1, 100_000, alternative='greater', random_state=42)
print(pval)
# 0.6731