如何解决试图在 python 中进行 Kruskall Wallis 事后测试,但统计数据不同?
我正在努力解决这个问题。我是来自 SPSS 背景的 Python 新手。基本上,一旦您完成了 Kruskal Wallis 检验并返回低 p 值,正确的程序是进行事后 Dunn 检验。我一直在努力计算数学,但我找到了这篇文章 (https://journals.sagepub.com/doi/pdf/10.1177/1536867X1501500117),我认为它已经说明了一切。
除了找出 P 值外,Python 似乎没有 Dunn 测试,但我希望有一个与您可以在 SPSS 中获得的成对比较测试类似的输出。这包括使用 Bonferroni 的 z-stat/检验统计量、标准差、标准差误差、p 值和调整后的 p 值。
现在我正在努力使测试统计正确,以便我可以完成其余的工作。我的数据是多个组,我已将其拆分为多个数据框。例如,我的数据如下所示:
df1 |因素 1 |因素 2 | | -------- | -------- | | 3.45 | 8.95 | | 5.69 | 2.35 | row_total=31 df2 |因素 1 |因素 2 | | -------- | -------- | | 5.45 | 7.95 | | 4.69 | 5.35 | row_total=75 等等等等
所以基本上我正在尝试测试 df1["Factor1"] 和 df2["Factor1]。 我目前拥有的是:
def dunn_test(df1,df2,colname):
##Equation is z= yI/Oi
##Where yi is the mean rankings of the two groups
## oi is the standard deviation of yi
#Data Needed
x=df1[colname]
z=df2[colname]
grouped = pd.concat([x,z])
N =len(grouped)
#calculating the Mean Rank of the Two Groups
rank1= stats.rankdata(x)
rank2=stats.rankdata(z)
Wa = rank1.sum()/len(x)
Wb = rank2.sum()/len(z)
#yi
y= Wa-Wb
#standard deviation of yi
#tied Ranks
ranks= stats.rankdata(grouped)
tied=pd.DataFrame([Counter(ranks)]).T
tied= tied.reset_index()
tied = tied.rename(columns={"index":"ranks",0:'ties'})
count_ties = tied[tied.ties >=2].count()
#standard Deviaton formula
t= tied["ties"]
for tied in t:
e = t**3-t
e = [i for i in e if i != 0]
oi=((N*(N+1)/2) - sum(e)/12*(N-1))*(1/len(x) + 1/len(z))
zstat=y/oi
return zstat
它输出 0.0630。我遇到的问题是,当我通过 SPSS 运行相同的测试时,数字是 -51.422。我不确定我做对了,有正确的方程式或我打算做什么。
任何帮助将不胜感激。
解决方法
我不得不做类似的事情。下面的代码应该适合你。它执行 Kruskal-Wallis 检验和 Dunn 检验。 Dunn 检验中的 p 值使用 Bonferroni 校正。数据需要在单列中进行结构化,并包含一些分层指标。 post_hoc_result_dict 依次返回变量名称、z 分数、p 值和更正后的 p 值。下面的代码应该对你有用。嗯。
函数调用:
f1 = df1['Factor 1'].to_frame(name='value')
f1['factor'] = 'Factor 1'
f2 = df1['Factor 1'].to_frame(name='value')
f2['factor'] = 'Factor 2'
correct_format = pd.concat([f1,f2])
k,p,post_hoc_result_dict = kw_test(correct_format,'factor','value')
功能:
def p_rounder(p_value):
if p_value < .0001:
p_value = '<.0001'
else:
p_value = str((round(p_value,4)))
return p_value
def bon_correct(p_value,k):
corrected_p = p_value * ((k *(k-1))/2)
return p_value,corrected_p
def kw_dunn_post_hoc(df,strat,comp_list,var):
post_hoc_result_dict = {}
N = df['rank'].count()
n_groups = df[strat].nunique()
for comp in comp_list:
m1 = df.loc[df[strat] == comp[0]]['rank'].mean()
n1 = df.loc[df[strat] == comp[0]]['rank'].count()
m2 = df.loc[df[strat] == comp[1]]['rank'].mean()
n2 = df.loc[df[strat] == comp[1]]['rank'].count()
Z = (m1 - m2)/sqrt(((N*(N+1))/12)*((1/n1)+(1/n2)))
Z = round(Z,4)
p = stats.norm.sf(abs(Z))
p,corrected_p = bon_correct(p,n_groups)
p = p_rounder(p)
corrected_p = p_rounder(corrected_p)
comparison = f'{comp[0]} vs. {comp[1]}'
post_hoc_result_dict[comparison] = [var,Z,corrected_p]
return post_hoc_result_dict
def kw_test(df,stratifier,var):
import sys
from math import sqrt
result_list = []
strat_list = []
comparison_list = []
counter = 0
temp_df = df[[stratifier,var]].copy()
temp_df['rank'] = temp_df[var].rank(method='average')
for strat in df[stratifier].unique():
result = df.loc[df[stratifier] == strat][var].values
result_list.append(result)
strat_list.append(strat)
for st in strat_list:
for st2 in strat_list:
if st != st2 and [st2,st] not in comparison_list:
comparison_list.append([st,st2])
post_hoc_result_dict = kw_dunn_post_hoc(temp_df,comparison_list,var)
if len(result_list) == 2:
k,p = stats.kruskal(result_list[0],result_list[1])
if len(result_list) == 3:
k,result_list[1],result_list[2])
elif len(result_list) == 4:
k,result_list[2],result_list[3])
elif len(result_list) == 5:
k,result_list[3],result_list[4])
else:
print('Stratifying levels greater than 5. Please modify code to accomodate.')
sys.exit()
k = round(k,4)
p = p_rounder(p)
return k,post_hoc_result_dict
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。