skbio.stats.power.subsample_power

skbio.stats.power.subsample_power(test, samples, draw_mode='ind', alpha_pwr=0.05, ratio=None, max_counts=50, counts_interval=10, min_counts=None, num_iter=500, num_runs=10)[源代码]

子样本数据迭代计算功率

状态：从0.4.0开始实验。

参数:

  • test (function) -- 接受值数组列表(样本ID或数值)并返回p值或p值的一维数组的统计测试。

  • samples (array_like) -- samples 可以是列表列表或数组列表，其中数组中的每个子列表或行对应于一个采样组。

  • draw_mode ({"ind", "matched"}, optional) -- 当样本中的观测值与其他组中的观测值有相应的观测值时，应使用“匹配”样本。例如，在处理以下情况下的回归数据时，这可能很有用 \(x_{1}, x_{2}, ..., x_{n}\) 映射到 \(y_{1}, y_{2}, ..., y_{n}\) 。在“匹配”模式下，样本向量的长度必须相同。如果样本之间没有倒数关系，则应使用“ind”模式。

  • alpha_pwr (float, optional) -- 用于计算功率的临界值。

  • ratio (1-D array, optional) -- 应分配给每组的样本计数的分数。如果这是一维数组，则其长度必须与 samples 。如果未提供任何值 (ratio 为无)，则将为每个样本抽取相同数量的观测值。在……里面 matched 模式，则该值将设置为1。

  • max_counts (positive int, optional) -- 为计算效果大小而绘制的每组的最大采样数。

  • counts_interval (positive int, optional) -- 每个子采样计数之间的差值。

  • min_counts (positive int, optional) -- 对于最小的子样本应该抽取多少个样本。如果这是NONE，则 counts_interval 将会被使用。

  • num_iter (positive int, optional) -- 要为曲线上的每个点生成的p值的数量。

  • num_runs (positive int, optional) -- 计算每条曲线的次数。

返回:

  • power ( array )--每次计数时为每个子样本计算的功率。该数组具有 num_runs 行的长度，元素数与 sample_counts 和深度，等于由 test 。如果 test 返回浮点数，则返回的数组将是二维的，而不是三维的。

  • sample_counts ( array )--每次幂计算时抽取的样本数。

抛出:

  • ValueError -- 如果 mode 为“匹配”，则将出现错误，如果 samples 是不一样的长度。

  • ValueError -- 如果采样数少于最小计数，则会出现ValueError。

  • ValueError -- 如果 counts_interval 大于样本开始值和最大值之间的差值，则该函数将引发ValueError。

  • ValueError -- 中没有相同数量的组 samples 在.中 ratios 。

  • TypeError -- test 不返回浮点数或一维NumPy数组。

示例

比方说，我们想要研究特定细菌的存在与 Gardnerella vaginalis 绝经前或绝经后妇女发生尿路感染(UTI)的可能性。健康女性在绝经前或绝经后参加了这项研究，并进行了8周的跟踪调查。参与者在研究开始时提交粪便样本，然后追踪尿路感染的临床症状。确诊的尿路感染是研究的终点。

利用现有文献和16S测序，确定了一组与尿路感染相关的候选分类群，包括 G. vaginalis 。在100名有尿路感染的妇女(绝经前和绝经后各50例)中，尿路感染的存在或不存在 G. vaginalis 经定量聚合酶链式反应证实。

我们可以对可检测到的概率进行建模 G. vaginalis 在这些样本中使用二项式模型发现。( Note that this is a simulation. )

>>> import numpy as np
>>> np.random.seed(25)
>>> pre_rate = np.random.binomial(1, 0.85, size=(50,))
>>> pre_rate.sum()
45
>>> pos_rate = np.random.binomial(1, 0.40, size=(50,))
>>> pos_rate.sum()
21

让我们建立一个测试函数，这样我们就可以测试在两组之间找到频率差异的概率。我们将使用 scipy.stats.chisquare 以寻找不同群体之间的频率差异。

>>> from scipy.stats import chisquare
>>> test = lambda x: chisquare(np.array([x[i].sum() for i in
...     range(len(x))]))[1]

让我们确保我们的两个发行版是不同的。

>>> print(round(test([pre_rate, pos_rate]), 3))
0.003

由于有偶数个样本，而我们没有足够的信息来尝试控制数据，因此我们将使用 skbio.stats.power.subsample_power 对两组进行比较。如果我们有其他风险因素的元数据，比如生育史、体重指数、吸烟情况，我们可能会想要使用 skbio.stats.power.subsample_paired_power 。我们还将使用“ind” draw_mode ，因为这两组样本之间没有关联。

>>> from skbio.stats.power import subsample_power
>>> pwr_est, counts = subsample_power(test=test,
...                                   samples=[pre_rate, pos_rate],
...                                   num_iter=100,
...                                   num_runs=5,
...                                   counts_interval=5)
>>> counts
array([ 5, 10, 15, 20, 25, 30, 35, 40, 45])
>>> np.nanmean(pwr_est, axis=0) 
array([ 0.056,  0.074,  0.226,  0.46 ,  0.61 ,  0.806,  0.952,  1.   ,
        1.   ])
>>> counts[np.nanmean(pwr_est, axis=0) > 0.8].min()
30

因此，我们可以估计，我们将看到显著不同的存在 G. vaginalis 如果我们每组至少有30个样本，在尿路感染前后妇女的粪便中。

如果我们想测试第二个候选分类群的关系，根据现有的文献，第二个候选分类群在种群中较为罕见，但可能具有类似的效果，我们可能还会从尝试在第二个候选分类群所在的组中识别30个样本开始。

现在，假设我们想要测试一种次生代谢物，只有在存在 G vaginalis 看看它是否也与尿路感染相关。我们可以将代谢物的丰度模拟为正态分布。

>>> met_pos = (np.random.randn(pre_rate.sum() + pos_rate.sum()) * 2000 +
...     2500)
>>> met_pos[met_pos < 0] = 0
>>> met_neg = met_neg = (np.random.randn(100 - (pre_rate.sum() +
...     pos_rate.sum())) * 2000 + 500)
>>> met_neg[met_neg < 0] = 0

让我们用Kruskal-Wallis检验来比较这些种群。从物理上讲，化学物质的浓度不可能是负的，所以我们把下限设为0。这意味着我们不能再假设我们的分布是正态的。

>>> from scipy.stats import kruskal
>>> def metabolite_test(x):
...     return kruskal(x[0], x[1])[1]
>>> print(round(metabolite_test([met_pos, met_neg]), 3))
0.005

当我们对所有数据进行统计检验时，你可能会注意到，来自患有 G. vaginalis 而不是那些没有的。当我们在测试力量时，解释这种差异可能是有意义的。因此，我们将设置 ratio 参数，这让我们可以从女性身上提取两倍的样本 G. vaginalis 。

>>> pwr_est2, counts2 = subsample_power(test=metabolite_test,
...                                     samples=[met_pos, met_neg],
...                                     counts_interval=5,
...                                     num_iter=100,
...                                     num_runs=5,
...                                     ratio=[2, 1])
>>> counts2
array([  5.,  10.,  15.,  20.,  25.,  30.])
>>> np.nanmean(pwr_est2, axis=0)
array([ 0.14 ,  0.272,  0.426,  0.646,  0.824,  0.996])
>>> counts2[np.nanmean(pwr_est2, axis=0) > 0.8].min()
25.0

当我们考虑功率分析中所需的每组样本的数量时，我们需要查看比率。分析说，我们需要25个样本，在最小的一组中，在这种情况下，没有 G. vaginalis 和50个样本来自女性 G. vaginalis 在80%的功率下，我们的次生代谢物的丰度有显著的不同。