선형 회귀 분석 (3)

 앞서의 예제에서 SSE/SSR/SST 를 구하는 것은 간단합니다. 그냥 공식에 대입하면 되는 것이죠. 공식에 대해서는 통계학 서적을 참조해 주시기 바랍니다. 

> yhat=beta0+beta*x

> SST=sum((y-ymean)^2)

> SSR=sum((yhat-ymean)^2)

> SSE=sum((y-yhat)^2)

> SST;SSR;SSE

[1] 9174.896

[1] 8314.219

[1] 860.6765

 SSE의 경우 더 간단하게 구할 수도 있습니다. 

> deviance(model) #SSE

[1] 860.6765

> 

 하지만 사실 SSE/SSR/SST 값은 통계 학습 시간에는 중요하게 다루더라도 실제 분석 결과로 나타내는 경우는 드물 것입니다. 주로 논문에서 결과값으로 표기하는 것은 회귀계수인 베타 값과 R^2 값, P 값일 것입니다. F 통계량을 포함 이 값들을 구하는 것 역시 공식에 대입해 쉽게 구할 수 있습니다. 

> R2=SSR/SST 

> Fval=SSR /(SSE/(length(y)-2)) 

> Pval=1-pf(Fval,1,(length(y)-2)) 

> c(R2,Fval,Pval)

[1]   0.9061922 463.6847041   0.0000000

 R2가 R^2 값, Fval 가 F 통계량, Pval는 P 값입니다. 여기서 주의할 점은 P값은 0이 될 수 없으므로 P<0 .001="" nbsp="" span="" summary="">

> set.seed(1234)

> x<-rnorm span="">

> set.seed(4567)

> b<-rnorm span="">

> y=3*x+b

> model=lm(y~x)

> summary(model)

Call:

lm(formula = y ~ x)

Residuals:

    Min      1Q  Median      3Q     Max 

-9.4838 -2.9000 -0.0445  2.4240  7.9616 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept)  26.6162    13.3737    1.99   0.0523 .  

x             2.9436     0.1367   21.53   <2e-16 span="">

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.234 on 48 degrees of freedom

Multiple R-squared:  0.9062, Adjusted R-squared:  0.9042 

F-statistic: 463.7 on 1 and 48 DF,  p-value: < 2.2e-16

 p-value: < 2.2e-16라는 표현은 P값이 매우 작다는 의미로 통상적인 논문에서는 P<0 -16="" .001="" 10="" e-16="" nbsp="" p="" r="" span="">

 summary로 주어지는 결과값 해석에서 중요한 것은 Std. Error 입니다. 모델에서 예측된 Y의 값은 실제 값과 근사할 뿐이지 동일하지는 않습니다. 모델에서 예측한 Y값의 분포는 정규 분포를 따르는 추정값입니다. 따라서 그 표준오차나 신뢰구간을 같이 표기해 주어야 합니다. 각각의 값을 보기 위해서는 R에서 predict 함수를 사용할 수 있습니다. predict(model, se.fit = TRUE)로 지정합니다. 

> predict(model, se.fit = TRUE)

$fit

       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18 

303.2101 325.0589 336.9365 286.4518 327.2915 328.4238 312.5167 312.9304 312.6681 307.8761 313.9524 306.2815 309.5508 321.9244 335.0975 319.3525 313.4547 307.5647 

      19       20       21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36 

308.6542 356.5319 322.9492 313.7538 314.4917 327.7399 310.7655 299.6610 329.4349 305.9095 320.7529 307.2004 337.1993 313.9759 310.5342 313.5982 296.9987 303.7907 

      37       38       39       40       41       42       43       44       45       46       47       48       49       50 

288.8899 301.2390 316.6443 314.1186 342.3093 305.2474 308.3864 316.8455 306.3410 306.7211 304.6782 302.5490 313.2660 313.6631 

$se.fit

 [1] 0.7900740 0.7796783 1.2095243 1.4255037 0.8498136 0.8878956 0.6045932 0.6022509 0.6036660 0.6691963 0.5990721 0.7053692 0.6382907 0.6955081 1.1361123 0.6430417

[17] 0.6001537 0.6757746 0.6538689 2.0502812 0.7208775 0.5993971 0.5989059 0.8647123 0.6210265 0.9062366 0.9231075 0.7146463 0.6694802 0.6837787 1.2201431 0.5990430

[33] 0.6239585 0.5997506 1.0023683 0.7727572 1.3236075 0.8526127 0.6085967 0.5989092 1.4316499 0.7318804 0.6589602 0.6103320 0.7039117 0.6947964 0.7473906 0.8104356

[49] 0.6007956 0.5995926

$df

[1] 48

$residual.scale

[1] 4.234473

 값이 많아서 눈에 잘 들어오지 않을 수 있지만, Yhat에 해당하는 값과 표준편차를 보여주고 있습니다. (예를 들어 303.2101와 0.7900740) 신뢰구간을 보고자 한다면 interval 인자를 주면 됩니다. lwr가 신뢰구간의 하한값, upr가 상한값이며 기본은 95%입니다. 

> predict(model, interval =  "confidence")

        fit      lwr      upr

1  303.2101 301.6216 304.7987

2  325.0589 323.4912 326.6265

3  336.9365 334.5046 339.3684

4  286.4518 283.5856 289.3179

5  327.2915 325.5829 329.0002

6  328.4238 326.6386 330.2090

7  312.5167 311.3011 313.7323

8  312.9304 311.7195 314.1413

9  312.6681 311.4543 313.8818

10 307.8761 306.5306 309.2216

11 313.9524 312.7479 315.1569

12 306.2815 304.8632 307.6997

13 309.5508 308.2674 310.8342

14 321.9244 320.5260 323.3228

15 335.0975 332.8132 337.3818

16 319.3525 318.0596 320.6454

17 313.4547 312.2480 314.6613

18 307.5647 306.2060 308.9235

19 308.6542 307.3395 309.9689

20 356.5319 352.4095 360.6543

21 322.9492 321.4998 324.3986

22 313.7538 312.5486 314.9590

23 314.4917 313.2875 315.6959

24 327.7399 326.0013 329.4785

25 310.7655 309.5169 312.0142

26 299.6610 297.8389 301.4832

27 329.4349 327.5789 331.2910

28 305.9095 304.4727 307.3464

29 320.7529 319.4068 322.0990

30 307.2004 305.8256 308.5752

31 337.1993 334.7460 339.6525

32 313.9759 312.7715 315.1804

33 310.5342 309.2796 311.7887

34 313.5982 312.3923 314.8041

35 296.9987 294.9833 299.0141

36 303.7907 302.2370 305.3444

37 288.8899 286.2286 291.5512

38 301.2390 299.5247 302.9533

39 316.6443 315.4206 317.8679

40 314.1186 312.9144 315.3228

41 342.3093 339.4308 345.1879

42 305.2474 303.7759 306.7190

43 308.3864 307.0615 309.7114

44 316.8455 315.6183 318.0726

45 306.3410 304.9257 307.7563

46 306.7211 305.3241 308.1181

47 304.6782 303.1755 306.1809

48 302.5490 300.9195 304.1785

49 313.2660 312.0580 314.4740

50 313.6631 312.4575 314.8686

 여기까지 내용은 좀 복잡할 수 있는데, 다행히 논문에 표기하는 것은 추정된 회귀계수와 표준 편차, P값 정도면 적당합니다. 그리고 모델의 예측력을 확인하기 위해 Adjusted R-squared를 추가로 기술하는 것인데, 이 내용은 나중에 예제를 통해 다시 알아보겠습니다. 

 아무튼 여기까지 보면 summary 기능이 왜 유용한지 알 수 있습니다. summary에 나오는 내용도 복잡하긴 한데, 그래도 위의 내용들을 모두 요약해서 보여주고 있기 때문이죠. 실제로 분석을 할 때는 summary를 중심으로 진행하게 되며 여기서 설명한 내용은 각각의 값을 눈으로 확인하고 싶을 때만 사용하게 될 것입니다. 

 그래도 여기까지는 사실 난이도가 높지 않은 부분입니다. 정말 곤란한 부분은 잔차부터 시작이라고 할 시 있을 것 같습니다.