Search code examples
pythonstatsmodels

Negative sample size by statsmodels.TTestIndPower().solve_power

I'm using solve_power method of TTestIndPower from statsmodels to solve required holdout size. Normally, we use it to solve treatment size. But here I need to solve the control size and the remaining of the population will be treatment. statsmodels seems don't work well in this case as it returns a negative sample size, while the required sample size should be 210.

from statsmodels.stats.power import TTestIndPower

effect_size = 0.2
alpha = 0.05
power = 0.8
total_size = 3617

ideal_holdout_size = math.ceil(total_size - 
                               TTestIndPower().solve_power(effect_size=effect_size, nobs1=None, alpha=alpha, power=power, 
                                                                             ratio=(total_size - nobs1) / nobs1, alternative='two-sided'))

print(f'Out of total {total_size} stores, the ideal holdout size is: {ideal_holdout_size}')

Here's the result of the above code.

enter image description here

Anyone can help fix this? Thank you!

Solution

  • Solving for nobs1 in this setting is not possible with solve_power because both nobs1 and nobs_ratio need to change in this case. solve_power can only search for the value of one keyword given the others.

    The function used for the rootfinding needs to be changed so that it solves for nobs1 for a given total size. There might not always exist a root (*) in this case or there could be multiple roots.

    In the following I check that power is monotonic over the relevant range of nobs1 and that it includes the desired power of 0.8 in this range of nobs1. Then we can use a scipy rootfinder to solve for nobs1.

    nobs1 = np.arange(100, 500, 10)
pow_ = TTestIndPower().power(effect_size=effect_size, nobs1=nobs1, alpha=alpha,
                             ratio=(total_size - nobs1) / nobs1, alternative='two-sided')
array([0.50469817, 0.54182567, 0.57681571, 0.60967519, 0.64043692,
       0.6691538 , 0.69589402, 0.72073695, 0.74376976, 0.76508462,
       0.78477648, 0.80294111, 0.81967375, 0.83506793, 0.84921453,
       0.86220121, 0.87411189, 0.88502644, 0.89502052, 0.90416541,
       0.91252805, 0.92017113, 0.92715306, 0.93352823, 0.93934712,
       0.94465644, 0.9494994 , 0.95391587, 0.95794252, 0.96161315,
       0.96495878, 0.96800784, 0.97078643, 0.97331846, 0.97562573,
       0.97772825, 0.97964426, 0.98139039, 0.98298186, 0.9844325 ])

from scipy import optimize
def power_func(nobs1):
    pow_ = TTestIndPower().power(effect_size=effect_size, nobs1=nobs1, alpha=alpha,
                                 ratio=(total_size - nobs1) / nobs1, alternative='two-sided')
    return pow_ - 0.8

optimize.brentq(power_func, 100, 500)
208.32446507391262

    (*) As a quick check whether the total sample size is large enough for the desired power, we can compute the required sample size at nobs_ratio equals to 1. Using equal sample sizes will have the largest power for fixed total size.

    In the example we need at least a total sample size of 787 to achieve the desired power of 0.8.

    nobs1 = TTestIndPower().solve_power(effect_size=effect_size, nobs1=None, alpha=alpha, power=power, 
                                   ratio=1, alternative='two-sided')
nobs1, 2 * nobs1
(393.4056989990322, 786.8113979980644)

    (Aside: I have never seen this usecase in the literature nor in other packages, so statsmodels sample size solvers were not designed for this case.)