sparklyr 에서 bucketing 하기 (기본 R과의 비교)

sparkly 에서는 spark 에서 실행하기 위한 함수들이 존재한다. 물론 그룹별 집계 등의 기본적인 데이터 처리의 경우, dplyr wrapper 를 통해 추가적인 학습 없이도 바로 사용할 수 있지만, 이외의 데이터 처리의 경우 sparklyr 의 고유 함수에 대한 학습이 필요하다. 

본 포스팅에서는 10분위수를 기준으로 변수를 10개의 카테고리로 분류하는 bucketing 작업을 sparklyr 을 통해 해보고 기본 R 과 비교해보았다. 

sparklyr 에서 bucketing 하기

get_spark_connection 의 경우 spark connection 을 잡는 함수인데, 관련해서 이 포스팅을 참고하길 바란다. 이 데이터처리에서 중요함수는 sdf_quantile 과 ft_bucketizer 이다. sdf 는 spark dataframe 의 약자이고, ft는 feature transformation 의 약자이다.

10분위수 구하기: sdf_quantile 을 사용하면 병렬적으로 변수의 quantile 을 구해서 리턴해준다. 이 작업은 데이터를 메모리로 로드후에 R 의 기본 quantile 을 실행하는 것보다 빠르다.

bucket 만들기: ft_bucketizer 는 연속형 변수에 대해 bucketizing 을 해준다. 이 때, 만약 10개의 bucket 으로 나눈다고 하면, splits 에 원소가 11개인 벡터를 전달해주어야한다. 10개의 bucket 으로 나눈다고 하면, 9개의 split point 가 필요한데, 여기에 최소값과 최대값을 추가로 더해 11개라고 볼 수 있다. quantile 함수에서 9개의 decile 과 최대값을 구했으므로, 최소값인 0을 추가로 더해 split 함수에 넣어준다. 

bucket 이름 만들기: ft_bucketizer 함수에는 아쉽게도 bucket 의 label 을 지정하는 방법이 없다. 단지 index 만 만들어줄 뿐이다. 따라서 index 와 label 을 매핑하는 데이터프레임을 따로 만들어서 해결해볼 수 있다. quantile 의 개수만큼 seq_along 을 통해 인덱스를 만들고, names 함수를 통해 벡터의 이름을 구해서 bucket_df 라는 데이터프레임을 만든 후에, 이를 데이터와 조인하여 bucket 의 label 을 만들 수 있다. 

library(tidyverse)
library(sparklyr)

sc <- get_spark_connection()

cars <- copy_to(sc, cars, "cars", overwrite=TRUE)

val_quantile <- sdf_quantile(cars %>% filter(speed > 0), column = "speed", probabilities = seq(0.1, 1, 0.1))
val_quantile <- val_quantile[!duplicated(val_quantile, fromLast = TRUE)]

val_quantile
#    10%       20%       30%       40%       50%       60%       70%       80%       90%      100% 
#   15900     27000     35000     49900     63800     86200    113800    158310    246000 608623010 

names(val_quantile)

bucket_df <- data.frame(bucket_index = seq_along(val_quantile) - 1, bucket = names(val_quantile))
colnames(bucket_df) <- c("speed_index", "speed_cat")

# bucketing: 10분위수 총 9개와 함께, 최솟값, 최댓값을 벡터로 만들어 전달. (총 원소가 11개인 벡터)
cars <- ft_bucketizer(cars,
                        input_col = "speed",
                        output_col = "speed_index",
                        splits = c(0,val_quantile))

# spark dataframe 으로 변경
bucket_df <- sdf_copy_to(sc, bucket_df, "bucket_df", overwrite = TRUE)
cars <- left_join(cars, bucket_df, by = "speed_index")
cars %>% select(speed, speed_index, speed_cat) %>% sample_n(50)

기본 R 에서 10분위 bucketing 하기

기본 R dml

val_quantile <- quantile(cars %>% filter(speed > 0) %>% pull(speed), probs=seq(0.1, 1, 0.1))
val_quantile
# 10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
#   8.9 11.0 12.7 14.0 15.0 17.0 18.3 20.0 23.1 25.0 


cars$speed_cat <- if_else(cars$speed == 0, "X", 
                          as.character(cut(cars$speed, breaks = c(0,  val_quantile), labels = names(val_quantile))))
cars %>% sample_n(50)

# speed dist speed_cat
# 1     14   60       40%
# 2     13   26       40%
# 3     19   36       80%
# 4     13   34       40%
# 5     12   24       30%
# 6      4   10       10%

R Standard evaluation 과 Tidy evaluation

R 에서 특정 데이터프레임(df) 를 특정 열(var1,var2) 로 정렬하기 위해서는 아래와 같이 할 수 있습니다.

df %>% arrange(var1, var2)

Standard evaluation 

그런데 변수 이름이 동적인 경우, 즉 열 이름이 항상 같지 않고 코드 런타임에 결정되는 경우에는 표준 평가(Standard Evaluation, SE) 방법을 사용해야 합니다.

var1 <- "col1"
var2 <- "col2"

df %>% arrange(!!sym(var1), !!sym(var2))

Tidy evaluation

tidy evaluation 의 경우 rlang 의 {{}} (curly-curly) operation 을 사용하는 방법입니다. 그런데, 이 방법은 함수 안에서만 사용할 수 있습니다. 즉, tidy evaluation 도 동일한 기능을 수행하지만, 아래와 같이 함수 안에서만 사용할 수 있다는 제한점이 있습니다. 

my_function <- function(df, var){
  df %>% arrange({{var}})
}

티스토리 블로그 하단 자동 광고 제거하기

글 작성일: 2024년 6월 14일 기준

언제부턴가 티스토리 블로그 하단에 앵커 광고가 뜨기 시작했습니다. 해당 광고는 위치가 고정되어 스크롤 업/다운을 하더라도 화면을 가려서 글을 읽는데 불편했는데요. 

현재 티스토리 관리자 페이지에서 해당 광고 제거하는 것은 불가능하고 구글 애드센스에 직접 접속해서 해당 광고는 제거하는 설정이 필요합니다.

1) 애드센스 접속후 왼쪽 사이드바 "광고" 메뉴 클릭

2) 수정할 블로그 수정 버튼 클릭

3) 자동 광고 토글 Off

이렇게 설정하고 완료를 누르면 최대 1시간 뒤에 반영된다는 팝업이 뜨는데요. 제 경우에는 바로 광고 제거가 적용되지 않았고, 1시간 뒤쯤 접속해서 확인해보니 앵커 광고가 사라졌습니다. 

기본적인 2x2 테이블 

이진 분류 모델의 최종 예측 값은 일반적으로 0~1사이로 나오게된다. 특정 임계치를 기준으로 테스트 양성과 테스트 음성을 분류한다. 예를 들어, 임계치를 0.5로 잡는다면, 0.5 이상인 경우를 테스트 양성, 0.5 미만인 경우를 테스트 음성으로 정의한다. 이 방법을 통해 아래와 같은 2x2 테이블을 만들 수 있다. 

• True Positive (TP): 50
• False Positive (FP): 10
• False Negative (FN): 5
• True Negative (TN): 100

Sensitivity = Recall (민감도)

sensitivity 는 실제 질병인 사람 중에 테스트 양성인 사람의 비율이다. 

-> 50/55 = 0.909

Specificity = Negative Recall (특이도)

specificity 는 실제 질병이 아닌 사람 중에 테스트 음성인 사람의 비율이다. 

-> 100/110 = 0.909

Positive Predictive Value = Precision (양성 예측도, PPV) 

ppv 는 양성으로 예측한 사람 중에 실제 질병인 사람의 비율이다.

-> 50/60 = 0.833

Negative Predictive Value (음성 예측도, NPV) 

npv 는 음성으로 예측한 사람 중에 실제 질병이 아닌 사람의 비율이다. 

-> 100/105 = 0.952 

ROC 커브와 AUC (Area under curve)  

임계치를 변화시키면서 1-specificity, sensitivity 그래프를 그린 것이 ROC 커브이다. 위 두 지표를 통해 그래프를 그리는 이유는 sensitivity 와 specificity 간에 트레이드오프관계가 있기 때문에, 이 관계를 시각적으로 표현하여 모델의 성능을 평가하기 위해서이다. 

1-specificity = False Positive Rate (FPR)

sensitivity = True Positive Rate (TPR)

임계치가 낮아지면 모델은 더 많은 사례를 양성으로 분류하게 되어 True Positive Rate와 False Positive Rate가 모두 증가한다. 반대로, 임계치가 높아지면 모델은 더 적은 사례를 양성으로 분류하게 되어 True Positive Rate와 False Positive Rate가 모두 감소한다. 

회귀 모델의 선택

1) 두 변수가 nested 관계에 있을 때

--> Likelihood ratio test 를 통해 유의미하게 좋은 모델을 선택한다. 

두변수가 nested 관계에 있다는 것은 full model, reduced model 관계에 있다는 것을 의미한다. 

$$ Model1: g(\pi_i) = \beta_0 + \beta_1x_{1i} $$

$$ Model2: g(\pi_i) = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} $$  

만약 모델 2에서 beta2가 0인 경우, 모델1 이 된다. 따라서 두 모델은 nested 관계에 있다. 

이 경우 변수가 많은 모델2가 무조건 likelihood 가 높게 된다. 

이 때, Likelihood ratio 가 카이제곱분포를 따르게된다. 자유도는 두 모델의 모수 개수의 차이가 된다. 여기서 L0가 간소한 reduced 모델이고, L1이 full model 이다. 이는 full model 과 reduce model 의 log likelihood 의 차이의 2배이다. 

$$ -2ln(\frac{L_0}{L_1}) = -2(lnL_0 - lnL_1) \sim \chi(1)$$

만약 위 통계량이 유의미한 카이제곱 값을 가지면, full model 이 reduced model 보다 좋은 것이다. 따라서 full model 을 채택한다. 만약 카이제곱값의 p-value 가 0.05보다 크다면, full model 이 reduced model 보다 좋지 않은 것이므로, reduced 모델을 채택한다. 

2) 두 변수가 unnested 관계 일 때

--> AIC 가 작은 모델을 선택한다.

$$ Model1: g(\pi_i) = \beta_0 + \beta_1x_{1i} $$

$$ Model3: g(\pi_i) = \beta_0 + \beta_1x_{2i} + \beta_2x_{3i} $$  

이러한 경우 AIC 를 통해 두 모델 중 어떤 모델이 좋은지를 판단할 수 있다. p는 변수의 숫자로 패널티텀이다. unnested 관계일 때, likelihood ratio test 를 적용할 수 없는 이유는 unnested 관계일 때는 likelihood ratio 가 카이제곱분포를 따르지 않기 때문이다. 

$$ AIC = -2(Log L_1 - p) $$ 

만약 모델1 의 로그 우도가 -120 이고, 모델3의 로그 우도가 -115 라고하자. 모델3의 로그 우도가 더 높다.  

모델1의 AIC = -2(-120 - 2) = 244

모델3의 AIC = -2(-115 - 3) = 236 

만약 모델3의 추정 파라미터 수가 2개였다면, AIC 는 234였을 것이다. 파라미터로 인한 패널티 2점이 들어갔음을 알 수 있다. 이 경우 패널티를 고려해도 모델3의 AIC가 낮기 때문에 모델3을 채택한다. 

R - 롱테일 분포의 히스토그램 그리기

실무를 하다보면 롱테일 분포를 많이 접하게 됩니다. 예를 들어서, 어떠한 이커머스 서비스에서 "구매 금액" 이라는 변수를 살펴보면, 대부분의 유저는 구매금액이 0~1만원 사이에 들어있지만, 일부 유저는 구매금액이 몇 백만원 심지어는 몇 억원에 이르는 경우를 심심찮게 볼 수 있습니다. 극심한 right-skewed 분포 (또는 롱테일 분포)의 예라고 볼 수 있습니다. 

이러한 롱테일 분포에 일반적인 히스토그램을 적용하게 되면 꼬리가 너무 길어져 가시성이 좋지 않습니다. 이런 경우에 특정 cutff 지점을 정해 따로 범주를 만들곤 합니다. 예를 들어, 구매금액이 백만원 이상인 유저는 '100만원 이상' 이라는 bucket 을 따로 만드는 것이죠. 꼬리 부분이 너무 길기 때문에 이 부분을 따로 모으는 것입니다. 

R 코드로는 다음과 같이 작성해볼 수 있습니다. 포인트는 raw 데이터에 적용하는 geom_hist 를 사용하는 것이 아니라, 집계 데이터를 먼저 만든 후, geom_bar 를 통해 히스토그램을 그리는 것입니다. 그리고, 집계 데이터를 만들기 위해 cut 함수를 사용합니다. 

R 코드

• anal_table 데이터 프레임의 value 컬럼이 histogram 을 그리고자하는 변수입니다.

top_1_percent <- quantile(anal_table$value, 0.99, na.rm=T) # 상위 1% 경계값 찾기

# bucket size 동적으로 설정
bucket_size <- 10^ceiling(log10(top_1_percent)) # 초기 bucket size

# while loop를 통해 bucket size 조정
while(TRUE) {
  breaks <- seq(0, top_1_percent, by = bucket_size) # 상위 1% 까지의 bucket  
  if(length(breaks) > 100) break # bucket 개수가 100개 이상이면 loop 탈출
  bucket_size <- bucket_size / 10 # bucket size 재조정
}

labels <- breaks
cutoff <- max(labels)+bucket_size
# 기본적으로 break 에서 좌측을 포함하지 않고 우측을 포함함(include lowest 를 통해 가장 좌측은 포함)
# right=FALSE 를 통해 우측을 포함하지 않게 지정
anal_table$bucket <- cut(anal_table$value, breaks = seq(0, cutoff, by = bucket_size), 
                         include.lowest = TRUE, 
                         right=FALSE,
                         labels = labels)

# bucket 이 없는 경우는, cutoff 이상인 경우로, 따로 만든 bucket 에 속하도록 바꾸어줌 
anal_table <- anal_table %>% mutate(bucket = if_else(is.na(bucket), as.character(ceiling(cutoff)), bucket))
anal_table$bucket <- factor(anal_table$bucket, levels = c(labels, ceiling(cutoff)))

summary_data <- anal_table %>% group_by(bucket) %>% count()
summary_data

summary_data <- summary_data %>% mutate(var_name = var_name)

val_quantile <- quantile((anal_table %>% select(value) %>% pull), probs=seq(0.1, 1, 0.1))

quantile_keys <- names(val_quantile)
quantile_values <- unname(val_quantile)

df_quantile <- data.frame(t(quantile_values))
colnames(df_quantile) <- quantile_keys

df_avg <- anal_table %>% summarize(avg = mean(value))
df_quantile <- cbind(df_quantile, df_avg)
df_quantile <- df_quantile %>% mutate(var_name = var_name)

total_ticks <- 10  

breaks <- pretty_breaks(n = total_ticks)(range(as.numeric(as.character(summary_data$bucket))))
ggplot(summary_data, aes(x = as.numeric(as.character(bucket)), y = n)) +     
  scale_y_continuous(labels = scales::label_comma()) +     
  geom_bar(stat = "identity", fill = "black") +    
  labs(x = "X", y = "Y") +    
  scale_x_continuous(breaks = breaks,  # breaks는 pretty_breaks를 사용해 계산된 값
                     labels = breaks) +   # labels도 breaks를 사용
  theme_bw(base_size = 10, base_family = "Kakao Regular") +    
  ggtitle("Histogram from Binned Data") +      
  theme(plot.margin = margin(0.5, 0.5, 0.5, 0.5, "cm")) +  
  geom_vline(aes(xintercept = df_quantile$avg), colour = "red") +  
  annotate("text", x = df_quantile$avg, y = max(summary_data$n),  
           label = paste("평균 =", round(df_quantile$avg, 2)),  
           vjust = 2, color = "black", size=3)

결과 히스토그램

• 위 코드를 통해 아래와 같이 지정된 bucket size 를 가지며, 상위 1% 이상은 하나의 bucket 으로 묶은 깔끔한 히스토그램을 그릴 수 있습니다.

위 코드에는 몇 가지 포인트가 있습니다. bucket size(bin)과 xtick 의 개수를 동적으로 결정한 부분인데요. 이 부분 코드를 좀 더 살펴보겠습니다. 

bucket size 를 동적으로 결정하기

• bucket 의 개수가 최소 100개가 되도록 하며, bucket size 가 1, 10, 100, 1000 처럼 10의 지수형태로 만드는 방법은 아래와 같습니다. 
• 또한 cut 함수의 labels 를 통해 label 을 이쁘게 만들어줍니다. 예를 들어 label 이 100이라고 하면, (100~199 사이의 bucket 을 의미하는 등)

# while loop를 통해 bucket size 조정
while(TRUE) {
  breaks <- seq(0, top_1_percent, by = bucket_size) # 상위 1% 까지의 bucket  
  if(length(breaks) > 100) break # bucket 개수가 100개 이상이면 loop 탈출
  bucket_size <- bucket_size / 10 # bucket size 재조정
}

labels <- breaks
cutoff <- max(labels)+bucket_size
# 기본적으로 break 에서 좌측을 포함하지 않고 우측을 포함함(include lowest 를 통해 가장 좌측은 포함)
# right=FALSE 를 통해 우측을 포함하지 않게 지정
anal_table$bucket <- cut(anal_table$value, breaks = seq(0, cutoff, by = bucket_size), 
                         include.lowest = TRUE, 
                         right=FALSE,
                         labels = labels)

동적 xtick 의 결정

• 아래 코드는 총 xtick 의 개수를 10개로 고정시키고, 변수에 따라 동적으로 xtick 간격을 조정하는 코드입니다. scalse 라이브러리의 pretty_breaks 라는 함수를 사용합니다. 

library(scales)
total_ticks <- 10  
breaks <- pretty_breaks(n = total_ticks)(range(as.numeric(as.character(summary_data$bucket))))

일반화 선형 모형의 개념 (Generalized Linear Model)

일반화 선형 모형의 식은 아래와 같다. 

$$ g(\mu) = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... \beta_px_p $$ 

x1, x2 ... 가 주어졌을 때, Y를 예측하고 싶다. 근데 특정 조건하에서 Y 는 정해진 값이 아니라 어떤 분포를 따른다고 가정하고, 그 평균을 예측하고 싶을 때, 일반화 선형 모형을 활용한다. 기본적인 회귀분석에서는 반응변수가 정규분포를 따른다고 가정하고 모델링하는데, 일반화 선형 모형은 Y가 다른 분포를 따르는 경우에도 활용할 수 있는 모델링 방법이라고 볼 수 있다. 

일반화 선형모형에서는 반응변수가 어떤 분포를 따른다고 가정하기 때문에 랜덤성분 (random component) 이라고 부르고, 반응변수의 평균을 설명하기 위한 설명 변수들의 함수 (위 식에서 우측부분) 를 체계적 성분 (systematic component) 이라고 부른다. 랜덤성분과 체계정 성분을 연결하는 함수를 연결함수(link function) 라고 부른다.

Y가 정규분포를 따른다면, 평균값이 -무한대~+무한대일 수 있고, Y 가 베르누이 분포를 따르면 Y의 평균이 0~1사이의 값이다. 따라서 적당한 연결함수를 통해 값의 범위를 변환하는 것이 필요하다. 

또한, 일반화 선형 모형에서는 Y 가 지수족 분포를 따른다고 가정한다. 지수족 분포에는 정규분포, 이항분포, 포아송분포, 감마분포 등이 있다. Y가 따른다고 가정한 분포에 따라 알맞는 연결함수를 적용해준다. GLM 에서 지수족 분포가 중요한 개념이지만, 다소 심플하게 내용을 설명하기 위해 지수족 관련 내용은 설명하지 않겠다. 

만약, Y가 정규분포를 따르는 경우에 가장 기본적으로 항등함수를 이용할 수 있다. 연결함수가 항등함수인 경우, 일반 선형 모형이라고 한다. (general liner model) (generalied linear model 과 다르다.). 

$$ \mu = \alpha + \beta x $$ 

연결함수가 항등함수인 경우 beta 값의 해석은 매우 쉽다. "X가 1단위 증가했을 때 반응 변수가 beta 만큼 증가한다" 고 해석한다. 

Y가 베르누이분포를 따르는 경우 0~1의 값을 무한대로 변환하는 연결함수로 여러가지를 이용할 수 있다. 가장 대표적인 것이 로짓함수이다. 로짓함수를 사용한 변수가 1개인 일반화 선형모형은 아래와 같이 정의된다. 이를 로지스틱 회귀분석 (logistic regression) 이라고 부른다.  

$$ log(\frac{\mu}{1-\mu}) = \alpha + \beta x $$ 

좌측을 살펴보면 log odds 라는 것을 알 수 있다. (=> log(성공확률/실패확률) 이므로) 즉, 로지스틱 회귀분석은 log odds 를 설명변수들의 조합으로 예측하는 것을 의미한다. odds 가 아닌 확률(평균) 의 관점에서 로지스틱 회귀분석은 아래와 같이 써볼 수 있다. 

$$ \mu = \frac{exp(\alpha + \beta x)}{1+exp(\alpha+\beta x)} $$

또한 로지스틱 회귀 분석에서 중요한 것은 beta 값의 해석이다. 만약 x가 연속형인 경우 x+1과 x의 odds 를 구해서 odds ratio 를 구해보자. 위 식에 넣어 계산해보면, OR = exp(beta) 가 나온다. 양변에 log 를 취해주면 log(OR) = beta 라는 것을 알 수 있다. 즉, x가 1단위 증가했을 때의 log(OR) 값이 beta 라는 것을 알 수 있다.

한편, Y가 베르누이 분포를 따르는 경우에 사용할 수 있는 다른 연결함수로는 프로빗 연결함수가 있다. 프로빗 연결함수를 사용한 일반화 선형 모형을 프로빗 모형이라고 부른다. 프로빗 모형은 표준정규분포의 누적분포함수의 역함수를 연결함수로 사용한다. 누적분포함수의 역함수를 연결함수로 사용한다는 의미가 무엇일까? 누적분포함수는 0~1사이의 값을 갖는다. 즉, 어떤 -무한대~무한대에 있는 X라고하는 값을 0~1 사이로 변환하는 함수이다. 이에 역함수이기 때문에 0~1사이의 값을 -무한대~무한대로 바꾸어주는 함수가 된다. 

도수 자료의 경우에는 일반화 선형모형중 포아송 회귀분석을 해볼 수 있다. 도수 자료란 반응 변수가 도수 (count)로 이루어진 자료를 의미한다 (예를 들어, 교통사고 수, 고장 수 등...). 도수자료는 양의 방향에서만 존재한다. 교통사고수가 마이너스일 수는 없다. 반면, 설명변수의 조합인 체계적 성분은 -무한대~무한대의 범위를 갖는다. 이를 변환하기 위해서, 포아송 회귀분석에서는 연결함수로 log 를 활용하여 좌변이 -무한대~무한대의 값을 갖도록 변환한다. 포아송 회귀분석 식은 아래와 같다. 

$$ log(\mu) = \alpha + \beta x $$ 

이는 평균의 관점에서는 아래와 같이 쓸 수 있다.

$$ \mu = exp(\alpha + \beta x) $$

x가 t일 때와 t+1일때의 mu 값을 비교해보자. 위 수식에 대입하면 x가 t+1 일 때의 mu 와 t 일 때의 mu 의 ratio 는 exp(beta) 가 됨을 알 수 있다. 즉, 포아송 회귀분석과 같은 log linear regression 에서 beta 를 해석하는 방법은 "x 가 1단위 증가했을 때 Y값의 평균이 exp(beta)배 증가한다." 이다. 

포아송 회귀 관련해서는 종종 이런 문제가 발생할 수 있다. 만약, X가 차량 사고수에 미치는 영향을 포아송 회귀로 모델링을 하려고하는데, 지역별로 데이터가 수집 되었고, 지역별로 기본적인 차량의 개수가 달라 사고수가 이에 영향을 받는다고 해보자. 이 때, "사고율" 을 반응 변수로해서 모델링할 수 있다. 차량의 개수를 t라고 하자.

$$ \log(\mu / t) = \alpha+\beta x $$

사고수의 관점에서 아래와 같은 수식으로 변환할 수 있다. 이 때, log(t) 를 offset 이라고 한다. 

$$ \log(\mu) = \log(t)+\alpha+\beta x $$

$$ \mu = texp(\alpha + \beta x) $$ 

선택편향

선택편향은 특정 그룹을 선택해서 분석했을 때, 다른 그룹 또는 전체를 대상으로 분석했을 때와 다른 결론이 나오는 것을 의미한다. 아래와 같이 왼쪽 그림에서는 X,Y 의 연관성이 없지만, X+Y가 1.2 이상인 그룹만 선택해서 봤을 때는 X,Y의 음의 상관성이 생기는 것을 알 수 있다. 이러한 선택 편향은 우리의 실생활에서도 많이 발생한다. 

Collider bias

Collider bias는 X와 Y가 모두 영향을 미치는 Z라고 하는 변수가 있을 때, Z를 고정시켜 놓고 보면, X (exposure) 과 Y (outcome) 에 연관성에 편향이 생기는 현상을 의미한다. 

왜 Collider bias 가 발생할까? 이에 대해 사고적으로 이해하는 방법에는 "explaining away" 라고 하는 개념이 있다. 예를 들어, X 를 통계학 실력이라고 하고, Y를 아첨 능력이라고 하자. 그리고 X,Y 가 모두 승진 (Z) 에 영향을 준다고 해보자. 이 때, 승진 대상자만을 놓고 통계학 실력과 아첨 능력의 관계를 보면 둘 사이에는 음의 상관성을 확인할 수 있다. (이는 정확히 위 selection bias 에서 설명하는 그림과 같다.) 

이처럼 실제로는 통계 실력과 아첨 능력에는 아무런 상관성이 없으며, 승진에 영향을 주는 원인 변수일 뿐인데, 승진 대상자를 놓고 봤을 때는 둘 사이에 연관성이 생긴다 (false association). 승진한 어떤 사람이 아첨능력이 매우 좋다고 했을 때, 이것이 승진의 이유를 explain 해주므로, 이 사람의 통계학 실력은 좋지 않을 것이라고 '추정' 할 수 있을 것이다. 또한, 어떤 사람이 통계 실력이 매우 좋지 않음에도 불구하고 승진했을 때, 이 사람은 아첨 능력이 뛰어날 것이라고 추정할 수 있다. 이처럼 둘 사이에 음의 상관성이 존재하는 것을 직관적으로 이해할 수 있다.

CMH 검정과 통계량 계산 방법

범주형 자료 분석에서 코크란-멘텔-헨젤(Cochran-Mantel-Haenszel) 검정의 목표는 Z 가 주어질 때, X와 Y가 조건부 독립인지를 검정하는 것이다. 즉, Z를 고려했을 때, X-Y의 연관성이 존재하는지를 판단하는 검정이라고 할 수 있다. 이는 인과추론에서 말하는 X,Y가 조건부 독립 (conditional independence) 인지를 확인하는 검정이라고 할 수 있다. 보통 Z는 confounder 로 설정하는 경우가 많다. 만약, conditional independence 가 아니라고 한다면, Z 를 고려함에도 X-Y 연관성이 존재하는 것이고, 이는 X,Y 의 인과성에 대해 조금 더 근거를 더해준다고 할 수 있다.  CMH 검정은 2 X 2 X K 표에 대해서 활용할 수 있다. (K 는 Z의 수준 개수)  

그룹 i 에서의 흡연과 폐암의 연관성

주요 지표

n = a+b+c+d

p1 = (a+b)/n (흡연X 비율)

p2 = (a+c)/n (폐암X 비율) 

m = n*p1*p2

CMH 통계량의 계산

그룹 i 에서의 CMH 통계량은 아래와 같다. 

$$ \frac{(a-m)^2}{m(1-p_1)(1-p_2)} $$

최종적인 CMH 통계량은 모든 그룹 i에서 위 값을 다 구해서 더한 것이다. 이 값은 자유도가 1인 카이제곱분포를 따른다는 것을 이용해 검정한다. 만약, 충분히 이 값이 큰 경우 그룹을 고려했을 때, 흡연과 폐암에 연관성이 있다고 결론을 낼 수 있다. 

위 수식에서 a-m 은 관측값에서 기대값 (평균) 을 빼준 것이고, 분모는 a의 분산을 의미한다. 이 분산은 초기하분포의 분산이다. 즉, cmh 통계량에서는 a가 초기하분포를 따른다고 가정한다. 즉, 수식은 a 에서 평균을 빼주고 표준편차로 나눈 값에 제곱이라고 할 수 있다. 

MH 공통 오즈비

그룹1

=> OR = 10*40 / 20*30 = 2/3

그룹2 

=> OR = 4*4 = 16 

1) 두 그룹의 공통 오즈비를 구하는 방법에는 단순히 두 그룹의 오즈비의 평균을 구하는 방법이 있을 수 있다. 이 경우 그룹2의 샘플수가 적음에도 불구하고 평균 오즈비는 8에 가깝게 높게 나온다. 

2) a*d 의 값을 모두 더한 값을 b*c 를 모두 더한 값으로 나누어주는 방법이 있다. 이러면 (10*40 + 4*4) / (20*30+1) = 0.69 가 나오게 된다. 이 값은 샘플수가 많은 그룹의 값으로 지나치게 치우친다. 

3) MH 공통 오즈비는 중도적인 방법으로 두 방법의 단점을 보완한다. 2) 방법에서 샘플수의 역수로 가중치를 줌으로써, 샘플수가 많은 그룹이 계산에 미치는 영향력을 의도적으로 줄여준다. 

(10*40/100 + 4*4/10) / (20*30/100 + 1/10) = 0.91 

즉, MH 공통 오즈비를 사용하면, 지나치게 그룹1에 치우치지 않으면서 적당한 공통 오즈비가 추정된다. 또한, 로그 MH 공통 오즈비의 분산을 계산할 수 있기 때문에, 공통 오즈비의 신뢰구간 및 오즈비가 유의미한지를 추론할 수 있다는 장점이 있따. 

예를 들어, 공통 오즈비가 0.91인 경우 로그 공통 오즈비는 -0.094이다. 그리고, 로그 공통 오즈비의 표준편차를 예를 들어 0.02라고 하자. 그러면 공통 오즈비의 95% 신뢰구간은 아래와 같이 계산된다. 

[exp(-0.094-1.96*0.02), exp(-0.094+1.96*0.02) ] = [0.88, 0.95] 

유전학에서의 딥러닝 활용이 정밀의학에 어떻게 기여하는가?

유전학 분야에서 딥러닝의 발전은 정밀 의학(personalized medicine) 에 구체적으로 어떻게 기여할 수 있을까?

1. 질병에 영향을 주는 유전적 변이 찾기 : 정밀 의학의 한가지 목적은 개인의 질병에 대한 위험도를 정밀하게 추정함으로써, 질병의 조기 발견 및 예방을 하고자하는 것이다. 그리고 그 중심에 있는 것이 과거엔 분석이 어려웠던 유전 정보라고 할 수 있다. 딥러닝 모델은 대규모 유전 데이터에서 유의미한 연관성을 발견하는 것에 기여한다. 예를 들어, 딥러닝을 활용하면 유전자 변이와 특정 질병 간의 관계를 더욱 잘 파악할 수 있다. 더욱 잘 파악한다는 것은 무슨 의미일까? 대표적으로 유전적 변이간의 교호작용 (interaction) 을 예로 들 수 있다. 교호작용이란 쉽게 말해 '시너지' 이다. 에를 들어, A 라는 유전변이가 질병 위험도에 3만큼 기여하고, B 라는 유전변이가 질병 위험도에 5만큼 기여한다고 하자. A,B변이가 모두 있는 사람이 질병 위험도가 30이 증가한다고 하면 기대치 8보다 22높은 값이다. 이런 경우 유전적 변이간에 교호작용 (gene-gene interaction) 이 있다고 한다. 일반적인 통계적인 방법으로도 이를 찾을 수 있지만, 경우의 수가 너무 많아 computational cost 도 크며, 실제 존재하는 interaction 을 잘 찾아내지 못할 가능성 (낮은 statistical power) 도 높다고 알려져 있다.  

딥러닝은 이러한 interaction 을 detection 하는데 더 효율적이라고 알려져 있다. 따라서, 개인의 유전 정보 기반 질병의 위험도 평가를 더욱 정확하게 할 수 있고, 이는 질병의 조기 발견 및 예방에 기여할 수 있다. 참고로, 딥러닝에서 유전자 변이와 질병간의 연관성을 파악할 때는, SNP 데이터에 feature engineering 방법 (예를 들면, PCA) 등을 적용해 차원 축소를 하고, 모델의 input 으로 넣는 방법이 많이 사용된다. 

2. DNA 의 전사 (Gene expression) 에 영향을 주는 유전적 변이 찾기: 어떠한 유전자 변이가 질병에 영향을 주는 대표적인 경로는 유전자 변이가 유전자 발현(gene expression)에 영향을 주고, 이 유전자 발현의 영향이 질병에 영향을 주는 것이다. 이에, 반응변수(y) 를 질병이 아닌 gene expression 등으로 두고, gene expression 에 영향을 주는 유전자 변이를 찾는 연구가 많이 이루어지고 있다. 보통 coding-variant 의 경우 해당 variant 가 gene expression 에 영향을 준다는 것을 비교적 쉽게 파악할 수 있다. 그러나 문제는 genome 에 대부분을 차지하는 non-coding region 에 위치한 variant 라고 할 수 있다.  딥러닝을 통해 non-coding variant 에 대한 정보(annotation) 을 쌓아, 이를 GWAS 의 결과를 해석하는데 사용할 수 있다.

보통 질병에 영향을 주는 유전적 변이를 찾는 과정에서는 SNP array 등을 많이 사용하는데, gene expression 에 영향을 주는 변이를 딥러닝을 통해 찾는 과정에서는 sequence data (ATGC.... 와 같은) 를 직접적으로 input 으로 넣는 경우가 많다. SNP array 를 사용했을 때와 비교하여 sequence data 를 사용하는 경우, 정보의 손실 (insertion/deletion 등)이 적기 때문일 것이다. 이는 질병에 인과적인 영향을 주는 causal variant 를 찾는 과정에 도움을 주기 때문에 유전적 리스크를 평가하는데 도움을 줌으로써 정밀 의학에 기여할 수 있다. 

3. 약물 반응 예측: 정밀 의학의 다른 목표 중 하나는 맞춤형 약물이라고 할 수 있다. 어떤 사람 A 에게는 잘 듣는 약물이 B 라는 사람에게는 잘 안들을 수 있다. 지금까지는 '평균적으로 잘 working 하는 약물' 을 모든 환자에게 투약하는 방식으로 치료 등이 이루어졌다면, 정밀의학 시대에서는 개인에게 잘 맞는 약물을 투약하는 것이 목표라고 할 수 있다.  딥러닝은 환자의 유전적 프로파일을 바탕으로 약물 반응성을 예측할 수 있다 따라서 특정 약물에 대한 환자의 반응을 예측하고, 부작용의 가능성을 최소화하는 데 도움을 줄 수 있다. 이 때의 input 은 genetic data (SNP array, sequence) 등이 될 것이다. 반응변수y는 약물 반응성이 될 것이다. 방법론적 측면에서 보자면 '질병 위험도 예측' 과 '약물 반응성 예측' 은 거의 비슷하다고 볼 수 있을 것이다.