是否有 python (scipy) 函数来确定获得目标功率所需的参数？答案

【问题标题】：Is there a python (scipy) function to determine parameters needed to obtain a target power?是否有 python (scipy) 函数来确定获得目标功率所需的参数？
【发布时间】：2013-03-04 14:39:47
【问题描述】：

在 R 中，有一个非常有用的函数可以帮助确定双边 t 检验的参数以获得目标统计功效。

函数调用power.prop.test。

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/power.prop.test.html

您可以使用以下方式调用它：

power.prop.test(p1 = .50, p2 = .75, power = .90)

它会告诉你获得这种功效所需的样本量。这对于阻止测试的样本量非常有用。

scipy包里有类似的功能吗？

【问题讨论】：

如果有的话，我想应该是here。
该函数也是用纯R 编写的，所以在没有() 的情况下调用它会显示源代码。如果 numpy 的端口尚不存在，它将是直截了当的。
感谢@Justin 这有助于创建以下内容。
感谢@Raufio 我使用您链接到的页面找到了下面的 isf 函数。

标签： python r numpy scipy

【解决方案1】：

我已经成功地使用下面的 n 公式和来自 scipy.stats 的逆生存函数 norm.isf 复制了该函数

from scipy.stats import norm, zscore

def sample_power_probtest(p1, p2, power=0.8, sig=0.05):
    z = norm.isf([sig/2]) #two-sided t test
    zp = -1 * norm.isf([power]) 
    d = (p1-p2)
    s =2*((p1+p2) /2)*(1-((p1+p2) /2))
    n = s * ((zp + z)**2) / (d**2)
    return int(round(n[0]))

def sample_power_difftest(d, s, power=0.8, sig=0.05):
    z = norm.isf([sig/2])
    zp = -1 * norm.isf([power])
    n = s * ((zp + z)**2) / (d**2)
    return int(round(n[0]))

if __name__ == '__main__':

    n = sample_power_probtest(0.1, 0.11, power=0.8, sig=0.05)
    print n  #14752

    n = sample_power_difftest(0.1, 0.5, power=0.8, sig=0.05)
    print n  #392

【讨论】：

你考虑过把它捐给 SciPy 吗？这肯定是一个有用的功能。
您需要在 GitHub 上注册，然后 fork their repo，将您的更改放入并提交拉取请求。（不幸的是，目前 SciPy 开发者文档有点乱……）
感谢@larsmans 我在 github 上，所以我会分叉并做到这一点。干杯
这看起来很有希望。您有机会在下面回答@erikwestlund 的答案吗？
我假设：（d = 均值差）和（s = 标准差）。但是 p1 和 p2 是什么？

【解决方案2】：

一些基本的功率计算现在可以在 statsmodels 中使用

http://statsmodels.sourceforge.net/devel/stats.html#power-and-sample-size-calculations http://jpktd.blogspot.ca/2013/03/statistical-power-in-statsmodels.html

博客文章尚未考虑对 statsmodels 代码的最新更改。另外，我还没有决定提供多少包装函数，因为许多功率计算只是减少到基本分布。

>>> import statsmodels.stats.api as sms
>>> es = sms.proportion_effectsize(0.5, 0.75)
>>> sms.NormalIndPower().solve_power(es, power=0.9, alpha=0.05, ratio=1)
76.652940372066908

在 R 统计中

> power.prop.test(p1 = .50, p2 = .75, power = .90)

     Two-sample comparison of proportions power calculation 

              n = 76.7069301141077
             p1 = 0.5
             p2 = 0.75
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

 NOTE: n is number in *each* group

使用 R 的 pwr 包

> library(pwr)
> h<-ES.h(0.5,0.75)
> pwr.2p.test(h=h, power=0.9, sig.level=0.05)

     Difference of proportion power calculation for binomial distribution (arcsine transformation) 

              h = 0.5235987755982985
              n = 76.6529406106181
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

 NOTE: same sample sizes

【讨论】：

这有问题，产生的答案取决于您使用的是 R 还是 Python，尤其是当您改变比率时。有什么想法有什么问题吗？
这是 R stats 和 Stata 与 R pwr 和 statsmodels。有关详细信息，请参阅github.com/statsmodels/statsmodels/issues/1197 和相关的邮件列表线程。我不记得 SAS 在哪里。

【解决方案3】：

Matt 获得所需 n（每组）的答案几乎是正确的，但有一个小错误。

给定 d（均值差）、s（标准差）、sig（显着性水平，通常为 0.05）和功效（通常为 0.80），计算每组观察数的公式为：

n= (2s^2 * ((z_(sig/2) + z_power)^2) / (d^2)

正如你在他的公式中看到的，他有

n = s * ((zp + z)**2) / (d**2)

“s”部分是错误的。重现 r 功能的正确函数是：

def sample_power_difftest(d, s, power=0.8, sig=0.05):
    z = norm.isf([sig/2]) 
    zp = -1 * norm.isf([power])
    n = (2*(s**2)) * ((zp + z)**2) / (d**2)
    return int(round(n[0]))

希望这会有所帮助。

【讨论】：

【解决方案4】：

你还有：

from statsmodels.stats.power import tt_ind_solve_power

并在您想要获得的值中输入“None”。例如，要获得 effect_size = 0.1、power = 0.8 等情况下的观察次数，您应该输入：

tt_ind_solve_power(effect_size=0.1, nobs1 = None, alpha=0.05, power=0.8, ratio=1, alternative='two-sided')

并获得：1570.7330663315456 作为所需的观察次数。或者，在固定其他值的情况下获得你可以获得的力量：

tt_ind_solve_power(effect_size= 0.2, nobs1 = 200, alpha=0.05, power=None, ratio=1, alternative='two-sided')

你得到：0.5140816347005553

【讨论】：