Mrthdrb

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.

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
t.test(A,B)$p.value
[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
 df 
17.98672 

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