使用`dplyr`保存残差

Question

我希望使用dplyr对数据框进行分组，拟合线性回归并将残差保存为原始未分组数据框中的一列。

以下是一个例子：

> iris %>%
   select(Sepal.Length, Sepal.Width) %>%
   group_by(Species) %>%
   do(mod = lm(Sepal.Length ~ Sepal.Width, data=.)) %>%

返回值：

     Species     mod
1     setosa <S3:lm>
2 versicolor <S3:lm>
3  virginica <S3:lm>

相反，我希望原始数据框中有一个包含残差的新列。

例如，

    Sepal.Length Sepal.Width  resid
1   5.1         3.5  0.04428474
2   4.9         3.0  0.18952960
3   4.7         3.2 -0.14856834
4   4.6         3.1 -0.17951937
5   5.0         3.6 -0.12476423
6   5.4         3.9  0.06808885

- Austin Richardson

Answer 1

我从http://jimhester.github.io/plyrToDplyr/中找到了一个示例，并进行了适应。

r <- iris %>%
  group_by(Species) %>%
  do(model = lm(Sepal.Length ~ Sepal.Width, data=.)) %>%
  do((function(mod) {
     data.frame(resid = residuals(mod$model))
  })(.))

corrected <- cbind(iris, r)

更新 另一种方法是使用包中的augment函数:

r <- iris %>%
  group_by(Species) %>%
  do(augment(lm(Sepal.Length ~ Sepal.Width, data=.))

这会返回：

Source: local data frame [150 x 10]
Groups: Species

   Species Sepal.Length Sepal.Width  .fitted    .se.fit      .resid       .hat
1   setosa          5.1         3.5 5.055715 0.03435031  0.04428474 0.02073628
2   setosa          4.9         3.0 4.710470 0.05117134  0.18952960 0.04601750
3   setosa          4.7         3.2 4.848568 0.03947370 -0.14856834 0.02738325
4   setosa          4.6         3.1 4.779519 0.04480537 -0.17951937 0.03528008
5   setosa          5.0         3.6 5.124764 0.03710984 -0.12476423 0.02420180
...

- Austin Richardson

---

我理解正在发生的事情，但我自己永远不会想到这一点。例如，为什么在第二个“do”中需要一个匿名函数，而在第一个“do”中不需要呢？ - Austin Richardson

Answer 2

一个比之前提出的解决方案更简单且更接近原始问题代码的解决方案是：

iris %>%
   group_by(Species) %>%
   do(data.frame(., resid = residuals(lm(Sepal.Length ~ Sepal.Width, data=.))))

结果：

# A tibble: 150 x 6
# Groups:   Species [3]
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species   resid
          <dbl>       <dbl>        <dbl>       <dbl> <fct>     <dbl>
 1          5.1         3.5          1.4         0.2 setosa   0.0443
 2          4.9         3            1.4         0.2 setosa   0.190 
 3          4.7         3.2          1.3         0.2 setosa  -0.149 
 4          4.6         3.1          1.5         0.2 setosa  -0.180 
 5          5           3.6          1.4         0.2 setosa  -0.125 
 6          5.4         3.9          1.7         0.4 setosa   0.0681
 7          4.6         3.4          1.4         0.3 setosa  -0.387 
 8          5           3.4          1.5         0.2 setosa   0.0133
 9          4.4         2.9          1.4         0.2 setosa  -0.241 
10          4.9         3.1          1.5         0.1 setosa   0.120 

- Gilles San Martin

Answer 3

由于您需要对每个组运行完全相同的回归分析，因此您可能会发现，事先将回归模型定义为一个function()，然后使用mutate为每个组执行它会更加简单。

model<- function(y,x){ 
  a<- y + x 
  if( length(which(!is.na(a))) <= 2  ){
    return( rep(NA, length(a)))
  } else {
    m<- lm( y ~ x, na.action = na.exclude)
    return( residuals(m))
    } 
}

请注意，此函数的第一部分是为了防止出现任何错误消息，以防您的回归在自由度小于零的组上运行（如果您有一个包含多个分组变量和许多级别或回归的大量独立变量的数据框架，例如lm(y~ x1 + x2)，并且不能检查每个变量是否具有足够的非NA观察值，则可能会出现这种情况）。请注意，保留HTML标记。

因此，您的示例可以重写如下：

iris %>% group_by(Species) %>% 
  mutate(resid = model(Sepal.Length,Sepal.Width) ) %>% 
  select(Sepal.Length,Sepal.Width,resid)

应该产生以下结果：

   Species Sepal.Length Sepal.Width       resid
    <fctr>        <dbl>       <dbl>       <dbl>
1   setosa          5.1         3.5  0.04428474
2   setosa          4.9         3.0  0.18952960
3   setosa          4.7         3.2 -0.14856834
4   setosa          4.6         3.1 -0.17951937
5   setosa          5.0         3.6 -0.12476423
6   setosa          5.4         3.9  0.06808885

这种方法在计算上与使用的方法没有太大区别。(我曾在包含几亿个观测值的数据集上使用了这两种方法，认为在速度方面与使用函数相比没有显著差异)。

另外，请注意省略，或者使用而不是，会导致从残差输出向量中排除具有NAs的行(在估计之前被删除)。因此，相应的向量将没有足够的来与数据集合并，可能会出现一些错误信息。

- steffchy