Различные результаты с взвешенным t-тестом между R и Python - PullRequest
Новая
       51

Различные результаты с взвешенным t-тестом между R и Python

3 голосов
/ 25 февраля 2020

Я запускаю взвешенный t-тест в Python и вижу другие результаты. Похоже, что проблема заключается в расчете степени свободы. Хотелось бы понять, почему я вижу разные выходы.

Вот пример кода.

В 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)

Значение P равно .06 в R, .02 в Python.
В исходном коде 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

Числитель здесь выглядит так же, но знаменатель выглядит совсем иначе.

1 Ответ

1 голос
/ 25 февраля 2020

Проблема, с которой вы столкнулись, заключается в том, что 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
...