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!!!
与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…