[회귀] Lack of fit test with R/ Simple Linear regression/ 적합 결여 검정 / 단순 선형 회귀

<적합 결여 검정>

 

설명변수와 반응변수의 관계를 적절히 설명하는 회귀식을 찾았다면, 이 회귀모형이 올바른 모형인지 확인 할 수 있어야 한다. 이에 대한 검정을 '적합 결여 검정' (lack of fit test) 라고 한다. 

 

앞선 포스팅 https://jangpiano-science.tistory.com/61 에서 살펴 보았듯이, MSE(Mean Square error)는 추정량의 분산과 편향(bias)의 제곱에 의해 결정된다. 따라서, MSE의 크기는 잘못된 회귀모형으로 인한 편향에 영향을 받는다. 이 특징으로 인해, MSE는 회귀모형이 관측값을 잘 설명하는 모형인지 검토하는데 유용하게 쓰이는 요소이다. 

 

적합 결여 검정에 대한 자세한 설명을 시작하겠다. 

앞선 포스팅에서도 언급하였듯이, 추정 회귀 모형이 실제 관측값(실제 모형)을 잘 설명할때, MSE는 분산의 비편향 추정량이 된다. 

 

즉, 주로 분산의 추정값으로 쓰이는 MSE는 모형의 적합도에 의존하여, 적합한 회귀식을 사용할때는 MSE와 분산의 차이가 크지 않고, 반면에 적합하지 못한 회귀식을 사용 할 경우, 높은 편향에 의해 MSE와 분산의 차이가 커진다. 

여기서, MSE를 비교하는 값이 모분산이 되면 좋겠지만, 실제로 회귀분석을 진행할때, 분산값이 알려져 있는 경우는 매우 드물다.

따라서, 분산이 알려져 있는 경우는 알려져있는 분산과 관측값을 이용하여 구한 MSE를 비교하여 추정회귀식이 적합한지 검정할 수 있다. 

보편적인 경우인, 분산이 알려져 있지 않는 경우에 대한 검정 방법도 밑에서 설명하도록 하겠다.

 

귀무가설과 대립가설을 다음과 같이 설정 하였으므로, 본격적으로 적합 결여 검정을 수행하여 보자. 

H0 (null hypothesis) : 주어진 모형은 관측값을 설명하는데 적합하다. 

H1 (alternative hypothesis) : 주어진 모형은 관측값을 설명하는데 적합하지 않다. 

 

본격적으로, 가설검정을 수행해보자. 귀무가설(null hypothesis) 이 '주어진 모형은 올바르다' 로 설정됨으로 인해, 귀무가설을 기각하게 되는 경우, 관측값을 설명하기에 적합한 새로운 모형을 만들어야 한다. 

 

<σ^2 이 알려져 있을때> 

 

SSE와 분산의 관계를 설명하는, 자유도가 n-p인 카이제곱 분포의 성질로 인해, 적합 결여 검정은 카이 제곱 검정(Chi -square test)에 의해 수행 될수 있다. (단순 회귀의 경우 모수의 수는 B0, B1 두개이므로 자유도가 n-2가 된다, 이에 대한 증명은 아래와 같다 ) 

 

 

 

위 단순 선형 회귀 모형을 일반화 시키자면, 이러하고, p 는 모수의 개수를 의미하고, n 은 관측값의 개수를 의미한다. 

즉, 모수가 p 개인 회귀모형을 검정할때는, 자유도가 n-p인 카이 제곱 검정을 필요로 한다. 

 

(모수가 p 개인 경우에 주로, B0, B1, ..., B(p-1) 과 B1, B2, ..., B(p) 모수의 개수를 셀때 헷갈리지 않도록 주의하세요.)

 

 

분산이 알려져 있는 경우, SSE/σ^2 값이 자유도를 n-p로 하는 카이제곱 분포의 기각영역에 들어가게 되면, H0을 기각하고, 그렇지 않은 경우, H0 을 기각하지 않는다. 

 

또한 , 카이제곱 검정통계량 (test statistics) 보다 극단적인 값을 얻게 될 확률인 p-값이 유의수준(significance level) α 보다 작게 나온 경우에도 H0을 기각한다. 

 

이렇듯, 카이 제곱 검정을 하기 위해서는, 모집단의 분산의 값을 알아야 한다. 하지만, 모집단의 분산의 값이 알려져 있는 경우는 매우 드물고, 따라서 분산이 알려지지 않은 경우를 대비해 회귀 모형에 의존하지 않는 분산의 추정량이 필요하다. 

*MSE는 회귀 모형의 적합도에 의존함. 

 

 

<σ^2 이 알려져 있지 않을때>

 

분산값이 알려져 있지 않는경우, 관측값이 주어져 있을때, 동일한 설명변수의 값에 대하여 반응변수의 값이 측정되었다면, 이러한 반응변수 값들을 이용해 분산을 구하는 방법이 보편화 되어 있다. 이는 모형에 의존하지 않는 분산 추정방법이다. 모수의 분산이 알려지지 않았고, 동일한 설명변수에 대한 반응변수의 값이 관측되지 않는다면, 모형에 의존하지 않는 분산을 추정값을 구할 수 없기 때문에, 단호하게 적합결여 검정이 불가능하다고 얘기할 수 있다. 

 

설명변수에 의해 구분되는 그룹 : m 

각각의 설명변수 값에서 측정되는 반응변수의 관찰값 : ni (i = 1, 2, ..., m)

 

즉 i 그룹 안에서의 표본분산을 모두 구하여 더한 값은 순수오차 (pure error) 로 SSpe (Sum of squares due to pure error)로 표기 된다.  이 순수오차값은 자유도로 나눈 값이 분산의 새로운 추정값이 된다. 

 

 

 

MSE값이 새롭게 추정된 분산값보다 크지 않으면, 가정한 모형은 적합성이 좋다고 해석할 수 있고, MSE와 분산의 추정값의 차이가 크다면, MSE의 편향값이 크다는 말이므로, 적합성이 의심된다. 

 

또다른, 검정법으로는, 순수 오차(SSpe : Sum of squares due to pure error) 와 '적합 결여에 따른 제곱합 (SSlof : Sum of square error due to lack of fit)' 을 비교 검정하는 F 검정에 의해 이루어질 수 있다. 적합 결여에 따른 제곱합이란, 오차제곱합(SSE: Sum of square error) 과 순수 오차(SSpe : Sum of squares due to pure error) 의 차이로 계산되는, 순수오차로 설명될 수 없는 SSE의 부분으로, 편향을 야기하는 부분이라고 할 수 있다.  

 

 

F 검정통계량 (F - test statistics) 는 다음과 같이 계산된며, m-p 과 n-m 를 df1, df2 로 가진다. 

 

 

따라서, F 검정 통계량 값이 자유도를 m-p, n-m 으로 하는 F 분포의 기각영역에 들어가게 되면, H0을 기각하고, 그렇지 않은 경우, H0 을 기각하지 않는다. 

 

또한 , F 검정통계량 (test statistics) 보다 극단적인 값을 얻게 될 확률인 p-값이 유의수준(significance level) α 보다 작게 나온 경우에도 H0을 기각한다. 

 

H0을 기각하게 되는 경우, 관측값을 적합하게 설명 가능한 새로운 회귀모델을 만들어야 하는데, 주로 설명변수나 반응변수를 변환시키거나, 새로운 설명변수를 추가함으로써 적합한 추정 모델을 찾을 수 있다. 

 

 

<R > 

대부분의 경우에, 모수의 분산이 알려져 있지 않은 경우가 많으므로, 모수의 분산이 알려지지 않은경우를 다루는 R 코드를 작성해 보았다.

 

밑의 예시는, 동일한 설명변수(ex, 1) 에 대해 반복적인 반응변수 (ex, 2,4) 가 나타나므로, 아래 예시를 가지고 분산을 추정할 수 있게 된다. 

 

> X = c(1,1, 3,3,3, 4, 5,5, 7,7,7,7, 8)

> Y = c(2,4, 2,5,3, 6, 8,7, 11,13,10,10, 15)

 

> data.frame(X,Y)

 X  Y

 1  2

 1  4

 3  2

 3  5

 3  3

 4  6

 5  8

 5  7

 7 11

 7 13

 7 10

 7 10

 8 15

 

> a = lm(Y ~ X)

> a

 

Call:

lm(formula = Y ~ X)

 

Coefficients:

(Intercept)            X  

    -0.2942       1.6365  

 

 

> b = lm(Y ~ factor(X))

 

> anova(a, b) 

Analysis of Variance Table

 

Model 1: Y ~ X

Model 2: Y ~ factor(X)

  Res.Df    RSS       Df     Sum of Sq      F          Pr(>F)

1     11        28.912                           

2      7       13.167    4    15.745          2.0927      0.1852     

               

          #SSpe 는 13.167 이므로 새로운 분산의 추정량은 13.167 / 7 이다. 

          # F - statistics = 2.0927이다. 

          # p-value 가 0.1852 로, 충분히 작지 않음으로, H0을 기각하지 않는다. 

 

 

 

SSlof 를 구하기 위해선, SSE 가 필요하다. 

 

> anova(a)

Analysis of Variance Table

 

Response: Y

                Df  Sum Sq    Mean Sq  F value      Pr(>F)    

X               1     184.165 184.165     70.069     0.000004235 ***

Residuals 11     28.912     2.628                        

---

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

 

                      #SSlof = SSE - SSpe =  28.912 - 13.167 = 15.745

 

<function 만들기> 

 

lack_of_fit_test = function(X, Y, a){

  

  t = lm(Y~X)

  

  estimator_B0 = t$coefficients[[1]]

  estimator_B1 = t$coefficients[[2]]

  SSE = anova(t)[2,2]

  

  print(paste("estimator for B0 : ", estimator_B0 , "estimator for B1 : ", estimator_B1))

  print(paste("estimator for Yi : ", estimator_B0,"+",estimator_B1,"x" ))

  print(paste("SSE : " , SSE))

  

  fit_test = lm(Y ~ factor (X))

  SSpe = anova(t, fit_test)[2,2]

  SSlof = SSE - SSpe

  print(paste("SSpe: ", SSpe ))

  print(paste("SSlof : ", SSlof ))

  

  estimator_var = anova(t ,fit_test)[2,2]/anova(t, fit_test)[2,1]

  print(paste("estimator for variance : ", estimator_var ))

  

  p_value = anova(t, fit_test)[[6]][2]

  

  if (p_value >= a){

    print("the model is a good fit model")

    }

  else{

    print("the model is not a good fit model")

    }

  

  }

 

> X = c(1,1, 3,3,3, 4, 5,5, 7,7,7,7, 8)

> Y = c(2,4, 2,5,3, 6, 8,7, 11,13,10,10, 15)

 

> lack_of_fit_test(X, Y, 0.05)

[1] "estimator for B0 :  -0.294183445190157 estimator for B1 :  1.63646532438479"

[1] "estimator for Yi :  -0.294183445190157 + 1.63646532438479 x"

[1] "SSE :  28.9116331096197"

[1] "SSpe:  13.1666666666667"

[1] "SSlof :  15.744966442953"

[1] "estimator for variance :  1.88095238095238"

[1] "the model is a good fit model"