programing

선형 회귀 분석에서 p-값 및 r-값 추출

newstyles 2023. 7. 8. 10:40

선형 회귀 분석에서 p-값 및 r-값 추출

간단한 선형 회귀 모형에서 p-값(단일 설명 변수가 0이 아닌 계수의 유의성)과 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-수정:요약 개체에서 직접 r 제곱 값을 반환할 수 있습니다.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-값을 반환하고 계수에 대한 p-값은 다음을 사용하여 추출할 수 있습니다.

summary(fit)$coefficients[,4]  

계수의 을 는다음 p-값 수 있습니다.anova(fit)위의 요약 개체와 유사한 방식으로 개체를 표시합니다.

:summary(fit)필요한 모든 정보가 포함된 개체를 생성합니다.베타, 집합, 집합 및 p 벡터가 저장됩니다.계수 행렬의 네 번째 열(요약 개체에 저장)을 선택하여 p-값을 가져옵니다.

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

ㅠㅠstr(summary(fit))이 개체에 포함된 모든 정보를 확인합니다.

편집: 기본적으로 제가 여기서 말하는 것에 어떻게 가는지를 알려주는 Chase의 답변을 잘못 읽었습니다.

반된개구체볼수있다습니조를의에 의해 반환된 를 볼 수 .summary()화로로 str(summary(fit))각 부품은 다음을 사용하여 액세스할 수 있습니다.$F 통계량에 대한 p-값은 반환된 개체로부터 더 쉽게 얻을 수 있습니다.anova.

간단히 말해, 다음과 같은 작업을 수행할 수 있습니다.

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

유사한 문제에 대한 제안된 해결책을 탐색하던 중 이 질문을 발견했습니다. 향후 참조를 위해 패키지를 활용하는 솔루션으로 사용 가능한 답변 목록을 업데이트하는 것이 가치가 있을 것으로 생각합니다.

샘플코드

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"

위의 두 가지 답변 모두 좋지만 개체의 일부를 추출하는 절차가 더 일반적입니다.

대부분의 경우 함수는 목록을 반환하고 개별 구성 요소는 다음을 사용하여 액세스할 수 있습니다.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]  
  • 위의 각 인스턴스에서 열 번호를 열 이름으로 바꿀 수 있습니다.

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

여전히 요약 표에서 값에 액세스하는 방법을 모르는 경우str()요약 표의 구조를 확인합니다.

str(summary(fit))

p-값을 끌어오는 가장 쉬운 방법은 다음과 같습니다.

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

저는 이 lmp 기능을 꽤 많이 사용했습니다.

데이터 분석을 강화하기 위해 새로운 기능을 추가하기로 결정했습니다.저는 R이나 통계량에 대한 전문가는 아니지만 사람들은 보통 선형 회귀 분석의 다른 정보를 봅니다.

  • p-값
  • a와 b
  • 그리고 물론 점 분포의 측면.

예를 들어 보겠습니다.여기 있습니다.

다음은 여러 변수를 사용한 재현 가능한 예제입니다.

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)

확실히 이 기능보다 빠른 해결책이 있지만 효과가 있습니다.

마지막에 표시되는 최종 p-값의 경우summary()사용하는 함수pf()에서 계산하면summary(fit)$fstatistic가치.

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

출처: [1], [2]

또 다른 옵션은 lm 대신 cor.test 기능을 사용하는 것입니다.

> 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

사용:

(summary(fit))$coefficients[***num***,4]

어디에num계수 행렬의 행을 나타내는 숫자입니다.모형에 있는 피쳐 수와 p-값을 추출할 피쳐 수에 따라 달라집니다.예를 들어, 변수가 하나만 있는 경우 절편에 대한 p-값은 [1,4]이고 실제 변수에 대한 다음 p-값은 [2,4]입니다.그럼 두 분이 되시는 거군요.

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

언급URL : https://stackoverflow.com/questions/5587676/pull-out-p-values-and-r-squared-from-a-linear-regression