我使用purrr在分组数据集的多列上运行一系列单一线性回归,但在不删除整个组的情况下排除没有数据的变量组时遇到了问题。
多亏了andrew_reece,我得到了如下的基本代码:
library(tidyverse)
ivs <- colnames(mtcars)[3:ncol(mtcars)]
names(ivs) <- ivs
mtcars %>%
group_by(cyl) %>%
group_modify(function(data, key) {
map_df(ivs, function(iv) {
frml <- as.formula(paste("mpg", "~", iv))
lm(frml, data = data) %>% broom::glance()
}, .id = "iv")
}) %>%
select(cyl, iv, r.squared, p.value)
它给出了这种格式的tibble:
cyl iv r.squared p.value
4 disp 0.6484051396 0.002782827
4 hp 0.2740558319 0.098398581
4 drat 0.180 0.193
4 wt 0.509 0.0137
4 qsec 0.0557 0.485
4 vs 0.00238 0.887
...
6 disp 0.0106260401 0.825929685
...
不幸的是,我的真实数据集很混乱,包含多个组变量组合,只有NA,或者实际值少于两个,lm无法处理。为了显示这一点,这里是mtcars,"disp"中的一些数据被NA替换。运行以上代码,mtna抛出一个NA错误。
#create mtcars dataset that will have a cyl group with entirely NA disp
mtna <- mtcars
mtna$disp[mtna$disp < 147] <- NA
test <- mtna %>% group_by(cyl) %>% summarize(mean = mean(disp))
我试图通过使lm具有条件来处理这个问题,并首先使用sum(!is.na(来检查是否有足够的实值来运行lm。这允许lm成功运行。
mtna %>%
group_by(cyl) %>%
group_modify(function(data, key) {
map_df(ivs, function(iv) {
tmpvar <- eval(parse(text = paste0("data$", iv)))
if(sum(!is.na(tmpvar)) < 3) {return(NA)} else {
frml <- as.formula(paste("mpg", "~", iv))
lm(frml, data = data) %>% broom::glance()
}}, .id = "iv")
}) %>%
select(cyl, iv, r.squared, p.value)
#which gives:
cyl iv r.squared p.value
1 4 NA NA NA
2 6 disp 0.0115 0.840
3 6 hp 0.0161 0.786
4 6 drat 0.0132 0.807
...
然而,当你查看结果时,你可以看到NA已经扩展到整个组,包括disp之外的变量(disp是唯一一个缺少值的变量(。现在根本没有与cyl = 4
相关的数据,即使在hp和drat这样没有缺失数据的组中也是如此。
我所希望的是:
cyl iv r.squared p.value
4 disp NA NA
4 hp 0.2740558319 0.098398581 # Currently missing
4 drat 0.1799791311 0.193450651 # Currently missing
4 wt 0.5086325963 0.013742782. # Currently missing
...
6 disp 0.0106260401 0.825929685
6 hp 0.0161462379 0.78602021
...
我怀疑这与数据格式有关——我想我正在将NA映射到该组的所有结果中,而不仅仅是一个变量。但我不知道如何解决这个问题。非常感谢您的帮助!
我想你差不多到了。看看下面的代码。对更改进行了注释。
library(tidyverse)
ivs <- colnames(mtcars)[3:ncol(mtcars)]
names(ivs) <- ivs
mtna <- mtcars
mtna$disp[mtna$disp < 147] <- NA
mtna %>%
group_by(cyl) %>%
group_modify(function(data, key) {
map_df(ivs, function(iv) {
tmpvar<-eval(parse(text = paste0("data$", iv)))
if(!is.na(sum(tmpvar))) { #only use complete data
frml <- as.formula(paste("mpg", "~", iv))
lm(frml, data = data) %>% broom::glance()
}}, .id = "iv")
}) %>%
select(cyl, iv, r.squared, p.value) %>%
right_join(.,expand.grid(cyl=unique(mtna$cyl),iv=ivs),
by=c("cyl","iv")) %>% #populating with NA for columns lost before
arrange(., cyl, iv) %>% #sort by cyl and iv
as.data.frame()
#which gives:
cyl iv r.squared p.value
1 4 am 0.2872892493 0.08921640
2 4 carb 0.0378466325 0.56650426
3 4 disp NA NA
4 4 drat 0.1799791311 0.19345065
5 4 gear 0.1146552225 0.30840242
6 4 hp 0.2740558319 0.09839858
7 4 qsec 0.0556742395 0.48487154
8 4 vs 0.0023819528 0.88668635
9 4 wt 0.5086325963 0.01374278
10 6 am 0.2810551424 0.22094408
11 6 carb 0.0000661434 0.98619369
12 6 disp NA NA
13 6 drat 0.0131597553 0.80653099
14 6 gear 0.0000901510 0.98388187
15 6 hp 0.0161462379 0.78602021
16 6 qsec 0.1753245893 0.34980138
17 6 vs 0.2810551424 0.22094408
18 6 wt 0.4645101506 0.09175766
19 8 am 0.0024647887 0.86615546
20 8 carb 0.1550637156 0.16359436
21 8 disp 0.2701577717 0.05677488
22 8 drat 0.0022975223 0.87074078
23 8 gear 0.0024647887 0.86615546
24 8 hp 0.0804491933 0.32575378
25 8 qsec 0.0108860059 0.72261712
26 8 vs 0.0000000000 NA
27 8 wt 0.4229655365 0.01179281