Обратные признаки тестовой статистики с использованием lm() по сравнению с t.test()

ProgramBox

Обратные признаки тестовой статистики с использованием lm() по сравнению с t.test()

Post author:admin
Запись опубликована:15 января, 2022
Post category:Вопросы по программированию

#r #linear-regression

Вопрос:

Почему значения тестовой статистики имеют обратные знаки при использовании lm() сравнения t.test() для одной 2-факторной переменной для линейного мореля? Например, ниже:

 sim_data = data.frame(
  groups = c(rep("A", n / 2), rep("B", n / 2)),
  values = rep(0, n))

sim_data$values = rnorm(n, mean = 42, sd = 3.5) 
summary(lm(values ~ groups, data = sim_data))
t.test(values ~ groups, data = sim_data, var.equal = TRUE)

Оба t.test() и lm() дают одну и ту же тестовую статистику, но с разными знаками: 0,02 и -0,02. Это потому, что, например, здесь R использует группу A в качестве опорного уровня, lm() но использует группу B в качестве опорного уровня t.test() ? Пожалуйста, просвети меня.

1. Да, просто эти две функции используют разную параметризацию. lm Принимает A в качестве элемента управления. t.test смотрит на разницу «между группой А и группой В» — то есть В является контролем. Статистика теста идентична, как и выводы, которые вы можете сделать.

Ответ №1:

Поскольку в справке по t.test это не обсуждается, давайте выпишем различные возможности, чтобы мы могли эмпирически определить, как он действует (или посмотрите ссылки на исходный код в комментариях).

Короткий ответ заключается в том, что мы можем заставить lm дать тот же результат, если мы используем contr.Контрасты SAS в lm (как в № 4 в разделе сравнения), которые совпадают с контрастами обработки по умолчанию, за исключением того, что последний уровень используется в качестве эталонного уровня вместо первого.

метод по умолчанию

Без метода формулы t.test можно записать как t.test(x, y, ...) , и статистика t, полученная в результате, представляет собой среднее значение x минус среднее значение y, деленное на стандартную ошибку, поэтому знак статистики определяется порядком, в котором x и y передаются в качестве аргументов.

 # the standard error does not depend on which level is reference level
# so just get it in most convenient way
out <- t.test(values ~ groups, sim_data, var.equal = TRUE)
out$statistic
       t 
1.118583 

stderr <- out$stderr
x <- with(sim_data, values[groups == "A"])
y <- with(sim_data, values[groups == "B"])

(mean(x) - mean(y)) / stderr
## [1] 1.118583

(mean(y) - mean(x)) / stderr
## [1] -1.118583

способ получения формулы

Если мы используем интерфейс формулы с символьным столбцом для групп, он эффективно преобразует группы в фактор с первой группой в алфавитном порядке в качестве первого аргумента и второй группой в алфавитном порядке в качестве второго аргумента.

 # A is passed as the first argument to t.test.default
gA <- factor(sim_data$groups, levels = c("A", "B"))
t.test(values ~ gA, sim_data, var.equal = TRUE)$statistic
##        t 
## 1.118583 

# B is passed as the first argument to t.test.default
gB <- factor(sim_data$groups, levels = c("B", "A"))
t.test(values ~ gB, sim_data, var.equal = TRUE)$statistic
##         t 
## -1.118583

лм

С помощью lm очевидно, что первый уровень эффективно передается в t.test.default в качестве второго аргумента и называется базовым или эталонным уровнем.

 # A is reference level
coef(summary(lm(values ~ gA, data = sim_data)))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 42.533345  0.8827383 48.183414 1.750263e-28
## gAB         -1.396417  1.2483805 -1.118583 2.728234e-01

# B is reference level
coef(summary(lm(values ~ gB, data = sim_data)))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 41.136929  0.8827383 46.601500 4.407622e-28
## gBA          1.396417  1.2483805  1.118583 2.728234e-01

или мы можем управлять базовым уровнем, также известным как базовый уровень, в lm, используя аргумент контрастов (или используя релевантность, как в другом ответе).:

 # A is reference level (default is first level)
coef(summary(lm(values ~ groups, data = sim_data)))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 42.533345  0.8827383 48.183414 1.750263e-28
## groupsB     -1.396417  1.2483805 -1.118583 2.728234e-01

# same - groups2 in output is groups with the first level, A, as base 
coef(summary(lm(values ~ groups, data = sim_data, 
  contrasts = list(groups = contr.treatment(2, base = 1)))))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 42.533345  0.8827383 48.183414 1.750263e-28
## groups2     -1.396417  1.2483805 -1.118583 2.728234e-01

# change reference level - groups1 in output is groups with 2nd level, B, as base
coef(summary(lm(values ~ groups, data = sim_data, 
  contrasts = list(groups = contr.treatment(2, base = 2)))))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 41.136929  0.8827383 46.601500 4.407622e-28
## groups1      1.396417  1.2483805  1.118583 2.728234e-01

# contr.SAS is like contr.treatment except last level is reference
coef(summary(lm(values ~ groups, data = sim_data, 
  contrasts = list(groups = contr.SAS))))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 41.136929  0.8827383 46.601500 4.407622e-28
## groups1      1.396417  1.2483805  1.118583 2.728234e-01

Сравнение

Возможно, будет легче сравнивать, если мы напишем похожие выражения вместе. Все они дают одну и ту же статистику t.

 attach(sim_data)   
out <- t.test(values ~ groups, var.equal = TRUE)
stderr <- out$stderr
x <- values[groups == "A"]
y <- values[groups == "B"]

# 1
(mean(x) - mean(y)) / stderr
## 1.118583 

# 2
t.test(x, y, var.equal = TRUE)$statistic
##        t 
## 1.118583 

# 3
coef(summary(lm(values ~ groups, 
  contrasts = list(groups = contr.treatment(2, base = 2)))))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 41.136929  0.8827383 46.601500 4.407622e-28
## groups1      1.396417  1.2483805  1.118583 2.728234e-01

# 4 - contr.SAS is like contr.treatment except reference level 
#  defaults to last, rather than first, level
coef(summary(lm(values ~ groups, 
  contrasts = list(groups = contr.SAS))))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 41.136929  0.8827383 46.601500 4.407622e-28
## groups1      1.396417  1.2483805  1.118583 2.728234e-01

и все они дают одну и ту же статистику t.

 # 1
(mean(y) - mean(x)) / stderr
## [1] -1.118583

# 2
t.test(y, x, var.equal = TRUE)$statistic
##         t 
## -1.118583 

# 3
coef(summary(lm(values ~ groups, 
  contrasts = list(groups = contr.treatment(2, base = 1)))))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 42.533345  0.8827383 48.183414 1.750263e-28
## groups2     -1.396417  1.2483805 -1.118583 2.728234e-01

# 4
coef(summary(lm(values ~ groups)))
##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 42.533345  0.8827383 48.183414 1.750263e-28
## groupsB     -1.396417  1.2483805 -1.118583 2.728234e-01

Примечание

Для воспроизводимости входных данных необходимо определить n и выдать set.seed из-за использования случайных чисел. Мы использовали:

 set.seed(123)
n <- 30
sim_data <- data.frame(
  groups = c(rep("A", n / 2), rep("B", n / 2)),
  values = rnorm(n, mean = 42, sd = 3.5))

1. Исправили предыдущую версию ответа.

2. Большое спасибо за ваш обстоятельный ответ! Возможно , подумайте о добавлении указателя на соответствующий код: например , t in t.test.default и split in t.test.formula , где порядок результата будет зависеть от уровней факторов, которые затем передаются t.test.default . Просто мысль, чтобы дополнить вашу превосходную эмпирическую демонстрацию.

Ответ №2:

вы можете задать справочную группу

 sim_data$groups <- relevel(factor(sim_data$groups), ref = "B")

Метки: Обратные признаки тестовой статистики с использованием lm() по сравнению с t.test()