7
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ Haven't dug into the calculations, but R gives df=9.94 whereas your manual calculation gives df=7.94. deparse(body(stats:::t.test.default))[73:77] will show you R's internal df calculations ... $\endgroup$
    – Ben Bolker
    Commented Nov 19, 2018 at 13:10

1 Answer 1

10
$\begingroup$

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)
$\endgroup$
5
  • $\begingroup$ at first glance, Wikipedia's Welch test df formula appears to be different. I wonder if there are different variants ... ??? $\endgroup$
    – Ben Bolker
    Commented Nov 19, 2018 at 13:16
  • $\begingroup$ I used wikipedia formula, and it exactly replicate the R result. $\endgroup$
    – Neeraj
    Commented Nov 19, 2018 at 13:17
  • $\begingroup$ This is Welch's t-test, which does not assume variance are equal for both the group. $\endgroup$
    – Neeraj
    Commented Nov 19, 2018 at 13:18
  • $\begingroup$ I knew it was Welch, but I may have been looking too quickly at the formula. $\endgroup$
    – Ben Bolker
    Commented Nov 19, 2018 at 13:22
  • $\begingroup$ @Ben There are indeed variants of Welch-Satterthwaite. $\endgroup$
    – Glen_b
    Commented Nov 20, 2018 at 1:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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