解释使用 ggplot 可视化的 R 中线性模型的代码

Question

我想了解在分析基因表达数据时如何使用线性模型来替代 t 检验。对于单个基因，我在第 1 组（n=10）和第 2 组（n=10）中共有 20 个基因表达值的数据框。

gexp = data.frame(expression = c(2.7,0.4,1.8,0.8,1.9,5.4,5.7,2.8,2.0,4.0,3.9,2.8,3.1,2.1,1.9,6.4,7.5,3.6,6.6,5.4),
                 group = c(rep(1, 10), rep(2, 10)))

可以使用 ggplot 绘制数据（方框），如下所示：

plot <- gexp %>%
  ggplot(aes(x = group, y = expression)) +
  geom_boxplot() +
  geom_point()
plot

我希望使用回归公式对第 1 组和第 2 组中的表达式建模： Y = Beta0 + (Beta1 x X) + e 其中 Y 是我要建模的表达式，X 代表分别编码为 0 和 1 的两个组。因此，组 1 中的表达式（当 x = 0 时）等于 Beta0；组 2 中的表达式（当 x = 1 时）等于 Beta0 + Beta1。如果这是建模：

mod1 <- lm(expression ~ group, data = gexp)
mod1

以上代码输出的截距为 2.75，斜率为 1.58。是我不明白的线性模型的可视化。如果您能清楚地解释以下代码，我将不胜感激：

plot +
  geom_point(data = data.frame(x = c(1, 2), y = c(2.75, 4.33)),
               aes(x = x, y = y),
               colour = "red", size = 5) +
    geom_abline(intercept = coefficients(mod1)[1] - coefficients(mod1)[2],
                slope = coefficients(mod1)[2])

我明白为什么选择 data.frame 值（4.33 的值是截距 Beta0 和斜率 Beta1 的总和），但它是 geom_abline 参数我不明白。为什么截距计算如图所示？在我使用它的文本中，'...我们需要在绘制线性模型时从截距中减去斜率，因为组 1 和组 2 在模型中被编码为 0 和 1，但绘制为图中的 1 和 2。' 我不理解这一点，请提供解释，但不要太专业。

Answer 1

如果 group 变量被编码为一个因子，我相信你的代码是正确的。

library(ggplot2)
gexp = data.frame(expression = c(2.7,0.4,1.8,0.8,1.9,5.4,5.7,2.8,2.0,4.0,3.9,2.8,3.1,2.1,1.9,6.4,7.5,3.6,6.6,5.4),
                  group = factor(c(rep(1, 10), rep(2, 10))))

plot <- 
  ggplot(gexp, aes(x = group, y = expression)) +
  geom_boxplot() +
  geom_point()

mod1 <- lm(expression ~ group, data = gexp)

plot +
  geom_point(data = data.frame(x = c(1, 2), y = c(2.75, 4.33)),
             aes(x = x, y = y),
             colour = "red", size = 5) +
  geom_abline(intercept = coefficients(mod1)[1] - coefficients(mod1)[2],
              slope = coefficients(mod1)[2])

^{由 reprex package (v2.0.1)}

于 2022-03-30 创建

要了解指定线性模型时因子和整数之间的区别，您可以查看模型矩阵。

model.matrix(y ~ f, data = data.frame(f = 1:3, y = 1))
#>   (Intercept) f
#> 1           1 1
#> 2           1 2
#> 3           1 3
#> attr(,"assign")
#> [1] 0 1

model.matrix(y ~ f, data = data.frame(f = factor(1:3), y = 1))
#>   (Intercept) f2 f3
#> 1           1  0  0
#> 2           1  1  0
#> 3           1  0  1
#> attr(,"assign")
#> [1] 0 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$f
#> [1] "contr.treatment"

^{由 reprex package (v2.0.1)}

于 2022-03-30 创建

在第一个模型矩阵中，您指定的就是您得到的：您将某些东西建模为截距和 f 变量的函数。在此模型中，您考虑 f = 2 是 f = 1.

的两倍

当 f 是一个因素时，这会有点不同。 k 级因子被拆分为 k-1 个虚拟变量，其中每个虚拟变量都用 1 或 0 编码，无论它是否偏离参考水平（第一个因子水平）。通过以这种方式建模，您不会认为第二个因子水平可能是第一个因子水平的两倍。

因为在 ggplot2 中，第一个因子水平显示在位置 = 1 而不是位置 = 0（它是如何建模的），你计算的截距是关闭的。您需要从计算的截距中减去 1 * slope 以使其在 ggplot2 中正确显示。

解释使用 ggplot 可视化的 R 中线性模型的代码

Explain the code underlying a linear model in R visualised with ggplot

r

linear-regression

ggplot2