R语言t.test函数执行t检验的底层数学公式与结果差异排查

问题 #

请问是否有人知晓R语言中t.test函数执行t检验所采用的具体数学公式？
我在互联网上检索到多种不同的t检验数学公式，曾尝试使用如下3种公式手动完成t检验计算（对应公式见2张附图）：

但将手动计算得到的p值与R中t.test函数的输出结果对比时，二者差异极大：

• 使用3种不同公式手动计算得到的p值均为不显著结果
• 调用R内置t.test函数计算得到的p值为显著结果

我的测试代码如下：

#Loading the data
library("lingpsych")
data("df_gibsonwu")

#Preparing the data
bysubj <- aggregate(rt~subj + type, mean, data=df_gibsonwu)
ORSR <- bysubjOR$rt - bysubjSR$rt
meanOR <- mean(bysubjOR$rt)
meanSR <- mean(bysubjSR$rt)

# Paired t-test by hand:
n <- 37

SE0 <- sqrt(((sd(bysubjOR$rt)^2) /n)+((sd(bysubjSR$rt)^2)/n))
SE1 <- sqrt(((n * sum(ORSR))^2 - (sum(ORSR)))^2 / n-1)
SE2 <- sd(bysubj$rt) / sqrt(n)

tvalue0 <- (meanOR - meanSR) / SE0
tvalue1 <- (sum(ORSR)) / SE1 #Which I think, is the correct formula for a paired-t-test?
tvalue2 <- (meanOR - meanSR) / SE2

critT <- qt(0.975,n-1)
pvalue0 <- 2*(pt(tvalue0, n-1, lower.tail=FALSE))
pvalue1 <- 2*(pt(tvalue1, n-1, lower.tail=FALSE))
pvalue2 <- 2*(pt(tvalue2, n-1, lower.tail=FALSE))

# Paired t-test using the R function:
t.test(bysubj$rt ~ bysubj$type, paired=TRUE)

运行得到的结果如下：
pvalue0 = 1.959，pvalue1 = 1.000，pvalue2 = 1.994
调用R的t.test函数得到的p值为：0.01248
提前感谢各位的解答！:)

回答 #

R中调用t.test()设置paired=TRUE时，执行的本质是配对差值的单样本t检验，具体数学逻辑非常明确：

1. 对n对配对观测，逐对计算差值$d_i = x_i - y_i$（必须保证x和y的顺序一一对应，不能打乱配对关系）
2. 计算差值的样本均值$\bar{d} = \frac{1}{n}\sum_{i=1}^n d_i$
3. 计算差值的样本标准差$s_d = \sqrt{\frac{1}{n-1}\sum_{i=1}^n (d_i - \bar{d})^2}$
4. 标准误$SE = \frac{s_d}{\sqrt{n}}$
5. t统计量$t = \frac{\bar{d}}{SE}$，自由度为$df = n-1$，再根据t分布计算对应p值即可。

手动计算结果和内置函数差异极大，核心是代码里有两个致命错误：

1. 数据预处理步骤缺失

代码里直接调用了bysubjOR、bysubjSR两个对象，但这两个对象从来没有被定义过，aggregate输出的bysubj是长格式数据，必须先按type拆分才能拿到两组配对的rt值，否则后续计算全是无效值。
正确的拆分代码：

bysubjOR <- bysubj[bysubj$type == "OR", ]
bysubjSR <- bysubj[bysubj$type == "SR", ]
# 要保证两个子集的subj顺序完全一致，差值才是配对的
ORSR <- bysubjOR$rt - bysubjSR$rt

2. 三个标准误的公式全部写错

• SE0是独立样本t检验的标准误计算方式，完全没考虑配对数据两组的相关性，不适用于配对t检验。
• SE1的公式结构完全错误，平方和的计算逻辑、运算优先级全不对，算出来的数值没有任何统计意义。
• SE2直接用了长格式全量数据的标准差算标准误，既没分组也没计算配对差值，结果必然错误。

可复现的正确手动计算代码

运行下面的代码，得到的结果会和R内置t.test()输出完全一致：

# 加载数据
library("lingpsych")
data("df_gibsonwu")

# 数据预处理
bysubj <- aggregate(rt~subj + type, mean, data=df_gibsonwu)
bysubjOR <- bysubj[bysubj$type == "OR", ]
bysubjSR <- bysubj[bysubj$type == "SR", ]
d <- bysubjOR$rt - bysubjSR$rt
n <- length(d)

# 按公式手动计算
d_mean <- mean(d)
d_sd <- sd(d)
SE <- d_sd / sqrt(n)
t_val <- d_mean / SE
p_val <- 2 * pt(abs(t_val), df = n-1, lower.tail = FALSE)

# 对比内置函数结果
t_val # 约2.62
p_val # 约0.01248
t.test(bysubj$rt ~ bysubj$type, paired=TRUE)

内容的提问来源于stack exchange，提问作者VicEllis