从线性回归中提取 p 值和 r 平方
Posted
技术标签:
【中文标题】从线性回归中提取 p 值和 r 平方【英文标题】:pull out p-values and r-squared from a linear regression 【发布时间】:2011-08-01 01:19:15 【问题描述】:如何从简单的线性回归模型中提取 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 平方值,但我希望能够将它们粘贴到其他变量中。
【问题讨论】:
如果你不将输出分配给一个对象,它只会显示值(例如r <- summary(lm(rnorm(10)~runif(10)))
不显示任何东西)。
【参考方案1】:
r-squared:您可以直接从摘要对象summary(fit)$r.squared
返回r-squared 值。请参阅names(summary(fit))
了解您可以直接提取的所有项目的列表。
模型p值:如果要获取整体回归模型的p值, this blog post 概述了一个返回 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]
或者,您可以采用与上述摘要对象类似的方式从 anova(fit)
对象中获取系数的 p 值。
【讨论】:
直接使用inherits
比直接使用class
好一点。也许你想要unname(pf(f[1],f[2],f[3],lower.tail=F))
?
如果您更喜欢单线:summary(fit)$fstatistic %>% unname(pf(.[1],.[2],.[3],lower.tail=F))
或作为无管道单衬管:with(summary(fit), pf(fstatistic[1],fstatistic[2],fstatistic[3],lower.tail=F))
【参考方案2】:
请注意,summary(fit)
会生成一个包含您需要的所有信息的对象。 beta、se、t 和 p 向量存储在其中。通过选择系数矩阵的第 4 列(存储在汇总对象中)来获取 p 值:
summary(fit)$coefficients[,4]
summary(fit)$r.squared
尝试str(summary(fit))
查看此对象包含的所有信息。
编辑:我误读了 Chase 的回答,它基本上告诉你如何得到我在这里给出的内容。
【讨论】:
注意:这是唯一可以让您轻松访问截距的 p 值以及其他预测变量的方法。到目前为止最好的。 这是正确的答案。评价最高的答案对我不起作用。 如果您想轻松访问 P 值,请使用此答案。为什么要编写多行函数或创建新对象(即方差分析输出),而您只需要更加努力地在摘要输出本身中找到 p 值。要隔离单个 p 值本身,您可以在文森特的答案中添加一个行号:例如,summary(fit)$coefficients[1,4]
用于拦截
注意:此方法适用于使用lm()
创建的模型,但不适用于gls()
模型。
Chase 的答案返回模型的 p 值,这个答案返回系数的 p 值。在简单回归的情况下,它们是相同的,但在具有多个预测变量的模型的情况下,它们就不一样了。因此,这两个答案都很有用,具体取决于您要提取的内容。【参考方案3】:
你可以通过调用str(summary(fit))
查看summary()
返回的对象的结构。每件作品都可以使用$
访问。 F 统计量的 p 值更容易从 anova
返回的对象中获得。
简而言之,您可以这样做:
rSquared <- summary(fit)$r.squared
pVal <- anova(fit)$'Pr(>F)'[1]
【讨论】:
这仅适用于回归的 p val 与预测变量相同的单变量回归【参考方案4】:我在探索类似问题的建议解决方案时遇到了这个问题;我认为为了将来的参考,可能值得使用broom
包的解决方案更新可用的答案列表。
示例代码
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"
【讨论】:
这是一个很好的答案!【参考方案5】:虽然上述两个答案都不错,但提取部分对象的过程更为通用。
在许多情况下,函数返回列表,并且可以使用str()
访问各个组件,这将打印组件及其名称。然后,您可以使用 $ 运算符访问它们,即myobject$componentname
。
对于 lm 对象,可以使用许多预定义的方法,例如 coef()
、resid()
、summary()
等,但您不会总是那么幸运。
【讨论】:
【参考方案6】:@Vincent 的answer 的扩展:
对于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]
注意,在上述每个实例中,您都可以将列号替换为列名:
summary(fit)$coefficients[1,"Pr(>|t|)"] ##lm
summary(fit)$tTable[1,"p-value"] ##gls
如果您仍然不确定如何访问汇总表中的值,请使用str()
找出汇总表的结构:
str(summary(fit))
【讨论】:
【参考方案7】:这是提取 p 值的最简单方法:
coef(summary(modelname))[, "Pr(>|t|)"]
【讨论】:
我试过这个方法,但是如果线性模型包含任何NA术语,它会失败【参考方案8】:我用过这个 lmp 函数很多次了。
有一次,我决定添加新功能来增强数据分析。我不是 R 或统计学方面的专家,但人们通常会查看线性回归的不同信息:
p 值 a 和 b r² 当然还有点分布方面让我们举个例子。你在这里
这是一个具有不同变量的可重现示例:
Ex<-structure(list(X1 = c(-36.8598, -37.1726, -36.4343, -36.8644,
-37.0599, -34.8818, -31.9907, -37.8304, -34.3367, -31.2984, -33.5731
), X2 = c(64.26, 63.085, 66.36, 61.08, 61.57, 65.04, 72.69, 63.83,
67.555, 76.06, 68.61), Y1 = c(493.81544, 493.81544, 494.54173,
494.61364, 494.61381, 494.38717, 494.64122, 493.73265, 494.04246,
494.92989, 494.98384), Y2 = c(489.704166, 489.704166, 490.710962,
490.653212, 490.710612, 489.822928, 488.160904, 489.747776, 490.600579,
488.946738, 490.398958), Y3 = c(-19L, -19L, -19L, -23L, -30L,
-43L, -43L, -2L, -58L, -47L, -61L)), .Names = c("X1", "X2", "Y1",
"Y2", "Y3"), row.names = c(NA, 11L), class = "data.frame")
library(reshape2)
library(ggplot2)
Ex2<-melt(Ex,id=c("X1","X2"))
colnames(Ex2)[3:4]<-c("Y","Yvalue")
Ex3<-melt(Ex2,id=c("Y","Yvalue"))
colnames(Ex3)[3:4]<-c("X","Xvalue")
ggplot(Ex3,aes(Xvalue,Yvalue))+
geom_smooth(method="lm",alpha=0.2,size=1,color="grey")+
geom_point(size=2)+
facet_grid(Y~X,scales='free')
#Use the lmp function
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)
# create function to extract different informations from lm
lmtable<-function (var1,var2,data,signi=NULL)
#var1= y data : colnames of data as.character, so "Y1" or c("Y1","Y2") for example
#var2= x data : colnames of data as.character, so "X1" or c("X1","X2") for example
#data= data in dataframe, variables in columns
# if signi TRUE, round p-value with 2 digits and add *** if <0.001, ** if < 0.01, * if < 0.05.
if (class(data) != "data.frame") stop("Not an object of class 'data.frame' ")
Tabtemp<-data.frame(matrix(NA,ncol=6,nrow=length(var1)*length(var2)))
for (i in 1:length(var2))
Tabtemp[((length(var1)*i)-(length(var1)-1)):(length(var1)*i),1]<-var1
Tabtemp[((length(var1)*i)-(length(var1)-1)):(length(var1)*i),2]<-var2[i]
colnames(Tabtemp)<-c("Var.y","Var.x","p-value","a","b","r^2")
for (n in 1:length(var1))
Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),3]<-lmp(lm(data[,var1[n]]~data[,var2[i]],data))
Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),4]<-coef(lm(data[,var1[n]]~data[,var2[i]],data))[1]
Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),5]<-coef(lm(data[,var1[n]]~data[,var2[i]],data))[2]
Tabtemp[(((length(var1)*i)-(length(var1)-1))+n-1),6]<-summary(lm(data[,var1[n]]~data[,var2[i]],data))$r.squared
signi2<-data.frame(matrix(NA,ncol=3,nrow=nrow(Tabtemp)))
signi2[,1]<-ifelse(Tabtemp[,3]<0.001,paste0("***"),ifelse(Tabtemp[,3]<0.01,paste0("**"),ifelse(Tabtemp[,3]<0.05,paste0("*"),paste0(""))))
signi2[,2]<-round(Tabtemp[,3],2)
signi2[,3]<-paste0(format(signi2[,2],digits=2),signi2[,1])
for (l in 1:nrow(Tabtemp))
Tabtemp$"p-value"[l]<-ifelse(is.null(signi),
Tabtemp$"p-value"[l],
ifelse(isTRUE(signi),
paste0(signi2[,3][l]),
Tabtemp$"p-value"[l]))
Tabtemp
# ------- EXAMPLES ------
lmtable("Y1","X1",Ex)
lmtable(c("Y1","Y2","Y3"),c("X1","X2"),Ex)
lmtable(c("Y1","Y2","Y3"),c("X1","X2"),Ex,signi=TRUE)
当然有比这个功能更快的解决方案,但它确实有效。
【讨论】:
【参考方案9】:对于在summary()
末尾显示的最终 p 值,该函数使用pf()
从summary(fit)$fstatistic
值进行计算。
fstat <- summary(fit)$fstatistic
pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE)
来源:[1],[2]
【讨论】:
【参考方案10】:另一种选择是使用 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
【讨论】:
【参考方案11】:x = cumsum(c(0, runif(100, -1, +1)))
y = cumsum(c(0, runif(100, -1, +1)))
fit = lm(y ~ x)
> names(summary(fit))
[1] "call" "terms"
[3] "residuals" "coefficients"
[5] "aliased" "sigma"
[7] "df" "r.squared"
[9] "adj.r.squared" "fstatistic"
[11] "cov.unscaled"
summary(fit)$r.squared
【讨论】:
愿意提供一个解释,即使是简短的,解释为什么这个代码有效? 这对现有答案(尤其是已接受的答案)有何改进?【参考方案12】:用途:
(summary(fit))$coefficients[***num***,4]
其中num
是一个数字,表示系数矩阵的行。这将取决于您的模型中有多少特征,以及您想为哪一个提取 p 值。例如,如果您只有一个变量,则截距的一个 p 值为 [1,4],而实际变量的下一个 p 值为 [2,4]。所以你的 num
将是 2。
【讨论】:
以上是关于从线性回归中提取 p 值和 r 平方的主要内容,如果未能解决你的问题,请参考以下文章