如何 运行 嵌套数据帧上的多个多元回归模型？
How to Run Several Multivariate Regression Models on Nested Dataframes?
编辑：我转而使用 mtcars 数据集，因为它更好地代表了我的真实数据的结构。忽略下面的虹膜示例：
我正在寻找 运行 多个具有多个因变量和自变量的回归模型。
假设我正在处理 iris 数据集。
我的因变量是 Petal.Length 和 Petal.Width。
我的自变量是 Sepal.Length 和 Sepal.Width.
我想运行分别为数据集中的每个物种建立以下模型（setosa、versicolor、virginica）：
lm(Petal.Length ~ Sepal.Length)
lm(Petal.Length ~ Sepal.Width)
lm(Petal.Width ~ Sepal.Length)
lm(Petal.Width ~ Sepal.Width)
我的真实数据集有更多的因变量和自变量，这就是为什么我想编写不那么重复的代码。
到目前为止，这是我的尝试，我在 Species 列中嵌套了 Tidyr，并改变了一个包含每个物种的回归数据的新列：
library(dplyr)
library(tidyr)
library(purrr)
iris %>%
nest(data = !Species) %>%
mutate(
fit = map(data, ~ lm (cbind(Petal.Length, Petal.Width) ~ ., data = .x))
)
这给出了每个物种的等价物：
lm(Petal.Length ~ Sepal.Length + Sepal.Width)
lm(Petal.Width ~ Sepal.Length + Sepal.Width)
让我们改用 mtcars 数据集：
mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
...
我的因变量是 disp 和 hp。
我的自变量是 drat、wt 和 qsec。
我想 运行 为每个 vs 和 am 分别输入以下模型：
lm(disp ~ drat)
lm(disp ~ wt)
lm(disp ~ qsec)
lm(hp ~ drat)
lm(hp ~ wt)
lm(hp ~ qsec)
我的真实数据集有更多的因变量和自变量，这就是为什么我想编写不那么重复的代码。
到目前为止，这是我的尝试，其中我用 Tidyr 嵌套了 vs 和 am 列，并改变了一个新列，其中包含每个 vs 和 am :
library(dplyr)
library(tidyr)
library(purrr)
mtcars %>%
group_by(vs,am) %>%
nest() %>%
mutate(
fit = map(data, ~ lm (cbind(disp, hp) ~ drat + wt + qsec, data = .x))
)
# A tibble: 4 x 4
# Groups: vs, am [4]
vs am data fit
<dbl> <dbl> <list> <list>
1 0 1 <tibble [6 x 9]> <mlm>
2 1 1 <tibble [7 x 9]> <mlm>
3 1 0 <tibble [7 x 9]> <mlm>
4 0 0 <tibble [12 x 9]> <mlm>
对于 vs 和 am 的每个组合给出以下公式的等价物：
lm(disp ~ drat + wt + qsec)
lm(hp ~ drat + wt + qsec)
但这不是我想要的。当我想要六个简单回归时，上面的代码为每个 vs 和 am 组合提供了两个多元回归。我不知道如何在 R 中生成它。
这是建模，在对数据进行整形后使用正确的公式，无需循环即可获得所需的结果。
new_data <- pivot_longer(iris, c(Sepal.Length, Sepal.Width),
names_to = 'Sepal', names_prefix = 'Sepal')
lm(cbind(Petal.Width, Petal.Length) ~ value/(Species:Sepal) - value +
Sepal/Species - Sepal + 0, new_data)
Call:
lm(formula = cbind(Petal.Width, Petal.Length) ~ value/(Species:Sepal) -
value + Sepal/Species - Sepal + 0, data = new_data)
Coefficients:
Petal.Width Petal.Length
Speciessetosa:Sepal.Length -0.17022 0.80305
Speciesversicolor:Sepal.Length 0.08326 0.18512
Speciesvirginica:Sepal.Length 1.22611 0.61047
Speciessetosa:Sepal.Width 0.02418 1.18292
Speciesversicolor:Sepal.Width 0.16691 1.93492
Speciesvirginica:Sepal.Width 0.66406 3.51090
value:Speciessetosa:Sepal.Length 0.08314 0.13163
value:Speciesversicolor:Sepal.Length 0.20936 0.68647
value:Speciesvirginica:Sepal.Length 0.12142 0.75008
value:Speciessetosa:Sepal.Width 0.06471 0.08141
value:Speciesversicolor:Sepal.Width 0.41845 0.83938
value:Speciesvirginica:Sepal.Width 0.45795 0.68632
如果你想运行 for 循环：
iris %>%
nest(data = !Species) %>%
summarise(model = map(c('Sepal.Length', 'Sepal.Width'),
~map(data,~lm(reformulate(.y, 'cbind(Petal.Width, Petal.Length)'),.x) %>%
coef()%>%
as_tibble(rownames = 'rn'), .y = .x)%>%
set_names(Species)%>%
bind_rows(.id='Species')))%>%
unnest(model)
# A tibble: 12 x 4
Species rn Petal.Width Petal.Length
<chr> <chr> <dbl> <dbl>
1 setosa (Intercept) -0.170 0.803
2 setosa Sepal.Length 0.0831 0.132
3 versicolor (Intercept) 0.0833 0.185
4 versicolor Sepal.Length 0.209 0.686
5 virginica (Intercept) 1.23 0.610
6 virginica Sepal.Length 0.121 0.750
7 setosa (Intercept) 0.0242 1.18
8 setosa Sepal.Width 0.0647 0.0814
9 versicolor (Intercept) 0.167 1.93
10 versicolor Sepal.Width 0.418 0.839
11 virginica (Intercept) 0.664 3.51
12 virginica Sepal.Width 0.458 0.686
编辑：
对于 mtcars 数据集，以下内容就足够了：
pivot_longer(mtcars, c(drat, wt, qsec)) %>%
unite('am_vs', c(am, vs)) %>%
lm(cbind(disp, hp)~value/(am_vs:name)-value + name/am_vs - name + 0, .)
运行 使用 map:
mtcars%>%
nest_by(vs, am)%>%
rowwise()%>%
summarise(vs, am, model = map(c('drat', 'wt', 'qsec'),
~lm(sprintf('cbind(disp, hp)~%s',.x), data)%>%
coef()%>%
as_tibble(rownames = 'rn')))%>%
unnest(model)
编辑：我转而使用 mtcars 数据集，因为它更好地代表了我的真实数据的结构。忽略下面的虹膜示例：
我正在寻找 运行 多个具有多个因变量和自变量的回归模型。
假设我正在处理 iris 数据集。
我的因变量是 Petal.Length 和 Petal.Width。
我的自变量是 Sepal.Length 和 Sepal.Width.
我想运行分别为数据集中的每个物种建立以下模型（setosa、versicolor、virginica）：
lm(Petal.Length ~ Sepal.Length)
lm(Petal.Length ~ Sepal.Width)
lm(Petal.Width ~ Sepal.Length)
lm(Petal.Width ~ Sepal.Width)
我的真实数据集有更多的因变量和自变量，这就是为什么我想编写不那么重复的代码。
到目前为止，这是我的尝试，我在 Species 列中嵌套了 Tidyr，并改变了一个包含每个物种的回归数据的新列：
library(dplyr)
library(tidyr)
library(purrr)
iris %>%
nest(data = !Species) %>%
mutate(
fit = map(data, ~ lm (cbind(Petal.Length, Petal.Width) ~ ., data = .x))
)
这给出了每个物种的等价物：
lm(Petal.Length ~ Sepal.Length + Sepal.Width)
lm(Petal.Width ~ Sepal.Length + Sepal.Width)
让我们改用 mtcars 数据集：
mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
...
我的因变量是 disp 和 hp。
我的自变量是 drat、wt 和 qsec。
我想 运行 为每个 vs 和 am 分别输入以下模型：
lm(disp ~ drat)
lm(disp ~ wt)
lm(disp ~ qsec)
lm(hp ~ drat)
lm(hp ~ wt)
lm(hp ~ qsec)
我的真实数据集有更多的因变量和自变量，这就是为什么我想编写不那么重复的代码。
到目前为止，这是我的尝试，其中我用 Tidyr 嵌套了 vs 和 am 列，并改变了一个新列，其中包含每个 vs 和 am :
library(dplyr)
library(tidyr)
library(purrr)
mtcars %>%
group_by(vs,am) %>%
nest() %>%
mutate(
fit = map(data, ~ lm (cbind(disp, hp) ~ drat + wt + qsec, data = .x))
)
# A tibble: 4 x 4
# Groups: vs, am [4]
vs am data fit
<dbl> <dbl> <list> <list>
1 0 1 <tibble [6 x 9]> <mlm>
2 1 1 <tibble [7 x 9]> <mlm>
3 1 0 <tibble [7 x 9]> <mlm>
4 0 0 <tibble [12 x 9]> <mlm>
对于 vs 和 am 的每个组合给出以下公式的等价物：
lm(disp ~ drat + wt + qsec)
lm(hp ~ drat + wt + qsec)
但这不是我想要的。当我想要六个简单回归时，上面的代码为每个 vs 和 am 组合提供了两个多元回归。我不知道如何在 R 中生成它。
这是建模，在对数据进行整形后使用正确的公式，无需循环即可获得所需的结果。
new_data <- pivot_longer(iris, c(Sepal.Length, Sepal.Width),
names_to = 'Sepal', names_prefix = 'Sepal')
lm(cbind(Petal.Width, Petal.Length) ~ value/(Species:Sepal) - value +
Sepal/Species - Sepal + 0, new_data)
Call:
lm(formula = cbind(Petal.Width, Petal.Length) ~ value/(Species:Sepal) -
value + Sepal/Species - Sepal + 0, data = new_data)
Coefficients:
Petal.Width Petal.Length
Speciessetosa:Sepal.Length -0.17022 0.80305
Speciesversicolor:Sepal.Length 0.08326 0.18512
Speciesvirginica:Sepal.Length 1.22611 0.61047
Speciessetosa:Sepal.Width 0.02418 1.18292
Speciesversicolor:Sepal.Width 0.16691 1.93492
Speciesvirginica:Sepal.Width 0.66406 3.51090
value:Speciessetosa:Sepal.Length 0.08314 0.13163
value:Speciesversicolor:Sepal.Length 0.20936 0.68647
value:Speciesvirginica:Sepal.Length 0.12142 0.75008
value:Speciessetosa:Sepal.Width 0.06471 0.08141
value:Speciesversicolor:Sepal.Width 0.41845 0.83938
value:Speciesvirginica:Sepal.Width 0.45795 0.68632
如果你想运行 for 循环：
iris %>%
nest(data = !Species) %>%
summarise(model = map(c('Sepal.Length', 'Sepal.Width'),
~map(data,~lm(reformulate(.y, 'cbind(Petal.Width, Petal.Length)'),.x) %>%
coef()%>%
as_tibble(rownames = 'rn'), .y = .x)%>%
set_names(Species)%>%
bind_rows(.id='Species')))%>%
unnest(model)
# A tibble: 12 x 4
Species rn Petal.Width Petal.Length
<chr> <chr> <dbl> <dbl>
1 setosa (Intercept) -0.170 0.803
2 setosa Sepal.Length 0.0831 0.132
3 versicolor (Intercept) 0.0833 0.185
4 versicolor Sepal.Length 0.209 0.686
5 virginica (Intercept) 1.23 0.610
6 virginica Sepal.Length 0.121 0.750
7 setosa (Intercept) 0.0242 1.18
8 setosa Sepal.Width 0.0647 0.0814
9 versicolor (Intercept) 0.167 1.93
10 versicolor Sepal.Width 0.418 0.839
11 virginica (Intercept) 0.664 3.51
12 virginica Sepal.Width 0.458 0.686
编辑： 对于 mtcars 数据集，以下内容就足够了：
pivot_longer(mtcars, c(drat, wt, qsec)) %>%
unite('am_vs', c(am, vs)) %>%
lm(cbind(disp, hp)~value/(am_vs:name)-value + name/am_vs - name + 0, .)
运行 使用 map:
mtcars%>%
nest_by(vs, am)%>%
rowwise()%>%
summarise(vs, am, model = map(c('drat', 'wt', 'qsec'),
~lm(sprintf('cbind(disp, hp)~%s',.x), data)%>%
coef()%>%
as_tibble(rownames = 'rn')))%>%
unnest(model)