尝试在 python 中进行 Kruskal Wallis 事后测试，但统计数据不同？

Question

我正在努力解决这个问题。我是来自 SPSS 背景的 python 新手。基本上，一旦您完成了 Kruskal Wallis 检验并返回低 p 值，正确的程序是进行事后 Dunn 检验。我一直在努力弄清楚数学，但我找到了这篇文章 (https://journals.sagepub.com/doi/pdf/10.1177/1536867X1501500117)，我认为它已经说明了一切。

除了计算 P 值之外，Python 似乎没有 Dunn 检验，但我希望得到与可以在 SPSS 中获得的成对比较检验类似的输出。这包括 z-stat/test 统计量、标准偏差、标准偏差误差、p 值和使用 Bonferroni 调整的 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。我不确定我做对了，有正确的方程式或我打算做什么。

任何帮助将不胜感激。

Answer 1

我不得不做类似的事情。下面的代码应该适合你。它执行 Kruskal-Wallis 检验和 Dunn 检验。 Dunn 检验的 p 值使用 Bonferroni 校正。数据需要在单列中结构化，并包含一些分层指标。 post_hoc_result_dict 以该顺序返回变量名称、z 分数、p 值和更正的 p 值。下面的代码应该按原样为您工作。 lmk。

函数调用：

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,p,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,stratifier,comparison_list,var)
    if len(result_list) == 2:
        k,p = stats.kruskal(result_list[0],result_list[1])
    if len(result_list) == 3:
        k,p = stats.kruskal(result_list[0],result_list[1],result_list[2])
    elif len(result_list) == 4:
        k,p = stats.kruskal(result_list[0],result_list[1],result_list[2],result_list[3])
    elif len(result_list) == 5:
        k,p = stats.kruskal(result_list[0],result_list[1],result_list[2],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, p, post_hoc_result_dict