从线性回归中提取p值和r平方
你如何从一个简单的线性回归模型中得出p值（单个解释variables系数的非显着性的意义）和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平方值，但是我想能够将这些值粘贴到其他variables中。
您可以直接从摘要对象summary(fit)$r.squared返回r平方值。 请参阅names(summary(fit))以获取可以直接提取的所有项目的列表。
这个博客文章概述了一个返回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
或者，您可以按照与上面的摘要对象类似的方式从anova(fit)对象中获取p值。
注意， summary(fit)生成一个包含所有需要的信息的对象。 beta，se，t和p向量存储在其中。 通过select系数matrix的第四列（存储在摘要对象中）获取p值：
summary(fit)$coefficients[,4] summary(fit)$r.squared
尝试str(summary(fit))来查看该对象包含的所有信息。
编辑：我误解了蔡斯的答案，基本上告诉你如何得到我在这里给的。
您可以通过调用str(summary(fit))来查看summary()返回的对象的结构。 每件作品都可以使用$进行访问。 F统计量的p值更容易从anova返回的对象中获得。
简而言之，你可以这样做：
rSquared <- summary(fit)$r.squared pVal <- anova(fit)$'Pr(>F)'[1]
虽然上面的两个答案都是好的，但是提取部分对象的过程更为一般。
在许多情况下，函数返回列表，并且可以使用str()来访问各个组件， str()将打印组件及其名称。 然后可以使用$运算符访问它们，即myobject$componentname 。
在lm对象的情况下，可以使用coef() ， resid() ， 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]
-
请注意，您可以使用以上每个实例中的列名称replace列号：
summary(fit)$coefficients[1,"Pr(>|t|)"] ##lm summary(fit)$tTable[1,"p-value"] ##gls
如果您仍然不确定如何访问汇总表中的值，请使用str()来计算汇总表的结构：
str(summary(fit))
我研究这个问题，同时探讨针对类似问题的build议解决scheme; 我认为，为了将来的参考，可能值得用利用broom包的解决scheme来更新可用的答案列表。
示例代码
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
旁注
我发现浏览function有用，因为它整齐地总结了有用的价值。 作为一个额外的好处，结果被存储为一个data.frame ，使得进一步的操作变得简单：
>> class(glance(fit)) [1] "data.frame"
这是拉取p值最简单的方法：
coef(summary(modelname))[, "Pr(>|t|)"]
我用了很多次这个lmp函数。
有一次，我决定添加新的function来增强数据分析。 我不是R或统计专家，但是人们通常会看到线性回归的不同信息：
- p值
- a和b
- R 2
- 当然还有点分布的方面
举个例子吧。 你在这里
这里有一个不同variables的可重复的例子：
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)
肯定有一个比这个function更快的解决scheme，但它的工作原理。
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
使用：
(summary(fit))$coefficients[***num***,4]
其中num是表示系数matrix的行的数字。 这将取决于您在模型中具有多lessfunction，以及您想要为哪个function提取p值。 例如，如果你只有一个variables，将会有一个P-值为[1,4]，而下一个为你的实际variables[2,4]。 所以你的num将是2。
另一种select是使用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
