본문 바로가기
공부/IT-R프로그래밍

R프로그래밍 강좌 - [15] [고급통계] 선형회귀분석, 선형회귀분석 검증, 다중 선형회귀분석

by 썸볼 2018. 1. 17.

1. 회귀분석 
(1) 선형회귀분석

- 회귀분석은 진화론에 해당하는 키가 큰 아버지의 자식은 점점 커질것인고 키가 작은 아버지의 자식은 점점 작아질 것이라는 가설을
  반박 하기 위해 연구되었다.

- 칼톤은 아들의 키는 아버지의 키에 영향을 받는다 하더라도 결국 평균으로 돌아가려는 현상이 있다는 것을 발견한다.
- 칼톤의 연구에 이어받은 피어슨은 신장자료를 이어받아 다음과 같은 함수식을 구한다.

         Y = 33.73 + 0.516X

   X는 아버지의 키이다.
   70을 X에 대입하면 69.85;         

   72을 X에 대입하면 72.45; 
   65을 X에 대입하면 67.27; 
   60을 X에 대입하면 64.69:의 값이 나온다.  값이 평균에서 크게 벗어나지 않는다. 즉, 평균으로 회귀하려는 모습을 볼 수 있다.

- 위에 내용에서 아버지 키와 아들의 키가 관련성이 있다는것을 알수 있다.

- 회귀분석은 변수간의 관계에 대해서 말하는 것이다. 
- 상관계수는 관계의 긴밀함을 수치로적으로 표현하는 것이며 회귀분석은 한 변수의 변화에 따른 다른 변수의 값을 파악 할 수 있게 한다.

- 아래는 사람들의 월소득과 월지출을 조사했다. 

   월 소득은 100만원인 사람의 월 지출은 50만원이다.
   
월 소득은 200만원인 사람의 월 지출은 100만원이다.
   월 소득은 300만원인 사람의 월 지출은 150만원이다.
     월 소득은 400만원인 사람의 월 지출은 200만원이다.

- 위내용에서 만약 500만원 소득인 사람의 월 지출은 250만원일 것이라 추측할 수 있다.
- 수학 함수식으로 표현하면 y=0.5x 이다. y=0.5x만 알고 있다면 어떤 값을 대입해도 결과를 알 수 있다.
- 위와 같은 식을 회귀모델이라 부르며 x축은 독립변수, y축은 종속변수라 한다.

- 회귀분석은 독립변수가 종속변수에 어떤 영향을 미치는지 알아보는 방법이다.
 

- 그러나 위 예시와 같이 극단적으로 데이터가 정확히 나오는 경우는 거의 없다. 월 소득이 600만원인데 지출이 331만원일 수 있고,
  월 소득이 300만원 인데 월 지출이 175만원일 수 있다.
- 이런 데이터라면 위와 같은 y= 0.5x와 같은 수학 함수식을 만들 수 없다. 그래서 회귀모델은 정확한 함수식이 아닌 모든 데이터를
  가장 만족시킬 만한 최선의 함수식을 만든다.


- y=a + bx는 직관적이고 쉬운데 통계에서 사용하지 않는다. a, b는 통계에서 잘 사용하지 않는 기호 이다.
  그래서 아래와 같이 사용한다.



- 모회귀식은 모집단이며, 표본회귀식은 표본이다. 표본에 의해 회귀모델이 구해지면 회귀모델은 모집단의 특징이 된다.

- 특징이 되었다는 것은 모집단의 데이터가 이 함수식에 의해서 설명되어야 한다 것을 의미한다. 하지만 회귀모델 선 위에
  없는 데이터가 대부분이다. 

- 선위에 없다는 것은 함수식에 의해 설명될 수 없다는 것을 의미 하며, 오차가 생겼다고 할 수 있다. 그 오차는 실험에서
  나와야 할 데이터에서 크게 벗어나지 않는다.



- R이 없었다면 수학적으로 미분을 해서 가장 적합한 함수식을 찾을 수 있지만 R을 이용하면 lm()함수만 실행하면 바로 알 수 있다.

- R에서 회귀모델을 구하는 방법은

 
> lm(mpg~hp, data=mtcars)

 
Call:
lm(formula = mpg ~ hp, data = mtcars)
 
Coefficients:
(Intercept)           hp  
   30.09886     -0.06823  
 


- hp밑의 -0.06823이 b이다.
- Intercept 가
 a이다.
- y= 30.09-0.06x라는 회귀모델식이 만들어진다.






- 위의 그림처럼 회귀모델(파란선)로 부터 데이터가 떨어진 값을 잔차(Residuals)라 한다. 표준편차도 평균으로 부터 떨어진 거리이다.
- 평균처럼 회귀모델도 무게중심으로 역할을 한다. 평균이 모든 값을 끌어안듯이 회귀모델도 모든 데이터를 끌어안아 자신을 나타낸다.
- 회귀모델을 유도 하는 식은 평균이 들어가며, 회귀모델은 평균으로 회귀한다는 의미를 가진다.


(2) R에서 선형회귀모델 실습
1) 안타와 홈런변수를 활용한 회귀분석 예제
- 2015년 프로야구 자료를 가지고 각 팀별 안타횟수에 따른 홈런의 수에 대해 회귀분석
- 자료는 한국 KBO 홈페이지에서 구한 자료이다.
- H-안타, HR-홈런

#데이터 읽기

> df <- read.csv("./data/example_kbo2015.csv")
> str(df)
 
'data.frame': 10 obs. of  26 variables:
 $ ranking: int  1 2 3 4 5 6 7 8 9 10
 $ team   : Factor w/ 10 levels "KIA","kt","LG",..: 9 6 7 4 8 5 10 2 3 1
 $ AVG    : num  0.3 0.3 0.291 0.288 0.277 0.276 0.27 0.27 0.259 0.256
 $ G      : int  102 102 99 101 104 99 102 102 103 100
 $ PA     : int  4091 4151 3950 3994 4082 3888 4037 3988 3984 3835
 $ AB     : int  3553 3620 3410 3485 3557 3373 3404 3498 3491 3362
 $ R      : int  634 658 570 591 545 476 504 466 458 469
 $ H      : int  1066 1085 991 1002 985 931 919 944 905 861
 $ X2B    : int  188 227 176 205 176 146 162 161 171 161
 $ X3B    : int  19 12 16 20 14 9 13 13 16 14
 $ HR     : int  127 155 98 110 130 92 82 87 83 95
 $ TB     : int  1673 1801 1493 1577 1579 1371 1353 1392 1357 1335
 $ RBI    : int  606 622 541 557 515 451 467 443 423 440
 $ SAC    : int  49 40 51 47 50 75 114 62 60 54
 $ SF     : int  44 37 45 36 23 17 27 26 29 23
 $ BB     : int  398 393 378 351 395 360 406 355 347 327
 $ IBB    : int  15 13 17 17 16 12 23 13 17 11
 $ HBP    : int  47 61 66 74 57 63 86 47 57 69
 $ SO     : int  650 815 544 717 866 720 780 799 803 778
 $ GDP    : int  78 79 95 65 105 73 78 75 68 68
 $ SLG    : num  0.471 0.498 0.438 0.453 0.444 0.406 0.397 0.398 0.389 0.397
 $ OBP    : num  0.374 0.374 0.368 0.362 0.356 0.355 0.36 0.343 0.334 0.332
 $ OPS    : num  0.845 0.872 0.806 0.815 0.8 0.761 0.757 0.741 0.723 0.729
 $ MH     : int  101 100 98 101 104 99 102 101 103 100
 $ RISP   : num  0.3 0.294 0.284 0.296 0.272 0.277 0.264 0.259 0.238 0.26
 $ PH.BA  : num  0.216 0.268 0.262 0.274 0.19 0.236 0.196 0.243 0.231 0.214


#안타와 홈런 변수 확인

> df$H
 [1] 1066 1085  991 1002  985  931  919
 [8]  944  905  861
> df$HR
 [1] 127 155  98 110 130  92  82  87  83  95

#두변수간의 상관계수를 구한다.
#꽤 상관관계가 있다고 볼수 있다.


> cor(df$H,df$HR)
[1] 0.8384116

# 두변수간 산점도 그린다.
> plot(HR~H,data=df,pch=20,col='grey',cex=1.5)



#회귀모델을 그린다


> lm<- lm(HR~H,data=df)
> lm
 
Call:
lm(formula = HR ~ H, data = df)
 
Coefficients:
(Intercept)            H  
  -172.3105       0.2871  
 
> abline(lm,lwd=2,col="red")



- 간단하게 두변수를 잘 설명해주는 직선의 방정식, 즉 회귀모델이 만들어졌다.
  그래프에 나타나는 회귀절면(Intercept) -172.3105의 의미는 안타를 한 번도 못쳤을 경우(hit=0)
  홈런을 -172.3105번 친다는 의미로 해석하지만, 안타가 0인 경우를 해석하기 쉽지 않기때문에 
  회귀모델에서 절편의 의미를 크게 해석하지 않는다.

- 회귀분석은 기울기에 조금 더 관심이 있다. 안타를 10번치면 2.8번은 혼런을 친다라는 사실을 알 수 있다.

(3) 선형회귀모델 검증하기
- b0(절편), b1(기울기)인 회귀계수를 검증한다 특히 기울기(b1)를 검증하는 것이 중요하다.
- b1이 0일 가능성은 없는지, b1의 질은 좋은지를 살펴보는 것이 회귀모델을 검증하는 방법이 된다.

1) 결정계수

- 결정계수 R²은 회귀모델이 질을 나타내는 중요한 값 중에 하나이다. 
R² = 1  : 모든 값이 회귀모델 위에 있다. R² 이 작아질수록 회귀모델에서 멀어진다.
R²는 상관계수를 제곱해서 얻은 값이다. 만약 상관계수가 0.7이면 R²은 0.49로 낮아진다.
  그만큼 회귀모델에 적합하려면 상관계수가 높아야 한다. 0.9일때 
R² 는 0.81가 된다.

2) R에서 결정계수 보기
- R에서는 summary()함수로 lm 객체를 보면 알수 있다.


> summary(lm)
 
Call:
lm(formula = HR ~ H, data = df)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-14.246  -8.874  -4.978  11.068  20.082 
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -172.3105    64.0968  -2.688  0.02757 * 
H              0.2871     0.0660   4.351  0.00244 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 14 on 8 degrees of freedom
Multiple R-squared:  0.7029, Adjusted R-squared:  0.6658 
F-statistic: 18.93 on 1 and 8 DF,  p-value: 0.002442
 
- 노란마킹부분이  R² 이다


3) 선형 회귀모델 가설 검증

- 모집단의 모회귀식  Y =  
βo + β ₁ χ 위에 없는 점은 모집단 관점에선 오차라 할 수있다.
  오차가 신뢰구간 안에 있는지 검증함으로써 모집단의 특징이라는 회귀모델을 검증한다.

- 결정계수가 0.45로 애매한 값이 나오면 F 통계량으로 회귀모델을 판단한다.
- F통계량으로 회귀모델을 검증할 때는 모집단을 검증하므로 
 β ₁ 이 0일 확률을 검증한다.
 
 β ₁이 0이면 βo만 남고 함수식은 상수가된다.
 
β ₁가 0이 아니어야 회귀모델은 의미를 가진다.

- F통계량은  
β ₁이 0일 가능성을 보여주는 값으로 0.05보다 낮으면 β ₁는 0일 가능성이 거의 없게 되므로 
  0.05보다 낮아야 한다.

summary(lm)의 맨아래 p-value값이 0.005보다 작으면 
β ₁이 0일 가능성이 낮아지면서 이 회귀모델은 적합하다고
  통계적으로 말할 수 있다. 위 푸른색 마킹된 부분이다.


β ₁의 신뢰구간을 알고 싶다면 아래와 같이 입력한다.
- 0.2871은 신뢰구간에 잘 들어가 있다.


> confint(lm)
                   2.5 %     97.5 %
(Intercept) -320.1180464 -24.502900
H              0.1349531   0.439328


4) 잔차 관련된 그래프를 그려 회귀적합성 검증하기
- 회귀적합성을 좀더 정확히 검증하는 방법은 '잔차와  관련된그래프를 그려 확인하는 것이다.
- R에서 plot(lm)으로 잔차와 관련된 그래프를 볼 수 있다.



- 좌측상단에 있는 결과물은 0을 기준으로 골고루 퍼져 있어야 좋다. 그러나 그래프가 +,-를 왔다갔다하거나
  특정 패턴있으면 이런 경우 회귀적합성이 없다고 봐도 된다.

- 우축상단에 있는 결과물은 잔차에 대한 정규분포의 분위수-분위수 그래프(QQplot)이라고 부르며 데이터의 분포가
  정규성을  만족하는지 못하는지 대략적으로 파악하기 위한 도구로서 사용된다. 정규분포 형태라면 점선위에 점들이 놓여있게된다.

- 좌측하단 그림은 잔차에서 개선된 표준화된 잔차를 제곱해 표현한 것으로 잔차의 정도를 파악하는데 조금 더 용이한 결과물이다.
  골고루 퍼져 있는 그래프가 좋다.

- 우측하단에 있는 그림은 영향점, 이상점 판별을 도와 주는 그래프로 cook's distance는 통상 1값이 넘어가면 그 관측치를 영향점
  으로 판별한다. 이상점은 등고선 안쪽에 동그라미가 위치하게되며 이상점이 있으면 회귀모델에 좋치 않은 영향을 미친다.


 




(4) 다중 선형회귀분석
- 종속변수에 영향을 미치는 독립변수가 한개 이상인 경우 다중선형회부 모델이라 한다. 
- R에서 lm()를 사용하며 '+'로 변수를 추가하면 된다.
  lm(HR~R+H+AO)

1) 2015 KBO 야구 데이터 다중선형회귀분석
- RBI(타점), XBH(장타), TB(루타), XR(추정득점), GW.RBI(결정타)등을 독립변수로 사용한다.

> df2 <- read.csv("./data/example_kbo2015_player.csv",stringsAsFactors = F,na='-')
> str(df2)
'data.frame': 337 obs. of  26 variables:
 $ 순위  : int  1 2 3 4 5 6 7 8 9 10 ...
 $ 선수명: chr  "유한준" "박병호" "홍성갑" "고종욱" ...
 $ 팀명  : chr  "넥센" "넥센" "넥센" "넥센" ...
 $ AVG   : num  0.364 0.349 0.333 0.324 0.318 0.316 0.311 0.3 0.295 0.287 ...
 $ G     : int  101 106 7 85 44 96 68 99 79 103 ...
 $ PA    : int  436 476 10 328 26 405 253 389 339 436 ...
 $ AB    : int  376 410 9 299 22 370 219 333 295 383 ...
 $ R     : int  82 100 0 59 17 53 40 52 62 68 ...
 $ H     : int  137 143 3 97 7 117 68 100 87 110 ...
 $ X2B   : int  35 29 0 22 2 21 10 22 21 28 ...
 $ X3B   : int  1 1 0 3 0 0 0 1 1 4 ...
 $ HR    : int  19 42 0 8 0 14 8 13 17 14 ...
 $ TB    : int  231 300 3 149 9 180 102 163 161 188 ...
 $ RBI   : int  84 111 2 40 3 65 21 67 53 57 ...
 $ SAC   : int  0 0 0 3 2 0 1 5 3 3 ...
 $ SF    : int  6 3 0 1 0 3 2 6 3 4 ...
 $ XBH   : int  55 72 0 33 2 35 18 36 39 46 ...
 $ GO    : int  97 70 2 73 6 95 60 96 39 80 ...
 $ AO    : int  95 74 0 67 2 105 57 84 75 106 ...
 $ GO.AO : num  1.02 0.95 NA 1.09 3 0.9 1.05 1.14 0.52 0.75 ...
 $ GW.RBI: int  7 10 0 6 0 2 1 5 1 3 ...
 $ BB.K  : num  1 0.42 0.25 0.3 0.29 0.45 0.78 0.63 0.31 0.46 ...
 $ P.PA  : num  4 4.07 4.7 3.89 4.08 3.98 4.19 3.85 4.14 4.22 ...
 $ ISOP  : num  0.25 0.383 0 0.174 0.091 0.17 0.155 0.189 0.251 0.204 ...
 $ XR    : num  88.4 113.3 1.3 47.9 4.2 ...
 $ GPA   : num  0.351 0.378 0.263 0.293 0.271 0.287 0.293 0.292 0.304 0.285 ...

#문자열 변수를 실수로 변환

 
> df2$AVG <- as.numeric(df2$AVG)
> df2$GO.AO <- as.numeric(df2$GO.AO)
> df2$BB.K <- as.numeric(df2$BB.K)
> df2$P.PA <- as.numeric(df2$P.PA)
> df2$ISOP <- as.numeric(df2$ISOP)

#홈런과 다른변수들간의 상관계수를 살펴본다
# use= "pairwise.complete.obs"
#
상관계수가 계산되는 변수들만을 대상으로 결측값이 있는 case 제거한  상관계수 계산



> cors <- cor(df2$HR, df2[,5:length(df2)],use='pairwise.complete.obs')
> cors
             G        PA        AB
[1,] 0.6573121 0.7607232 0.7519703
             R         H       X2B
[1,] 0.8012273 0.7815465 0.7756991
           X3B HR        TB       RBI
[1,] 0.3089791  1 0.8828075 0.9217014
            SAC        SF       XBH
[1,] 0.02882778 0.6766359 0.9190164
            GO        AO      GO.AO
[1,] 0.6096195 0.7477335 -0.1618067
        GW.RBI      BB.K       P.PA
[1,] 0.7889489 0.2839297 0.03248409
          ISOP        XR       GPA
[1,] 0.7813991 0.8815765 0.5505212

#정렬해서 상관관계를 다시 본다.

 
> cors <- cors[,order(cors)]
> cors
      GO.AO         SAC        P.PA 
-0.16180671  0.02882778  0.03248409 
       BB.K         X3B         GPA 
 0.28392966  0.30897909  0.55052123 
         GO           G          SF 
 0.60961949  0.65731209  0.67663588 
         AO          AB          PA 
 0.74773350  0.75197032  0.76072321 
        X2B        ISOP           H 
 0.77569907  0.78139909  0.78154650 
     GW.RBI           R          XR 
 0.78894892  0.80122734  0.88157647 
         TB         XBH         RBI 
 0.88280745  0.91901637  0.92170137 
         HR 
 1.00000000 

# 독립변수의 값이 0인 관측치를 제거한다.
> df2$HR[df2$HR==0]<-NA
> df2$RBI[df2$RBI==0]<-NA
> df2$XBH[df2$XBH==0]<-NA
> df2$TB[df2$TB==0]<-NA
> df2$XR[df2$XR==0]<-NA
> df2$GW.RBI[df2$GW.RBI==0]<-NA
 
> df2 <-na.omit(df2) #결측치 제거

#여러개의 독립변수로 회귀분석을 한다.

 
> lm2<-lm(HR~RBI+XBH+TB+XR+R, data=df2)
> lm2
 
Call:
lm(formula = HR ~ RBI + XBH + TB + XR + R, data = df2)
 
Coefficients:
(Intercept)          RBI          XBH  
   -1.08756      0.15104      0.51513  
         TB           XR            R  
   -0.07797      0.14409     -0.12632  

 
#회귀모델의 검증을 한다.
#Pr(>|t|) 의 값이 0.05보다 작으므로 0일 가능성이 매우 낮다.
#p-value, 결정계수등 모두 좋다.

 
> summary(lm2)
 
Call:
lm(formula = HR ~ RBI + XBH + TB + XR + R, data = df2)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-6.8941 -1.5712  0.2196  1.5315 10.5936 
 
Coefficients:
            Estimate Std. Error t value
(Intercept) -1.08756    0.62254  -1.747
RBI          0.15104    0.03954   3.820
XBH          0.51513    0.08172   6.304
TB          -0.07797    0.03510  -2.221
XR           0.14409    0.07046   2.045
R           -0.12632    0.04309  -2.932
            Pr(>|t|)    
(Intercept) 0.083267 .  
RBI         0.000215 ***
XBH         5.32e-09 ***
TB          0.028260 *  
XR          0.043094 *  
R           0.004057 ** 
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 3.027 on 117 degrees of freedom
Multiple R-squared:  0.8637, Adjusted R-squared:  0.8579 
F-statistic: 148.3 on 5 and 117 DF,  p-value: < 2.2e-16

댓글