如何解决从具有多个模型的 lm lapply 循环中提取最终 p 值统计
我有以下代码可以在我的自变量 (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 |
这些是我用来获得上述相关性的实际值
解决方法
我认为没有存储 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
,
好的,我按以下方式编辑了我的代码:
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 个解决方案。认为学习处理问题的替代方法可能没问题,并感谢您的提问,非常好。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。