정규성 검정 해석 - jeong-gyuseong geomjeong haeseog

2022년 08월 29일 버전입니다. 변수는 계속하여 추가될 예정입니다. 다음 카테고리에 있는 글에서 이용된 데이터입니다. - 기술 통계 - 통계 프로그램 사용 방법 1) 엑셀 파일 2) CSV 파일 3) 코드북

medistat.tistory.com

먼저 데이터를 불러온다

(데이터 불러오는 방법은 다음 링크 참고:

2022.08.04 - [통계 프로그램 사용 방법/SPSS] - [SPSS] 데이터 불러오기 및 저장하기)

정규성 검정 방법

1) 분석 (A)>기술통계량 (E) > 데이터 탐색 (E)

2) 분석하고자 하는 변수를 "종속변수(D)" 쪽으로 이동

3) 도표(T)를 누르고, "히스토그램(H)", "검정과 함께 정규성도표(O)" 체크박스를 클릭. 다 되었으면 "계속(C)" 클릭

4) "확인" 을 눌러 결과 확인

1) Q-Q plot

R에서 그린 그림과 비교해보면 약간 다른데, x축과 y축이 바뀌어 있기 때문이다.

R로 그린 Q-Q plot은 다음 링크에서 확인할 수 있다: 2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (1) : QQplot - qqnorm(), qqline() 

해석방법

데이터가 직선상에 있음: 정규성 따름

데이터가 직선에서 벗어나 있음: 정규성 따르지 않음

따라서, "대부분의 데이터들이 일직선 상에 있으므로 정규성을 따른다고 할 수 있다."

Q-Q plot 이론은 다음 링크에서 확인할 수 있다.

2022.08.12 - [통계 이론] - [이론] Q-Q Plot (Quantile-Quantile Plot)

[이론] Q-Q Plot (Quantile-Quantile Plot)

[이론] Q-Q Plot (Quantile-Quantile Plot) 정규성을 검정할 때 Q-Q Plot을 쓰곤 한다. 그런데 이런 궁금증이 들 수 있다. 왜 Q-Q Plot이 직선에 가까운 것이 정규성을 따른다는 뜻인가? 이에 대해 조목조목..

medistat.tistory.com

2) 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov test) 

해석방법

초급자: p-value>0.05이므로 정규성을 따른다.

중급자:

귀무 가설과 대립 가설은 다음과 같다.

-귀무 가설: ALT는 정규성을 따른다.

-대립 가설: ALT는 정규성을 따르지 않는다.

이때 제 1종 오류는 다음과 같다.

제 1종 오류

=귀무 가설이 참인데도 기각하고 대립 가설을 택함

=ALT는 정규성을 따르는데도 따르지 않는다고 결론 내림

p-value는 제 1종 오류를 범할 확률을 의미하므로 다음과 같이 결론내릴 수 있다.

결론: 정규성을 따르는 것이 진실인데, 따르지 않는다고 결론내렸을 확률이 0.05보다 크다.

=정규성을 따르는 것이 진실인데, 따르지 않는다고 결론내린 것은 잘못일 수 있다.

=정규성을 따르지 않는다고는 할 수 없다.

고급자:

-(비록 본 데이터는 n수가 충분하지만) n수가 굉장히 적다면 정규성 검정을 잘 통과하므로 결과를 맹신해서는 안 된다.

-본 데이터의 n수는 1,000명으로 2,000명 미만이므로 Shapiro-Wilk test의 결과를 인용한다. (SAS 기준, 2,000명까지는 Shapiro-Wilk test을 계산해주고, Kolmogorov-Smirnov test는 2,000명 이상이 필요하다.)

-또한, QQplot과 히스토그램의 결과도 고려하여 정규성 여부를 판단해야 한다.

R에서 시행한 Kolmogorov-Smirnov test 결과와 약간 다르다. R의 Kolmogorov-Smirnov test결과는 다음 링크에서 확인할 수 있다.2022.08.11 - [기술 통계/R] - [R] 정규성 검정 (3) : 정량적 검정 (Shapiro-Wilk, Kolmogorov-Smirnov) - shapiro.test(), ks.test()

다르게 나타나는 이유는 두 가지가 있다.

1) SPSS 결과의 p-value "0.200*" 이란 ">0.2"를 의미한다. 즉 적어도 0.2보다는 크다는 뜻을 의미하며, R에서의 결과인 "0.9717"과는 그렇게 배치되는 이야기가 아니다. 

2) SPSS는 Lilliefors가 수정한 방식의 유의확률을 계산해준다. 초창기에 만든 Kolmogorov-Smirnov test는 모종의 문제가 있어 Lilliefors가 p-value를 수정하게 되었다. SPSS는 이를 제시해주고 있다. R에서 Lilliefors가 수정한 Kolmogorov-Smirnov test를 하는 방법은 다음 링크에서 확인할 수 있다.

논문 자료 분석할 때 초보자들이 하기 쉬운 가장 흔한 실수중의 하나가 탐색적 자료분석 (Explanatory Data Analysis)에 들이는 시간을 아까워 한다는 점이다. 언뜻 보기에는 화면 전체에 복잡한 숫자와 p-value들이 난무하는 다변수 분석이 뭔가 하는 것 같고 멋있어 보이겠지만, 대부분의 결과는 그 이전에 이미 결정되어 있는 경우가 많다.

1. 정규성 검정이란?

정규분포 자료는 위의 그림과 같이 종모양이며 평균을 중심으로 양쪽으로 예쁘게 분포되어 있는 자료를 말한다. 정규분포가 중요한 이유는, 우리가 다루는 대부분의 통계적 분석법은 자료(특히 결과변수)가 정규분포를 만족한다고 가정하고 진행하는 것이기 때문이다.

2. 중심극한정리 (Central Limit Theorem)

확률론과 통계학에서, 중심 극한 정리(中心 極限 定理, 영어: central limit theorem, 약자 CLT)는 동일한 확률분포를 가진 독립 확률 변수 n개의 평균의 분포는 n이 적당히 크다면 정규분포에 가까워진다는 정리이다. -by Wiki

뭔 소린지 잘 모르는게 정상일 수도 있다 하지만 중요한 것은

라는 것이다. 중심극한정리를 아무 상황에서나 막 들이대면 안된다.

그럼 사전 설명은 이 정도로 마치고 실제 자료를 이용해서 정규성 검정을 해 보자.

3. Hadley Wicham의 다이아몬드 자료

Hadley Wicham은 R을 사용하는 사람이라면 누구나 한 번쯤 들어 봤을 인물이다. 지금 이용할 자료는 Hadley Wicham이 ggplot2 에 넣은 다이아몬드 데이터(50000개 짜리)를 2캐럿 미만의 500개 자료로 축약한 버전이다.

이번에 다룰 변수는 carat과 price이다. 위의 기술통계량 결과를 보면 여러 결과치들이 보이는 데 거기서 놓치지 말아야 할 부분이 최소값과 최대값 그리고 결측치의 수이다. 최소값과 최대값은 자료의 성질을 알 수 있게 해주며, 자료 분석의 과정을 결정하기도 한다 (당장 아래 보면 자료에 0이 포함되었는지 아닌지에 따라 뭔가 많이 달라진다). 거기에 이상치( 특히 입력 오차)를 알 수 있게 해주기도 한다. 예를 들어 키의 최대값이 1705cm라면 170.5cm의 입략 오류가 아닌가 확인할 필요가 있다 (아주 흔한 일이다. 자료를 입력해 주는 인턴이나 아랫년차가 그렇게 사명감을 가졌다고 생각하지 않는게 좋다. 더욱이 사람은 누구나 실수를 한다.).

상식적으로 다이아몬드 가격을 결정하는 가장 큰 요소는 다이아의 크기이다. 이를 확인하기 위해 산점도를 확인해 보면,

위 그림을 보면 다이아의 크기가 클 수록 가격이 비싸지는 듯 한데… 왠지 선형회귀 적합선 (빨간 점선)이 경과와 그리 잘 들어 맞지는 않는 것 같다. 0.2캐럿 미만의 다이아는 가격이 너무 낮게, 0.5~1캐럿 사이는 좀 더 비싸게 예측한다. 더구나 1캐럿을 넘어가면 점들의 분포가 기하급수적으로 커진다. 이렇게 뭔가 이상한 부분들이 보이는 것은 두 자료의 관계를 보기 전에 했어야 하는 과정들 (탐색적 자료분석) 과정이 생략되었기 때문에 발생한 것이다. 자료를 하나 하나 보자.

분석 > 기술통계량 > 데이터 탐색 ('도표' 항목에서 '검정과 함께 정규성 도표' 항목 체크)

SPSS에서 전통적으로 제공하는 정규성 검정은 Kolmogrov-Smirov 검정이다. 그 뒤의 Shapiro-Wilk 검정은 자료의 갯수가 5000미만일 때만 분석이 가능한데, 그 외에는 K-S 검정에 비해서 우수하다.(근데 자료의 수가 늘면 정규성 검사를 통과하기가 점점 힘들어 지기 때문에 정규성 검정 결과 자체를 신뢰할 수 없게 된다.) K-S 검정의 결과는 0.001 미만이다. 따라서 carat 변수는 정규분포를 만족하지 못한다. 왜도 (좌우 비뚤림) 0.690, 첨도 (상하 비뚤림) -0.277 으로 측정되었고, Q-Q plot에서도 이탈한 관측치들이 보인다.

p-value가 0.000으로 표현되는 경우는 p-value가 진짜 0 이어서가 아니라, 정확히 말하면 0.0005 미만인 경우이다. SPSS가 반올림 후 소수점 3자리 까지만 보여주기 때문에 이렇게 보이는 건데, 논문 표에 진짜 0.000 이라고 적으면 안된다. <0.001로 적는게 제대로 된 표현법이다.

히스토그램 상으로도 전형적인 정규분포모양 (종모양이며 양쪽으로 예쁘게 퍼진)이 아니라 우측으로 치우친 상태 (Right skewed, Positively skewed, 우리말 어감과는 다르다. 왼쪽으로 치우친 것 같은데…)이다.

4. 비뚤림의 해소

1. Tukey Ladder of Power (\(x'=x^\lambda\))

예로부터 내려오는 말중에 “뭔가 시원치 않으면 로그 씌워라.”라는 말이 있다. 지금 생각해보면 참 명언이라는 생각이 든다. 의학자료분석에서 상당수의 (아니 거의 대부분의) 자료들은 우로 치우친 형태이며 로그 씌우면 해결된다. 다만 여기서 너무 심했다 싶으면 루트 씌워보고, 너무 약하다 싶으면 \(-1/\sqrt{x}\)이나 \(-1/x\) 해보면 된다. 변환 시 가장 중요한 부분 중 하나가 변환하려면 자료가 0을 초과하는 양수여야 한다는 점이다.(\(ln(0)\)은 과연 얼마일까요?)

1. 장점: 해보기 쉽고 만만하다는 점. 그리고 나중에 계산하기 쉽다는 점이다. (이 부분은 나중에 설명)

2. 단점: 계산을 여러번 해야 하며, 노가다 같아서 뭔가 그럴 듯한 측면이 없다는 점이다.

2. Box-Cox transformation

if \(\lambda \not= 0, x>0\) \[x'=\frac{(x^\lambda-1)}{\lambda}\] if \(\lambda = 0, x>0\) \[x'=log(x)\]

Tukey법은 \(x^\lambda\)로만 변환하는데 비해서, Box-Cox법은 \(\frac{(x^\lambda-1)}{\lambda}\)를 해 줌으로써 계산은 복잡해 지지만, 최대치와 최소치는 변환전과 큰 차이가 없도록 만든다. 거기다 최적의 \(\lambda\) 수치를 계산할 수 있는 방법을 제공한다. (Box-Cox는 \(\lambda\)가 0일때는 log변환, 물론 변환 변수(x)는 0을 초과하는 양수여야 함. 나머지 값에 대해서는 공식에 입력해서 해결)

carat 변수에 대해 \(\lambda\)를 0, 0.25, -0.25 라고 가정한 후 Box-Cox 변환을 시행하였다.

변환 > 변수계산

이 세변수를 기술통계항목을 이용해서 정규성 검정을 하게 되면

로그변환을 기준으로 0.25와 -0.25 값을 시도해 본 결과 다 비슷하긴 한데 첨도와 왜도 기준으로 로그 변환이 제일 좋은 결과를 보였다 (첨도와 왜도 모두 0이 이상적인 정규분포 값이며 \(\pm\) 2 이내가 납득 가능한 수준). K-S 검정 결과는 세 가지 변환 모두 정규분포를 만족하지 못하는 것으로 나왔지만 변환 후의 변수가 정규성을 만족하는 경우는 그리 많지 않다. 따라서 목표를 한쪽으로 쏠리지 않고 평평하게 정도로 잡는 것이 현실적이다.

히스토그램 상으로도 로그 변환이 제일 한 쪽으로 쏠리지 않고 골고루 펼쳐져 있어 보인다. 별 차이 없으면 로그 변환이 제일 대중적 이기도 하고… 필요하면 다른 값을 넣어서 비교해 봐도 된다. 만약 \(\lambda\) -0.25가 제일 좋으면 -0.5라던가 -1 등 추가적인 수치에 대해 시도해 봐야 한다. (R에서는 최적의 \(\lambda\) 값을 계산할 수 있다. SPSS에서는 그게 안 되니 이런 식으로 수동으로 하는 방법이 최선이다. 아니면 우로 치우친 자료가 대부분이니 복잡하게 생각 할 필요 없이 로그 씌우는게 정답일 수도 있다.) 그리고 의학논문에서 명확하게 입력오류로 밝혀지지 않은 이상치를 찾고 제거하는 일은 극히 드문일이지만, 이상치를 찾아 보려면 정규변환 후에 해야만 한다.

3. Yeo-Johnson transformation

if \(\lambda \not= 0, x\ge0\) \[x'=\frac{((x+1)^\lambda-1)}{\lambda}\] if \(\lambda = 0, x\ge0\) \[x'=log(x+1)\]

흔히 Johnson transformation으로 불리는 변환법이다. 근데 앞의 Yeo가 숙명여대에 계시는 여인권교수님인거 알고 나서는 왠만하면 정식으로 이름을 다 부른다. 이 변환법의 의미는 변환할 변수 x가 0을 포함할 때의 해결법을 제시해 준다는 것이다.

4. price 변수 변환

price의 경우도 로그 변환하는 것이 적절하다.

5. 로그변환 전과 후의 차이

뭐가 달라졌는지는 길게 언급하지 않아도 될 것 같다.

회귀방정식의 R-squared 수치도 0.820에서 0.927로 증가하였다. 로그 변환 후의 회귀방정식 결과는 \[log(price)=1.672\times log(carat)+8.455\] \[ price= e^{1.672 \times log(carat)+8.455}\] 으로 표현할 수 있다. (다변수 분석 결과는 더 복잡해지며 Box-Cox 사용하면 더욱 혼잡해 진다. 먼저 언급했던 계산하기 쉽고 어려운게 이 의미이다.

사실 다이아몬드 예제는 로그 변환 후 모형결과가 굉장히 많이 좋아지는 예제이다. 그럼에도 불구하고 사용된 두 가지 변수 모두 변환후에도 정규성을 만족하지는 못했다. 우선 carat의 경우를 보면 봉우리가 3개 있는 것이 문제의 하나로 보인다. 그 지점을 보면 0.5, 1, 1.5 캐럿의 지점이다. 다이아 반지를 살 때 대개 5부, 1캐럿 반지를 많이 하지 0.6캐럿, 1.2캐럿 반지는 잘 못 들어 본 것 같다. 다이아를 자를 때 0.5, 1캐럿에 맞게 자르는 경우가 많을 것 같다. 예를 들어 0.6캐럿 다이아 라면 좀 흠집 있어 보이는 부분을 잘라내서 높은 수준의 0.5캐럿 반지를 만든다던가… 가격도 보면 5000불 에서 평행한 부분이 보이는데 로그변환 그래프를 보면 그 부분이 좀 더 두드러져 보인다 (\(e^{8.5}=4914\)). 가격의 심리적 저항선을 의미하는 것 같다.

결국 이렇게 봉우리가 생기는 것은 서로 다른 여러 집단이 혼재되어 있다는 의미이며, 정확한 결과를 얻기 위해서는 이를 분리해서 분석하는 것이 합리적이다. 나 같으면 0.5캐럿과 1캐럿 다이아만 따로 분리해서 분석한 후 이해가 완벽해지면 자료를 통합해서 분석하는 방법을 택할 것 같다. 의학논문의 특성 상 대상자가 사람이라서 대상자 한 사람이라도 아까워서 일부 대상자에 대해서만 분석한다는 생각을 하기 힘든데, 그게 전체적으로 봐서는 오히려 손해일 수도 있다. 

어쩔 수 없는 경우 히스토그램 상 종 모양 까지는 아니더라도 어느 정도 even하게 펼쳐져 있고, skew와 kurtosis 지표가 \(\pm 2\) 이내인 경우, 이 변수를 변환한 채로 사용할 지 말지의 결정을 다변수 분석 후 까지 미룬다. 그리고 원칙적으로 정규성을 만족해야 하는 것은 결과변수뿐이다. 원인변수는 그냥 even 하게 펼쳐져 있기만 하면 분석에 문제가 없다.

5. 노인 우울증 자료

65세 이상 단독가구 노인 361명의 우울지수(PHQ-9)에 대한 가상자료이다. ainc(가구 월수입)과 PHQ9을 중심으로 보자.

월수입은 우로 치우쳐 있으며, 노인 연금 때문인지 최소값은 17만원이라서 (0을 포함하지 않아서) 로그변환에 문제는 없어 보인다. 한 달에 1500만원 번다고 응답한 사람도 있어서 이를 이상값으로 고려해야 할까 고민할 수도 있겠지만, 먼저도 언급했듯이 정규변환 이전까지는 이상치에 대한 고민을 하지 말라고 언급했었다. 현 상황에서는 옛 어른들의 말씀대로 무작정 로그변환하고 상황을 봐도 되지만 적절한 \(\lambda\)값을 찾기 위해 -0.5,-0.25,0,0.25,0.5값을 입력해서 결과를 확인했다.

최적의 \(\lambda\)가 0도 아니고 -0.25에 가장 가까워 보인다. 이럴때는 그냥 0에 더 가깝기도 하고 옛 어른들의 발씀을 따를 겸 로그 씌우는 방법이 나을 수도 있겠지만, 월소득 변수는 원인변수이니 결과변수 일 때처럼 계산에 강박적이지 않아도 된다 (논문 중간에 정규성 확보 위해서 -0.25승 했다면 소득에 따른 우울감이 논문의 주제가 아니라면 적당히 넘어가는 게 가능하다). 하지만, 그 결정은 연구자가 해야한다. 내 논문에 소득이 중요하고 그 결과가 수식으로 필요하다면 로그변환이 나을 것이지만, 단순히 covariate로 보정하기 위해서라면 -0.25 Box-Cox 결과가 더 좋을 것이다.

마지막은 PHQ-9에 대해 볼 것이다.

PHQ-9 역시 우로 치우쳐 있으며, 0을 포함하기 때문에 만약 변환이 필요한 경우 Yeo-Johnson 변환을 이용해야 한다. 그 첫 단계는 변수가 0을 초과할 수 있도록 1을 더해주는 것이다. ainc 변수에서와 동일한 노가다법으로 얻은 최적의 \(\lambda\) 값은 -0.25이며, 이를 이용해서 Yeo-Johnson 변환을 시행하면,

PHQ-9의 경우에는 적절한 변환을 했음에도 불구하고 제대로 된 정규분포를 만족하지 않는다. 그 이유는 0점이 너무 많아서 이며, 만약 PHQ-9을 결과변수로 사용하려 했다면 이런 식으로 변환된 결과를 사용해서는 안된다.

PHQ-9을 결과변수로 논문을 계획했다면,

1. 적절한 지점을 기준으로 2군으로 잘라서 로지스틱 회귀분석을 돌리거나,

2. PHQ-9 자체가 자연수인 정수 값을 가지니 poission, zero-inflated poission, negative-binomial, zero-inflated negative-binomial 4가지 분석을 시행해서 최적 결과를 도출 (SPSS에서는 포아송 분석과 음이항 분석은 제공하지만 0과잉 모형에 대해서는 분석을 제공하지 않는다.)

둘 중 하나의 방법으로 진행해야 할 것이다.

6. Take Home Message

지금까지 여러 변수들의 정규성을 확인하고, 정규성을 만족하게 하기 위한 수정과정에 대해 알아 보았다. 로그 변환으로 쉽게 정규성을 회복하는 경우도 있었지만 그렇지 않은 경우도 있었다. 납득할 만한 결과를 얻기 위해 꾸준히 고민하는 것이 좋은 결과를 내는 가장 좋은 방법일 것이다.

SPSS는 약간의 지식만 가지면 통계에 대해 지식이 많지않은 비전공자라도 쉽게 결과를 얻을 수 있다는 장점이 있지만, 아무래도 부족한 부분이 존재한다. 이 게시물과 동일한 내용을 가지고 R을 이용해 분석한 게시물 이 있으니 한 번 살펴보는 것도 나쁘지 않을 것 같다.