I am trying to compute welch's t-test statistic and p-value manually for educational purpose.

I have written a function for welch's t-test as follows,

def welch_t_test(vec_a, vec_b):
    mean_a = vec_a.mean()
    mean_b = vec_b.mean()
    var_a = vec_a.var()
    var_b = vec_b.var()

    T = (mean_a - mean_b)/((var_a-vec_a.size) + (var_b-vec_b.size))
    return T

This function returns T-test statistic. To find p-value, I am using pt function from R-3.3.3.

My example inputs:

A = [6,3,0,7,7,0,5,7,7,6]
B = [10,18,14,18,12,13,15,18,16,18]

My test statisitc is:

Test statistic: -7.065217391304347

When I try to compute p-value using pt of R like pt(-7.06521,19) I get 5.038877e-07 whereas My Rs defualt t.test function result yields test statistic of -8.1321 and p-value of 1.95e07. I see the qualitatively both are similar but there exists a quantitative difference. Could someone point out why is that so? And also df in R shows some fractional number, in my case df of R is 17.987 which I suppose 19 following n-1 formula for degrees of freedom.

  • 1
    $\begingroup$ At first glance--In your first block of code, I think your denominator for the Welch t statistic is incorrect. I will give it a try. It's always a good idea to make sure software is giving answers that agree with known formulas. $\endgroup$ – BruceET Apr 19 at 8:07

Welch t statistic:

Using R as a calculator:

A = c(6,3,0,7,7,0,5,7,7,6)
B = c(10,18,14,18,12,13,15,18,16,18)
m = length(A);  n = length(B)
t = ( mean(A) - mean(B) ) / sqrt( var(A)/m + var(B)/n ); t
[1] -8.132062

From t.test in R, where Welch is the default 2-sample t test.

A = c(6,3,0,7,7,0,5,7,7,6)
B = c(10,18,14,18,12,13,15,18,16,18)
t.test(A, B)

        Welch Two Sample t-test

data:  A and B
t = -8.1321, df = 17.987, p-value = 1.95e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -13.086987  -7.713013
sample estimates:
mean of x mean of y 
      4.8      15.2 

So the two versions agree that $T = -8.1321.$

Welch degrees of freedom and P-values:

Traditionally, t distributions had integer degrees of freedom, but the approximate DF for the Welch t test gives non-integer DF and there is no problem extending the formula for the density function of t to accommodate non-integer DF.

The (moderately messy) formula for DF of the Welch test can be found in many statistics texts and in Wikipedia. This formula gives degrees of freedom $\nu$ with $$ \min(m-1, n-1) \le \nu \le m + n - 2.$$ Roughly speaking $\nu$ is near the upper limit when sample sizes and sample variances are nearly equal.

Two-sided alternative. For the default 2-sided t test, the P-value in R is obtained using the CDF function pt. The result agrees with the rounded value given in the printout. For closer comparison, the unrounded P-value can be obtained using $-notation. (I think the very slight difference is that R does not print the fractional degrees of freedom to more decimal places.)

2 * pt(-8.1321, 17.987)
[1] 1.949439e-07
[1] 1.949789e-07

Left-sided alternative. For a left-sided test, the P-values from pt and t.test are again in substantial agreement:

pt(-8.1321, 17.987)
[1] 9.747197e-08
t.test(A,B, alt="less")$p.value
[1] 9.748945e-08       # P-value captured from 't.test' procedure

t.test(A,B, alt="less")$parameter

pt(t.test(A,B, alt="less")$statistic, t.test(A,B, alt="less")$parameter)
9.748945e-08           # exact agreement with exact DF

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.