[R] 카이제곱 적합도 검정(Chi-squared goodness of fit test)으로 특정 분포인지 확인하기

이번 포스팅에서는 카이제곱 적합도 검정(Chi-squared goodness of fit test)를 활용하여 관측치가 특정 분포를 따르는지 아닌지를 검정하는 방법을 소개하겠습니다.

 

다음과 같은 문제가 있다고 해봅시다.

 

[문제]  시간이 지남에 따라 일어난 어떤 특정 사건의 발생건수에 대한 결과가 다음과 같이 주어졌다고 하자. 그 때 특정 사건의 발생 회수가 포아송(λ) 분포를 따르는지에 대한 검정을 유의수준 5%에서 실시하시오.

 

발생건수	0회	1회	2회	3회	4회	5회	6회	7회	8회	9회 이상
관찰회수	6	20	36	50	52	40	25	12	4	5

 

 

(풀이)

특정 관측치들이 어떤 이론적 분포를 따르고 있는지를 검정하는 데에는 적합도 검정(Chi-squared goodness of fit test) 를 이용합니다.

 

검정해야 할 귀무가설(H0)와 대립가설(H1)을 세워보면,

• 귀무가설 : 특정 사건의 발생 회수가 포아송 분포(Poisson(λ))를 따른다.
                   (즉, 관측값의 도수와 가정한 분포로부터의 기대관측도수가 같다)
• 대립가설 : 특정 사건의 발생 회수가 포아송 분포를 따르지 않는다. (Not )

 

이며,

카이제곱 검정 통계량은  입니다. (이 문제에서는 i=0, 1, ..., 9)

이때 k 는 범주(class)의 개수를 나타내며, 는 도수분포의 각 구간에 있는 관측도수,

는 도수분포의 각 구간의 기대도수를 나타냅니다.

 

 

주어진 데이터로부터 가정된 포아송 분포의 평균(λ)의 추정치을 구하면, =3.844가 됩니다.

 

= (0*6+1*20+2*36+3*50*4*52+5*40+6*25+7*12+8*4+9*5)/250 = 3.844

 

R 로 위의 관측치 데이터셋을 입력하고 가정된 포아송 분포의 평균 추정치를 구해보도록 하겠습니다.

 

> #=========================================== > # Chi-squared Goodness-of-fit Test > # Poisson Distribution, 5% significant level > #=========================================== > rm(list=ls()) # clear all > > # Class of Events > classnum <- c(0:9) > > # Observed Frequency > obsnum <- c(6, 20, 36, 50, 52, 40, 25, 12, 4, 5) > > # Sum of Observed Frequency > obssum <- sum(obsnum); obssum [1] 250 > > # Calculate Lambda > lamb <- sum(classnum*obsnum) / sum(obsnum); lamb [1] 3.844

 

 

포아송분포의 모수 λ에 대한 추정치로 3.844를 사용하여 각 class 별 가설적인 확률 (i=0, 1, ..., 9)를 구해보겠습니다.

 

에 i=0, 1, ..., 9 를 대입해서 풀어보면 아래의 표와 같이 가설적인 확률 를 구할 수 있습니다.

 

R 의 dpois() 함수를 사용하여 확률 구해봤는데요, sum(hyp_prob)=0.994로서 전체 합이 1이 안되고 0.994 이네요.

 

> # Calculate the hypothesised probabilities > hyp_prob = round(dpois(classnum, lamb), 3); hyp_prob [1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.011 > sum(hyp_prob) # not 1 [1] 0.994

 

확률의 공리적 정의 상 표본공간의 P(S)=1 이어야 하므로 모든 class별 확률을 전부 더했을 때 1이 되도록 보정을 해줄 필요가 있습니다. 만약 보정을 안한 상태에서 바로 카이제곱 검정을 실시하면 아래와 같이 "확률값들의 총계는 반드시 1 이어야 합니다"와 같은 에러 메시지가 뜹니다.

 

> # Chi-squared goodness-of-fit test > chisq.gof.test <- chisq.test(obsnum, p=hyp_prob) # error Error in chisq.test(obsnum, p = hyp_prob) : 확률값들의 총계는 반드시 1 이어야 합니다

 

 

R의 indexing을 사용하여 가장 마지막인 '9회 이상' class의 확률 를 조정하여 전체 class별 확률의 합이 1이 되도록 조정 해보겠습니다.

 

> # Change the last hypothesised probability to make sum(hyp_prob)=1 > hyp_prob[length(hyp_prob)] <- 1-sum(hyp_prob[1:(length(hyp_prob)-1)]); hyp_prob [1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.017 > sum(hyp_prob) # 1 [1] 1

 

 

발생확률 Pi를 구했으므로 전체 발생회수 합인 250과 각 class별 가정된 포아송분포로 부터의 확률을 곱한 250*로 부터 기대도수(Expected frequency)를 구할 수 있습니다.

 

발생건수	0회	1회	2회	3회	4회	5회	6회	7회	8회	9회 이상	계
발생확률	0.021	0.082	0.158	0.203	0.195	0.150	0.096	0.053	0.025	0.017	1
기대도수 (Ei=Pi*250)	5.25	20.50	39.50	50.75	48.75	37.50	24.00	13.25	6.25	4.25	250
관측도수 (Oi)	6	20	36	50	52	40	25	12	4	5	250

 

 

이제 카이제곱 적합도 검정(Chi-squared goodness of fit test)을 할 준비가 되었네요. 가정된 포아송분포(λ=3.844)로 부터의 기대도수(Expected frequency)와 250개 샘플 데이터로 부터 실제 관측한 도수(Observed frequency) 간에 차이가 있다고 볼 수 있는 것인지, 관측치가 포아송분포를 따른다고 볼 수 있는 것인지 가설을 5% 유의수준 하에서 검정해보겠습니다.

 

R의 chisq.test() 함수를 사용하여 카이제곱 적합도 검정을 수행한 결과는 아래와 같습니다.

 

> # Chi-squared goodness-of-fit test > chisq.gof.test <- chisq.test(obsnum, p=hyp_prob) Warning message: In chisq.test(obsnum, p = hyp_prob) : 카이제곱 approximation은 정확하지 않을수도 있습니다 > chisq.gof.test Chi-squared test for given probabilities data:  obsnum X-squared = 1.9258, df = 9, p-value = 0.9926

 

검정통계량인 X-squared 통계량이 1.9258이고 P-value 가 0.9926 이므로 유의수준 5% 하에서 관측값이 포아송 분포를 따른다는 귀무가설을 채택(Accept Null Hyperthesis)합니다.

 

 

R로 카이제곱 적합도 검정을 한 결과를 chisq.gof.test 라는 이름의 객체로 저장을 했는데요, names() 함수로 확인을 해보면 chisq.gof.test 객체 리스트 안에 여러 분석 결과 통계량들이 들어있음을 알 수 있습니다.

 

> names(chisq.gof.test) [1] "statistic" "parameter" "p.value" "method" "data.name" [6] "observed" "expected" "residuals" "stdres"

 

 

이들 중에서 index을 해서 관측도수(Observed frequency)와 기대도수(Expected frequency)를 확인해보겠습니다.

 

> chisq.gof.test$observed [1]  6 20 36 50 52 40 25 12  4  5 > chisq.gof.test$expected [1]  5.25 20.50 39.50 50.75 48.75 37.50 24.00 13.25  6.25  4.25 > round(hyp_prob, 3) [1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.017

 

 

카이제곱 검정을 할 때 각 class의 기대도수가 너무 작을 경우 검정이 유효하지 않을 수가 있습니다. 위에 R분석 결과에서 빨간색으로 "카이제곱 approximation은 정확하지 않을 수도 있습니다"라는 warning message가 나온 이유가 '9회 이상' class의 기대도수가 4.25로서 작기 때문입니다. 최소 기대도수에 대해서는 일반적인 합의는 없으며, 보통 3, 4, 5개가 종종 사용되며, 3~5개 이하일 경우 인접 class와 병합해서 재분석을 하기도 합니다.

 

이번 문제에서는 최소 기대도수가 4.25이고, 또 P-value가 0.9926으로서 매우 높은 수준으로 귀무가설을 채택하였으므로 '9회 이상' class를 '8회'와 합쳐서 재분석하는 것은 필요없어 보입니다. (병합해서 재분석해도 결과는 동일하게 귀무가설 채택함)

 

마지막으로, 각 class별 관측도수와 기대도수를 그래프로 그려서 시각화로 비교를 해보겠습니다.

 

> # Plot of Poisson Distribution: Observed vs. Expected > # Observed Frequency > plot(c(0:9), chisq.gof.test$observed, + main="Poisson Distribution: Observed vs. Expected", + type='b', + pch=0, + col='blue', + xlab="Number of Events", + ylab="Frequency", + ylim=c(0,55)) > > # Dual Y-axes plot (secondary Y-axis) > par(new=T) > > # Expected frequency > plot(c(0:9), chisq.gof.test$expected, + type='b', + pch=1, + col='red', + xlab="", + ylab="", + ylim=c(0,55)) > > legend(x=6.5, y=50, + c("Observed", "Expected"), + pch=c(0,1), + col=c('blue', 'red'))

 

 

Poisson(λ=3.844) 로부터의 각 class별 기대도수(Expected frequency)와 250개 샘플로 부터의 관측도수(Observed frequency)가 매우 유사한 분포를 띠고 있음을 눈으로도 확인할 수 있습니다.

 

많은 도움이 되었기를 바랍니다.