【问题标题】:Create tibble with all possible models使用所有可能的模型创建 tibble
【发布时间】:2017-12-09 18:03:26
【问题描述】:

我创建了一个循环,它为给定的数据集创建所有可能的模型组合。有 63 个可能的模型,我需要将它们放入一个带有模型编号、解释变量子集、模型公式和结果(特别是 r 平方值)的小标题中。

Cols <- names(finalprojectdata3)
Cols <- Cols[! Cols %in% 'debt']
n <- length(Cols)

id <- unlist(
  lapply(1:n,
         function(i)combn(1:n, i, simplify = FALSE)
         ),
  recursive = FALSE)

Formulas <- sapply(id, function(i)
  paste('debt~', paste(Cols[i],collapse="+")))

models <- lapply(Formulas, function(i)
summary(lm(as.formula(i), data = finalprojectdata3)))

models

输出是每个模型的摘要,但我需要它在一个易于阅读的 tibble 中。

【问题讨论】:

  • 也许你应该试试broom。 (Introduction to broom 小插图。)
  • 我一直在尝试,但似乎没有任何效果。
  • 您在这里有一个大部分可重复的问题,但了解Cols 的样子可能会有所帮助(我们真的不需要知道finalprojectdata3 的样子,因为您只是做的是创建一个因子设计)。 (这也可能受益于tidyr。)
  • Cols is "distance" "scholarship" "parents" "car" "nhouse" "major2" 有没有办法上传文件?我可以这样做..
  • 制作一个玩具示例,或者在您的数据集上使用dput(),或者(如果非常大)在您的数据集的一部分上使用dput()。见这里:stackoverflow.com/questions/5963269/…

标签: r


【解决方案1】:

目前还不清楚您希望如何输出,但这里有一个建议的路径,遵循tidyr 嵌套对象。

按照上述方法制作一些虚假数据:

dat <- mtcars[,1:5]
Cols <- names(dat)
Cols <- Cols[! Cols %in% 'mpg']
n <- length(Cols)

id <- unlist(
  lapply(1:n,
         function(i)combn(1:n, i, simplify = FALSE)
         ),
  recursive = FALSE)
str(id)
# List of 15
#  $ : int 1
#  $ : int 2
#  $ : int 3
#  $ : int 4
#  $ : int [1:2] 1 2
#  $ : int [1:2] 1 3
#  $ : int [1:2] 1 4
#  $ : int [1:2] 2 3
#  $ : int [1:2] 2 4
#  $ : int [1:2] 3 4
#  $ : int [1:3] 1 2 3
#  $ : int [1:3] 1 2 4
#  $ : int [1:3] 1 3 4
#  $ : int [1:3] 2 3 4
#  $ : int [1:4] 1 2 3 4

Formulas <- sapply(id, function(i)
  paste('mpg ~', paste(Cols[i], collapse=" + ")))
head(Formulas)
# [1] "mpg ~ cyl"        "mpg ~ disp"       "mpg ~ hp"         "mpg ~ drat"      
# [5] "mpg ~ cyl + disp" "mpg ~ cyl + hp"  

这就是我与你的道路不同的地方。

library(dplyr)
library(tidyr)
library(purrr)

x <- data_frame(Formulas) %>%
  mutate(
    lms = map(Formulas, ~ lm(as.formula(.), data = dat)),
    summaries = map(lms, ~ summary(.)),
    coefs = map(summaries, ~ as.data.frame(coef(.)))
  )
x
# # A tibble: 15 × 4
#                        Formulas      lms        summaries                coefs
#                           <chr>   <list>           <list>               <list>
# 1                     mpg ~ cyl <S3: lm> <S3: summary.lm> <data.frame [2 × 4]>
# 2                    mpg ~ disp <S3: lm> <S3: summary.lm> <data.frame [2 × 4]>
# 3                      mpg ~ hp <S3: lm> <S3: summary.lm> <data.frame [2 × 4]>
# 4                    mpg ~ drat <S3: lm> <S3: summary.lm> <data.frame [2 × 4]>
# 5              mpg ~ cyl + disp <S3: lm> <S3: summary.lm> <data.frame [3 × 4]>
# 6                mpg ~ cyl + hp <S3: lm> <S3: summary.lm> <data.frame [3 × 4]>
# 7              mpg ~ cyl + drat <S3: lm> <S3: summary.lm> <data.frame [3 × 4]>
# 8               mpg ~ disp + hp <S3: lm> <S3: summary.lm> <data.frame [3 × 4]>
# 9             mpg ~ disp + drat <S3: lm> <S3: summary.lm> <data.frame [3 × 4]>
# 10              mpg ~ hp + drat <S3: lm> <S3: summary.lm> <data.frame [3 × 4]>
# 11        mpg ~ cyl + disp + hp <S3: lm> <S3: summary.lm> <data.frame [4 × 4]>
# 12      mpg ~ cyl + disp + drat <S3: lm> <S3: summary.lm> <data.frame [4 × 4]>
# 13        mpg ~ cyl + hp + drat <S3: lm> <S3: summary.lm> <data.frame [4 × 4]>
# 14       mpg ~ disp + hp + drat <S3: lm> <S3: summary.lm> <data.frame [4 × 4]>
# 15 mpg ~ cyl + disp + hp + drat <S3: lm> <S3: summary.lm> <data.frame [5 × 4]>

我是分段进行的,保留模型摘要,主要用于演示和以防您重复使用lm(可能用于predict)。如果您知道永远不需要保留原始的 lm 输出,则可以将它们组合成一个函数调用。

我相信您要求的是系数的 data.frame,在这种情况下:

x$summaries[[1]]
# Call:
# lm(formula = as.formula(.), data = dat)
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -4.9814 -2.1185  0.2217  1.0717  7.5186 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  37.8846     2.0738   18.27  < 2e-16 ***
# cyl          -2.8758     0.3224   -8.92 6.11e-10 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Residual standard error: 3.206 on 30 degrees of freedom
# Multiple R-squared:  0.7262,  Adjusted R-squared:  0.7171 
# F-statistic: 79.56 on 1 and 30 DF,  p-value: 6.113e-10
coef(x$summaries[[1]])
#             Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) 37.88458  2.0738436 18.267808 8.369155e-18
# cyl         -2.87579  0.3224089 -8.919699 6.112687e-10

不幸的是,如果您尝试将所有这些系数汇总合并到一个 data.frame 中,则行名称会在 dplyr::bind_rows 中丢失:

bind_rows(map(x$summaries[1:2], ~ as.data.frame(coef(.))))
#      Estimate  Std. Error   t value     Pr(>|t|)
# 1 37.88457649 2.073843606 18.267808 8.369155e-18
# 2 -2.87579014 0.322408883 -8.919699 6.112687e-10
# 3 29.59985476 1.229719515 24.070411 3.576586e-21
# 4 -0.04121512 0.004711833 -8.747152 9.380327e-10

尽管您缺少“哪个模型”组件,但始终可以使用基础 R:

do.call(rbind.data.frame, map(x$summaries[1:2], ~ as.data.frame(coef(.))))
#                 Estimate  Std. Error   t value     Pr(>|t|)
# (Intercept)  37.88457649 2.073843606 18.267808 8.369155e-18
# cyl          -2.87579014 0.322408883 -8.919699 6.112687e-10
# (Intercept)1 29.59985476 1.229719515 24.070411 3.576586e-21
# disp         -0.04121512 0.004711833 -8.747152 9.380327e-1

我们可以通过在原始管道中使用tibble::rownames_to_column 重新引入:

x <- data_frame(Formulas) %>%
  mutate(
    lms = map(Formulas, ~ lm(as.formula(.), data = dat)),
    summaries = map(lms, ~ summary(.)),
    coefs = map(summaries, ~ tibble::rownames_to_column(as.data.frame(coef(.))))
  )
select(x, Formulas, coefs) %>% unnest()
# # A tibble: 47 × 6
#            Formulas     rowname    Estimate `Std. Error` `t value`   `Pr(>|t|)`
#               <chr>       <chr>       <dbl>        <dbl>     <dbl>        <dbl>
# 1         mpg ~ cyl (Intercept) 37.88457649  2.073843606 18.267808 8.369155e-18
# 2         mpg ~ cyl         cyl -2.87579014  0.322408883 -8.919699 6.112687e-10
# 3        mpg ~ disp (Intercept) 29.59985476  1.229719515 24.070411 3.576586e-21
# 4        mpg ~ disp        disp -0.04121512  0.004711833 -8.747152 9.380327e-10
# 5          mpg ~ hp (Intercept) 30.09886054  1.633920950 18.421246 6.642736e-18
# 6          mpg ~ hp          hp -0.06822828  0.010119304 -6.742389 1.787835e-07
# 7        mpg ~ drat (Intercept) -7.52461844  5.476662574 -1.373942 1.796391e-01
# 8        mpg ~ drat        drat  7.67823260  1.506705108  5.096042 1.776240e-05
# 9  mpg ~ cyl + disp (Intercept) 34.66099474  2.547003876 13.608536 4.022869e-14
# 10 mpg ~ cyl + disp         cyl -1.58727681  0.711844271 -2.229809 3.366495e-02
# # ... with 37 more rows

【讨论】:

    【解决方案2】:

    考虑通过调整您的最后一次 lapply 调用以返回数据帧来保持在基础 R 中:

    df_list <- lapply(seq_along(Formulas), function(i) { 
         mod <- summary(lm(as.formula(Formulas[[i]]), data = finalprojectdata3))
    
         data.frame(model_num = i,
                    formula = Formulas[[i]],
                    r2 = mod$r.squared,
                    adjr2 = mod$adj.r.squared
         )       
    })
    
    final_df <- do.call(rbind, df_list)
    
    final_tibble <- as_data_frame(finaldf)          # requires tidyverse
    

    使用 mtcars(借鉴 @r2evans 的可重现示例)

    final_tibble
    # A tibble: 15 x 4
    #    model_num                      formula        r2     adjr2
    #  *     <int>                       <fctr>     <dbl>     <dbl>
    #  1         1                    mpg ~ cyl 0.7261800 0.7170527
    #  2         2                   mpg ~ disp 0.7183433 0.7089548
    #  3         3                     mpg ~ hp 0.6024373 0.5891853
    #  4         4                   mpg ~ drat 0.4639952 0.4461283
    #  5         5             mpg ~ cyl + disp 0.7595658 0.7429841
    #  6         6               mpg ~ cyl + hp 0.7407084 0.7228263
    #  7         7             mpg ~ cyl + drat 0.7402482 0.7223343
    #  8         8              mpg ~ disp + hp 0.7482402 0.7308774
    #  9         9            mpg ~ disp + drat 0.7310094 0.7124583
    # 10        10              mpg ~ hp + drat 0.7411716 0.7233214
    # 11        11        mpg ~ cyl + disp + hp 0.7678877 0.7430186
    # 12        12      mpg ~ cyl + disp + drat 0.7650941 0.7399256
    # 13        13        mpg ~ cyl + hp + drat 0.7693992 0.7446920
    # 14        14       mpg ~ disp + hp + drat 0.7750131 0.7509073
    # 15        15 mpg ~ cyl + disp + hp + drat 0.7825119 0.7502914
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2019-09-06
      • 1970-01-01
      • 2021-12-05
      • 2021-08-19
      • 2023-04-09
      • 1970-01-01
      相关资源
      最近更新 更多