时间戳

首页 > 解决方案 > R split/map from purrr 结合 tidy from broom 按组获取线性模型统计

r - R split/map from purrr 结合 tidy from broom 按组获取线性模型统计

问题描述

我在 tidyverse.org 的这个链接上在线找到了这段代码：

mtcars %>%
  split(.$cyl) %>%
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%
  map_dbl("r.squared")

在此处输入图像描述

该代码按预期工作。我现在正在使用相同的结构进行练习，但使用的是长数据框。可以看代码；大部分是一样的。首先，我转换为 tibble，为汽车添加行名，选择数值变量，并使数据框成为长数据框。

mtcars <- as_tibble(mtcars, rownames = 'car')

mtcars_numeric <- mtcars %>%
  select(car, mpg, disp, hp, drat, wt, qsec) 

mtcars_long_numeric <- pivot_longer(mtcars_numeric, names_to = 'names', values_to = 'values', 3:7)


mtcars_long_numeric %>%
  split(.$names) %>%
  map(~ lm(mpg ~ values, data = .)) %>%
  map(summary) %>%
  map_df("r.squared") %>%
  pivot_longer(., names_to = 'explanatory_variable_to_mpg', values_to = 'r_squared', 1:5) %>%
  arrange(desc(r_squared))

在此处输入图像描述

但是其他模型统计数据，比如 p 值呢？我该如何提取它？如果我只是将“r.squared”更改为“p.value”，它就不起作用。我尝试了其他变体，例如“p_value”和“pvalue”，但它不起作用。我也不知道如何为这些对象找到正确的名称。

在此处输入图像描述

我可以创建一个线性模型对象并查看摘要中的 r.squared 并获得正确的值。

mtcars_linear_model <- lm(mpg ~ wt, mtcars)

summary(mtcars_linear_model)$r.squared

在此处输入图像描述

...但是在这个小插曲之外，我不知道我怎么会知道 r.squared 存在于线性模型的摘要中。如果我只是在 summary(lm) 之后键入美元符号，我会得到不存在的值。（这是一个错误吗？）

在此处输入图像描述

然后我尝试了不同的策略。我可以看到，如果我使用扫帚并整理线性模型对象，我还有其他统计信息：

broom::tidy(mtcars_linear_model)

在此处输入图像描述

有没有办法将 broom::tidy 函数添加到这些涉及 purrr:map 的数据帧中？目的是弄清楚如何提取其他模型统计数据，如 p 值。另外，我如何找到可以从线性模型对象摘要（lm）$'的摘要中提取的项目的综合列表？

以下代码不起作用。我尝试了一些变体，例如 %>% tidy() 或者像这样将 tidy 包裹在 map(summary) 周围： tidy(map(summary)) 但它不起作用。

mtcars_long_numeric %>%
  split(.$names) %>%
  map(~ lm(mpg ~ values, data = .)) %>%
  map(summary) %>%
  tidy() %>% #### ????????
  map_df("r.squared") %>%
  pivot_longer(., names_to = 'explanatory_variable_to_mpg', values_to = 'r_squared', 1:5) %>%
  arrange(desc(r_squared))

标签： rlinear-regressionpurrrtidybroom

解决方案

这个？。您需要使用glance而不是tidy模型统计信息。

mtcars_long_numeric %>% 
  nest_by(names) %>% 
  mutate(model = list(lm(mpg ~ values, data = data))) %>% 
  summarise(glance(model))
`summarise()` has grouped output by 'names'. You can override using the `.groups` argument.
# A tibble: 5 × 13
# Groups:   names [5]
  names r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC deviance df.residual  nobs
  <chr>     <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1 disp      0.718         0.709  3.25     76.5  9.38e-10     1  -82.1  170.  175.     317.          30    32
2 drat      0.464         0.446  4.49     26.0  1.78e- 5     1  -92.4  191.  195.     604.          30    32
3 hp        0.602         0.589  3.86     45.5  1.79e- 7     1  -87.6  181.  186.     448.          30    32
4 qsec      0.175         0.148  5.56      6.38 1.71e- 2     1  -99.3  205.  209.     929.          30    32
5 wt        0.753         0.745  3.05     91.4  1.29e-10     1  -80.0  166.  170.     278.          30    32

推荐阅读