R 和 Python 之间加权 t 检验的不同结果

Question

我正在 Python 中运行加权 t 测试，但看到了不同的结果。看来问题是自由度计算。想了解为什么我看到不同的输出。

这是一些示例代码。

在 R 中：

library(weights)
x <- c(373,398,245,272,238,241,134,410,158,125,198,252,577,272,208,260)
y <- c(411,471,320,364,311,390,163,424,228,144,246,371,680,384,279,303)

weightsa = c(rep(1,8), rep(2,8))
weightsb = c(rep(2,8), rep(1,8))

wtd.t.test(x = x,
           y = y,
           weight = weightsa,
           weighty = weightsb, samedata=F)

$test [1] "Two Sample Weighted T-Test (Welch)"

$coefficients
    t.value          df     p.value 
-1.88907197 29.93637837  0.06860382 

$additional Difference     Mean.x     Mean.y   Std. Err   -80.50000 
267.12500  347.62500   42.61352

在 Python 中：

import numpy as np
from statsmodels.stats.weightstats import ttest_ind
x = np.asarray([373,398,245,272,238,241,134,410,158,125,198,252,577,272,208,260])
y = np.asarray([411,471,320,364,311,390,163,424,228,144,246,371,680,384,279,303])

weightsa = [1] * 8 + [2] * 8
weightsb = [2] * 8 + [1] * 8

ttest_ind(x, y, usevar='unequal', weights=(weightsa, weightsb))

(-2.3391969704691085, 0.023733058922455107, 45.90244683439944)

R 中的 P 值为 .06，Python 中为 .02。
R 源代码使用 Satterthwaite 公式计算自由度：

df <- (((vx/n) + (vy/n2))^2)/((((vx/n)^2)/(n - 1)) + 
      ((vy/n2)^2/(n2 - 1)))

Python 函数源代码也声称使用此公式：

def dof_satt(self):
        '''degrees of freedom of Satterthwaite for unequal variance
        '''
        d1 = self.d1
        d2 = self.d2
        #this follows blindly the SPSS manual
        #except I use  ``_var`` which has ddof=0
        sem1 = d1._var / (d1.nobs-1)
        sem2 = d2._var / (d2.nobs-1)
        semsum = sem1 + sem2
        z1 = (sem1 / semsum)**2 / (d1.nobs - 1)
        z2 = (sem2 / semsum)**2 / (d2.nobs - 1)
        dof = 1. / (z1 + z2)
        return dof

这里的分子看起来一样，但分母却大不相同。

Answer 1

你在这里遇到的问题是weights::wtd.t.test()有一个（对我来说，奇怪的）默认参数mean1 = TRUE，它控制“是否应该强制权重的平均值为 1”（来自help("wtd.t.test")）。

如果我们使用mean1 = FALSE，我们会得到与以下相同的行为ttest_ind()：

wtd.t.test(x = x,
           y = y,
           weight = weightsa,
           weighty = weightsb,
           samedata = FALSE,
           mean1 = FALSE)

$test
[1] "Two Sample Weighted T-Test (Welch)"

$coefficients
    t.value          df     p.value 
-2.33919697 45.90244683  0.02373306 

$additional
Difference     Mean.x     Mean.y   Std. Err 
 -80.50000  267.12500  347.62500   34.41352