jerome-white/sakai-book

Statistical methods from Tetsuya Sakai's IR evaluation book

sakai-book's Issues

Inaccurate power calculations for this F-test

There are discrepancies between the values statsmodels F-test power calculations produce
and those from R.

Consider the following R program (specifically taken from 2waynorep):

library(pwr)
fstat <- 0.63
m <- 4
n <- 17
alpha <- 0.05
pwr.f2.test(u=m-1, v=(m-1)*(n-1), f2=fstat, sig.level=alpha)$power

Produces 0.9984019. However, the "semantically" equivalent Python program:

from statsmodels.stats.power import FTestPower

(fstat, m, n, alpha) = (0.63, 4, 17, 0.05)
kwargs = {
    'effect_size': fstat / (n - 1),
    'df_num': m - 1,
    'df_denom': (m - 1) * (n - 1),
    'alpha': 0.05,
}

fpow = FTestPower()
print(fpow.solve_power(**kwargs))

Produces 0.05011698892318561.

According to blog post by Josef Perktold (section titled, "F-Test for linear Model"), the correct way to use solve_power is to interpret the arguments differently. When applying the logic laid out in the corresponding section of that post:

import math
from statsmodels.stats.power import FTestPower

(fstat, m, n, alpha) = (0.63, 4, 17, 0.05)
kwargs = {
    'effect_size': math.sqrt(fstat / (n - fstat)),
    'df_denom': m - 1,
    'df_num': (m - 1) * (n - 1),
    'alpha': 0.05,
}

fpow = FTestPower()
print(fpow.solve_power(**kwargs))

Produces 0.18002339044585902, which, although not exactly what is seen in R, is much closer.