Расчет доверительных интервалов в R

Question

Я пытаюсь вычислить доверительные интервалы из t-теста в R вручную, и я подозреваю, что способ их вычисления отключен.

Вот как я сейчас рассчитываю доверительные интервалы вручную

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 
tvalue <-a1$statistic

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])
#formula for standard error
stan_error=sqrt((sd1/n1)+(sd2/n2))

тогда я беру с этой страницы формулу о расчете доверительных интервалов http://onlinestatbook.com/2/estimation/difference_means.html

Нижний доверительный интервал, который я вычисляю вот так

lower=mean_diff - (tvalue * stan_error)'

и получается результат -4,147333

Но доверительные интервалы

t.test(mpg ~ am, mtcars)

являются

95 percent confidence interval:
 -11.280194  -3.209684

Есть идеи?

Answer 1

Во-первых, критическое значение для t неверно.

tvalue <- a1$statistic необходимо заменить на tvalue <- abs(qt(0.05/2, 30)).

Обратите внимание, что это не 32, потому что мы теряем 2 степени свободы.

И вам не хватает ^2 (то есть в степени двух) в формуле для стандартной ошибки.То, что у вас есть в sd1 и sd2, является стандартными ошибками, поэтому вам нужно преобразовать их в дисперсии.Итак, правильная формула:

stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

Таким образом, новый код становится:

library(broom)
data("mtcars")
a1=tidy(t.test(mpg ~ am, mtcars))
mean_diff<-a1$estimate 

t_cv<- abs(qt(0.05/2, 30))

#standard error 
sd1=sd(mtcars$mpg[mtcars$am==0])
sd2=sd(mtcars$mpg[mtcars$am==1])
n1=length(mtcars$mpg[mtcars$am==0])
n2=length(mtcars$mpg[mtcars$am==1])

#formula for standard error
stan_error = sqrt((sd1^2 / n1) + (sd2^2 / n2))

lower=mean_diff - (t_cv* stan_error)

lower
[1] -11.17264

Но это все равно не соответствует доверительному интервалу при использовании функции t.test, поскольку t.test использует Уэлчаt-тест (https://en.wikipedia.org/wiki/Welch%27s_t-test) Таким образом, ваше t-критическое значение по t-критерию Уэлча должно быть

# Welch's t test degrees of freedom
welch_df <- (sd1^2/n1 + sd2^2/n2)^2 / (sd1^4/(n1^2*(n1-1)) + sd2^4/(n2^2*(n2-1)))
t_cv <- abs(qt(0.05/2, welch_df))  

# Recalculate lower confidence interval
lower= mean_diff - (t_cv* stan_error)
lower
[1] -11.28019 # this matches confidence interval in t.test