[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.


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

Created on Sat Jan 12 21:48:06 2013

Author: Josef Perktold

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

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
        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,
#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


More information about the SciPy-User mailing list