Python中加权Wilcoxon秩和检验实现的正确性验证问询

用户问题 #

我拥有一组有序数据，希望在两组（例如男性与女性）之间进行Wilcoxon秩和检验。但由于其他类别（如年龄）的样本与总体比例存在不平衡，因此我想为每个数据点赋予权重。我以R中的加权Wilcoxon测试实现为参考，编写了一段Python代码（假设已获取x组和y组的权重），代码如下：

import numpy as np
from scipy.special import erfc

def weighted_ranksum_test(x: np.ndarray, y: np.ndarray, wx: np.ndarray, wy: np.ndarray):
    U = 0
    for iy, weight_y in zip(y, wy):
        smaller = x < iy
        equal = x == iy

        sum_smaller = np.sum(wx[smaller] * weight_y)
        sum_equal = np.sum(wx[equal] * weight_y / 2)
        sum_tot = sum_smaller + sum_equal

        U += sum_tot

    nY = np.sum(wy)
    nX = np.sum(wx)

    mU = nY * nX / 2
    sigU = np.sqrt((nY * nX * (1 + nY + nX)) / 12)
    zU = (U - mU) / sigU

    pU = erfc(zU / np.sqrt(2)) / 2

    return pU

我已在等权重假设下将该代码与scipy.stats.ranksums的结果对比，二者p值一致，但不确定此实现是否适用于我的样本不平衡场景，特咨询该实现是否正确。

核心逻辑的合理性 #

你的代码在等权重场景下与scipy.stats.ranksums结果一致，说明基础的U统计量计算、均值（mU）的推导完全符合标准Wilcoxon秩和检验的逻辑。针对加权场景，你实现的U统计量计算方式——对每个y样本的权重，累加x中小于y的样本权重乘积，加上等于y的样本权重乘积的一半——这是加权Wilcoxon检验中U统计量的标准计算思路，和R中同类实现的核心逻辑对齐，这部分是正确的。

样本不平衡场景的适配性 #

你提到的“其他类别样本与总体比例不平衡”属于抽样权重的应用场景，此时加权的目的是让样本回归总体结构，你的代码在这种场景下的适配需要注意两个关键点：

• 权重定义的正确性：只要wx和wy是正确的抽样权重（即每个样本代表总体中对应数量的个体），代码中用np.sum(wx)和np.sum(wy)替代原始样本量来计算有效样本量的方式，完全符合加权检验的处理原则，这部分没问题。
• 方差计算的局限性：当前代码中的sigU公式是标准Wilcoxon检验的方差公式，没有考虑加权场景下的结（tie）修正。如果你的数据中有较多相同取值（结），在加权场景下必须对方差进行修正，否则会导致z值偏差，进而影响p值的准确性。标准加权Wilcoxon检验的方差修正逻辑需要纳入所有结的权重和，示例如下：

  # 补充结修正的方差计算逻辑
  combined_data = np.concatenate([x, y])
  combined_weights = np.concatenate([wx, wy])
  unique_values = np.unique(combined_data)
  
  tie_correction = 0
  for val in unique_values:
      weight_sum = np.sum(combined_weights[combined_data == val])
      tie_correction += weight_sum ** 3 - weight_sum
  
  total_weight = nX + nY
  sigU_corrected = np.sqrt(
      (nY * nX / 12) * 
      ( (total_weight + 1) - tie_correction / (total_weight * (total_weight - 1)) )
  )
  

  如果数据中结很少，当前的方差计算可以近似使用；但结较多时，必须加入修正才能保证结果可靠。

结论 #

你的代码在无结或结极少的加权抽样场景下是正确的，能够适配样本结构不平衡的需求；但如果数据存在较多结，需要补充方差的结修正逻辑，否则p值可能会出现偏差。

内容的提问来源于stack exchange，提问作者Pedro H. Forli