数据框每列上的多个线性模型

Question

我有一个包含 8 列 (y,x1,...,x7) 的通用 csv 文件，其中 y 代表响应变量，x 代表潜在预测变量，每个都有 100 个观察值。

我的目标是使用 R 来...

1 - 创建 7 个 y 散点图与每个 x 的散点图。

2 - 用每个 x 创建一个 y 的线性模型。

目前我只是重复输入所有内容...

p9 <- ggplot(generic, aes(x1,y)) + geom_point() + labs(title = "y vs x1")
p10 <- ggplot(generic, aes(x2,y)) + geom_point() + labs(title = "y vs x2")
p11 <- ggplot(generic, aes(x3,y)) + geom_point() + labs(title = "y vs x3")
p12 <- ggplot(generic, aes(x4,y)) + geom_point() + labs(title = "y vs x4")
p13 <- ggplot(generic, aes(x5,y)) + geom_point() + labs(title = "y vs x5")
p14 <- ggplot(generic, aes(x6,y)) + geom_point() + labs(title = "y vs x6")
p15 <- ggplot(generic, aes(x7,y)) + geom_point() + labs(title = "y vs x7")
grid.arrange(p9,p10,p11,p12,p13,p14,p15,ncol = 3)

和

x1mod <- lm(generic$y~generic$x1)
x2mod <- lm(generic$y~generic$x2)
x3mod <- lm(generic$y~generic$x3)
x4mod <- lm(generic$y~generic$x4)
x5mod <- lm(generic$y~generic$x5)
x6mod <- lm(generic$y~generic$x6)
x7mod <- lm(generic$y~generic$x7)
summary(x1mod)
summary(x2mod)
summary(x3mod)
summary(x4mod)
summary(x5mod)
summary(x6mod)
summary(x7mod)

我想减少重复性。 我试图用一个 for 循环来完成这个，但它变得有点混乱。 我还阅读了使用 purr 将函数映射到数据的内容，但我无法完全弄清楚如何使它适合我的情况，因为我没有尝试将数据按任何因素划分。 我对 R 比较陌生，所以如果我的问题简单得可笑，我深表歉意。

Answer 1

您可以在列级别使用apply() 来构建模型，并使用ggplot2 结合一些tidyverse 函数来达到您的结果。这是使用虚拟数据的代码：

library(ggplot2)
library(tidyverse)
#Randomness
set.seed(123)
#Data
df <- data.frame(y=rpois(15,0.8),
                 x1=rnorm(15,0,1),
                 x2=rnorm(15,1,1),
                 x3=rnorm(15,3,1),
                 x4=rnorm(15,4,1),
                 x5=rnorm(15,5,1),
                 x6=rnorm(15,6,1),
                 x7=rnorm(15,7,1))
#Models
Lmods <- apply(df[,-1],2,function(x) lm(y~x,data = df))
lapply(Lmods, summary)
#Plot
df %>% pivot_longer(-y) %>%
  mutate(name=paste0('y vs. ',name)) %>%
  ggplot(aes(x=value,y=y,color=name))+
  geom_point()+
  facet_wrap(name~.,scales = 'free')+
  theme(legend.position = 'none')

绘图的输出：

Answer 2

这里是创建 lm 的示例。

您可以开始定义哪个是您的 y 和您的 x。

require(tidyverse)

# Get all col names
cols <- names(mtcars)

y <- "mpg"

# Exclude y from cols
cols <- cols[cols != y]

您在这里定义一个函数，该函数将在输入中接收代表 x、y 和数据集的字符串。

# Create a function that return a lm for each x you want
make_lm <- function(col, y, dataset){
  formula <- paste0(y, "~", col)
  return(lm(formula, data = dataset))
}

# This will return a list with a model for each col
l_model <- cols %>% map(make_lm, y, mtcars)

# Here you can get the summary for each model
l_model %>% map(summary)

输出将是一个列表，其中包含 cols 中每个 col 的模型。

我建议使用 tidyverse 来处理这类事情。

Answer 3

给定数据：

set.seed(1)
df <- as.data.frame(replicate(8, rnorm(100)))
names(df) <- c("y", paste0("x", 1:7))

---

您可以在 tidyverse 工作流程中非常整洁地完成所有工作，这非常酷。

# libraries
library(dplyr)
library(tidyr)
library(ggplot2)

df %>%
 tibble::rowid_to_column() %>% 
 pivot_longer(-c(rowid,y)) %>% 
 nest_by(name) %>% 
 summarise(plot = list(ggplot(data, aes(value,y)) + geom_point() + labs(title = paste("y vs", name))),
           lm = list(summary(lm(y ~ value, data))))

#> `summarise()` regrouping output by 'name' (override with `.groups` argument)
#> # A tibble: 7 x 3
#> # Groups:   name [7]
#>   name  plot   lm        
#>   <chr> <list> <list>    
#> 1 x1    <gg>   <smmry.lm>
#> 2 x2    <gg>   <smmry.lm>
#> 3 x3    <gg>   <smmry.lm>
#> 4 x4    <gg>   <smmry.lm>
#> 5 x5    <gg>   <smmry.lm>
#> 6 x6    <gg>   <smmry.lm>
#> 7 x7    <gg>   <smmry.lm>

现在您在数据框中拥有了所有图表和摘要，您可以按照自己的方式处理它们。

---

或者更传统的解决方案：

# libraries
library(dplyr)
library(tidyr)
library(ggplot2)

# reorganize data
df_lng <- df %>%
 tibble::rowid_to_column() %>% 
 pivot_longer(-c(rowid,y))

# lm results
df_lng %>% 
 nest_by(name) %>% 
 summarise(broom::tidy(lm(y ~ value, data)))

#> `summarise()` regrouping output by 'name' (override with `.groups` argument)
#> # A tibble: 14 x 6
#> # Groups:   name [7]
#>    name  term         estimate std.error statistic p.value
#>    <chr> <chr>           <dbl>     <dbl>     <dbl>   <dbl>
#>  1 x1    (Intercept)  0.109       0.0903   1.20      0.231
#>  2 x1    value       -0.000932    0.0947  -0.00984   0.992
#>  3 x2    (Intercept)  0.108       0.0903   1.20      0.233
#>  4 x2    value        0.0160      0.0877   0.182     0.856
#>  5 x3    (Intercept)  0.111       0.0903   1.23      0.221
#>  6 x3    value       -0.0457      0.0914  -0.500     0.618
#>  7 x4    (Intercept)  0.113       0.0893   1.27      0.208
#>  8 x4    value        0.113       0.0767   1.47      0.144
#>  9 x5    (Intercept)  0.104       0.0898   1.16      0.248
#> 10 x5    value       -0.102       0.0933  -1.10      0.276
#> 11 x6    (Intercept)  0.123       0.0915   1.35      0.181
#> 12 x6    value        0.0717      0.0836   0.857     0.393
#> 13 x7    (Intercept)  0.109       0.0902   1.21      0.230
#> 14 x7    value        0.0293      0.0826   0.355     0.723


# plots
ggplot(df_lng, aes(x = value, y = y, colour = name)) +
 geom_point() + 
 labs(title = "y vs x") +
 facet_wrap("name")

Answer 4

reformulate() 是您应该知道的基本包中包含的基本功能。您可以使用它轻松创建模型拟合公式，特别适用于使用lapply 的多种方法：

xes <- names(dat)[-8]  ## storing x variable names

## fitting all possible bivariate models and store them in a list
res <- setNames(lapply(xes, function(x) {
  fo <- reformulate(x, "y")  ##dynamically create formula
  do.call("lm", list(fo, quote(dat)))
}), xes)

^{归功于 @G.Grothendieck 这个简洁的 do.call 方法，它将干净的调用属性分配给结果！}

## calculate summary accessing model X1
summary(res$X1)
# Call:
# lm(formula = y ~ X1, data = dat)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.52761 -0.43187 -0.06834  0.49291  1.55059 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)   2.8910     0.1321  21.889  < 2e-16 ***
# X1            1.0301     0.2185   4.714 8.04e-06 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.6565 on 98 degrees of freedom
# Multiple R-squared:  0.1848,  Adjusted R-squared:  0.1765 
# F-statistic: 22.22 on 1 and 98 DF,  p-value: 8.042e-06

您可以访问存储在模型调用中的公式并直接使用它们进行绘图。

op <- par(mfrow=c(2, 4))  ## sets pars
lapply(xes, function(x) {
  plot(res[[x]]$call$formula, dat, main=paste("y vs", x))
  abline(res[[x]])  ## optional for regression line
 })
par(op)  ## restores pars

注意：如果公式变得更复杂，我们可以使用 as.formula()（即对于具有随机效应或工具变量的模型，在 | 之后包含另一个项）。

res2 <- setNames(lapply(xes, function(x) {
  fo <- as.formula(paste("y ~ ", x))
  do.call("lm", list(fo, quote(dat)))
}), xes)

stopifnot(all.equal(res1, res2))

---

数据：

m <- 100;n <- 7
set.seed(42)
dat <- data.frame(matrix(runif(m*n), m, n))
dat <- transform(dat, y=X1 + X2 + X3 + X4 + X5 + X6 + X7)