我使用pylr包运行了大约 30,000 个多元回归模型。
fit = plyr::dlply(data, "id", function(x) lm(t ~ d1 + d2 + d3, data = x))
我们有超过 30,000 个 ID，因此它推出了超过 30,000 个模型。我无法一一阅读模型。是否有任何有效的方法可以根据条件查看模型。我想根据以下情况查看模型
top 5 p-value (from the smallest)
top r squred
all the models with d3 significant
此刻，我只能用l_ply(fit, summary, .print = TRUE)逐个看模型
我尝试在broom包中使用tidy函数，以便将结果更改为清理格式，以便我可以导出.csv文件，但它不起作用。
test = broom::tidy(fit)
我收到错误消息，告诉我Error in tidy.list(fit) : No tidying method recognized for this list
任何建议将不胜感激！
以下是添加的一些数据：
df <- data.frame(id=sample(LETTERS[1:10], replace= T, 10),
d1=rnorm(20), d2=rnorm(20), d3=rnorm(20), t= rnorm(20), stringsAsFactors = F)
fit1 = plyr::dlply(df, "id", function(x) lm(t ~ d1 + d2 + d3, data = x))
下面是基于多模型一文的工作流，在嵌套数据帧上使用broom函数。请注意，我增加了数据框的大小，只是为了有更多的观测值可供使用。
使用此方法，可以按 ID 对数据进行分组，嵌套数据，并通过映射嵌套数据来创建lm对象。broom::glance为您提供模型的摘要输出：
library(tidyverse)
library(broom)
set.seed(124)
df <- data.frame(id=sample(LETTERS[1:10], size = 10, replace = T),
d1=rnorm(100), d2=rnorm(100), d3=rnorm(100), t= rnorm(100), stringsAsFactors = F)
models_df <- df %>%
group_by(id) %>%
nest() %>%
mutate(mod = map(data, function(x) lm(t ~ d1 + d2 + d3, data = x))) %>%
mutate(glance = map(mod, glance)) %>%
unnest(glance)
models_df
#> # A tibble: 6 x 14
#> id data mod r.squared adj.r.squared sigma statistic p.value df
#> <chr> <list> <lis> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 A <tibb… <S3:… 0.495 0.243 0.955 1.96 0.221 4
#> 2 E <tibb… <S3:… 0.142 -0.0189 0.907 0.883 0.471 4
#> 3 F <tibb… <S3:… 0.179 0.0247 1.05 1.16 0.355 4
#> 4 D <tibb… <S3:… 0.818 0.727 0.522 8.97 0.0123 4
#> 5 C <tibb… <S3:… 0.0514 -0.0580 0.994 0.470 0.706 4
#> 6 J <tibb… <S3:… 0.202 -0.197 0.933 0.506 0.693 4
#> # ... with 5 more variables: logLik <dbl>, AIC <dbl>, BIC <dbl>,
#> # deviance <dbl>, df.residual <int>
现在，您已经拥有了一个包含汇总统计数据(如 R 平方(的数据框，您可以使用常规dplyr命令(如top_n(来获取排名：
# smallest 5 p-value
models_df %>% top_n(-5, p.value)
#> # A tibble: 5 x 14
#> id data mod r.squared adj.r.squared sigma statistic p.value df
#> <chr> <list> <lis> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 A <tibb… <S3:… 0.495 0.243 0.955 1.96 0.221 4
#> 2 E <tibb… <S3:… 0.142 -0.0189 0.907 0.883 0.471 4
#> 3 F <tibb… <S3:… 0.179 0.0247 1.05 1.16 0.355 4
#> 4 D <tibb… <S3:… 0.818 0.727 0.522 8.97 0.0123 4
#> 5 J <tibb… <S3:… 0.202 -0.197 0.933 0.506 0.693 4
#> # ... with 5 more variables: logLik <dbl>, AIC <dbl>, BIC <dbl>,
#> # deviance <dbl>, df.residual <int>
# top r-sq
models_df %>% top_n(1, r.squared)
#> # A tibble: 1 x 14
#> id data mod r.squared adj.r.squared sigma statistic p.value df
#> <chr> <list> <lis> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 D <tibb… <S3:… 0.818 0.727 0.522 8.97 0.0123 4
#> # ... with 5 more variables: logLik <dbl>, AIC <dbl>, BIC <dbl>,
#> # deviance <dbl>, df.residual <int>
对于重要的d3项，我返回到模型数据框并在模型之间映射broom::tidy，它给出了每个项的摘要统计信息，然后过滤以研究d3项。
# significant d3--assuming alpha = 0.05
models_df %>%
mutate(tidied = map(mod, tidy)) %>%
unnest(tidied, .drop = T) %>%
filter(term == "d3") %>%
filter(p.value < 0.05)
#> # A tibble: 1 x 17
#> id r.squared adj.r.squared sigma statistic p.value df logLik AIC
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
#> 1 D 0.818 0.727 0.522 8.97 0.0123 4 -5.13 20.3
#> # ... with 8 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
#> # term <chr>, estimate <dbl>, std.error <dbl>, statistic1 <dbl>,
#> # p.value1 <dbl>
创建于 2018-06-20 由 reprex 软件包 (v0.2.0(.
下面是提取 3 个最高 R^2 模型的示例。
r2 = unlist(lapply(fit, function(x) summary(x)$r.squared))
fit[order(r2, decreasing = T)[1:3]]
其他的都是模拟的：如果你想按这些值排名，请参阅这篇文章 https://stackoverflow.com/a/5587781/2761575 了解如何提取p值