在R的多元线性回归分析中到底发生了什么？

Question

我正在做一个用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文件，并为每个组分别生成线性模型。将它们放入自己的文件中，还可以更容易地将周星驰测试作为一种事后分析。当我再次对它们进行回归时，我得到了完全不同的系数。例如，下面是第一组的摘要：

> 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和Group。“组”是根据分组分配的“一”、“二”、“三”和“四”等字的一栏。

ptmod2的总结是否适合于整个数据集的线性回归？我真的不知道发生了什么。

非常感谢你所能提供的洞察力。也许我按组比较回归的代码是不正确的。

Answer 1

这在某种程度上是编程和统计问题之间的分歧。它可能更适合于crossvalidation。然而，这个问题很简单，足以理解。

你的问题可以分为以下几个子问题：

我在ptmod2

• How中的完整(整体)数据集上拟合一个模型，我在分组数据集中估计多个模型吗？

• ，如何正确地分析这种情况的系数？

我是否在ptmod2的完整(完整)数据集上拟合模型？

长和短都是“是”。在R和statistics中，当向数据集添加“组”变量时，这不等于将数据集拆分为多个组。相反，它正在添加指示符变量(0或1)，指示特定组，包括引用级别。因此，在您的例子中，您有4个组，从1到4，并且添加了一个指示符，指示某人是在group 1、group 2、group 3还是(参考级别) group 4中。这是一种衡量各个组之间的截距有多大差异的度量。例如：这些变量的解释如下：

，如果模型有一个共同的斜率，avg.pulse在特定组所解释的avg.time中有显着性差异吗？

您只看到3个组而不是4个组的原因是，通过将所有其他组设置为FALSE来解释第四个组。例如：如果你不在第1组、第2组或第3组，你是第4组的一部分。所以，第4组的“效果”是而不是存在于第1、2或3组(本例中)。

研究这个问题的一种方法，我的许多学生似乎觉得很有帮助，那就是学习一个小版本的model.matrix，例如：

data(mtcars)
mtcars$cyl <- as.factor(mtcars$cyl)
model.matrix(mpg ~ hp + cyl, data = mtcars)

您可以非常积极地看到，(intercept)、hp有一个列，cyl6和cyl8有两个列( 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个组，这意味着它们是“avg.time在组X中的影响”，而在后一组中，您将再次拥有一个引用级别。注意，这两个方法之间的一个明显区别是，在后一个方法中，所有group都具有相同的截获，而在后者中，所有组都可以具有不同的截获。

在一般统计中，后者是首选方法，除非您有一个--非常好的理由--不期望每个组具有不同的平均值。这与经验法则非常相似:除非你有一个很好的理由，否则不要测试你的拦截，即使这样你也不应该“。基本上是因为遵循这些规则有很多逻辑上的意义(尽管要花几天的时间才能意识到原因)。

分析这种情况的系数的正确方法是什么？

如果您坚持使用后两种方法中的一种，则分析类似于正常的回归分析。系数可以使用t-tests、anova等进行测试(使用summary/drop1和anova)，如果您有理由可以使用标准测试测试组合并(尽管如果它们不重要，则很少有理由以任何方式合并它们)。整个技巧变成了“我如何解释系数”。

对于方法1，这是非常明显的。“第一组的avg.pulse效应如此之大”等等。对于方法2，它稍微微妙一些。如果系数为正(负)，则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 (作为因素)作为分组变量。置信区间被移除，因为它们对插图不重要。重要的是要注意两个模型是如何不同的 (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