在不丢失 dplyr/broom 中的额外列的情况下创建模型和扩充数据
Creating models and augmenting data without losing additional columns in dplyr/broom
考虑以下数据/示例。每个数据集包含大量样本，其中包含一项观察值和一项估计值：
library(tidyverse)
library(broom)
data = read.table(text = '
dataset sample_id observation estimate
A A1 4.8 4.7
A A2 4.3 4.5
A A3 3.1 2.9
A A4 2.1 2
A A5 1.1 1
B B1 4.5 4.3
B B2 3.9 4.1
B B3 2.9 3
B B4 1.8 2
B B5 1 1.2
', header = TRUE)
我想为每个数据集计算一个线性模型以消除观察和估计之间的任何线性偏差，并获得与原始值相邻的拟合值：
data %>%
group_by(dataset) %>%
do(lm(observation ~ estimate, data = .) %>% augment)
然而，这样做的目的是删除 sample_id 列，我需要保留该列以便继续基于该唯一 ID 使用此数据集进行计算：
# A tibble: 10 x 10
# Groups: dataset [2]
dataset observation estimate .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 A 4.80 4.70 4.68 0.107 0.115 0.478 0.152 0.491 1.04
2 A 4.30 4.50 4.49 0.0996 -0.193 0.416 0.0609 0.957 -1.64
3 A 3.10 2.90 2.97 0.0693 0.135 0.201 0.156 0.120 0.976
4 A 2.10 2.00 2.11 0.0849 -0.00583 0.303 0.189 0.000444 -0.0452
5 A 1.10 1.00 1.15 0.120 -0.0508 0.602 0.180 0.206 -0.521
6 B 4.50 4.30 4.31 0.109 0.191 0.468 0.0597 1.20 1.65
7 B 3.90 4.10 4.09 0.100 -0.193 0.396 0.0844 0.798 -1.56
8 B 2.90 3.00 2.91 0.0713 -0.00630 0.201 0.195 0.000247 -0.0443
9 B 1.80 2.00 1.83 0.0898 -0.0275 0.319 0.193 0.0103 -0.210
10 B 1.00 1.20 0.964 0.125 0.0355 0.616 0.191 0.104 0.360
如何保留原始数据集中的附加列？
我以前看过this answer用nest来折叠数据，但我仍然只用这种方法得到模型参数。我想我可以提取每个数据集的参数：
data %>%
group_by(dataset) %>%
nest() %>%
mutate(
mod = map(data, linear_adj_model),
pars = map(mod, tidy)
) %>%
unnest(pars) %>%
select(dataset, term, estimate) %>%
spread(term, estimate)
... 这给了我这个：
# A tibble: 2 x 3
dataset `(Intercept)` estimate
* <fct> <dbl> <dbl>
1 A 0.196 0.955
2 B -0.330 1.08
… 然后将其与原始数据左连接，然后 mutate 每个 estimate 得到线性调整的数据，但这似乎太复杂了。
我发现的另一个丑陋的技巧是将列作为虚拟变量添加到模型中：
data %>%
group_by(dataset) %>%
do(lm(observation ~ estimate + 0 * sample_id, data = .) %>% augment)
是否有更简单（整洁）的解决方案，无需手动指定我要保留的变量？
这个怎么样：
DF %>%
group_by(dataset) %>%
do(cbind(sample_id = .$sample_id, lm(observation ~ estimate, data = .) %>% augment))
还是太丑了？
您可以使用 broom::augment_columns 而不是 augment。我们需要的函数的两个参数是 x -- "a model" -- 和 data -- "original data onto which columns should be added".
library(tidyverse)
library(broom)
split(data, data$dataset) %>%
map(., ~lm(formula = observation ~ estimate, data = .)) %>%
map2(.x = ., .y = split(data, f = data$dataset), .f = ~augment_columns(x = .x, data = .y)) %>%
bind_rows() %>%
select(-.rownames)
# dataset sample_id observation estimate .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
#1 A A1 4.8 4.7 4.6845093 0.10675590 0.115490737 0.4781238 0.15157780 0.4911635931 1.03547990
#2 A A2 4.3 4.5 4.4934963 0.09956065 -0.193496255 0.4158455 0.06089193 0.9570799385 -1.63978525
#3 A A3 3.1 2.9 2.9653922 0.06929022 0.134607804 0.2014190 0.15623754 0.1200409795 0.97563873
#4 A A4 2.1 2.0 2.1058337 0.08491818 -0.005833662 0.3025227 0.18902495 0.0004439221 -0.04524332
#5 A A5 1.1 1.0 1.1507686 0.11979870 -0.050768624 0.6020891 0.18032220 0.2055920869 -0.52129162
#6 B B1 4.5 4.3 4.3087226 0.10879087 0.191277434 0.4679235 0.05965705 1.1954021471 1.64881395
#7 B B2 3.9 4.1 4.0929657 0.10006757 -0.192965672 0.3958920 0.08438937 0.7984863377 -1.56105324
#8 B B3 2.9 3.0 2.9063028 0.07128455 -0.006302757 0.2009004 0.19471901 0.0002470587 -0.04433279
#9 B B4 1.8 2.0 1.8275183 0.08983650 -0.027518289 0.3190771 0.19335019 0.0103015495 -0.20968503
#10 B B5 1.0 1.2 0.9644907 0.12484420 0.035509285 0.6162071 0.19051943 0.1042741368 0.36040302
想法是 split 数据集的数据，将模型拟合到列表的每个组件，然后使用 map2 迭代模型和用于模型的（完整）数据构建，即并行 split(data, f = data$dataset) 的结果。
augment_columns 添加了一个 .rownames 列，因此在最后一行添加了 select。
编辑
相同的解决方案，但希望更容易阅读。
data_split <- split(data, data$dataset)
models <- map(data_split, ~lm(formula = observation ~ estimate, data = .))
map2(.x = models, .y = data_split, .f = ~augment_columns(x = .x, data = .y)) %>%
bind_rows() %>%
select(-.rownames)
作为具有四个参数的函数的第一个代码块：df、split_var、dependend_var 和 explanatory_var。
augment_df <- function(df, split_var, dependend_var, explanatory_var) {
require(tidyverse)
require(broom)
split(df, df[split_var]) %>%
map(., ~lm(formula = as.formula(paste0(dependend_var, " ~ ", explanatory_var)), data = .)) %>%
map2(.x = ., .y = split(df, df[split_var]), .f = ~augment_columns(x = .x, data = .y)) %>%
bind_rows() %>%
select(-.rownames)
}
augment_df(df = data, split_var = "dataset", dependend_var = "observation", explanatory_var = "estimate")
numbered_data <-
data %>%
mutate(row = row_number())
numbered_data %>%
group_by(dataset) %>%
do(augment(lm(observation ~ estimate + 0*row, data = .))) %>%
left_join(numbered_data %>% select(-observation, -estimate), by=c('dataset', 'row')) %>%
select(-row)
这种方法也依赖于在模型中使用虚拟变量，但它以与列无关的方式进行。通过定义新的虚拟列 row，您可以 left_join() 将原始数据返回到 augment() 的结果中，无需手动指定即可恢复任意数量的列。
我发现它比其他解决方案更具可读性，但它仍然有点老套。从 left_join 中删除重复的列有点乏味。您可能不希望在输出中包含 observation.x 和 estimate.y 之类的列，除非您执行 select(-observation, -estimate) 部分，否则您将拥有这些列。
这与 Markus 的回答基本相同，但可能更清晰一些。
library(tidyverse)
library(broom)
data = read.table(text = '
dataset sample_id observation estimate
A A1 4.8 4.7
A A2 4.3 4.5
A A3 3.1 2.9
A A4 2.1 2
A A5 1.1 1
B B1 4.5 4.3
B B2 3.9 4.1
B B3 2.9 3
B B4 1.8 2
B B5 1 1.2
', header = TRUE)
data %>%
group_by(dataset) %>%
nest() %>%
mutate(mod = map(data, ~lm(observation ~ estimate, data = .)),
aug = map2(mod, data, ~augment_columns(.x, .y))) %>%
unnest(aug)
考虑以下数据/示例。每个数据集包含大量样本，其中包含一项观察值和一项估计值：
library(tidyverse)
library(broom)
data = read.table(text = '
dataset sample_id observation estimate
A A1 4.8 4.7
A A2 4.3 4.5
A A3 3.1 2.9
A A4 2.1 2
A A5 1.1 1
B B1 4.5 4.3
B B2 3.9 4.1
B B3 2.9 3
B B4 1.8 2
B B5 1 1.2
', header = TRUE)
我想为每个数据集计算一个线性模型以消除观察和估计之间的任何线性偏差，并获得与原始值相邻的拟合值：
data %>%
group_by(dataset) %>%
do(lm(observation ~ estimate, data = .) %>% augment)
然而，这样做的目的是删除 sample_id 列，我需要保留该列以便继续基于该唯一 ID 使用此数据集进行计算：
# A tibble: 10 x 10
# Groups: dataset [2]
dataset observation estimate .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 A 4.80 4.70 4.68 0.107 0.115 0.478 0.152 0.491 1.04
2 A 4.30 4.50 4.49 0.0996 -0.193 0.416 0.0609 0.957 -1.64
3 A 3.10 2.90 2.97 0.0693 0.135 0.201 0.156 0.120 0.976
4 A 2.10 2.00 2.11 0.0849 -0.00583 0.303 0.189 0.000444 -0.0452
5 A 1.10 1.00 1.15 0.120 -0.0508 0.602 0.180 0.206 -0.521
6 B 4.50 4.30 4.31 0.109 0.191 0.468 0.0597 1.20 1.65
7 B 3.90 4.10 4.09 0.100 -0.193 0.396 0.0844 0.798 -1.56
8 B 2.90 3.00 2.91 0.0713 -0.00630 0.201 0.195 0.000247 -0.0443
9 B 1.80 2.00 1.83 0.0898 -0.0275 0.319 0.193 0.0103 -0.210
10 B 1.00 1.20 0.964 0.125 0.0355 0.616 0.191 0.104 0.360
如何保留原始数据集中的附加列？
我以前看过this answer用nest来折叠数据，但我仍然只用这种方法得到模型参数。我想我可以提取每个数据集的参数：
data %>%
group_by(dataset) %>%
nest() %>%
mutate(
mod = map(data, linear_adj_model),
pars = map(mod, tidy)
) %>%
unnest(pars) %>%
select(dataset, term, estimate) %>%
spread(term, estimate)
... 这给了我这个：
# A tibble: 2 x 3
dataset `(Intercept)` estimate
* <fct> <dbl> <dbl>
1 A 0.196 0.955
2 B -0.330 1.08
… 然后将其与原始数据左连接，然后 mutate 每个 estimate 得到线性调整的数据，但这似乎太复杂了。
我发现的另一个丑陋的技巧是将列作为虚拟变量添加到模型中：
data %>%
group_by(dataset) %>%
do(lm(observation ~ estimate + 0 * sample_id, data = .) %>% augment)
是否有更简单（整洁）的解决方案，无需手动指定我要保留的变量？
这个怎么样：
DF %>%
group_by(dataset) %>%
do(cbind(sample_id = .$sample_id, lm(observation ~ estimate, data = .) %>% augment))
还是太丑了？
您可以使用 broom::augment_columns 而不是 augment。我们需要的函数的两个参数是 x -- "a model" -- 和 data -- "original data onto which columns should be added".
library(tidyverse)
library(broom)
split(data, data$dataset) %>%
map(., ~lm(formula = observation ~ estimate, data = .)) %>%
map2(.x = ., .y = split(data, f = data$dataset), .f = ~augment_columns(x = .x, data = .y)) %>%
bind_rows() %>%
select(-.rownames)
# dataset sample_id observation estimate .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
#1 A A1 4.8 4.7 4.6845093 0.10675590 0.115490737 0.4781238 0.15157780 0.4911635931 1.03547990
#2 A A2 4.3 4.5 4.4934963 0.09956065 -0.193496255 0.4158455 0.06089193 0.9570799385 -1.63978525
#3 A A3 3.1 2.9 2.9653922 0.06929022 0.134607804 0.2014190 0.15623754 0.1200409795 0.97563873
#4 A A4 2.1 2.0 2.1058337 0.08491818 -0.005833662 0.3025227 0.18902495 0.0004439221 -0.04524332
#5 A A5 1.1 1.0 1.1507686 0.11979870 -0.050768624 0.6020891 0.18032220 0.2055920869 -0.52129162
#6 B B1 4.5 4.3 4.3087226 0.10879087 0.191277434 0.4679235 0.05965705 1.1954021471 1.64881395
#7 B B2 3.9 4.1 4.0929657 0.10006757 -0.192965672 0.3958920 0.08438937 0.7984863377 -1.56105324
#8 B B3 2.9 3.0 2.9063028 0.07128455 -0.006302757 0.2009004 0.19471901 0.0002470587 -0.04433279
#9 B B4 1.8 2.0 1.8275183 0.08983650 -0.027518289 0.3190771 0.19335019 0.0103015495 -0.20968503
#10 B B5 1.0 1.2 0.9644907 0.12484420 0.035509285 0.6162071 0.19051943 0.1042741368 0.36040302
想法是 split 数据集的数据，将模型拟合到列表的每个组件，然后使用 map2 迭代模型和用于模型的（完整）数据构建，即并行 split(data, f = data$dataset) 的结果。
augment_columns 添加了一个 .rownames 列，因此在最后一行添加了 select。
编辑
相同的解决方案，但希望更容易阅读。
data_split <- split(data, data$dataset)
models <- map(data_split, ~lm(formula = observation ~ estimate, data = .))
map2(.x = models, .y = data_split, .f = ~augment_columns(x = .x, data = .y)) %>%
bind_rows() %>%
select(-.rownames)
作为具有四个参数的函数的第一个代码块：df、split_var、dependend_var 和 explanatory_var。
augment_df <- function(df, split_var, dependend_var, explanatory_var) {
require(tidyverse)
require(broom)
split(df, df[split_var]) %>%
map(., ~lm(formula = as.formula(paste0(dependend_var, " ~ ", explanatory_var)), data = .)) %>%
map2(.x = ., .y = split(df, df[split_var]), .f = ~augment_columns(x = .x, data = .y)) %>%
bind_rows() %>%
select(-.rownames)
}
augment_df(df = data, split_var = "dataset", dependend_var = "observation", explanatory_var = "estimate")
numbered_data <-
data %>%
mutate(row = row_number())
numbered_data %>%
group_by(dataset) %>%
do(augment(lm(observation ~ estimate + 0*row, data = .))) %>%
left_join(numbered_data %>% select(-observation, -estimate), by=c('dataset', 'row')) %>%
select(-row)
这种方法也依赖于在模型中使用虚拟变量，但它以与列无关的方式进行。通过定义新的虚拟列 row，您可以 left_join() 将原始数据返回到 augment() 的结果中，无需手动指定即可恢复任意数量的列。
我发现它比其他解决方案更具可读性，但它仍然有点老套。从 left_join 中删除重复的列有点乏味。您可能不希望在输出中包含 observation.x 和 estimate.y 之类的列，除非您执行 select(-observation, -estimate) 部分，否则您将拥有这些列。
这与 Markus 的回答基本相同，但可能更清晰一些。
library(tidyverse)
library(broom)
data = read.table(text = '
dataset sample_id observation estimate
A A1 4.8 4.7
A A2 4.3 4.5
A A3 3.1 2.9
A A4 2.1 2
A A5 1.1 1
B B1 4.5 4.3
B B2 3.9 4.1
B B3 2.9 3
B B4 1.8 2
B B5 1 1.2
', header = TRUE)
data %>%
group_by(dataset) %>%
nest() %>%
mutate(mod = map(data, ~lm(observation ~ estimate, data = .)),
aug = map2(mod, data, ~augment_columns(.x, .y))) %>%
unnest(aug)