Ytukyg

I'm trying to reproduce how a t-test works using the example from page 211 in The Art Of Computer System Performance Analysis by Raj Jain.

The calculation is as follows:

# system A sample and statistics

a <- c(5.36, 16.57, 0.62, 1.41, 0.64, 7.26)

x_a <- mean(a)

s2_a <- var(a)

n_a <- length(a)



# system B sample and statistics

b <- c(19.12, 3.52, 3.38, 2.5, 3.6, 1.74)

x_b <- mean(b)

s2_b <- var(b)

n_b <- length(b)



# computation of t-test

s <- (s2_a/n_a + s2_b/n_b)^(1/2)

v <- ((s2_a/n_a + s2_b/n_b)^2)/((1/(n_a - 1))*((s2_a/n_a)^2) + (1/(n_b - 1))*((s2_b/n_b)^2)) - 2



# 90% confidence interval

(x_a - x_b) + c(1, -1) * qt(c(0.95), v) * s

The output of the last computation is (6.55, -7.22) which matches the result given in the book (see erratas here: https://www.cse.wustl.edu/~jain/books/ftp/errors_all.pdf).

However, the built-in t-test gives a different result:

> t.test(a, b, conf.level = 0.9)



    Welch Two Sample t-test



data:  a and b

t = -0.09015, df = 9.9434, p-value = 0.93

alternative hypothesis: true difference in means is not equal to 0

90 percent confidence interval:

 -7.038828  6.372161

sample estimates:

mean of x mean of y 

 5.310000  5.643333 

The built-in function gives a confidence interval of (-7.04, 6.37). I'm unable to reproduce this interval. What is the cause of the difference? Is my calculation wrong?

Update: The different result is due to a different value for $v$, indicating the degrees of freedom. As, Ben Bolker points out, R uses 9.94 whereas the manual calculation in the book uses $9.94 - 2 = 7.94$ from it. The book does not say why it subtracts 2. It only assumes that the observations are independent and unpaired but is silent about the variance.

The answer by Neeraj reproduces the calculation of the degrees of freedom that is given in the Wikipedia article to the Welch's t-test.

You are making mistake in calculating your degree of freedom.

Here is the code, that exactly reproduce the R t.test results.

a <- c(5.36, 16.57, 0.62, 1.41, 0.64, 7.26)

b <- c(19.12, 3.52, 3.38, 2.5, 3.6, 1.74)



v1 <- var(a)

v2 <- var(b)



n1 <- length(a)

n2 <- length(b)



se <- sqrt(v1/n1 + v2/n2)



nu <- se^4 / ((v1^2 /(n1^2*(n1 -1))) + (v2^2/(n2^2*(n2-1)))) #degree of freedom



#Confidence Interval

mean(a) - mean(b) + c(1, -1)* qt(.95, nu)*se



> 6.372161 -7.038828

It exactly matches with t.test results

t.test(a, b, conf.level = 0.9)