difference in test statistic and p-value R and Python script

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.

• 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. – 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
 -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.949439e-07 t.test(A,B)$p.value
 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)
 9.747197e-08
t.test(A,B, alt="less")$p.value  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