5

我知道有几种方法可以比较回归模型。创建模型（从线性到多个）并比较 R2、Adjusted R2 等的一种方法：

Mod1: y=b0+b1
Mod2: y=b0+b1+b2
Mod3: y=b0+b1+b2+b3 (etc)

我知道有些包可以执行逐步回归，但我试图用 purrr来分析它。我可以创建几个简单的线性模型（感谢这里的这篇文章），现在我想知道如何创建回归模型，将特定的 IV 添加到方程：

可重现的代码

data(mtcars)
library(tidyverse)
library(purrr)
library(broom)
iv_vars <- c("cyl", "disp", "hp")
make_model <- function(nm) lm(mtcars[c("mpg", nm)])
fits <- Map(make_model, iv_vars)
glance_tidy <- function(x) c(unlist(glance(x)), unlist(tidy(x)[, -1]))
t(iv_vars %>% Map(f = make_model) %>% sapply(glance_tidy))

输出 线性模型的输出

我想要的是：

Mod1: mpg ~cyl
Mod2: mpg ~cly + disp
Mod3: mpg ~ cly + disp + hp

非常感谢。

4

2 回答 2

5

我将首先创建一个存储您的公式的列表小标题。然后将模型映射到公式上，并在模型上映射一瞥。

library(tidyverse)
library(broom)

mtcars %>% as_tibble()

formula <- c(mpg ~ cyl, mpg ~ cyl + disp)

output <-
  tibble(formula) %>% 
  mutate(model = map(formula, ~lm(formula = .x, data = mtcars)),
         glance = map(model, glance))

output$glance

output %>% unnest(glance)
于 2017-11-29T16:47:29.970 回答
3

您可以累积粘贴到您的向量上id_vars以获得您想要的组合。我使用此答案中的代码来执行此操作。

我使用加号作为变量之间的分隔符，以便为lm.

cumpaste = function(x, .sep = " ") {
     Reduce(function(x1, x2) paste(x1, x2, sep = .sep), x, accumulate = TRUE)
}

( iv_vars_cum = cumpaste(iv_vars, " + ") )

[1] "cyl"             "cyl + disp"      "cyl + disp + hp"

然后切换make_model函数以使用公式和数据集。由加号分隔的解释变量在公式中的波浪号之后传递给函数。一切都粘贴在一起，lm方便地解释为公式。

make_model = function(nm) {
     lm(paste0("mpg ~", nm), data = mtcars)
}

我们可以看到它按预期工作，返回一个具有两个解释变量的模型。

make_model("cyl + disp")

Call:
lm(formula = as.formula(paste0("mpg ~", nm)), data = mtcars)

Coefficients:
(Intercept)          cyl         disp  
   34.66099     -1.58728     -0.02058

您可能需要重新考虑如何将信息组合在一起，因为由于系数数量的增加，您现在将看到不同的列数。

一个可能的选项是添加dplyr::bind_rows到您的glance_tidy函数中，然后使用map_dfrfrom purrr作为最终输出。

glance_tidy = function(x) {
     dplyr::bind_rows( c( unlist(glance(x)), unlist(tidy(x)[, -1]) ) )
}

iv_vars_cum %>% 
     Map(f = make_model) %>% 
     map_dfr(glance_tidy, .id = "model")

# A tibble: 3 x 28

            model r.squared adj.r.squared    sigma statistic      p.value    df    logLik      AIC
            <chr>     <dbl>         <dbl>    <dbl>     <dbl>        <dbl> <dbl>     <dbl>    <dbl>
1             cyl 0.7261800     0.7170527 3.205902  79.56103 6.112687e-10     2 -81.65321 169.3064
2      cyl + disp 0.7595658     0.7429841 3.055466  45.80755 1.057904e-09     3 -79.57282 167.1456
3 cyl + disp + hp 0.7678877     0.7430186 3.055261  30.87710 5.053802e-09     4 -79.00921 168.0184 ...
于 2017-11-29T16:32:51.963 回答