我正在按照 R for Data Science 一书中的多模型示例生成单个多元线性回归模型。 如下面的可重现示例所示，我的数据框由三列组成：id, val1, val2, val3,. 我能够使用哈德利书中详述的map函数拟合线性模型。 但是，我一直在努力从每个模型中提取系数并将这些值作为列添加到my.list对象的数据帧中。 模型系数当前存储的方式使得调用代码的其他部分变得困难/繁琐。
到目前为止，我想出的最好的方法是通过迭代my.list中的每个数据帧来列出my.list长度并提取系数：Name1, Name2, Name3。 这意味着现在我的全局环境中有另一个列表，并且coef.list不再包含来自my.listName1, Name2, Name3的因素;这些现在已替换为[[1]], [[2]], [[3]].
在使用多个模型时，谁能建议一种"更干净"的方法来提取模型系数？ 我的首选输出只是为每个系数创建一个列：intercept, val1, val2. 这些列将位于现有数据框中Name1, Name2, Name3my.list内，以便我可以直接在数据框中使用mutate：
# reproducible example
set.seed(1363)
d1 <- data.frame(id=c("Name1", "Name2", "Name3"),
val1=c(rnorm(n=15, mean=5)),
val2=c(rnorm(n=15, mean=3)),
val3=c(rnorm(n=15, mean=8)))
# linear model function
lm.fun <- function(df){
lm(val3 ~ val1+val2, data = df)
}
# map lm function
d1 <- d1 %>%
group_by(id) %>%
nest() %>%
mutate(model = map(data, lm.fun)) %>%
unnest(data, .drop = FALSE)
#split data frame by 'id' and convert into list
my.list <- split(d1, d1$id)
# make list of coefficients
coef.list <- list(length(my.list))
for (i in seq_along(my.list)) {
coef.list[[i]] <- my.list[[i]][["model"]][[1]][["coefficients"]]
}
>head(coef.list, n=1)
[[1]]
(Intercept) val1 val2
9.03278337 -0.07096932 0.02119088
期望的输出
my.list$Name1
id val1 val2 val3 intc coef1 coef2
Name1 1 2 3 9.03 -.070 .021
Name1 3 1 5 9.03 -.070 .021
Name1 2 6 8 9.03 -.070 .021
根据描述，如果我们需要创建一些带有系数的列，一种选择是使用do与group_by. 按"id"分组后，将系数提取为do内的list并rename列(如果需要(
library(tidyverse)
d1 %>%
group_by(id) %>%
do(data.frame(., as.list(coef(lm(val3 ~ val1 + val2, data = .))))) %>%
rename_at(5:7, ~c("intc", "coef1", "coef2"))
# A tibble: 15 x 7
# Groups: id [3]
# id val1 val2 val3 intc coef1 coef2
# <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 Name1 6.76 2.64 9.85 9.03 -0.0710 0.0212
# 2 Name1 6.78 1.52 6.94 9.03 -0.0710 0.0212
# 3 Name1 4.14 4.31 8.30 9.03 -0.0710 0.0212
# 4 Name1 5.55 2.16 9.97 9.03 -0.0710 0.0212
# 5 Name1 6.18 3.64 8.32 9.03 -0.0710 0.0212
# 6 Name2 6.08 1.12 9.96 7.33 0.468 -0.488
# 7 Name2 3.54 4.71 6.43 7.33 0.468 -0.488
# 8 Name2 5.66 4.54 8.69 7.33 0.468 -0.488
# 9 Name2 6.88 4.15 7.79 7.33 0.468 -0.488
#10 Name2 4.89 1.27 8.72 7.33 0.468 -0.488
#11 Name3 6.41 4.38 6.22 20.1 -2.15 0.118
#12 Name3 5.06 3.28 9.42 20.1 -2.15 0.118
#13 Name3 6.25 3.16 8.15 20.1 -2.15 0.118
#14 Name3 6.03 3.63 7.78 20.1 -2.15 0.118
#15 Name3 6.51 1.46 5.90 20.1 -2.15 0.118
或者我们可以使用broom包函数。tidy、glance提取感兴趣的列(以及更多(等函数。nest分组数据集后，构建模型，使用tidy提取组件，select感兴趣的列并unnest
library(broom)
d1 %>%
group_by(id) %>%
nest %>%
mutate(model1 = map(data, ~
lm(val3 ~ val1 + val2, data = .) %>%
tidy %>%
dplyr::select(term, estimate) %>%
spread(term, estimate) %>%
rename_all(~ c("intc", paste0("coef", 1:2))))) %>%
unnest(model1) %>%
unnest(data)
# A tibble: 15 x 7
# id intc coef1 coef2 val1 val2 val3
# <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 Name1 9.03 -0.0710 0.0212 6.76 2.64 9.85
# 2 Name1 9.03 -0.0710 0.0212 6.78 1.52 6.94
# 3 Name1 9.03 -0.0710 0.0212 4.14 4.31 8.30
# 4 Name1 9.03 -0.0710 0.0212 5.55 2.16 9.97
# 5 Name1 9.03 -0.0710 0.0212 6.18 3.64 8.32
# 6 Name2 7.33 0.468 -0.488 6.08 1.12 9.96
# 7 Name2 7.33 0.468 -0.488 3.54 4.71 6.43
# 8 Name2 7.33 0.468 -0.488 5.66 4.54 8.69
# 9 Name2 7.33 0.468 -0.488 6.88 4.15 7.79
#10 Name2 7.33 0.468 -0.488 4.89 1.27 8.72
#11 Name3 20.1 -2.15 0.118 6.41 4.38 6.22
#12 Name3 20.1 -2.15 0.118 5.06 3.28 9.42
#13 Name3 20.1 -2.15 0.118 6.25 3.16 8.15
#14 Name3 20.1 -2.15 0.118 6.03 3.63 7.78
#15 Name3 20.1 -2.15 0.118 6.51 1.46 5.90
或者不使用broom
d1 %>%
group_by(id) %>%
nest %>%
mutate(model = map(data, ~ lm(val3 ~ val1 + val2, data = .) %>%
coef %>%
as.list %>%
as_tibble)) %>%
unnest(model) %>%
unnest(data)
如果我们只需要每个"id"的摘要输出
d1 %>%
group_by(id) %>%
nest %>%
mutate(model1 = map(data, ~
lm(val3 ~ val1 + val2, data = .) %>%
tidy %>%
dplyr::select(term, estimate) %>%
spread(term, estimate) %>%
rename_all(~ c("intc", paste0("coef", 1:2))))) %>%
dplyr::select(-data) %>%
unnest
或者使用data.table，我们可以更简洁地做到这一点，在按"id"分组并通过提取模型的coef作为list来分配(:=(新列
library(data.table)
setDT(d1)[, c('intc', 'coef1', 'coef2') :=
as.list(coef(lm(val3 ~ val1 + val2))), id]
d1[order(id)]
# id val1 val2 val3 intc coef1 coef2
# 1: Name1 6.755964 2.642874 9.849828 9.032783 -0.07096932 0.02119088
# 2: Name1 6.776666 1.522431 6.937053 9.032783 -0.07096932 0.02119088
# 3: Name1 4.141883 4.307537 8.301940 9.032783 -0.07096932 0.02119088
# 4: Name1 5.551850 2.163882 9.971588 9.032783 -0.07096932 0.02119088
# 5: Name1 6.179506 3.635832 8.319042 9.032783 -0.07096932 0.02119088
# 6: Name2 6.083243 1.116293 9.960934 7.325156 0.46840770 -0.48806159
# 7: Name2 3.536476 4.708967 6.427627 7.325156 0.46840770 -0.48806159
# 8: Name2 5.663909 4.541081 8.691523 7.325156 0.46840770 -0.48806159
# 9: Name2 6.883746 4.150780 7.791050 7.325156 0.46840770 -0.48806159
#10: Name2 4.890291 1.269559 8.723792 7.325156 0.46840770 -0.48806159
#11: Name3 6.414915 4.383609 6.220188 20.106581 -2.14601530 0.11770877
#12: Name3 5.059774 3.276510 9.421862 20.106581 -2.14601530 0.11770877
#13: Name3 6.251416 3.157157 8.147720 20.106581 -2.14601530 0.11770877
#14: Name3 6.028100 3.630858 7.783118 20.106581 -2.14601530 0.11770877
#15: Name3 6.505153 1.460564 5.895564 20.106581 -2.14601530 0.11770877
或者，如果我们不想更新初始数据，请执行join
setDT(d1)[, as.list(coef(lm(val3 ~ val1 + val2))), id][d1, on = .(id)]
注意："d1"是初始数据集
使用sapply，给定数据d1(包含model列(：
coeff <- sapply(d1$model, function(x) return(coefficients(x)))
library(dplyr)
bind_cols(d1, data.frame(t(coeff))) %>%
rename_at(6:8, ~ c("intc", "coef1", "coef2")) %>%
distinct(id, .keep_all = TRUE)
library(tidyverse)
d1%>%
group_by(id)%>%
mutate(s=toString(coef(lm(val3~val1+val2))))%>%
separate(s,c("intercept","coef1","coef2"),sep=",",convert = T)%>%
arrange(id)
# A tibble: 15 x 7
# Groups: id [3]
id val1 val2 val3 intercept coef1 coef2
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Name1 6.76 2.64 9.85 9.03 -0.0710 0.0212
2 Name1 6.78 1.52 6.94 9.03 -0.0710 0.0212
3 Name1 4.14 4.31 8.30 9.03 -0.0710 0.0212
4 Name1 5.55 2.16 9.97 9.03 -0.0710 0.0212
5 Name1 6.18 3.64 8.32 9.03 -0.0710 0.0212
6 Name2 6.08 1.12 9.96 7.33 0.468 -0.488
7 Name2 3.54 4.71 6.43 7.33 0.468 -0.488
8 Name2 5.66 4.54 8.69 7.33 0.468 -0.488
9 Name2 6.88 4.15 7.79 7.33 0.468 -0.488
10 Name2 4.89 1.27 8.72 7.33 0.468 -0.488
11 Name3 6.41 4.38 6.22 20.1 -2.15 0.118
12 Name3 5.06 3.28 9.42 20.1 -2.15 0.118
13 Name3 6.25 3.16 8.15 20.1 -2.15 0.118
14 Name3 6.03 3.63 7.78 20.1 -2.15 0.118
15 Name3 6.51 1.46 5.90 20.1 -2.15 0.118