使用 lm() 与 t.test() 反转测试统计的符号
Reversed signs of test statistic using lm() vs. t.test()
为什么当使用 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()？请赐教。
你可以设置一个参考组
sim_data$groups <- relevel(factor(sim_data$groups), ref = "B")
由于 t.test 的帮助没有讨论这个，让我们写出各种可能性，以便我们可以凭经验确定它是如何工作的（或查看评论中的源代码链接）。
简短的回答是，如果我们在 lm 中使用 contr.SAS 对比（如比较部分的#4），我们可以使 lm 给出相同的结果，这与默认的 contr.treatment 对比只是最后一个水平被用作参考水平而不是第一个。
默认方法
不用公式方法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
使用 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 中使用 contrasts 参数（或使用其他答案中的 relevel）控制参考水平，也称为基准水平：
# 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))
为什么当使用 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()？请赐教。
你可以设置一个参考组
sim_data$groups <- relevel(factor(sim_data$groups), ref = "B")
由于 t.test 的帮助没有讨论这个，让我们写出各种可能性，以便我们可以凭经验确定它是如何工作的（或查看评论中的源代码链接）。
简短的回答是，如果我们在 lm 中使用 contr.SAS 对比（如比较部分的#4），我们可以使 lm 给出相同的结果，这与默认的 contr.treatment 对比只是最后一个水平被用作参考水平而不是第一个。
默认方法
不用公式方法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
使用 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 中使用 contrasts 参数（或使用其他答案中的 relevel）控制参考水平，也称为基准水平：
# 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))