我希望使用dplyr对data.frame进行分组,拟合线性回归,并将残差保存为原始未分组的data.frame中的一列。
下面是一个例子
> 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>
相反,我希望原始的data.frame包含一个包含残差的新列。
例如,
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
发布于 2020-01-08 13:29:38
一个似乎比迄今提出的方案更容易并更接近原问题守则的解决办法是:
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
发布于 2014-12-12 13:52:26
我从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)
更新另一种方法是在broom包中使用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
...
发布于 2017-09-20 08:13:55
由于您正在对每个组运行完全相同的回归,因此您可能会发现,只需事先将回归模型定义为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))
}
}
请注意,此函数的第一部分是为了防止出现的任何错误消息,以防您的回归在一个小于零自由度的组上运行(如果您有一个具有多个dataframe
的分组变量和多个levels
,或者有许多用于回归的独立变量(例如lm(y~ x1 + x2)
),并且无法对每个变量进行足够的非NA观测,则可能是这种情况)。
因此,您的示例可以重写如下:
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
这种方法在计算上不应该与使用augment()
的方法有太大的不同。(我不得不在包含数亿个观测值的数据集中使用这两种方法,并且认为与使用do()
函数相比,在速度方面没有明显的差异)。
另外,请注意,省略na.action = na.exclude
或使用m$residuals
而不是residuals(m)
将导致将具有NAs (在估计之前被删除)的行排除在残差的输出向量中。因此,对应的向量将没有足够的length()
以便与数据集合并,并且可能会出现一些错误消息。
https://stackoverflow.com/questions/27452491
复制