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는 약물 반응성이 될 것이다. 방법론적 측면에서 보자면 '질병 위험도 예측' 과 '약물 반응성 예측' 은 거의 비슷하다고 볼 수 있을 것이다. 

rstudio server 에서 github copilot 사용하기

rstudio server 환경에서 github copilot 을 사용하려고 보니 아래처럼

Github Copilot integration has been disabled by the administrator 라는 문구가 떴다. 

알아보니, Rstudio Server나 Posit Workbench에서는 관리자 설정 이후 사용이 가능하다고 한다. 

관리자 설정 하는 방법은 아래와 같다. rsession.conf 파일을 열고 copilot-enabled=1 문구를 추가해주면 된다!

cd /etc/rstudio/
sudo vim rsession.conf

이렇게 Enable github copilot 옵션이 잘 나온다. 

폰트 다운로드 후 위치 변경

sudo mv [폰트 파일 이름].ttf /usr/share/fonts/truetype/

extrafont 를 이용해 R 환경에 폰트 설치

library(extrafont)
font_import(prompt=FALSE)  # 폰트 설치 디렉토리를 검색하여 사용 가능한 폰트를 가져옵니다.
fonts()  # 사용 가능한 폰트 목록을 출력합니다.

교차표에서 효과를 추정하는 방법은 아래 3가지가 있다.

교차표를 본다는 것은 범주형으로 이루어진 X,Y 의 변수간의 연관성을 파악하고 싶다는 의미이다. 

1) 비율의 차를 이용한 방법

2) 비율의 비를 이용한 방법

3) 오즈비를 이용한 방법 

흡연과 폐암의 관계

1) 비율의 차를 이용하는 방법

흡연자중 폐암 비율 = 90/1000 = 0.09

비흡연자중 폐암 비율 = 10/1000 = 0.01 

비율의 차 = 0.08 (risk difference)

2) 비율의 비를 이용하는 방법

0.09 / 0.01 = 9 (relative risk)

3) 오즈비를 이용하는 방법

(90/100) / (10/990) =  9.79 (odds ratio)

또는 p1 / (1-p1) / p2 / (1-p2) = ((90/1000) / (1-90/1000)) / ((10/1000) / (990/1000)) = 9.79 (odds ratio)

각 효과 추정 방법의 특징 및 장단점

비율의 비와 오즈비는 어떤 treatment의 효과(effect) 를 설명할 때 좋다. 비율의 비는 특히, 설명할 때 좋다. 

- 예를 들면, 어떤 위험인자가 질병에 미치는 영향이 있는지를 설명할 때는 비율의 비를 활용하는 것이 좋다. 

- 흡연의 reltavie risk 가 3이라는 말은 흡연을 하면 폐암 발생 위험을 3배 높인다고 해석할 수 있다. 

비율의 차는 전체 모수에서의 impact 를 설명할 때 좋다.

- 어떤 요인 A의 risk difference 는 10% 인데 relative risk 는 2라고 하자.

- 어떤 요인 B 의 risk difference 는 1%인데 relative risk 는 10이라고 하자.

- 이 때, 요인 B 의 effect size는 더 크지만, 실제 요인의 중요도는 A가 더 클 수 있다.  

오즈비는 y=1의 비율이 적을 때, 상대위험도와 값이 유사하다. 

- 만약의 y가 폐암과 같이 질병인 경우, P(Y=1) 은 유병률이다. 

- 즉, 유병률이 작은 질병의 경우 오즈비를 relative risk 처럼 해석할 수 있다.  

오즈비는 샘플이 불균형하게 추출한 경우에도 사용할 수 있는 지표이다. 

- 비율의 차 또는 비율의 비는 샘플링 바이어스의 영향을 받는다. 

- 만약, 흡연자 100명, 비흡연자100명을 선정해서 폐암여부를 비교할 때 비율의 비(relative risk) 에는 bias 가 생긴다. 이는 모집단에서 계산한 값과 차이가 생긴다는 의미이다. 

- 그러나, 오즈비의 경우 모집단에서의 값과 오즈비와 비교하여 bias 가 없게 된다. 

- 왜 오즈비는 샘플링 영향이 없는지 관련해서는 이 포스팅을 참고할 수 있다. 

ratio 에 로그를 취한 값은 유용하다.  

- 비율의 비 또는 오즈비는 매우 skew 된 값이다. 

- ratio는 0에서 무한대의 값을 갖는다.

- ratio 에 log 를 취해주면 -무한대~ +무한대의 값을 갖게 된다. 

- 만약, 어떤 A약의 효과가 B약의 효과보다 1.5배 있다 라는 것을 반대로 말하면 B 약의 효과가 비 약의 효과보다 1/1.5배 = 0.67배 있다라는 것이다. 그러나, 1.5배와 0.67배가 한눈에 역수 관계에 있다는 것을 알기 어렵다. 1.5배는 1로부터 0.5 떨어져 있고, 0.67은 1로부터 0.33 떨어져있다. 만약, 1.5와 0.67에 log를 취해주면, 각각 0.405, -0.405로 나오게 되어, 역관계에 있다는 것을 바로 확인할 수 있다. 

효과가 유의미한지 보려면 어떻게 할까? 

-> 신뢰구간을 보고, 이 값이 0을 포함하지 않으면 유의미하다고 판단할 수 있다.

비율의 차의 신뢰구간

$$ p_1 = n_{11}/n_1 $$

$$ p_2 = n_{21}/n_2 $$

p1,p2비율의 차의 standard error (s.e) 는 아래와 같다.  

$$ \sigma = \sqrt{\frac{p_1(1-p_1)}{n_1} + \frac{p_2(1-p_2)}{n_2}} $$

따라서 95% 신뢰구간은 1.96*±s.e 이다. 

오즈비의 신뢰구간

오즈비의 신뢰구간을 구하기 위해서는 로그 오즈비를 통해서 구하는 것이 좋다. 

로그 오즈비의 standard errer (s.e)는 아래와 같다. 

$$ \sigma = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}} }  $$ 

따라서 로그 오즈비의 95% 신뢰구간은 1.96*±s.e이 이며, 이를 오즈비의 신뢰구간으로 변환하기 위해서는 exponential 을 취해주면 된다. 따라서 오즈비의 95% 신뢰구간은 exp(1.96*±s.e) 이다.

어떤 컬럼의 값이 아래와 같은 문자열로 저장되어있을 때

["2021_12","2022_3","2022_1","2022_12","2023_4"....] 

해당 문자열 컬럼을 벡터컬럼으로 바꾸고 해당값을 unnest 하는 예시 

하나의 컬럼 값이 벡터형테인 경우 nested 라고 하고, 이를 row 로 변경하는 것을 unnest 라고 한다. 

# 문자열 parsing하여 year와 month로 분리하고 각 row로 만들기

df$dates <- lapply(df$month_ids, function(x) {
  unlist(fromJSON(x, simplifyVector = TRUE))
})
df<- df%>% 
  mutate(month_id = map(dates, str_split, pattern = ",")) %>%
  unnest(month_id)

X-> Y 의 인과적 관계 파악을 위해, 간단한게 심플 회귀 분석을 진행할 수 있다. 만약 X 와 correlation 이 있고, Y 의 determinants 인 Z 라고 하는 변수가 보정되지 않는다면, omitted variable bias 가 발생한다. 이러한 상황에서 omitted variable bias 의 방향은 다음과 같이 알 수 있다. 

1) Z->X 에 영향을 주는 방향 

2) Z->Y 에 영향을 주는 방향

1) 2) 를 곱하면 이것이 bias 의 방향이 된다. 

예를 들어, 소득(X)이 의료비 지출(Y)에 주는 영향을 파악하려고 한다. 이 때, 건강 상태(이를 개인이 갖고 있는 질병의 갯수라고 하자) 를 보정하지 않으면, omitted variable bias 가 발생하게 된다. 

질병의 개수는 소득에 negative effect 이다. 질병의 개수가 증가할 수록 소득은 감소한다.

질병의 개수는 의료비 지출에 positive effect 이다. 질병의 개수가 증가할 수록 의료비 지출은 증가한다. 

1) 2) 를 곱하면 negative 가 되기 때문에 bias 의 방향은 negative 가 된다. 따라서 건강 상태를 변수로 포함하지 않고 소득과 의료비 지출의 관계를 파악하여 나온 회귀 계수는 underestimate 이 되었다고 볼 수 있다. 만약, 동일한 건강상태에 있는 사람들만을 대상으로 소득과 의료비 지출의 연관성은 더욱 강하게 측정될 것이다.