在 R 中完成的多元线性回归分析中实际发生了什么？
What is actually occurring in this multiple linear regression analysis done in R?
我正在做一个用 R 语言进行数据分析的项目。
我想要做的是确定一个数据集是否可以用线性回归模型来描述，然后尝试测试该数据集的某些子组是否比整个数据集具有更强的相关性。更具体地说，我正在比较学生记录他们的脉搏和时间估计的数据集，并检查数据子组中是否存在更强的相关性，其中未发现学生与变量和子组有每日节律计算出学生在时间估计和心率方面都有每日节律。我使用的值是时间估计和心率的每日平均值。
我运行整个数据集的线性模型：
> summary(ptmod1)
Call:
lm(formula = avg.time ~ avg.pulse, data = pulsetime)
Residuals:
Min 1Q Median 3Q Max
-11.7310 -1.6725 -0.0162 2.0134 9.8548
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.82047 2.99244 22.998 <2e-16 ***
avg.pulse -0.10449 0.04115 -2.539 0.0125 *
并且还尝试运行每个子组的线性回归
> summary(ptmod2)
Call:
lm(formula = avg.time ~ avg.pulse + Group, data = pulsetime)
Residuals:
Min 1Q Median 3Q Max
-12.9884 -1.7723 -0.1873 2.4900 8.7424
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.45350 2.92287 23.420 < 2e-16 ***
avg.pulse -0.08566 0.03985 -2.149 0.03388 *
GroupOne -1.22325 0.91444 -1.338 0.18386
GroupThree 0.11062 0.97666 0.113 0.91003
GroupTwo -3.09096 0.95446 -3.238 0.00161 **
但是，我想确保我所看到的是正确的，因为我真的没想到这么多组具有显着系数。因此，我将这些组分成各自的 .csv 文件，并分别为每个组生成线性模型。将它们切割成自己的文件也使得 运行 Chow 测试更容易作为 post-hoc 分析。当我 运行 再次对它们进行回归时，我得到了完全不同的系数。
例如，这是第一组的摘要：
> summary(mod1)
Call:
lm(formula = avg.time ~ avg.pulse, data = group1)
Residuals:
Min 1Q Median 3Q Max
-7.1048 -1.6529 -0.7279 1.4063 5.6574
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 67.41445 4.15917 16.209 8.99e-15 ***
avg.pulse -0.08916 0.05657 -1.576 0.128
这让我质疑我从 ptmod2 摘要中得出的结果究竟意味着什么？我不确定如何为按单个子组排序的线性模型设置 R 代码，所以我的代码是
> ptmod2<-lm(avg.time~avg.pulse+Group, data=pulsetime)
在我的电子表格文件中，我有三列：avg.pulse、avg.time 和组。 “组”是根据子组分配的单词“一”、“二”、“三”和“四”的列。
ptmod2 的摘要是否适合整个数据集的线性回归？我真的不知道发生了什么。
非常感谢您提供的任何见解。也许我按组比较回归的代码不正确。
这在某种程度上是编程问题和统计问题之间的分离。它可能更适合 crossvalidation。但是问题很简单，可以理解。
您的问题可以分为以下子问题：
- 我是否在
ptmod2 中的完整（整个）数据集上拟合模型？
- 如何跨分组数据集估计多个模型？
- 这种情况下系数的正确分析方法是什么？
我是否在 ptmod2 中的完整（整个）数据集上拟合模型？
多空是“是”。在 R 和统计中，当您向数据集添加“组”变量时，这并不等同于将数据集拆分为多个组。相反，它添加了指示特定组的指示变量（0 或 1），包括参考水平。因此，在您的情况下，您有 4 个组，从 1 到 4，并且您正在添加一个指标来判断某人是否处于 group 1、group 2、group 3 或（参考水平）group 4.这是对组间截距差异程度的度量。例如。这些变量有解释：
If the models share a common slope avg.pulse are there a significant difference in the avg.time explained by the specific group?
您只看到 3 个组而不是 4 个的原因是第四组是通过将所有其他组设置为等于 FALSE 来解释的。例如。如果您不在第 1、2 或 3 组中，那么您属于第 4 组。因此，在第 4 组中的“效果”是 not 在第 1、2 组中的效果或 3（在本例中）。
研究这个的方法，似乎我的很多学生都觉得有用，就是研究 model.matrix 的一个小版本，例如：
data(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)
model.matrix(mpg ~ hp + cyl, data = mtcars)
您可以非常积极地看到 (intercept)、hp 有一个列，cyl6 和 cyl8 有 2 个列（没有 cyl4 这是参考）。将 cyl6 和 cyl8 中的索引与 mtcars 中的值相匹配说明 cyl6 列中的 1 表示 cyl == 6.
如何跨分组数据集估计多个模型？
根据所寻求的问题，有多种方法可以做到这一点。在您的情况下，您似乎对“avg.pulse 的效果是否对每个组有显着差异？”这个问题感兴趣。例如，您想估计每个组的 avg.pulse 系数。一个是像你后来做的那样，并估计每个组的模型
groups <- split(pulsetime, pulsetime$Group)
models <- lapply(groups, function(df)lm(avg.time ~ avg.pulse, data = df))
lapply(models, summary)
这给出了估计。这里的问题是“如何比较这些”。有这样做的方法，通过比较每个模型之间参数之间的协方差，这被称为“多元统计分析”或多元回归模型。然而，这过于复杂，因为模型共享一个共同的结果。
一种更简单的方法是通过使用指示变量为每个组添加“额外”效果来合并不同的估计。这类似于添加组变量，但不是单独添加它（表示在组 X 中的效果），我们使用
之一将它乘以相关变量
# Let each group have their own `avg.pulse` variable
ptmod2_1 <- lm(formula = avg.time ~ avg.pulse : Group, data = pulsetime)
# Let each group have their own `avg.pulse` variable and account for the effect of `Group`
ptmod2_2 <- lm(formula = avg.time ~ avg.pulse * Group, data = pulsetime)
在前面你会看到 avg.time:GroupX，对于所有 4 个组，这意味着这些是“X 组中 avg.time 的效果”，而在后面你会再次看到参考水平。请注意，这两种方法之间的明显区别在于，在后者中，所有 group 具有相同的截距，而在后者中，所有组都可以具有不同的截距。
在一般统计中，后者是首选方法，除非您有 非常好的理由 不期望每个组有不同的平均值。它与经验法则非常相似：除非你有充分的理由，否则不要测试你的拦截，即使那样你也可能不应该”。基本上是因为遵循这些规则在逻辑上很有意义（尽管可能需要几天的反思才能意识到原因）。
这种情况下系数的正确分析方法是什么？
如果您坚持使用后两种方法中的一种，则该分析类似于正态回归分析。可以使用 t-tests、anova 等（使用 summary/drop1 和 anova）来测试系数，如果你有理由可以测试组也使用标准测试进行合并（尽管如果它们无关紧要，则很少有理由以任何一种方式合并它们）。整个技巧变成了“我如何解释系数”。
对于方法 1，这是显而易见的。 “第 1 组有 avg.pulse 的影响”等等。对于方法 2，它稍微更微妙一些。 The effect of avg.pulse in group 1 is avg.pulse + avg.pulse:GroupOne. Because you have to note that avg.pulsedoes **not** disappear when you change group. It is the reference level, and every other effect is the **additional** effect onavg.pulseof going from being in group X to being in group Y. Visually yourslope` 变化图变成如果系数为正（负），则更陡（更平）。
我在下面使用 mtcars 数据集进行了可视化，使用 mpg 表示结果，hp 表示数字变量，cyl（作为一个因素）作为一个分组变量。置信区间被删除，因为它们对于说明并不重要。重要的是要注意 2 个模型有多么不同（cyl == 4 一个是正的，另一个是负的！）。这进一步说明了为什么方法 2 通常比之前的方法“更正确”。
再现性代码
下面是我用于插图和示例的代码
data(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)
model.matrix(mpg ~ hp + cyl, data = mtcars)
#split-fit across groups and run summary
groups <- split(mtcars, mtcars$cyl)
models <- lapply(groups, function(df)lm(mpg ~ hp, data = df))
lapply(models, summary)
#Fit using interaction effects
fit_11 <- lm(mpg ~ hp:cyl , mtcars)
fit_12 <- lm(mpg ~ hp*cyl , mtcars)
summary(fit_11)
summary(fit_12)
#Illustrate interaction effects
library(sjPlot)
library(sjmisc)
library(ggplot2)
library(patchwork)
theme_set(theme_sjplot())
p1 <- plot_model(fit_11, type = "pred", terms = c("hp","cyl"), ci.lvl = 0) + ggtitle("Same intercept different slope") +
geom_point(aes(x = hp, y = mpg, col = cyl, fill = cyl), mtcars, inherit.aes = FALSE)
p2 <- plot_model(fit_12, type = "pred", terms = c("hp", "cyl"), ci.lvl = 0) + ggtitle("Different intercept and slope") +
geom_point(aes(x = hp, y = mpg, col = cyl, fill = cyl), mtcars, inherit.aes = FALSE)
p1 / p2
我正在做一个用 R 语言进行数据分析的项目。
我想要做的是确定一个数据集是否可以用线性回归模型来描述，然后尝试测试该数据集的某些子组是否比整个数据集具有更强的相关性。更具体地说，我正在比较学生记录他们的脉搏和时间估计的数据集，并检查数据子组中是否存在更强的相关性，其中未发现学生与变量和子组有每日节律计算出学生在时间估计和心率方面都有每日节律。我使用的值是时间估计和心率的每日平均值。
我运行整个数据集的线性模型：
> summary(ptmod1)
Call:
lm(formula = avg.time ~ avg.pulse, data = pulsetime)
Residuals:
Min 1Q Median 3Q Max
-11.7310 -1.6725 -0.0162 2.0134 9.8548
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.82047 2.99244 22.998 <2e-16 ***
avg.pulse -0.10449 0.04115 -2.539 0.0125 *
并且还尝试运行每个子组的线性回归
> summary(ptmod2)
Call:
lm(formula = avg.time ~ avg.pulse + Group, data = pulsetime)
Residuals:
Min 1Q Median 3Q Max
-12.9884 -1.7723 -0.1873 2.4900 8.7424
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.45350 2.92287 23.420 < 2e-16 ***
avg.pulse -0.08566 0.03985 -2.149 0.03388 *
GroupOne -1.22325 0.91444 -1.338 0.18386
GroupThree 0.11062 0.97666 0.113 0.91003
GroupTwo -3.09096 0.95446 -3.238 0.00161 **
但是，我想确保我所看到的是正确的，因为我真的没想到这么多组具有显着系数。因此，我将这些组分成各自的 .csv 文件，并分别为每个组生成线性模型。将它们切割成自己的文件也使得 运行 Chow 测试更容易作为 post-hoc 分析。当我 运行 再次对它们进行回归时，我得到了完全不同的系数。 例如，这是第一组的摘要：
> summary(mod1)
Call:
lm(formula = avg.time ~ avg.pulse, data = group1)
Residuals:
Min 1Q Median 3Q Max
-7.1048 -1.6529 -0.7279 1.4063 5.6574
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 67.41445 4.15917 16.209 8.99e-15 ***
avg.pulse -0.08916 0.05657 -1.576 0.128
这让我质疑我从 ptmod2 摘要中得出的结果究竟意味着什么？我不确定如何为按单个子组排序的线性模型设置 R 代码，所以我的代码是
> ptmod2<-lm(avg.time~avg.pulse+Group, data=pulsetime)
在我的电子表格文件中，我有三列：avg.pulse、avg.time 和组。 “组”是根据子组分配的单词“一”、“二”、“三”和“四”的列。
ptmod2 的摘要是否适合整个数据集的线性回归？我真的不知道发生了什么。
非常感谢您提供的任何见解。也许我按组比较回归的代码不正确。
这在某种程度上是编程问题和统计问题之间的分离。它可能更适合 crossvalidation。但是问题很简单，可以理解。
您的问题可以分为以下子问题：
- 我是否在
ptmod2中的完整（整个）数据集上拟合模型？ - 如何跨分组数据集估计多个模型？
- 这种情况下系数的正确分析方法是什么？
我是否在 ptmod2 中的完整（整个）数据集上拟合模型？
多空是“是”。在 R 和统计中，当您向数据集添加“组”变量时，这并不等同于将数据集拆分为多个组。相反，它添加了指示特定组的指示变量（0 或 1），包括参考水平。因此，在您的情况下，您有 4 个组，从 1 到 4，并且您正在添加一个指标来判断某人是否处于 group 1、group 2、group 3 或（参考水平）group 4.这是对组间截距差异程度的度量。例如。这些变量有解释：
If the models share a common slope
avg.pulseare there a significant difference in theavg.timeexplained by the specific group?
您只看到 3 个组而不是 4 个的原因是第四组是通过将所有其他组设置为等于 FALSE 来解释的。例如。如果您不在第 1、2 或 3 组中，那么您属于第 4 组。因此，在第 4 组中的“效果”是 not 在第 1、2 组中的效果或 3（在本例中）。
研究这个的方法，似乎我的很多学生都觉得有用，就是研究 model.matrix 的一个小版本，例如：
data(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)
model.matrix(mpg ~ hp + cyl, data = mtcars)
您可以非常积极地看到 (intercept)、hp 有一个列，cyl6 和 cyl8 有 2 个列（没有 cyl4 这是参考）。将 cyl6 和 cyl8 中的索引与 mtcars 中的值相匹配说明 cyl6 列中的 1 表示 cyl == 6.
如何跨分组数据集估计多个模型？
根据所寻求的问题，有多种方法可以做到这一点。在您的情况下，您似乎对“avg.pulse 的效果是否对每个组有显着差异？”这个问题感兴趣。例如，您想估计每个组的 avg.pulse 系数。一个是像你后来做的那样，并估计每个组的模型
groups <- split(pulsetime, pulsetime$Group)
models <- lapply(groups, function(df)lm(avg.time ~ avg.pulse, data = df))
lapply(models, summary)
这给出了估计。这里的问题是“如何比较这些”。有这样做的方法，通过比较每个模型之间参数之间的协方差，这被称为“多元统计分析”或多元回归模型。然而，这过于复杂，因为模型共享一个共同的结果。
一种更简单的方法是通过使用指示变量为每个组添加“额外”效果来合并不同的估计。这类似于添加组变量，但不是单独添加它（表示在组 X 中的效果），我们使用
之一将它乘以相关变量# Let each group have their own `avg.pulse` variable
ptmod2_1 <- lm(formula = avg.time ~ avg.pulse : Group, data = pulsetime)
# Let each group have their own `avg.pulse` variable and account for the effect of `Group`
ptmod2_2 <- lm(formula = avg.time ~ avg.pulse * Group, data = pulsetime)
在前面你会看到 avg.time:GroupX，对于所有 4 个组，这意味着这些是“X 组中 avg.time 的效果”，而在后面你会再次看到参考水平。请注意，这两种方法之间的明显区别在于，在后者中，所有 group 具有相同的截距，而在后者中，所有组都可以具有不同的截距。
在一般统计中，后者是首选方法，除非您有 非常好的理由 不期望每个组有不同的平均值。它与经验法则非常相似：除非你有充分的理由，否则不要测试你的拦截，即使那样你也可能不应该”。基本上是因为遵循这些规则在逻辑上很有意义（尽管可能需要几天的反思才能意识到原因）。
这种情况下系数的正确分析方法是什么？
如果您坚持使用后两种方法中的一种，则该分析类似于正态回归分析。可以使用 t-tests、anova 等（使用 summary/drop1 和 anova）来测试系数，如果你有理由可以测试组也使用标准测试进行合并（尽管如果它们无关紧要，则很少有理由以任何一种方式合并它们）。整个技巧变成了“我如何解释系数”。
对于方法 1，这是显而易见的。 “第 1 组有 avg.pulse 的影响”等等。对于方法 2，它稍微更微妙一些。 The effect of avg.pulse in group 1 is avg.pulse + avg.pulse:GroupOne. Because you have to note that avg.pulsedoes **not** disappear when you change group. It is the reference level, and every other effect is the **additional** effect onavg.pulseof going from being in group X to being in group Y. Visually yourslope` 变化图变成如果系数为正（负），则更陡（更平）。
我在下面使用 mtcars 数据集进行了可视化，使用 mpg 表示结果，hp 表示数字变量，cyl（作为一个因素）作为一个分组变量。置信区间被删除，因为它们对于说明并不重要。重要的是要注意 2 个模型有多么不同（cyl == 4 一个是正的，另一个是负的！）。这进一步说明了为什么方法 2 通常比之前的方法“更正确”。
再现性代码
下面是我用于插图和示例的代码
data(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)
model.matrix(mpg ~ hp + cyl, data = mtcars)
#split-fit across groups and run summary
groups <- split(mtcars, mtcars$cyl)
models <- lapply(groups, function(df)lm(mpg ~ hp, data = df))
lapply(models, summary)
#Fit using interaction effects
fit_11 <- lm(mpg ~ hp:cyl , mtcars)
fit_12 <- lm(mpg ~ hp*cyl , mtcars)
summary(fit_11)
summary(fit_12)
#Illustrate interaction effects
library(sjPlot)
library(sjmisc)
library(ggplot2)
library(patchwork)
theme_set(theme_sjplot())
p1 <- plot_model(fit_11, type = "pred", terms = c("hp","cyl"), ci.lvl = 0) + ggtitle("Same intercept different slope") +
geom_point(aes(x = hp, y = mpg, col = cyl, fill = cyl), mtcars, inherit.aes = FALSE)
p2 <- plot_model(fit_12, type = "pred", terms = c("hp", "cyl"), ci.lvl = 0) + ggtitle("Different intercept and slope") +
geom_point(aes(x = hp, y = mpg, col = cyl, fill = cyl), mtcars, inherit.aes = FALSE)
p1 / p2