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.
Anyone can help fix this? Thank you!
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.)