预测多个独立组的线性回归
Predict linear regression with multiple separate groups
我想预测单个数据框中多个组的线性回归值。
我发现以下博客文章几乎可以满足我的所有需求：https://www.r-bloggers.com/2016/09/running-a-model-on-separate-groups/
但是，我无法将其与带有新数据的 predict() 函数结合使用。
对于一组，我使用以下内容：
m <- lm(y ~ x, df)
new_df <- data.frame(x=c(5))
predict(m, new_df)
这给出了 y 在 x=5 时的预测值。
当我的 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
使用注释中可重复显示的输入，因为我们只需要拟合值，所以不需要使用 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
这是使用 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 class 和名字。
您所描述的是具有不同截距和斜率的估计。
与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，即因子变量和因子与变量之间的交互作用。
因此，每组（此处为 species）将有一个截距和一个系数与 x1 相关联。
然后可以像往常一样使用预测方法，这将导致请求的结果。这是在不需要循环或 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 每个物种有一个系数（变化的斜率）。
我想预测单个数据框中多个组的线性回归值。 我发现以下博客文章几乎可以满足我的所有需求：https://www.r-bloggers.com/2016/09/running-a-model-on-separate-groups/
但是，我无法将其与带有新数据的 predict() 函数结合使用。 对于一组，我使用以下内容：
m <- lm(y ~ x, df)
new_df <- data.frame(x=c(5))
predict(m, new_df)
这给出了 y 在 x=5 时的预测值。
当我的 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
使用注释中可重复显示的输入，因为我们只需要拟合值，所以不需要使用 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
这是使用 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 class 和名字。
您所描述的是具有不同截距和斜率的估计。
与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，即因子变量和因子与变量之间的交互作用。
因此，每组（此处为 species）将有一个截距和一个系数与 x1 相关联。
然后可以像往常一样使用预测方法，这将导致请求的结果。这是在不需要循环或 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每个物种有一个系数（变化的斜率）。