BUG: Inconsistent behaviour in power calculation

[ryanward-io]

Describe the bug

I don't understand why I'm getting different results between statsmodels and manually doing it in scipy. I presume I'm doing something wrong but would love if someone would help me understand.

Code Sample, a copy-pastable example if possible

Manually calculating power with scipy

import scipy.stats as st

mu1 = 50 # H0. Two tailed.
mu2 = 55 # Proposed critical parameter for beta calculation
sd = 20 # Assume same for both
n = 100
se = sd/(n)**0.5

# Upper value at alpha=0.05
critical_value = mu1 + st.norm.ppf(0.95) * se

# Calculating z-score for critical parameter
z = (critical_value - mu2)/se
beta = st.norm.cdf(z)
power = 1 - beta
print(power)
# 0.7054139024424575

Calculating power with statsmodels

import statsmodels.stats.power as smp

effect_size = (mu2-mu1)/sd

smp_power = smp.ttest_power(effect_size, nobs=n, alpha=0.05)
print(smp_power)
# 0.696980269099517

There is only a 0.4 difference but I don't understand why there is any difference at all.

Thanks in advance for your help!

Output of import statsmodels.api as sm; sm.show_versions()

INSTALLED VERSIONS ------------------ Python: 3.10.5.final.0 OS: Linux 5.10.16.3-microsoft-standard-WSL2 #1 SMP Fri Apr 2 22:23:49 UTC 2021 x86_64 byteorder: little LC_ALL: None LANG: C.UTF-8

statsmodels

Installed: 0.13.5 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/statsmodels)

Required Dependencies

cython: Not installed
numpy: 1.23.1 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/numpy)
scipy: 1.8.1 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/scipy)
pandas: 1.4.3 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/pandas)
dateutil: 2.8.2 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/dateutil)
patsy: 0.5.2 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/patsy)

Optional Dependencies

matplotlib: 3.5.2 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/matplotlib)
backend: module://matplotlib_inline.backend_inline
cvxopt: Not installed
joblib: 1.1.0 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/joblib)

Developer Tools

IPython: 8.4.0 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/IPython)
jinja2: 3.0.3 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/jinja2)
sphinx: Not installed
pygments: 2.11.2 (/home/rward/downloads/ENTER/envs/stats/lib/python3.10/site-packages/pygments)
pytest: Not installed
virtualenv: Not installed

[josef-pkt]

I don't manage to make the difference go away completely, but the difference becomes smaller if

• in your calculation you use the z-test not the t-test
• one-sided: you computation checks one tail at 0.05, but your result is for tail prob 0.025
  your power only looks at one tail, but as far as I can see the opposite tail prob is very small

with one sided alternative and t-test, I get

at a 0.025 tail

smp_power = smp.ttest_power(effect_size, nobs=n, alpha=0.025, alternative="larger")
print(smp_power)
0.6969756881165668

ddof = 0
se = sd/(n)**0.5
# Upper value at alpha=0.05
critical_value = mu1 + st.t.ppf(0.975, n-ddof) * se
​
# Calculating z-score for critical parameter
z = (critical_value - mu2)/se
beta = st.t.cdf(z, n - ddof)
power = 1 - beta
print(power)
0.6965131749922717

and at a 0.05 tail

smp_power = smp.ttest_power(effect_size, nobs=n, alpha=0.05, alternative="larger")
print(smp_power)
0.7989854657792202

ddof = 0
se = sd/(n)**0.5
# Upper value at alpha=0.05
critical_value = mu1 + st.t.ppf(0.95, n-ddof) * se
​
# Calculating z-score for critical parameter
z = (critical_value - mu2)/se
beta = st.t.cdf(z, n - ddof)
power = 1 - beta
print(power)

0.798478502107568

I currently don't see where the remaining difference comes from.
I thought there should be a degrees of freedom correction, but that didn't work.

[josef-pkt]

using the ztest power replicates your results
(so the problem with t-test is still likely some degrees of freedom handling.)

smp.zt_ind_solve_power(effect_size, nobs1=n, alpha=0.05, alternative="larger", ratio=0)
0.8037649400154937

smp.zt_ind_solve_power(effect_size, nobs1=n, alpha=0.025, alternative="larger", ratio=0)  # 0.7054139024424575
0.705413902442457

smp.zt_ind_solve_power(effect_size, nobs1=n, alpha=0.05, alternative="two-sided", ratio=0)
0.7054180011138002

[ryanward-io]

Thanks for investigating - really helpful. I think this can be closed then