配对数据比较：回归还是配对 t 检验？

机器算法验证 回归 t检验 配对数据
2022-03-27 14:33:23

问题
我总是使用配对 t 检验或 wilcoxon 符号秩检验（当然取决于数据集）来检查两种方法（平均而言）是否产生相同的结果。在了解了有关回归的更多信息之后，我认为这也适用于回归，但是我不明白在哪种情况下哪个会“更好”？

示例
让我们以这个（太小）示例数据集为例，并假设它是正态分布的。

data <- read.table(text = "  sample methodx methody
                   1      1    0.52    0.53
                   2      2    0.50    0.51
                   3      3    0.48    0.48
                   4      4    0.40    0.41
                   5      5    0.36    0.36
                   6      6    0.30    0.32
                   7      7    0.28    0.30
                   8      8    0.28    0.29", header = T)

# Regression analysis
model <- lm(data$methodx ~ data$methody)
summary(model)

# Residuals:
#   Min        1Q     Median        3Q       Max 
# -0.007317 -0.004931 -0.002012  0.004596  0.011341 
#
# Coefficients:
#             Estimate Std. Error    t value  Pr(>|t|)    
#  (Intercept)  -0.02341    0.01181  -1.983   0.0946 .  
#  data$methody  1.03354    0.02879  35.900   3.11e-08 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.007374 on 6 degrees of freedom
# Multiple R-squared:  0.9954,  Adjusted R-squared:  0.9946 
# F-statistic:  1289 on 1 and 6 DF,  p-value: 3.115e-08

# Paired t-test
t.test(data$methodx, data$methody, paired = TRUE)

# Paired t-test
#
# data:  data$methodx and data$methody
# t = -3.7417, df = 7, p-value = 0.007247
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
#   -0.016319724 -0.003680276
# sample estimates:
#   mean of the differences 
#                     -0.01

在查看回归时：我看到一个高相关性 ( 0.9954)，这似乎是线性的，因为线的 rc 是1.03354。配对 t 检验告诉我拒绝 H0，可能是因为这个数据集太小了。但总的来说，两者似乎都能告诉我这些方法是否平均给出了相同的结果。那么在比较两种方法时，什么时候选择线性回归，什么时候选择配对t检验呢​​？

3个回答

您以“相关样式”设置回归的方式确实回答了一个不同的问题（@olooney 很好地解释了这一点）。但是，您可以使用回归来测试均值差异。最后，这就是我们所说的 AN(C)OVA。同样，均值差的t检验可以看作是回归斜率检验的空间案例。

为了估计和之间的平均差异methodx，我们首先将数据转换为长格式。methodylm()

library(reshape2)

d_long <- melt(data, id.vars = "sample", variable.name = "method")
d_long$method <- relevel(d_long$method, ref = "methody")

执行以下回归会产生正确的均值差 (-0.01)，它由斜率估计值表示。截距的估计对应于参考组的平均值，即methody。

summary(lm(value ~ method, d_long))
[...]
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.40000    0.03485  11.479 1.65e-08 ***
methodmethodx -0.01000    0.04928  -0.203    0.842 
[...]

但是，标准误差/自由度是错误的，因为该模型将每个样本视为独立观察。因此，该模型为您提供了一个对应于非配对t检验的p值。 

t.test(data$methodx, data$methody, var.equal = TRUE) 
    Two Sample t-test

data:  data$methodx and data$methody
t = -0.20292, df = 14, p-value = 0.8421
[...]
sample estimates:
mean of x mean of y 
     0.39      0.40

为了合并依赖结构，我们可以使用多级回归/线性混合模型，并使用 Kenward-Roger 方法近似自由度。此分析和配对t检验密切相关（例如，请参阅配对 t 检验作为线性混合效应建模的特例），在这种情况下，它为您提供与配对t检验完全相同的结果：

library(lmerTest)

summary(lmer(value ~ method + (1|sample), d_long), ddf = "Kenward-Roger")
[...]
Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)    0.400000   0.034847  7.020619  11.479 8.37e-06 ***
methodmethodx -0.010000   0.002673  7.000000  -3.742  0.00725 ** 
[...]

您的 t 检验正在回答您想要的问题，即（用您自己的话）“检查两种方法（平均而言）是否产生相同的结果”，并且在这方面您的分析看起来是正确的。鉴于您的样本量较小 n=7，这很简单、正确且合适。

但是，您的回归模型并未设置为回答相同的问题，因此直接比较这两种方法是没有意义的。让我们先看看你指定的模型实际上在做什么。

给出的模型只是methodx从预测methody。这个模型总是形式的methodx = slope * methody + Intercept。

模型图

如果 H0 为真，我们预计slope为 1 和Intercept0。事实上，尤其是在 H0 下，我们预计相关性非常高，因此您看到它的事实告诉我们很少。不仅如此，如果 Intercept 发生变化，相关性实际上根本不会改变！所以我会完全忽略R2默认回归分析中的数字 - 它们告诉我们这个特定问题没有什么有趣的。

为了证伪 H0，我们必须做出并捍卫以下陈述之一：

  1. 截距明显不同于 0。
  2. 斜率与 1 显着不同。
  3. 斜率和截距的组合与 1/0 显着不同。

声明 #1 可以（几乎）从回归诊断中读取：t 值和 p 值用于 Intercept 不同于零的假设检验。但请注意，与进行直接 t 检验时相比，t 值要低得多，而 p 值要高得多。事实上，我们必须将单个系数的 t 检验理解为控制模型中的所有其他变量：在这种情况下是斜率。这是一个与我们最初提出的问题根本不同的问题。

语句#2不能直接从回归诊断中读取，因为默认情况下它将斜率与 0 进行比较，而您希望将其与斜率 1 进行比较，因为这就是您的零假设下的斜率。您可以通过从拟合斜率系数中减去一个，除以第二列中斜率参数的给定“标准误差”并将 t 检验应用于生成的 t 统计量，将其作为练习自己做。但是，我们仍然存在将斜率和截距混合在一起的问题，而且这不是您要寻找的“均值差异”，因此我也不建议这样做。

对于声明#3，我们希望两者都做。我们可以用 ANOVA 做到这一点：

首先，我们构建了一个强制斜率为 1 且截距为 0 的（有些微不足道的）模型。该模型实际上没有自由参数，但 R 很乐意lm为我们构建一个类对象，这是我们想要的。

null_model <- lm(methodx ~ offset(1*methody) -1, data=data)

然后，您可以将其与上面适合的模型进行比较：

anova(model, null_model)

当我为您的数据执行此操作时，我得到：

方差表分析

Model 1: methodx ~ methody
Model 2: methodx ~ offset(1 * methody) - 1
  Res.Df        RSS Df   Sum of Sq      F  Pr(>F)  
1      6 0.00032622                                
2      8 0.00120000 -2 -0.00087378 8.0355 0.02009 *

所以综合 p 值为 0.02。这基本上是在底层使用 F 检验并提出一个问题，“当使用更复杂的模型时，我是否会解释更多额外的 RSS，而不是纯粹的机会解释？”

如果您不确定与 methodx 和 methody 相关的真实斜率是否恰好为 1 ，您可能想要使用这种综合检验。这是我们在简单地对组间均值差异应用 t 检验时获得的唯一好处。请注意，这意味着对于我们观察到的任何斜率差异，我们都在“接受”（在报告较低、更显着的 p 值的意义上）。根据实验，这实际上可能是完全倒退的——例如，不同于 1 的斜率可能表明两次测量之间的方法不一致。

两种方法回答的问题略有不同。t 检验是关于均值的，正如您使用的那样，回归是关于找到最佳线性关系。当然，如果截距是0并且斜率是1，那么一切都很容易。但是，如果不是呢？

x = runif(30, -5, 10)
y = jitter(1.2*x)

summary(lm(y~x))

截距是0，斜率不是。你能从回归结果中看出均值是否相同吗？你不能。将有不同的方法取决于您是否开始

x1 = runif(30, -5, 10)
x2 = runif(30, -5, 5)
x3 = runif(30, -5, 0)

你无法从回归的结果中看出这一点。所以不：两种方法的报告都不一样。

其它你可能感兴趣的问题
上一篇生存分析与否？ 下一篇在 R 包 mgcv 中，对两个连续变量进行平滑随机效应是否有效？