有没有python（scipy）函数可以确定获得目标功率所需的参数？

Question

在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)

它会告诉你为了达到这个功效所需的样本大小n。这在确定测试的样本大小时非常有用。

在scipy包中有没有类似的函数呢？

Answer 1

Matt的回答关于每组需要多少个样本几乎是对的，但有一个小错误。

给定d（均值差）、s（标准差）、sig（显著性水平，通常是0.05）和power（通常是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]))

希望这对你有帮助。

Answer 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 

Answer 3

我已经成功用下面的公式和来自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