탐색적 자료분석 1. 정규성 검정 (R)

normal.Rdata

0.01MB

 

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

 

1. 정규성 검정이란?

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

 

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

 

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

 

 

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

 

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

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

 

 

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

Hadley Wicham은 R을 사용하는 사람이라면 누구나 한 번쯤 들어봤어야 하는 인물이다. ggplot2, dplyr … 의 기라성 같은 패키지의 제작자이며 대인배 기질을 지닌 프로그램 교육 사이트인 데이터 캠프의 주요 멤버이다. 지금 이용할 자료는 ggplot2 에 포함된 자료(50000개 짜리)를 2캐럿 미만의 500개 자료로 축약한 버전이다.

load("D:/Data/normal.Rdata")
library(dplyr)
library(ggplot2)
library(car)
library(nortest)
library(psych)

자료와 필요한 패키지를 불러온 후,

head(diamond)

##   carat       cut color clarity price
## 1  0.46     Ideal     I     VS1   997
## 2  0.90 Very Good     H     VS1  4021
## 3  0.60     Ideal     G     VS1  2109
## 4  1.04   Premium     I     SI1  4290
## 5  0.50     Ideal     J     VS2  1188
## 6  1.54   Premium     G    VVS2 15983

str(diamond)

## 'data.frame':    500 obs. of  5 variables:
##  $ carat  : num  0.46 0.9 0.6 1.04 0.5 1.54 1.01 0.33 1.22 0.9 ...
##  $ cut    : Factor w/ 5 levels "Fair","Good",..: 5 3 5 4 5 4 3 4 5 3 ...
##  $ color  : Factor w/ 7 levels "D","E","F","G",..: 6 5 4 6 7 4 3 5 7 1 ...
##  $ clarity: Factor w/ 8 levels "I1","SI2","SI1",..: 5 5 5 3 4 6 2 8 3 3 ...
##  $ price  : int  997 4021 2109 4290 1188 15983 3709 868 5586 4523 ...

이번에 다룰 변수는 carat과 price이다.

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

diamond %>% ggplot(aes(carat, price))+geom_point()+geom_smooth(method = "lm",se=F, col="red", lty=2)+theme_classic()

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

 

예전에 R공부할 때 너무 힘들었던 것 중 하나가 책이나 교재마다 쓰는 명령어가 달랐다는 점이다. 어찌 보면 하나의 결과를 얻기 위한 여러 방법이 있다는게 R의 장점이긴 한데, 익숙해지기까지 너무 힘들었다. 이럴 때 가장 좋은 방법은 모르는 부분이 있어도 아무 고민하지 말고 직접 입력해 가면서 익숙해지는 것 같다. 방법이 많다고 해봐야 몇 개 안되고, 익숙해지면 그게 그거다. 더욱이 혹시 새 방법이 나한테 더 잘 맞을지도…

 

shapiro.test(diamond$carat)

## 
##  Shapiro-Wilk normality test
## 
## data:  diamond$carat
## W = 0.91947, p-value = 1.089e-15

describe(diamond$carat)

##    vars   n mean   sd median trimmed  mad  min max range skew kurtosis
## X1    1 500 0.76 0.41    0.7    0.72 0.47 0.23   2  1.77 0.69     -0.3
##      se
## X1 0.02

diamond %>% ggplot(aes(carat))+geom_density(col="red", fill="red", alpha=0.05)+theme_classic()

qqPlot(diamond$carat) # car package

## [1] 185 266

R에서 기본으로 지원하는 정규성 검정은 Shapiro-Wilk 검정이다. Sahpiro test는 자료의 갯수가 5000미만일 때만 분석이 가능한데, 그 외에는 다른 검정에 비해서 모든 부분에서 우수하다 (근데 자료의 수가 많아지면 정규성 검정 자체가 의미 없어진다). 자료의 수가 많아지면 사용 가능한게 Anderson-Darling test이고, Lilliefors(Kolmogorov-Smirov) test는 SPSS 사용자 들에게 익숙한 검사이다.

Shapiro 검정의 결과는 p-value 1.089e-15로 0.001 미만이다. 따라서 carat 변수는 정규분포를 만족하지 못한다. 아래 psych package의 describe 명령의 결과물을 확인 해 봐도 skew (좌우 비뚤림) 0.69, kurtosis (상하 비뚤림) -0.3으로 측정되었고, 마지막의 density plot을 봐도 우측(Right skewed, Positively skewed, 우리말 어감과는 반대이다. 왼쪽으로 치우친 것 같은데…) 으로 치우쳐 있다. 이렇게 한 쪽으로 치우쳐 있는 것이 정규분포를 하지 않는다는 가장 큰 증거다. Q-Q plot에서도 이탈한 관측치들이 보인다.

 

4. 추가적인 정규성 검사

nortest package를 사용하면 원하는 정규성 검사를 원 없이 할 수 있다. 개인적으로는 자료의 분포가 중요하지 정규성 검사 결과 자체는 큰 의미 없다고 믿는 쪽이지만 그래도 제목이 정규성 검정이니 언급은 하는게 맞는 것 같다.

# R base function
shapiro.test(diamond$carat)

## 
##  Shapiro-Wilk normality test
## 
## data:  diamond$carat
## W = 0.91947, p-value = 1.089e-15

# nortest function
ad.test(diamond$carat)

## 
##  Anderson-Darling normality test
## 
## data:  diamond$carat
## A = 12.139, p-value < 2.2e-16

lillie.test(diamond$carat)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  diamond$carat
## D = 0.13856, p-value < 2.2e-16

 

5. 비뚤림의 해소

1. Tukey Ladder of Power (\(x'=x^\lambda\), 직접 몸빵법)

\(\lambda\)	변수상태	변환법
3	좌로 치우침	\(x^3\)
2	\(\uparrow\)	\(x^2\)
1	정상 상태	변환 필요없음
1/2	\(\downarrow\)	\(\sqrt{x}\)
1/3	\(\downdownarrows\)	\(x^{1/3}\)
0	우로 치우침	\(ln(x)\)
-1/2	\(\downarrow\)	\(-1/\sqrt{x}\)
-1	\(\downarrow\)	\(-1/x\)
-2	\(\downdownarrows\)	\(-1/x^2\)

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

 

• 장점은 해보기 쉽고 만만하다는 점. 그리고 나중에 계산하기 쉽다는 점이다. (이 부분은 나중에 설명)
• 단점은 계산을 여러번 해야 하며, 노가다 같아서 뭔가 그럴 듯한 측면이 없다는 점이다.

 

2. Box-Cox transformation

 

\(x'\)=\(\frac{(x^\lambda-1)}{\lambda}\) if \(\lambda \not= 0\)

\(x'\)=\(log(x)\) if \(\lambda = 0\)

 

Tukey법은 \(x^\lambda\)로만 변환하는데 비해서, Box-Cox법은 \(\frac{(x^\lambda-1)}{\lambda}\)를 해 줌으로써 계산은 복잡해 지지만, 최대치와 최소치는 변환전과 큰 차이가 없도록 만든다. 거기다 최적의 \(\lambda\) 수치를 계산할 수 있는 방법을 제공한다.

boxCox(diamond$carat~1)

정확한 \(\lambda\)의 최대치는 0 근처이다.(정확히는 95% 신뢰구간 내에 0이 포함된다.

만약 \(\lambda\)의 최대치가 0.15가 나왔다 해보자. \(2^{0.15}\) 계산 가능한가?

) 그럼 그 결과를 믿고 log 씌우면 된다. (Box-Cox는 \(\lambda\)가 0일때는 log변환, 물론 변환 변수는 0을 초과하는 양수여야 함. 나머지 값에 대해서는 공식에 입력해서 해결)

 

3. Yeo-Johnson transformation

 

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

\(x'=log(x+1)\) if \(\lambda = 0, x\ge0\)

 

흔히 Johnson transformation으로 불리는 변환법이다. 근데 앞의 Yeo가 숙명여대에 계시는 여인권교수님인거 알고 나서는 왠만하면 정식으로 이름을 다 부른다. 이 변환법의 의미는 변환변수 x가 0을 포함할 때 의 해결법을 제시해 준다는 것이다. 전체 공식은 Box-Cox로 해결이 되지 않는 x가 음의 값인 경우 및 0을 포함하는 경우일 때의 해결법을 제시해 주니 궁금하면 한 번 검색해 보는 것도 좋겠다.

 

4. 실제변환

diamond$log_carat <- log(diamond$carat)
shapiro.test(diamond$log_carat)

## 
##  Shapiro-Wilk normality test
## 
## data:  diamond$log_carat
## W = 0.94598, p-value = 1.574e-12

describe(diamond$log_carat)

##    vars   n  mean   sd median trimmed  mad   min  max range  skew kurtosis
## X1    1 500 -0.42 0.55  -0.36   -0.42 0.66 -1.47 0.69  2.16 -0.04    -1.24
##      se
## X1 0.02

diamond %>% ggplot(aes(log_carat))+geom_density(col="red", fill="red", alpha=0.05)+theme_classic()

qqPlot(diamond$log_carat)

## [1] 185 266

Shapiro 검정결과는 여전히 만족스럽지 못하지만(p-value<0.001), skweness와 kurtosis 수치가 좋아졌고, density plot결과와 Q-Q plot 결과도 나아졌다. p-value로 표현되는 정규성 검정 결과만을 맏지 말고, 여기서 만족하고 중단해야 한다. 그리고 의학논문에서 명확하게 입력오류로 밝혀지지 않은 이상치를 찾고 제거하는 일은 극히 드문일이지만, 이상치를 찾아 보려면 정규변환 후에 해야만 한다.

 

5. price 변수 변환

shapiro.test(diamond$price)

## 
##  Shapiro-Wilk normality test
## 
## data:  diamond$price
## W = 0.82824, p-value < 2.2e-16

describe(diamond$price)

##    vars   n    mean      sd median trimmed     mad min   max range skew
## X1    1 500 3597.16 3351.67 2409.5 3005.58 2482.61 374 17761 17387 1.55
##    kurtosis     se
## X1     2.35 149.89

diamond %>% ggplot(aes(price))+geom_density(col="blue", fill="blue", alpha=0.05)+theme_classic()

qqPlot(diamond$price)

## [1] 492 456

boxCox(diamond$carat~1)

price의 경우도 똑같이 로그 변환해주면 될 것 같다.

diamond$log_price <- log(diamond$price)
shapiro.test(diamond$log_price)

## 
##  Shapiro-Wilk normality test
## 
## data:  diamond$log_price
## W = 0.96363, p-value = 8.746e-10

describe(diamond$log_price)

##    vars   n mean   sd median trimmed mad  min  max range skew kurtosis
## X1    1 500 7.76 0.96   7.79    7.75 1.2 5.92 9.78  3.86 0.02    -1.12
##      se
## X1 0.04

diamond %>% ggplot(aes(log_price))+geom_density(col="blue", fill="blue", alpha=0.05)+theme_classic()

qqPlot(diamond$log_price)

## [1] 492 456

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

diamond %>% ggplot(aes(carat, price))+geom_point()+geom_smooth(method = "lm",se=F, col="red", lty=2)+theme_classic()

diamond %>% ggplot(aes(log_carat, log_price))+geom_point()+geom_smooth(method = "lm",se=F, col="red", lty=2)+theme_classic()

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

summary(lm(price~carat, data=diamond))

## 
## Call:
## lm(formula = price ~ carat, data = diamond)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6311.3  -754.9   -84.0   469.2  9274.4 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2109.3      135.7  -15.54   <2e-16 ***
## carat         7470.8      156.9   47.62   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1424 on 498 degrees of freedom
## Multiple R-squared:  0.8199, Adjusted R-squared:  0.8196 
## F-statistic:  2268 on 1 and 498 DF,  p-value: < 2.2e-16

summary(lm(log_price~log_carat, data=diamond))

## 
## Call:
## lm(formula = log_price ~ log_carat, data = diamond)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.98687 -0.16119 -0.01754  0.16639  1.08327 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.45525    0.01452  582.32   <2e-16 ***
## log_carat    1.67177    0.02104   79.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2592 on 498 degrees of freedom
## Multiple R-squared:  0.9269, Adjusted R-squared:  0.9268 
## F-statistic:  6316 on 1 and 498 DF,  p-value: < 2.2e-16

 

회귀방정식의 R-squared 수치도 0.8199에서 0.9269로 증가하였다. 로그 변환 후의 회귀방정식 결과는 \[log(price)=1.672\times log(carat)+8.455\]

price= exp(1.672*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캐럿 다이아만 따로 분리해서 분석한 후 이해가 완벽해지면 자료를 통합해서 분석하는 방법을 택할 것 같다. 의학논문의 특성 상 대상자가 사람이라서 대상자 한 사람이라도 아까워서 일부 대상자에 대해서만 분석한다는 생각을 하기 힘든데, 그게 전체적으로 봐서는 오히려 손해일 수도 있다.

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

 

6. 노인 우울증 자료

head(dep)

##   Women age edu      ainc health sick work PHQ9
## 1     1  75   2  57.91667      5    1    2   19
## 2     1  68   1  71.00000      1    2    1    0
## 3     0  65   1 100.00000      3    2    2    3
## 4     1  70   1  40.00000      5    1    2   18
## 5     1  68   3 175.00000      3    2    2   11
## 6     1  78   1  34.16667      3    1    1    0

str(dep)

## 'data.frame':    318 obs. of  8 variables:
##  $ Women : num  1 1 0 1 1 1 0 1 0 1 ...
##  $ age   : num  75 68 65 70 68 78 67 73 65 70 ...
##  $ edu   : Factor w/ 4 levels "1","2","3","4": 2 1 1 1 3 1 1 1 1 3 ...
##  $ ainc  : num  57.9 71 100 40 175 ...
##  $ health: Factor w/ 5 levels "1","2","3","4",..: 5 1 3 5 3 3 2 1 4 3 ...
##  $ sick  : Factor w/ 2 levels "1","2": 1 2 2 1 2 1 2 2 2 2 ...
##  $ work  : Factor w/ 2 levels "1","2": 2 1 2 2 2 1 2 2 2 1 ...
##  $ PHQ9  : num  19 0 3 18 11 0 0 0 2 1 ...
##  - attr(*, "na.action")= 'omit' Named int  8 14 28 32 33 63 96 97 106 110 ...
##   ..- attr(*, "names")= chr  "8" "14" "28" "32" ...

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

shapiro.test(dep$ainc)

## 
##  Shapiro-Wilk normality test
## 
## data:  dep$ainc
## W = 0.39498, p-value < 2.2e-16

describe(dep$ainc)

##    vars   n  mean    sd median trimmed   mad min  max range skew kurtosis
## X1    1 318 75.03 99.96  54.58   60.39 34.96  17 1500  1483 9.76    129.3
##      se
## X1 5.61

dep %>% ggplot(aes(ainc))+geom_density(col="red", fill="red", alpha=0.05)+theme_classic()

qqPlot(dep$ainc)

## [1] 21 46

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

boxCox(dep$ainc~1)

최적의 \(\lambda\)가 0도 아니고 -0.5도 아니게 어중간하게 되었다. 이럴때는 그냥 0에 더 가깝기도 하고 옛 어른들의 발씀을 따를 겸 로그 씌우는 방법이 나을 수도 있겠지만, 월소득 변수는 원인변수이니 결과변수 일 때처럼 계산에 강박적이지 않아도 된다 (논문 중간에 정규성 확보 위해서 -0.25승 했다면 소득에 따른 우울감이 논문의 주제가 아니라면 적당히 넘어가는 게 가능하다). 하여간 그림을 확대해 보자.

boxCox(dep$ainc~1, lambda = seq(-0.3,-0.1, 0.01))

로그변환과 -0.25승 Box-Cox 변환을 둘 다 해서 비교해 보면

dep$log_ainc <- log(dep$ainc)
dep$ainc2 <- ((dep$ainc^-0.25)-1)/(-0.25)
shapiro.test(dep$log_ainc)

## 
##  Shapiro-Wilk normality test
## 
## data:  dep$log_ainc
## W = 0.97047, p-value = 4.267e-06

shapiro.test(dep$ainc2)

## 
##  Shapiro-Wilk normality test
## 
## data:  dep$ainc2
## W = 0.98271, p-value = 0.000709

describe(dep$log_ainc)

##    vars   n mean  sd median trimmed  mad  min  max range skew kurtosis
## X1    1 318 4.02 0.7      4       4 0.66 2.83 7.31  4.48 0.56     1.06
##      se
## X1 0.04

describe(dep$ainc2)

##    vars   n mean   sd median trimmed  mad  min  max range  skew kurtosis
## X1    1 318 2.52 0.25   2.53    2.52 0.23 2.03 3.36  1.33 -0.02    -0.26
##      se
## X1 0.01

dep %>% ggplot(aes(log_ainc))+geom_density(col="red", fill="red", alpha=0.05)+theme_classic()

dep %>% ggplot(aes(ainc2))+geom_density(col="blue", fill="blue", alpha=0.05)+theme_classic()

qqPlot(dep$log_ainc, col.lines = "red")

## [1] 21 46

qqPlot(dep$ainc2, col.lines = "blue")

## [1] 21 46

Box-Cox 변환 결과가 좀 더 나아보이기는 하지만, 그 결정은 연구자가 해야한다. 내 논문에 소득이 중요하고 그 결과가 수식으로 필요하다면 로그변환이 나을 것이지만, 단순히 covariate로 보정하기 위해서라면 Box-Cox 결과가 더 좋을 것이다.

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

shapiro.test(dep$PHQ9)

## 
##  Shapiro-Wilk normality test
## 
## data:  dep$PHQ9
## W = 0.75248, p-value < 2.2e-16

describe(dep$PHQ9)

##    vars   n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 318 3.78 5.11      2    2.73 2.97   0  23    23 1.65     2.01 0.29

dep %>% ggplot(aes(PHQ9))+geom_density(col="red", fill="red", alpha=0.05)+theme_classic()

qqPlot(dep$PHQ9)

## [1] 137  41

PHQ-9 역시 우로 치우쳐 있으며, 0을 포함하기 때문에 만약 변환이 필요한 경우 Yeo-Johnson 변환을 이용해야 한다. 그 첫 단계는 변수가 0을 초과할 수 있도록 1을 더해주는 것이다.

boxCox((dep$PHQ9+1)~1)

boxCox((dep$PHQ9+1)~1, lambda = seq(-0.4,-0.1,0.01))

dep$PHQ2 <- (((dep$PHQ+1)^-0.25)-1)/(-0.25)
shapiro.test(dep$PHQ2)

## 
##  Shapiro-Wilk normality test
## 
## data:  dep$PHQ2
## W = 0.87244, p-value = 1.45e-15

describe(dep$PHQ2)

##    vars   n mean   sd median trimmed  mad min  max range skew kurtosis
## X1    1 318 0.84 0.74   0.96     0.8 1.17   0 2.19  2.19 0.15    -1.39
##      se
## X1 0.04

dep %>% ggplot(aes(PHQ2))+geom_density(col="blue", fill="blue", alpha=0.05)+theme_classic()

qqPlot(dep$PHQ2, col.lines = "blue")

## [1] 137  41

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

 

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

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

2. PHQ-9 자체가 자연수인 정수 값을 가지니 poission, zero-inflated poission, negative-binomial, zero-inflated negative-binomial 4가지 분석을 시행해서 최적 결과를 도출

 

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

 

7. Take Home Message

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