logo标签列表

使用dplyr / broom创建模型和增强数据时不丢失额外列

rdplyrtidyversebroom
10
考虑以下数据/示例。每个数据集包含多个样本，每个样本都有一个观察值和一个估计值：
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

如何保留原始数据集的额外列？

我看过这个答案，其中使用了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

...然后与原始数据进行左连接，然后对每个estimate进行mutate以获得线性调整值，但这看起来过于复杂。

我发现的另一种丑陋的方法是将列作为虚拟变量添加到模型中：

data %>%
  group_by(dataset) %>% 
  do(lm(observation ~ estimate + 0 * sample_id, data = .) %>% augment)

有没有更简单（简洁）的解决方案，不需要手动指定我想要保留的变量？
1
另一种hack-ish方法是将数据join()回所得到的数据框中。如果将do(lm(observation ~ estimate, data = .) %>% augment)的结果管道传输到left_join(data %>% select(dataset, observation, sample_id), by=c('dataset', 'observation'))中，你就可以得到想要的结果。但这取决于你的数据观测值是否唯一（每个数据集）。 - Curt F.
我考虑过这个问题，但不幸的是 - 这些观察结果并不唯一。不过，在 do(augment) 调用之后，我想我可以指望所有内容都按照相同的顺序排列？ - slhck
看我的回答。添加一个行号虚拟列就可以解决唯一性问题了... - Curt F.
4个回答
6
您可以使用broom :: augment_columns代替augment。我们需要的函数的两个参数是x——“模型”——和data——“应添加列的原始数据”。
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")
3
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()的结果上，恢复任意数量的列而无需手动指定。

我认为这比其他解决方案更易读，但仍然有点hackish。从left_join中去除重复列有点繁琐。您可能不希望在输出中出现像observation.x和estimate.y这样的列，除非您执行了select(-observation, -estimate)部分。

我在4年后回到这个答案，发现使用broom::tidy提取的线性拟合项（当使用0*row技巧时）不再包括截距。有什么办法可以避免这种情况吗？ - slhck
1
这基本上与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)
0
怎么样，这个怎么处理？
DF %>%
  group_by(dataset) %>% 
  do(cbind(sample_id = .$sample_id, lm(observation ~ estimate, data = .) %>% augment))

还是太丑了吗？

1
这不是保留多列的通用解决方案。 - David Arenburg
是的，就像@DavidArenburg所说的那样，这将是我使用这种方法的主要问题。（虽然公平地说，问题中从未指定第二个变量。） - slhck
重新思考一下，我们可以使用cbind函数将原始数据框中除了observation和estimate列之外的所有内容进行合并。 - slhck
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接