如何按组将一列拟合值添加到数据框？答案

【问题标题】：How to add a column of fitted values to a data frame by group?如何按组将一列拟合值添加到数据框？
【发布时间】：2015-10-22 10:58:10
【问题描述】：

假设我有一个这样的数据框：

X <- data_frame(
  x = rep(seq(from = 1, to = 10, by = 1), 3),
  y = 2*x + rnorm(length(x), sd = 0.5),
  g = rep(LETTERS[1:3], each = length(x)/3))

如何拟合按变量 g 分组的回归 y~x 并将来自 fitted 和 resid 泛型方法的值添加到数据框？

我知道我能做到：

A <- X[X$g == "A",]
mA <- with(A, lm(y ~ x))
A$fit <- fitted(mA)
A$res <- resid(mA)

B <- X[X$g == "B",]
mB <- with(B, lm(y ~ x))
B$fit <- fitted(mB)
B$res <- resid(mB)

C <- X[X$g == "C",]
mC <- with(B, lm(y ~ x))
C$fit <- fitted(mC)
C$res <- resid(mC)

然后是rbind(A, B, C)。然而，在现实生活中我没有使用lm（我在quantreg 包中使用rqss）。该方法偶尔会失败，所以我需要错误处理，我想将所有失败的行放在NA 的位置。另外，有超过 3 个组，所以我不想为每个组继续复制和粘贴代码。

我尝试将dplyr 与do 一起使用，但没有取得任何进展。我在想它可能是这样的：

make_qfits <- function(data) {
  data %>%
    group_by(g) %>%
    do(failwith(NULL, rqss), formula = y ~ qss(x, lambda = 3))
}

用这种方法会很容易做到吗？在基础 R 中还有其他方法吗？

【问题讨论】：

我想到了 broom 包中的 augment.columns 函数。您在do 中尝试过哪些不起作用的方法？
@aosmith 我在想它可能类似于我刚刚添加的代码。

标签： r error-handling dplyr regression quantreg

【解决方案1】：

您可以在此任务的分组数据上使用do，在do 中拟合每个组中的模型，并将模型残差和拟合值放入data.frame。要将这些添加到原始数据中，只需在输出data.frame 中包含代表进入do 的数据的.。

在你的简单情况下，这看起来像这样：

X %>%
    group_by(g) %>%
    do({model = rqss(y ~ qss(x, lambda = 3), data = .)
        data.frame(., residuals = resid.rqss(model), fitted = fitted(model))
            })

Source: local data frame [30 x 5]
Groups: g

    x         y g     residuals    fitted
1   1  1.509760 A -1.368963e-08  1.509760
2   2  3.576973 A -8.915993e-02  3.666133
3   3  6.239950 A  4.174453e-01  5.822505
4   4  7.978878 A  4.130033e-09  7.978878
5   5 10.588367 A  4.833475e-01 10.105020
6   6 11.786445 A -3.807876e-01 12.167232
7   7 14.646221 A  4.167763e-01 14.229445
8   8 15.938253 A -3.534045e-01 16.291658
9   9 19.114927 A  7.610560e-01 18.353871
10 10 19.574449 A -8.416343e-01 20.416083
.. ..       ... .           ...       ...

如果您需要捕获错误，事情看起来会更复杂。这是使用try 并用NA 填充残差和拟合列的样子，如果对组的拟合尝试导致错误。

X[9:30,] %>%
    group_by(g) %>%
    do({catch = try(rqss(y ~ qss(x, lambda = 3), data = .))
    if(class(catch) == "try-error"){
        data.frame(., residuals = NA, fitted = NA)
    }
    else{
        model = rqss(y ~ qss(x, lambda = 3), data = .)
        data.frame(., residuals = resid.rqss(model), fitted = fitted(model))
        }
    })
Source: local data frame [22 x 5]
Groups: g

    x         y g     residuals    fitted
1   9 19.114927 A            NA        NA
2  10 19.574449 A            NA        NA
3   1  2.026199 B -4.618675e-01  2.488066
4   2  4.399768 B  1.520739e-11  4.399768
5   3  6.167690 B -1.437800e-01  6.311470
6   4  8.642481 B  4.193089e-01  8.223172
7   5 10.255790 B  1.209160e-01 10.134874
8   6 12.875674 B  8.290981e-01 12.046576
9   7 13.958278 B -4.803891e-10 13.958278
10  8 15.691032 B -1.789479e-01 15.869980
.. ..       ... .           ...       ...

【讨论】：

【解决方案2】：

对于lm 型号，您可以尝试

library(nlme)     # lmList to do lm by group
library(ggplot2)  # fortify to get out the fitted/resid data
do.call(rbind, lapply(lmList(y ~ x | g, data=X), fortify))

这会为您提供“.resid”和“.fitted”列中的残差和拟合数据以及一堆其他拟合数据。默认情况下，行名将以g 中的字母作为前缀。

rqss 模型可能会失败

do.call(rbind, lapply(split(X, X$g), function(z) {
    fit <- tryCatch({
        rqss(y ~ x, data=z)
    }, error=function(e) NULL)
    if (is.null(fit)) data.frame(resid=numeric(0), fitted=numeric(0))
    else data.frame(resid=fit$resid, fitted=fitted(fit))
}))

【讨论】：

谢谢，但这不会迫使我在幕后使用lme 吗？我正在尝试使用rqss 进行非参数拟合，因为我的真实数据是非线性的。
@wdkrnls 好的，fortify 不能与 rqss 一起使用，所以我添加了拆分应用组合解决方案
@nongkrog 谢谢！这看起来像是我可以使用的东西。但偶尔会因idObject(.Object) : invalid class "matrix.csr" object: ra and ja don't have equal lengths 而失败。如果rqss 失败，你将如何处理错误？
@wdkrnls 我会添加一个tryCatch 子句，尝试编辑，我没有失败的数据，所以我不确定
嗯，我尝试将列绑定在一起，发现行数不同。 bind_cols(X, Y) > 错误：行数不兼容（20553，预期为 21561）`

【解决方案3】：

这是一个适用于基础 R 的版本：

modelit <- function(df) {
    mB <- with(df, lm(y ~ x, na.action = na.exclude))
    df$fit <- fitted(mB)
    df$res <- resid(mB)
    return(df)
}

dfs.with.preds <- lapply(split(X, as.factor(X$g)), modelit)
output <- Reduce(function(x, y) { rbind(x, y) }, dfs.with.preds)

【讨论】：