从具有多个模型的 lm lapply 循环中提取最终 p 值统计量

Question

我有以下代码在我的自变量 (Kpl) 和我的所有其他因变量 (Y1, Y2, ...., Yi) 之间自动执行 lm：

linear_summary <- lapply(testdata[,-1], function(x) summary(lm(Kpl ~ x)))

这个的输出是

Call:
lm(formula = Kpl ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.37567 -0.52392  0.04236  0.67444  0.81316 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.7282     0.3456   5.001 0.000402 ***
x            -0.1550     0.2712  -0.571 0.579196    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.772 on 11 degrees of freedom
Multiple R-squared:  0.02883,   Adjusted R-squared:  -0.05946 
F-statistic: 0.3265 on 1 and 11 DF,  p-value: 0.5792


$Y2

Call:
lm(formula = Kpl ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2472 -0.4236 -0.2057  0.7140  1.0348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.6900     0.9010   0.766    0.460
x             0.8832     0.8767   1.007    0.335

Residual standard error: 0.7495 on 11 degrees of freedom
Multiple R-squared:  0.08447,   Adjusted R-squared:  0.001238 
F-statistic: 1.015 on 1 and 11 DF,  p-value: 0.3354

等等。 （我只截断了前 2 个相关性）

我想为每个实例提取整个模型的最终 p 值（在这两种情况下分别为 0.5792 和 0.3354）。理想情况下，这将以某种表格形式出现，并带有相关的相关变量，即 Y1=0.5792 Y2=0.3354。

我能找到的大部分信息要么似乎只适用于单个相关性（而不是具有多个相关性的 sapply），要么我似乎无法让它工作，这可能是我原始代码的问题。

对于刚开始使用 R 的人有什么建议可以解决这个问题？

编辑：数据看起来像这样

|    X     |     Y1      |     Y2      |     Y3      |     Y4      |
| -------- | ------------|-------------|-------------|-------------|
| 0.33767  | 2.33063062  | 1.013212308 | 1.277996888 | 1.373238355 |
| 0.33767  | 0.095967324 | 0.508830529 | 0.789257027 | 0.815877121 |
| 1.010474 | 2.344657045 | 0.842490752 | 1.240582283 | 1.262360905 |
| 1.010474 | 0.08135992  | 0.912535398 | 0.384427466 | 0.409817599 |
| 1.183276 | 0.135626937 | 0.967877981 | 0.505801442 | 0.576288093 |
| 1.536974 | 1.507146148 | 1.428839993 | 1.316569449 | 1.392022619 |
| 1.536974 | 1.255210981 | 1.191822955 | 1.395769591 | 1.41903939  |
| 2.017965 | 1.410299711 | 1.121560244 | 1.369835675 | 1.385143026 |
| 2.017965 | 1.032587109 | 1.372235121 | 1.390878783 | 1.42741762  |
| 2.3436   | 1.275999998 | 0.930400789 | 1.19877482  | 1.217540034 |
| 2.3436   | 1.250513383 | 1.063880146 | 1.206719195 | 1.23325973  |
| 2.387598 | 0.182866909 | 0.89588293  | 0.416923749 | 0.45364797  |
| 2.387598 | 0.097133916 | 0.750430855 | 0.506463633 | 0.03434754  |

这些是我用来获得上述相关性的实际值

Answer 1

我认为p值没有存储，你需要从fstatistics中计算出来，大概是这样的：

set.seed(111)
testdata = data.frame(Kpl = rnorm(100), Y1 = rnorm(100),
                      Y2 = rnorm(100), Y3 = rnorm(100))

IV = colnames(testdata)[-1]
DV = "Kpl"

linear_summary <- lapply(IV,function(x){
         summary(lm(reformulate(response=DV,termlabels=x),data=testdata))
                         })

names(linear_summary) = IV

tab = lapply(IV,function(x){
  p = with(
       linear_summary[[x]],
       pf(fstatistic[1],fstatistic[2],fstatistic[3],lower.tail=FALSE)
          )
  data.frame(IV = x, p = p)
})

do.call(rbind,tab)

       IV         p
value  Y1 0.5757187
value1 Y2 0.4922582
value2 Y3 0.4009439

检查例如第一个摘要：

linear_summary[[1]]

Call:
lm(formula = reformulate(response = DV, termlabels = x), data = testdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.94515 -0.73325  0.05448  0.57901  2.76026 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.01382    0.10747  -0.129    0.898
Y1          -0.05950    0.10597  -0.562    0.576

Residual standard error: 1.075 on 98 degrees of freedom
Multiple R-squared:  0.003207,  Adjusted R-squared:  -0.006964 
F-statistic: 0.3153 on 1 and 98 DF,  p-value: 0.5757

Answer 2

好的，我按以下方式编辑了我的代码：

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

df %>%    # Solution 1
  pivot_longer(-X) %>%
  group_split(name) %>%
  set_names(nm = map(., ~ first(.x$name))) %>%
  map(~ tidy(lm(X ~ value, data = .))) %>%
  bind_rows(.id = "var") %>%
  filter(term == "value")

# A tibble: 4 x 6
  var   term  estimate std.error statistic p.value
  <chr> <chr>    <dbl>     <dbl>     <dbl>   <dbl>
1 Y1    value  -0.155      0.271   -0.571    0.579
2 Y2    value   0.883      0.877    1.01     0.335
3 Y3    value   0.0341     0.552    0.0618   0.952
4 Y4    value  -0.158      0.469   -0.337    0.743

或者你可以使用这个：

df %>%    # Solution 2
  pivot_longer(Y1:Y4) %>%
  group_by(name) %>%
  arrange(.by_group = TRUE) %>% 
  nest() %>%
  mutate(models = map(data, ~ lm(X ~ value, data = .)),
         glance = map(models, glance)) %>%
  unnest(glance)

# A tibble: 4 x 15
# Groups:   name [4]
  name  data    models r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
  <chr> <list>  <list>     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
1 Y1    <tibbl~ <lm>    0.0288        -0.0595  0.772   0.327     0.579     1  -14.0  34.0  35.7
2 Y2    <tibbl~ <lm>    0.0845         0.00124 0.750   1.01      0.335     1  -13.6  33.2  34.9
3 Y3    <tibbl~ <lm>    0.000348      -0.0905  0.783   0.00382   0.952     1  -14.2  34.4  36.1
4 Y4    <tibbl~ <lm>    0.0102        -0.0798  0.779   0.113     0.743     1  -14.1  34.2  35.9
# ... with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

我知道您已经得到了答案，但在这里我提出了另外 2 个解决方案。认为学习处理问题的替代方法可能没问题，感谢您的问题，非常好。