Reproducir la prueba t en R da un resultado diferente al de la función incorporada

Question

Estoy intentando reproducir cómo funciona una prueba t utilizando el ejemplo de la página 211 de The Art Of Computer System Performance Analysis de Raj Jain.

El cálculo es el siguiente:

# 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

El resultado del último cálculo es (6,55, -7,22), que coincide con el resultado dado en el libro (ver erratas aquí: https://www.cse.wustl.edu/~jain/libros/ftp/errores_todos.pdf ).

Sin embargo, la prueba t incorporada da un resultado diferente:

> 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 

La función incorporada da un intervalo de confianza de (-7,04, 6,37). No consigo reproducir este intervalo. ¿Cuál es la causa de la diferencia? ¿Mi cálculo es erróneo?

Actualización : El resultado diferente se debe a un valor distinto de $v$ que indica los grados de libertad. Como señala Ben Bolker, R utiliza 9,94 mientras que el cálculo manual del libro utiliza $9.94 - 2 = 7.94$ de él. El libro no dice por qué resta 2. Sólo asume que las observaciones son independientes y no apareadas, pero no dice nada sobre la varianza.

La respuesta de Neeraj reproduce el cálculo de los grados de libertad que se da en el artículo de Wikipedia al Prueba t de Welch .

Answer 1

Estás cometiendo un error al calcular tu grado de libertad.

Aquí está el código, que reproduce exactamente los resultados de R t.test.

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

Coincide exactamente con los resultados de la prueba t

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