python - power calculation by hand for the z-test (binary variable)

Question

I am trying to calculate the power of a recently published study, in which a z-test was conducted to proof that immunosuppresiv medication significantly increases the prelimenary case fatality rate, after being infected with sars-cov-2. However, zt_ind_solve_power from statsmodels.stats.power gives back different results and I dont understand why! Here is what I did! I calculated the expected standard error in the population, based on group sizes and standard deviations and then computed the area under the aternative curve (here: approx. distribution, based on drawn sample) right from the critical value (here: critical value based on reported p in the paper):

import numpy as np
import scipy.stats as stats
import scipy.integrate as integrate
from statsmodels.stats.power import zt_ind_solve_power
n1 = 29
n2 = 37

mean1 = 0.28 # mortaility rate in control
mean2 = 0.7 # mortaility rate in treatment

x1 = n1*mean1
x2 = n2*mean2

p = (x1+x2)/(n1+n2)

sd1 = n1*p*(1-p)
sd2 = n2*p*(1-p)

sdp = np.sqrt(((n1-1)*sd1**2+(n2-1)*sd2**2)/n1+n2-2)
d = (mean2-mean1)/sdp

mean_null = 0 # H0: mortality rate equal in both groups
mean_sample = mean2-mean1

se = np.sqrt(p*(1-p)*((1/n1)+(1/n2)))
x = np.linspace(-0.5,1,num=10000)
null_dist = stats.norm.pdf(x=x, loc=mean_null, scale=se)
sample_dist = stats.norm.pdf(x=x, loc=mean_sample, scale=se)
p_val = 0.0013
critical_one = stats.norm.sf(x=p_val, loc=mean_null, scale=se)# assuming one-tailed test!

for i in range(x.size):
    if x[i] > critical_one and x[i-1] <= critical_one:
        print(x[i]-critical_one)
        index_one = i
    if x[i] > critical_two and x[i-1] <= critical_two:
        index_two = i
        print(x[i]-critical_two)
power_one = integrate.simps(y=sample_dist[index_one:x.size+1],x=x[index_one:x.size+1])

I am trying to solve for power by hand as part of my bachelor's thesis (I am analyzing various researcher degrees of freedom). I know that otherwise, there is no need to do that, really... What made me wonder, is that the solve_power function gives back power=1.9% for p=0.0013 and power=6.6%, which appears to be way too low to be true, here. I got power=27% at p=1.3%.

zt_ind_solve_power(effect_size = d, nobs1 = n1, alpha=p_val, ratio=n2/n1, alternative='larger')

I cannot get my head around what I am doing wrong! I dont expect zt_ind_solve_power to make the mistake, so... Can you tell me what I am missing?

Thank you so much!!!

Answer 1

等待大神答复