대다수의 통계적 방법들은 특정한 가정을 전제한다. 그리고 그 가정이 성립할 때 잘 동작한다. 대표적 예시로는 관측한 data의 상관관계에 대한 가정 (독립적이거나, 특정한 covariance 형태를 가지거나) 혹은 data가 특정 형태의 분포를 가져야 한다는 것이다. 각각의 통계적 방법들은 적절한 가정 아래에서 이론적/실질적 당위성을 가진다. 따라서 우리는 분석하고자 하는 data를 자세히 살펴보면서, 사용하고 싶은 방법의 가정이 그 data에 적절한지 고민해야 한다.
수많은 "통계적 가정" 중에서, 실제 마주치는 data와 가장 거리가 먼 것은 무엇일까? 나는 data가 특정 분포를 가진다는 가정, 그 중에서도 Normal(Gaussian) distribution 에 대한 가정인 것 같다. 안타깝지만 대부분의 data는 Normal 하지 않거나, 전혀 Normal 하지 않다. 혹은 Normal한 data는 분석자에게 흥미롭지 않은 data일 수 있다.
그럼에도 불구하고 정규분포가 중요한 이유는 무엇인가? 꽤 어려워 보이는 확률분포함수에도 불구하고, 정규분포는 다루기 쉽다. 분포의 평균과 분산 두 가지 정보만 있다면, 우리는 정규분포의 모든 특성을 기술 할 수 있다. 또한, 우리가 자주 사용하는 통계량, 예를들어 평균, 의 분포는 정규분포이거나, asymptotic 하게 정규분포를 따르는 경우가 많다. 이러한 간단함과 범용성 때문에 우리 data가 정규분포를 따르는지 확인하는 작업이 필요하다.
일반적으로 data가 정규분포를 따르지 않는다고 해서 분석 결과가 무의미하진 않다. Linear regression은 data error가 평균이 0인 정규분포를 따른다는 가정을 한다. 하지만 (심지어 Linear regression 에서도) 그럼에도, 분포 shape를 조절하거나 data값을 적당히 조절해 분석 결과를 개선할 수 있다면, Non-Normal data를 Normal하게 변환하는 것이 의미있을 것이다.
이번 Post에서는, 주어진 data를 정규분포를 따르는 data로 변환하는 방법에 대해 정리해보자. Box-Cox (1964), Yeo-Johnson (2000) 모두 parametrized transform function 으로 주어진 data로 변환하고, likelihood 관점에서 변환된 data가 정규분포와 유사하게 하는 parameter 값을 추정하는 방법이다.
1. Transform for normality
1.1 Box-Cox Transform
Box-Cox 변환은 0보다 큰 실수를 가지는 original data \( X \)를 아래 수식을 사용해 \( y(\lambda) \)로 변환시킨다. Box-Cox 변환에서 \( \lambda \)값에 따라 \( y(\lambda) \)가 어떻게 결정되는지 알아보자.
\[ y(\lambda) = \begin{cases} \begin{split} \cfrac{ x^{\lambda}-1 }{ \lambda } &, & if \lambda \ne 0 \\ \ln x &, & if \lambda = 0 \end{split} \end{cases} \]
## Box Cox
bc = function(x, lambda){
if(lambda == 0){
return(log(x))
}
else{
((x^lambda) - 1) / lambda
}
}
## Original value to B-C Transformation
x_grid = seq(0.1, 2, by=0.01)
lambda_grid = seq(-1,2, by=0.2)
result = data.frame(x=NA,lambda=NA,bc=NA)
for(l in lambda_grid){
bc = data.frame(x=x_grid, lambda=l, bc=bc(x_grid, l))
result = rbind(result, bc)
}
library(ggplot2)
ggplot(result, aes(x=x, y=bc, color=lambda)) +
geom_point(size=0.5) +
scale_color_gradient(colours = rainbow(5))
1.2 Yeo-Johnson Transform
Yeo-Johnson는 Box-Cox을 좀 더 일반화 한 변환법이다. Box-Cox 보다 변환방법이 복잡하지만, Yeo-Johnson Transform은 실수 전체 구간에서 정의된 확률변수에 적용할 수 있다.
\[ x(\lambda) = \begin{cases} \begin{split} \cfrac{(y+1)^{\lambda} -1}{\lambda} &,& & & & & if \lambda \ne 0 &,& y \ge 0 \\ \ln (y+1) &,& & & & & if \lambda = 0 &,& y \ge 0 \\ \cfrac{ (1-y)^{2-\lambda} -1}{\lambda -2 } &,& & & & & if \lambda \ne 2 &,& y < 0 \\ - \ln(1-y) &,& & & & & if \lambda = 2 &,& y < 0 \end{split} \end{cases} \]
yj = function(x, lambda){
if(x >= 0){
if(lambda == 0){y = log(x+1) }
else{ y = ((x+1)^lambda-1)/lambda }
}
else{
if(lambda == 2){y = -1 * log((1-x))}
else{y = -1 * ( (1-x)^(2-lambda) - 1) / (2-lambda) }
}
return(y)
}
yj_list = function(x, l){ return(sapply(x, yj, lambda=l)) }
x_grid = seq(-2, 2, by=0.01)
lambda_grid = seq(-1,2, by=0.2)
result = data.frame(x=NA,lambda=NA,yj=NA)
for(l in lambda_grid){
yj_transform = data.frame(x=x_grid, lambda=l, yj = sapply(x_grid, yj, lambda=l))
result = rbind(result, yj_transform)
}
ggplot(result, aes(x=x, y=yj_value, color=lambda)) +
geom_point(size=0.5) +
scale_color_gradientn(colours = rainbow(5))
2. Parameter \( \lambda \) 추정
Box-Cox과 Yeo-Johnson Transform 에서 사용하는 변환식은 모두 parametric function이기 때문에, 적절한 \( \lambda \)를 설정해야 한다. 주어진 data를 사용해 \( \lambda \) 를 결정하는 방법을 알아보자. 대표적인 방법은 변환된 data \( Y \)의 Likelihood를 최대화 하는 \( \lambda \)값을 찾는 것이다. 여기에 Cross-validation을 응용해 parameter 추정시 과적합을 방지할 수 있다. 이번에는 Cross-valdiation 없이 전체 data에서 (Normal distribution에서) Likelihood를 최대화 하는 \( \lambda \)를 찾는 방법을 알아보자.
변환된 data \( Y \)가 정규분포를 따른다고 할때, Likelihood를 아래처럼 전개할 수 있다.
\[ L ( \lambda , \beta , \sigma^2 \mid X, y) = \cfrac{ \exp(- 1/2\sigma^2) (y ( \lambda) - X\beta )^{'}(y ( \lambda) - X\beta ) }{(2\pi\sigma^2 )^{n/2}} J( \lambda, y ) \]
이때 \( J( \lambda, y ) = \cfrac{d}{dy} y(\lambda) \) 는 \( X \rightarrow Y \) 로의 변수 변환에 따른 Jacobian matrix 이다. 두 가지 변환 방법별 log Jacobian은 아래와 같다.
Box-Cox Transform : \( (\lambda - 1) \sum_{i} ln(y_i) \)
Yeo-Johnson Transform : \( (\lambda -1) \sum_{i} sign(y_i) ln ( |y_i| + 1) \)
\( \beta, \sigma^{2} \)의 MLE는 \( \lambda \) 추정값과 무관하기 때문에 \( \hat{\beta}^{MLE} \), \( \hat{\sigma^2}^{MLE} \)를 먼저 구한 후, 두 추정값이 주어졌을때 \( \lambda^{MLE} \) 를 구한다. 두 추정값은 Linear model \( y(\lambda, g) \sim N(x\beta, \sigma^{2}) \) 에서 \( \beta, \sigma^{2} \)의 MLE와 동일하다.
\( \hat{\beta}^{MLE} (\lambda) = \hat{\beta} = (XX)^{-1}X^{'}y(\lambda) \)
\( \hat{\sigma^2}^{MLE} (\lambda) = \hat{\sigma^2} = y(\lambda) (I - X(XX)^{-1}X^{'}) y(\lambda)/n \)
따라서, 두 parameter의 MLE가 주어졌을때 \( \lambda \) 의 log Profile Likelihood 는 다음과 같다.
\( l( \lambda | \hat{\beta}, \hat{\sigma^2}) = - \cfrac{ (y(\lambda) - X\hat{\beta} )^{'}(y(\lambda) - X\hat{\beta}) }{2\hat{\sigma^2}} + ln(J(\lambda, y)) -\cfrac{n}{2} ln(2\pi \hat{\sigma^2}) \)
\( = - \cfrac{ y(\lambda)^{'}(I - H)y(\lambda) }{2\hat{\sigma^2}} + ln(J(\lambda, y)) -\cfrac{n}{2} ln(2\pi \hat{\sigma^2}) \)
\( = -\cfrac{1}{2} + ln(J(\lambda, y)) -\cfrac{n}{2} ln(2\pi \hat{\sigma^2}) \)
앞서 구한 두 가지 변환방법별 (log) Jacobian값을 적용하면,
(Box-Cox) = \( -\cfrac{1}{2} + (\lambda - 1) \sum_{i}^{n} ln(y_{i}) -\cfrac{n}{2} ln(2\pi \hat{\sigma^2}) \)
\( \propto \cfrac{n}{2} \sum_{i=1}^{n} \cfrac{2}{n} ln(y_{i}^{\lambda-1}) - \cfrac{n}{2} ln(\hat{\sigma^2}) \)
= \( - \cfrac{n}{2} ( - ln(\prod_{i=1}^{n} 2(y_{i}^{1/n})^{\lambda-1}) + ln(\hat{\sigma^2}) ) \)
= \( - \cfrac{n}{2} ( - ln[ \cfrac{\hat{\sigma^2}}{g^{2(\lambda - 1)}}]) \) where \( g = \prod_{i=1}^{n} y_{i}^{1/n}\)
= \( - \cfrac{n}{2} ( - ln[ \cfrac{y(\lambda, g)^{'} (I-H)y(\lambda, g)}{n}] ) \)
(Yeo-Johnson) = \( -\cfrac{1}{2} + (\lambda - 1) \sum_{i}^{n} sign(y_i) ln(y_{i}) -\cfrac{n}{2} ln(2\pi \hat{\sigma^2}) \)
\( \propto (\lambda - 1) \sum_{i}^{n} sign(y_i) ln(y_{i}) - \cfrac{n}{2} ln(\hat{\sigma^2}) \) ...(2)
= \( -\cfrac{n}{2} ln (\cfrac{y(\lambda)^{'}(I-H) y(\lambda)}{g^{\lambda-1} n g^{\lambda-1}}) -n \times ln(g^{\lambda - 1} ) \)
= \( -\cfrac{n}{2} ln( S^{2}_{\lambda}) - n(\lambda - 1)ln((\prod^{n}_{i} (y_i))^{1/n} ) \)
이때 \( g = \prod_{i=1}^{n} y_{i}^{1/n}\)는 \( y\)의 geometric mean 이며, \( y(\lambda, g) = y(\lambda)/g^{\lambda - 1}\)로 정의한다. 즉, Box-Cox의 Profile Log-likelihood에 대한 MLE를 찾는 문제는 Linear model \( y(\lambda, g) ~ N(x\beta, \sigma^{2}I) \) 의 Residual Sum of Sqaures(SSR) \( S_{\lambda}^{2} \)을 최소화 하는 \( \lambda \)를 찾는것과 동일하다.
다만 Python Scipy의 경우, MLE 추정시 식 (2) 을 목적함수로 하여 최적해를 구한다.
