Существует ли функция python (scipy) для определения параметров, необходимых для получения целевой мощности?

В 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?

Ответы

Ответ 1

Мне удалось реплицировать функцию, используя приведенную ниже формулу для n и функцию выживания norm.isf от scipy.stats

enter image description here

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

Ответ 2

Некоторые базовые вычисления мощности теперь доступны в статистических моделях

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

In R stats

> 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

Ответ 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]))

Надеюсь, что это поможет.