从线性回归中提取P值和R平方。

Question

如何从简单线性回归模型中获取 p 值（用于衡量单个自变量系数是否显著）和 R-squared 值？例如...

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-squared值，但我希望能够把这些信息放到其他变量中。

- grautur

---

只有当你不将输出分配给一个对象时才会显示值（例如， r <- summary(lm(rnorm(10)~runif(10))) 不会显示任何内容）。 - Joshua Ulrich

Answer 1

确定系数（r-squared）：您可以直接从摘要对象summary(fit)中返回确定系数值summary(fit)$r.squared。查看names(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]  

或者，你可以以类似的方式从anova(fit)对象中获取系数的p值，就像上面的摘要对象一样。

- Chase

---

使用inherits比直接使用class稍微好一点。也许你想要使用unname(pf(f[1],f[2],f[3],lower.tail=F))？ - hadley

如果您喜欢一行代码：summary(fit)$fstatistic %>% {unname(pf(.[1],.[2],.[3],lower.tail=F))} - Comfort Eagle

或者作为无管道的一行代码：with(summary(fit), pf(fstatistic[1],fstatistic[2],fstatistic[3],lower.tail=F)) - Comfort Eagle

Answer 2

请注意，summary(fit) 生成了一个包含您需要的所有信息的对象。beta、se、t 和 p 向量存储在其中。通过选择系数矩阵（存储在摘要对象中）的第4列来获取 p 值：

summary(fit)$coefficients[,4] 
summary(fit)$r.squared

尝试 str(summary(fit)) 查看此对象包含的所有信息。

编辑：我误读了Chase的答案，他基本上告诉您如何到达我在此处提供的内容。

- Vincent

---

注意：这是唯一的方法，可以轻松访问截距项和其他预测变量的P值。它是上述方法中最好的。 - Daniel Egan

这是正确的答案。排名第一的答案对我没有用。 - Chris

如果您想轻松获取P值，请使用此答案。为什么要编写多行函数或创建新对象（例如，ANOVA输出），当您只需更仔细地查找摘要输出中的P值即可。要单独分离一个 P 值本身，您需要将行号添加到 Vincent 的答案中：例如，对于截距，summary(fit)$coefficients[1,4]。 - theforestecologist

注意：此方法适用于使用lm()创建的模型，但不适用于gls()模型。 - theforestecologist

请查看下面的答案，以获取适用于 lm() 和 gls() 的扩展内容。 - theforestecologist

Chase的回答返回模型的p值，而这个回答返回系数的p值。在简单回归的情况下，它们是相同的，但在有多个预测变量的模型中，它们是不同的。因此，根据您想要提取的内容，这两个答案都很有用。 - Jeromy Anglim

Answer 3

您可以通过调用str(summary(fit))来查看summary()返回的对象的结构。每个部分都可以使用$访问。F统计量的p值可以更容易地从anova返回的对象中获得。

简而言之，您可以这样做：

rSquared <- summary(fit)$r.squared
pVal <- anova(fit)$'Pr(>F)'[1]

- jberg

---

此翻译仅适用于单变量回归分析，其中回归的p值与预测变量相同。 - Bakaburg

由于某些原因，当使用aov的summary对象时，此答案无法正常工作。 - seeker_after_truth

Answer 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"

- Konrad

---

这是一个很棒的回答！ - Andrew Brēza

Answer 5

虽然上面两个答案都很好，但提取对象部分的过程更为通用。

在许多情况下，函数返回列表，可以使用 str() 打印组件及其名称。然后可以使用 $ 运算符访问它们，即 myobject$componentname。

对于 lm 对象，有许多预定义的方法可供使用，例如 coef()、resid()、summary() 等，但你并不总是那么幸运。

- richiemorrisroe

Answer 6

在@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]  

• Note, you can replace the column number with the column name in each of the above instances:

  summary(fit)$coefficients[1,"Pr(>|t|)"]  ##lm 
  summary(fit)$tTable[1,"p-value"]         ##gls 
  

如果你还不确定如何从汇总表中获取一个值，请使用str()来确定汇总表的结构：

str(summary(fit))

- theforestecologist

Answer 7

这是提取p值最简单的方法：

coef(summary(modelname))[, "Pr(>|t|)"]

- RTrain3k

---

我尝试了这种方法，但如果线性模型包含任何NA项，它将失败。 - j_v_wow_d

Answer 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)

肯定有比这个函数更快的解决方案，但它能够工作。

- Dorian Grv

Answer 9

在summary()的最终p值中，该函数使用pf()函数从summary(fit)$fstatistic值计算得出。

fstat <- summary(fit)$fstatistic
pf(fstat[1], fstat[2], fstat[3], lower.tail=FALSE)

Source: [1], [2]

- Saftever

Answer 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

- dalloliogm