skbio.stats.power.subsample_paired_power

skbio.stats.power.subsample_paired_power(test, meta, cat, control_cats, order=None, strict_match=True, alpha_pwr=0.05, max_counts=50, counts_interval=10, min_counts=None, num_iter=500, num_runs=10)[源代码]

使用具有匹配元数据的样本迭代地估计幂

状态：从0.4.0开始实验。

参数:

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

  • meta (pandas.DataFrame) -- 与样本关联的元数据。

  • cat (str) -- 不同样本的元数据类别是不同的。

  • control_cats (list) -- 要用作控件的元数据类别。例如，如果您想改变年龄 (cat =“年龄”)，您可能希望控制性别和健康状况(即 control_cats = ["SEX", "HEALTHY"] )。

  • order (list, optional) -- 类别中组的顺序。这可用于限制选定的组。例如，如果有一个类别包含A组、B组和C组，而您只想查看A和B， order 将设置为 ['A', 'B'] 。

  • strict_match (bool, optional) -- 这决定了如何使用将数据分组 control_cats 。如果其中的样本 meta 中的任何列都具有未定义的值(NaN control_cats ，则该样本将不被视为具有匹配项，并且在以下情况下将被忽略 strict_match 是真的。如果 strict_match 为FALSE，则 control_cats 可以被视为匹配。

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

  • 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 -- 如果采样数少于最小计数，则会出现ValueError。

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

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

示例

假设你对髓系细胞中蛋白质易位的特定细胞因子的作用感兴趣。你可以培养两种巨噬细胞系(骨髓来源的吞噬细胞和腹膜来源的巨噬细胞)。由于情况不妙，您的成长媒体必须从多个来源(实验室、A公司、B公司)采购。同样不幸的是，您必须使用劳动密集型低通量分析。您有一些初步的测量结果，并且您想要预测需要分析多少(更多)细胞才能达到80%的功率。

您有关于60个细胞的信息，我们将在下面模拟这些信息。请注意，我们设置的是随机种子值以保持一致性。

>>> import numpy as np
>>> import pandas as pd
>>> np.random.seed(25)
>>> data = pd.DataFrame.from_dict({
...     'CELL_LINE': np.random.binomial(1, 0.5, size=(60,)),
...     'SOURCE': np.random.binomial(2, 0.33, size=(60,)),
...     'TREATMENT': np.hstack((np.zeros((30)), np.ones((30)))),
...     'INCUBATOR': np.random.binomial(1, 0.2, size=(60,))})
>>> data['OUTCOME'] = (0.25 + data.TREATMENT * 0.25) + \
...     np.random.randn(60) * (0.1 + data.SOURCE/10 + data.CELL_LINE/5)
>>> data.loc[data.OUTCOME < 0, 'OUTCOME'] = 0
>>> data.loc[data.OUTCOME > 1, 'OUTCOME'] = 1

我们将通过假设结果的分布不是正态分布来处理这一问题，并应用Kruskal-Wallis检验来比较细胞因子处理和未处理的细胞。

>>> from scipy.stats import kruskal
>>> f = lambda x: kruskal(*[data.loc[i, 'OUTCOME'] for i in x])[1]

让我们来看看细胞因子治疗是否对所有细胞都有显著影响。

>>> treatment_stat = [g for g in data.groupby('TREATMENT').groups.values()]
>>> round(f(treatment_stat), 17)
0.00193863362662502

现在，让我们选择控制类别。似乎有理由认为细胞系对治疗结果可能有影响，这可能归因于受体表达的差异。也有可能是由于细胞因子来源的不同。在整个实验过程中，培养箱保持在相同的条件下，在任何给定的时间，在温差1度内，在相同的二氧化碳水平内。因此，至少在一开始，让我们忽略由于孵化器而产生的差异。

尽管进一步迭代可能鼓励考虑效应大小与目标变量相似或大于目标变量的变量，但作为第一遍分析，建议基于什么可能与系统生物学相关的想法来选择控制变量。

>>> control_cats = ['SOURCE', 'CELL_LINE']
>>> from skbio.stats.power import subsample_paired_power
>>> pwr, cnt = subsample_paired_power(test=f,
...                                   meta=data,
...                                   cat='TREATMENT',
...                                   control_cats=control_cats,
...                                   counts_interval=5,
...                                   num_iter=25,
...                                   num_runs=5)
>>> cnt
array([  5.,  10.,  15.,  20.])
>>> pwr.mean(0)
array([ 0.24 ,  0.528,  0.68 ,  0.88 ])
>>> pwr.std(0).round(3)
array([ 0.088,  0.127,  0.168,  0.08 ])

从功率曲线估计，看起来每组20个细胞可能为这个实验提供足够的功率，尽管功率的巨大差异可能会建议延长曲线或增加每组的样本数量。