预测具有多个单独组的线性回归

Question

我想从单个数据框中的多个组的线性回归中预测值。 我发现以下博客文章几乎可以满足我的所有需求：https://www.r-bloggers.com/2016/09/running-a-model-on-separate-groups/

但是，我不能将它与 predict() 函数与 newdata 结合使用。 对于一组，我使用以下内容：

m <- lm(y ~ x, df)
new_df <- data.frame(x=c(5))
predict(m, new_df)

这给了我在 x=5 时 y 的预测值。

当我的 df 中有多个组时，我该怎么做？这是我尝试过的：

df %>%
    nest(-group) %>%
    mutate(fit = map(data, ~ lm(.$y ~ .$x)),
           results = map(fit, predict)) %>%
    unnest(results)

当我尝试使用 results = map(fit, predict(new_df)) 时，我只得到一个错误。有没有办法可以将我的 x 值（在本例中为 5）传递到上面的代码中？

理想情况下，我会得到一个新的 data.frame，其中包含两列、组和预测的 y 值。

这是一个示例 data.frame：

group   x   y
g1  1   2
g1  1.5 3
g1  2   4
g1  2.3 4.4
g1  3   6
g1  3.4 6.2
g1  4.11    7
g1  4.8 7.9
g1  5   8
g1  5.3 8.2
g2  2   5
g2  2.3 4
g2  4   2.2
g2  4.4 1.9
g2  7   0.3

编辑：

使用ggplot2绘制样本数据，得到如下图：

ggplot(df, aes(x,y,colour=group)) +
 geom_point() +
 stat_smooth(method="lm", se=FALSE)

使用下面的代码，我得到了预测的 y 值：

predict(lm(y ~ x, df[df$group =="g1", ]), new_df)
       1 
8.180285 

predict(lm(y ~ x, df[df$group =="g2", ]), new_df)
       1 
1.732136 

我想生成一个新的数据框，它应该看起来像这样并包含 x=5 处的预测 y 值：

group   y_predict  
g1  8.180285  
g2  1.732136

Answer 1

使用注释中可重复显示的输入，由于我们只需要拟合值，我们不需要使用nest，而可以只使用mutate：

library(dplyr)

df %>%
  group_by(group) %>%
  mutate(pred = fitted(lm(y ~ x))) %>%
  ungroup %>%
  select(group, pred)

给予：

# A tibble: 15 x 2
   group    pred
   <chr>   <dbl>
 1 g1     2.47  
 2 g1     3.19  
 3 g1     3.90  
 4 g1     4.33  
 5 g1     5.33  
 6 g1     5.90  
 7 g1     6.91  
 8 g1     7.89  
 9 g1     8.18  
10 g1     8.61  
11 g2     4.41  
12 g2     4.15  
13 g2     2.63  
14 g2     2.27  
15 g2    -0.0563

这也可以这样做：

library(dplyr)

df %>%
  mutate(pred = fitted(lm(y ~ x*group + 0, df))) %>%
  select(group, pred)

或者像这样只使用base R：

transform(df, pred = fitted(lm(y ~ x*group + 0, df)))[c("group", "pred")]

或使用来自 nlme 的 lmList（R 附带，因此不必安装）：

library(dplyr)
library(nlme)

df %>%
  mutate(pred = fitted(lmList(y ~ x | group, df))) %>%
  select(group, pred)

或使用不带 dplyr 的 lmList：

library(nlme)

transform(df, pred = fitted(lmList(y ~ x | group, df)))[c("group", "pred")]

注意

Lines <- "
group   x   y
g1  1   2
g1  1.5 3
g1  2   4
g1  2.3 4.4
g1  3   6
g1  3.4 6.2
g1  4.11    7
g1  4.8 7.9
g1  5   8
g1  5.3 8.2
g2  2   5
g2  2.3 4
g2  4   2.2
g2  4.4 1.9
g2  7   0.3"
df <- read.table(text = Lines, header = TRUE)

添加

关于注释，此代码按组生成 x = 5 的预测：

df %>%
  group_by(group) %>%
  summarize(pred = predict(lm(y ~ x), list(x = 5)), .groups = "drop") %>%
  select(group, pred)
## # A tibble: 2 x 2
##   group  pred
##   <chr> <dbl>
## 1 g1     8.18
## 2 g2     1.73

Answer 2

这是使用lapply 函数的完美案例。试试这个：

linear_model <- function(x) lm(y ~ x, x)
m <- lapply(split(df,df$group),linear_model)

现在，您有一个 list 或 linear models。让我们用它来预测所有模型的 new_df 的 y 值：

new_df <- data.frame(x=c(5))
my_predict <- function(m) predict(m,new_df)
sapply(m,my_predict)

输出：

#     g1.1     g2.1 
# 8.180285 1.732136

输出是带有名称的numeric 类。

Answer 3

您所描述的是具有不同截距和斜率的估计。

与lm

您可以直接使用lm：

base = iris
names(base) = c("y", "x1", "x2", "x3", "species")

newdata = data.frame(x1 = 5, species = c("setosa", "versicolor", "virginica"))
res_1 = lm(y ~ species/x1, base)
newdata$y = predict(res_1, newdata)
newdata
#>   x1    species        y
#> 1  5     setosa 6.091450
#> 2  5 versicolor 7.865123
#> 3  5  virginica 8.414509

快捷方式species/x1表示species + species:x1，即因子变量以及因子与变量的交互作用。 因此，每个组将有一个截距和一个与x1 关联的系数（此处为species）。

然后可以像往常一样使用 predict 方法，这将导致请求的结果。这不需要循环也不需要lapply。

替代方法

另一种方法是使用专门的包来估计这种模型，例如fixest。由于它专门用于固定效应估计，因此对于大型数据集，运行时间将大大缩短。

library(fixest)

# Using variables with varying slopes
res_2 = feols(y ~ 1 | species[x1], base)
predict(res_2, newdata)
#>        1        2        3 
#> 6.091450 7.865123 8.414509

一些解释：

• 您的group 是变量species。
• feols 等效于 lm，但您可以在管道之后定义固定效果。
• species[x1] 表示 species 固定效应（即每个物种一个截距）+ x1 每个物种有一个系数（不同的斜率）。