从 lm 函数中获取分组数据的 p 值答案

【问题标题】：Get the p-values from the lm function for grouped data从 lm 函数中获取分组数据的 p 值
【发布时间】：2016-07-04 02:09:48
【问题描述】：

我正在尝试使用lm() 函数和plyr 包为我的数据中的每个段拟合一个模型，因为我的数据是按一个键分组的。

我已经设法运行模型并获得系数以及 R^2 和 adj r 平方，但我正在努力处理 p 值。

library("plyr")
#Sample data
test_data <- data.frame(key = c("a","a","a","a","a","b","b","b","b","b"),
y = c(100,180,120,60,140,200,220,240,260,280),
x1 = c(50,60,79,85,90,133,140,120,160,170),
x2 =  c(20,18,47,16,15,25,30,25,20,15))

#model
model_1 <- dlply(test_data, .(key), 
    function(test_data) lm(y ~ x1 +     x2,data = test_data))

#coefficients
ldply(model_1, coef)

#adj r-squared
ldply(model_1, function(x) summary(x)$r.squared)

我已经尝试过这个，它得到了我的键和 p 值，但它没有变量的名称，我需要稍后将输出与模型中的系数合并。

#p-values but missing the variable names
ldply(model_1, function(x) summary(x)$coefficients)[,c(1,5)]

我尝试使用 Do 和 tidy 从 dplyr 包中拟合模型，这适用于小型数据集，因为它实际上返回了我需要的所有内容，但我的实际数据包含超过 1,000 种不同段和 RStudio 最终崩溃。

【问题讨论】：

nlme::lmList 可能对您来说也很有趣。看看这个！除此之外，我很惊讶地听到dplyr::do 方法与broom::tidy 结合崩溃。您是否有机会重新创建错误？
@coeffeinjunky - 抱歉，“崩溃”是指 RStudio 变得无响应，因此没有“错误”。当我在我的数据集上运行 dplyr 时，它在大约 4 秒内通过了 1000 多个模型，但是当我尝试“做”和“整理”时，它运行了将近 10 分钟，而 RStudio 变得没有响应。

标签： r dplyr plyr lm

【解决方案1】：

我正在使用“dplyr”包来格式化输出。在你在“dlply”函数中使用的函数中，你应该使用summary()到lm()，所以当你调用“coef”时，它也会包含p.values。

test_data <- data.frame(key = c("a","a","a","a","a","b","b","b","b","b"),
                        y = c(100,180,120,60,140,200,220,240,260,280),
                        x1 = c(50,60,79,85,90,133,140,120,160,170),
                        x2 =  c(20,18,47,16,15,25,30,25,20,15))

model<-by(test_data,test_data$key,function(x)summary(lm(y~x1+x2,x)))

R2<-t(data.frame(lapply(model,function(x)x$adj.r.squared)));colnames(R2)<-"R2_adj";R2

      R2_adj
a -0.8939647
b  0.4292186

Co<-as.data.frame(t(data.frame(lapply(model,function(x)x$coef))))

colnames(Co)<-c("intercept","x1","x2")

library(dplyr)

Co%>%
        mutate(key=substr(rownames(Co),1,1),
               variable=substr(rownames(Co),3,12))%>%
        select(key,variable,intercept,x1,x2)

  key   variable   intercept         x1          x2
1   a   Estimate 162.1822438 -0.6037364  0.07628315
2   a Std..Error 141.3436897  1.8054132  2.29385395
3   a    t.value   1.1474318 -0.3344035  0.03325545
4   a   Pr...t..   0.3699423  0.7698867  0.97649134
5   b   Estimate 271.0532276  0.3624009 -3.62853907
6   b Std..Error 196.2769562  0.9166979  3.25911570
7   b    t.value   1.3809733  0.3953330 -1.11335080
8   b   Pr...t..   0.3013515  0.7307786  0.38142882

【讨论】：

【解决方案2】：

不需要plyr 我认为sapply 就可以了。

sapply(model_1, function(x) summary(x)$coefficients[, 4])

                    a         b
(Intercept) 0.3699423 0.3013515
x1          0.7698867 0.7307786
x2          0.9764913 0.3814288

t() 将获得与您的估计相同的配置。

顺便说一句，您可能想看看multidplyr 包，毕竟与tidy 和dplyr::do 有关。

【讨论】：