# [SciPy-User] Power

josef.pktd@gmai... josef.pktd@gmai...
Sat Jan 12 21:38:04 CST 2013

```Oh, it's just power of statistical hypothesis tests

experiment: how to solve for root with respect to each argument of a function

first steps in power and sample size calculations.
coming (eventually) to statsmodels.stats

Dumping a script, in case anyone is interested.

Josef

-----------------
# -*- coding: utf-8 -*-
"""Statistical power, solving for nobs, ... - trial version

Created on Sat Jan 12 21:48:06 2013

Author: Josef Perktold

Example
roundtrip - root with respect to all variables

calculated, desired
nobs   33.367204205 33.367204205
effect 0.5 0.5
alpha  0.05 0.05
beta   0.8 0.8

"""

import numpy as np
from scipy import stats

#power
d = 0.68  #effect size
nobs = 20  #observations
alpha = 0.05
alpha_two = alpha / 2.
crit = stats.t.isf(alpha_two, nobs-1) #critical value at df=nobs
pow_t = stats.nct(nobs-1, d*np.sqrt(nobs)).sf(stats.t.isf(alpha_two, nobs-1))
pow_f = stats.ncf(1, nobs-1, d**2*nobs).sf(stats.f.isf(alpha, 1, nobs-1))

def ttest_power(effect_size, nobs, alpha, df=None, alternative='two-sided'):
'''Calculate power of a ttest
'''
d = effect_size
if df is None:
df = nobs - 1
if alternative == 'two-sided':
alpha_ = alpha / 2.  #no inplace changes, doesn't work
else:
alpha_ = alpha
pow_ = stats.nct(df, d*np.sqrt(nobs)).sf(stats.t.isf(alpha_, df))
return pow_

from scipy import optimize

def solve_power_nobs(effect_size, alpha, beta):
'''solve for required sample size for t-test to obtain power beta
'''

pow_n = lambda nobs: ttest_power(effect_size, nobs, alpha) - beta
return optimize.fsolve(pow_n, 10).item()

effect_size, alpha, beta = 0.5, 0.05, 0.8
nobs_p = solve_power_nobs(effect_size, alpha, beta)
print nobs_p
print ttest_power(effect_size, nobs_p, alpha), beta

#generic implementation
#----------------------

#to be reused for other tests

def ttest_power_id(effect_size=None, nobs=None, alpha=None, beta=None,
df=None, alternative='two-sided'):
'''identity for power calculation, should be zero

'''
return ttest_power(effect_size=effect_size, nobs=nobs, alpha=alpha,
df=df, alternative=alternative) - beta

#import functools as ft   #too limited

start_ttp = dict(effect_size=0.01, nobs=10., alpha=0.15,
beta=0.6)
#possible rootfinding problem for effect_size, starting small seems to work

def solve_ttest_power(**kwds):
'''solve for any one of the parameters of a t-test

for t-test the keywords are:
effect_size, nobs, alpha, beta

exactly one needs to be `None`, all others need numeric values

'''
#TODO: maybe use explicit kwds,
#      nicer but requires inspect? and not generic across tests
#TODO: use explicit calculation for beta=None
key = [k for k,v in kwds.iteritems() if v is None]
#print kwds, key;
if len(key) != 1:
raise ValueError('need exactly one keyword that is None')
key = key[0]

def func(x):
kwds[key] = x
return ttest_power_id(**kwds)

return optimize.fsolve(func, start_ttp[key]).item() #scalar

print '\nroundtrip - root with respect to all variables'
print '\n       calculated, desired'

print 'nobs  ', solve_ttest_power(effect_size=effect_size, nobs=None,
alpha=alpha, beta=beta), nobs_p
print 'effect', solve_ttest_power(effect_size=None, nobs=nobs_p,
alpha=alpha, beta=beta), effect_size

print 'alpha ', solve_ttest_power(effect_size=effect_size,
nobs=nobs_p, alpha=None, beta=beta), alpha
print 'beta  ', solve_ttest_power(effect_size=effect_size,
nobs=nobs_p, alpha=alpha, beta=None), beta

--------------
```