如何从一个简单的线性回归模型中提取p值(单个解释变量的系数的显著性为非零)和r平方值?例如……
x = cumsum(c(0, runif(100, -1, +1)))
y = cumsum(c(0, runif(100, -1, +1)))
fit = lm(y ~ x)
summary(fit)
我知道summary(fit)显示了p值和r平方值,但我希望能够将它们插入到其他变量中。
如何从一个简单的线性回归模型中提取p值(单个解释变量的系数的显著性为非零)和r平方值?例如……
x = cumsum(c(0, runif(100, -1, +1)))
y = cumsum(c(0, runif(100, -1, +1)))
fit = lm(y ~ x)
summary(fit)
我知道summary(fit)显示了p值和r平方值,但我希望能够将它们插入到其他变量中。
当前回答
对于summary()末尾显示的最终p值,该函数使用pf()从summary(fit)$fstatistic值中计算。
fstat <- summary(fit)$fstatistic
pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE)
来源:[1],[2]
其他回答
r-squared:您可以直接从摘要对象summary(fit)$r.squared返回r-squared值。有关可以直接提取的所有项的列表,请参阅名称(summary(fit))。
模型p值:如果要得到整体回归模型的p值, 这篇博文概述了一个返回p值的函数:
lmp <- function (modelobject) {
if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
f <- summary(modelobject)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=F)
attributes(p) <- NULL
return(p)
}
> lmp(fit)
[1] 1.622665e-05
在只有一个预测因子的简单回归的情况下,模型p值和系数p值将是相同的。
系数p值:如果你有一个以上的预测器,那么上面将返回模型p值,系数p值可以使用以下方法提取:
summary(fit)$coefficients[,4]
或者,您可以以类似于上面的摘要对象的方式从方差分析(fit)对象中获取系数的p值。
虽然上面的两个答案都很好,但提取对象部分的过程更一般。
在许多情况下,函数返回列表,并且可以使用str()访问各个组件,str()将打印组件及其名称。然后您可以使用$操作符访问它们,即myobject$componentname。
在lm对象的情况下,有许多预定义的方法可以使用,如coef()、remainder()、summary()等,但你不会总是那么幸运。
@Vincent回答的延伸:
对于lm()生成的模型:
summary(fit)$coefficients[,4] ##P-values
summary(fit)$r.squared ##R squared values
对于gls()生成的模型:
summary(fit)$tTable[,4] ##P-values
##R-squared values are not generated b/c gls uses max-likelihood not Sums of Squares
为了隔离单独的p值本身,你需要在代码中添加一个行号:
例如,要访问两个模型摘要中截距的p值:
summary(fit)$coefficients[1,4]
summary(fit)$tTable[1,4]
注意,你可以在上面的每个实例中用列名替换列号: 总结(适合)$系数[1,“公关(> | t |)“# # lm 总结(适合)$ tTable[1,“假定值 "] ## gl
如果你仍然不确定如何从汇总表中访问一个值,使用str()来找出汇总表的结构:
str(summary(fit))
我在探索类似问题的建议解决方案时遇到了这个问题;我认为,为了将来的参考,它可能是值得更新可用的答案列表与解决方案利用扫帚包。
示例代码
x = cumsum(c(0, runif(100, -1, +1)))
y = cumsum(c(0, runif(100, -1, +1)))
fit = lm(y ~ x)
require(broom)
glance(fit)
结果
>> glance(fit)
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
1 0.5442762 0.5396729 1.502943 118.2368 1.3719e-18 2 -183.4527 372.9055 380.7508 223.6251 99
一边笔记
我发现glance函数很有用,因为它简洁地总结了关键值。结果被存储为data.frame,这使得进一步的操作很容易:
>> class(glance(fit))
[1] "data.frame"
另一个选择是使用cor.test函数,而不是lm:
> x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
> y <- c( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8)
> mycor = cor.test(x,y)
> mylm = lm(x~y)
# r and rsquared:
> cor.test(x,y)$estimate ** 2
cor
0.3262484
> summary(lm(x~y))$r.squared
[1] 0.3262484
# P.value
> lmp(lm(x~y)) # Using the lmp function defined in Chase's answer
[1] 0.1081731
> cor.test(x,y)$p.value
[1] 0.1081731