선택편향

선택편향은 특정 그룹을 선택해서 분석했을 때, 다른 그룹 또는 전체를 대상으로 분석했을 때와 다른 결론이 나오는 것을 의미한다. 아래와 같이 왼쪽 그림에서는 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 이 되었다고 볼 수 있다. 만약, 동일한 건강상태에 있는 사람들만을 대상으로 소득과 의료비 지출의 연관성은 더욱 강하게 측정될 것이다. 

충분 통계량 (sufficient statistics) 

충분 통계량 (sufficient statistics) 란 계수를 설명할 수 있는 충분한 정보를 가지고 있는 통계량입니다. 예를 들어, 평균이 theta 이고 분산이 1인 정규분포에서 X1,X2,X3...Xn을 관찰하였고, X의 표본 평균값을 구해서 3이 나왔습니다. 그러면 X의 평균값 이외에 X1,X2,X3,...Xn 각각에 대한 데이터를 추가적으로 봄으로써 theta의 추정에 도움이 되는 다른 정보를 얻을 수 있을까요? 만약 추가적인 정보가 없다면, X의 표본 평균은 theta 를 설명할 수 있는 충분한 통계량이라고 볼 수 있습니다. 좀더 포멀하게는 아래와 같이 충분 통계량을 정의할 수 있습니다.

$$f(X|Y;\theta) \ne g(\theta)$$

위 식을 만족하면 Y는 theta 에 대한 충분통계량입니다. 위 식의 의미는 Y가 주어졌을 때, X의 pdf 가 theta 에 의존적이지 않다라는 의미입니다. 즉, Y가 theta 에 대한 모든 정보를 갖고 있다는 의미입니다.  

충분 통계량의 예시

예를 들어, 베르누이를 2회 시행해서 X1,X2 를 구했다고 합니다. X는 0 또는 1의 값을 갖습니다. 

이 때, X1 + X2 는 베르누이의 모수 p 에 대한 충분 통계량일까요? 

만약 X1+X2가 0이라면, X1,X2가 가질 수 있는 경우의 수는 (0,0) 입니다.

만약 X1+X2가 1이라면, (0,1), (1,0) 입니다.

만약 X1+X2가 2라면 (1,1) 입니다. 

즉, X1+X2 가 주어졌을 때, (X1,X2) 의 분포 (pmf) 는 theta 와 독립적입니다. 따라서 X1+X2 는 충분 통계량이라고 볼 수 있습니다. 이 의미는 만약에 X1+X2 가 1이라고 한다면, X1이 1이냐, X2가 1이냐는 적어도 계수를 추정함에 있어서 중요하지 않다는 것입니다. 위 예시에서는 베르누이 2회 시행에서 예시를 들었으나, N번의 베르누이 시행에서도 X1+...+Xn은 충분 통계량입니다. 만약 100번의 베르누이 시행에서 53번 성공했다라고 한다면, 몇 번째 시행에서 성공했는지는 계수 추정에 있어서 추가적인 정보를 주지는 않습니다. 

Maximum likelihood estimation

MLE는 관측한 데이터를 가장 잘 설명할 수 있는 계수(parameter) 를 찾는 것이 목적이라고 할 수 있습니다. 데이터 분석을 하다보면 데이터를 관측하고, 이를 통해 모델의 계수를 추정해야할 때가 있습니다. 우선 likelihood 를 정의해보면, 특정 계수에서 데이터를 관찰할 가능성을 의미합니다. 

예를 들어, 계수가 0.4인 베르누이분포를 100번 시행해서 45번의 성공이 발생했다고 합시다. likelihood 란 계수의 함수이며, 이 경우에서는 아래와 같이 나타낼 수 있습니다. 

$$L(\theta) = {100 \choose 45}\theta^{45}(1-\theta)^{55}$$

이 그래프는 계수가 0.1일 때 데이터를 관찰할 확률, 0.2 일때 관찰할 확률 등등.. 그 값을 모두 구해 계수와 likelihood 의 관계를 나타내는 그래프라고 볼 수 있습니다. 이 때, likelihood 는 계수가 0.45 일 때 최댓값을 갖습니다. 따라서 MLE 는 0.45 입니다. 이 경우 MLE 와 실제 계수에 차이가 있습니다. 하지만, 관측하는 데이터의 갯수가 증가할 수록, MLE 는 실제 계수인 0.4로 수렴해갈 것입니다 (하지만 모든 분포에서 실제 계수에 수렴하게 되는 것은 아닙니다).  

Probability 와 Likelihood 는 다른 개념입니다. Likelihood 는 데이터가 주어져있고, 이를 통해 계수를 추정하기 위한 것이라고 볼 수 있습니다. Probability 란, 어떤 모델의 계수가 주어졌을 때, 특정 outcome 을 예측하기 위한 것입니다. 

Maximum likelihood 를 만드는 지점을 찾기

그러면 likelihood function 의 값을 최대화하는 계수값을 어떻게 구할 수 있을까요? 일반적으로는 likelihood function 을 미분해서 값이 0인 지점을 찾으면 됩니다. 또는 계산상의 이점을 위해 log likelihood 를 미분해서 값이 0이 되는 지점을 찾습니다. likelihood 를 미분해서 0이 되는 지점이나, log likelihood 를 미분해서 0이 되는 지점이나 값은 같기 때문입니다. 

예를 들어, 앞서 살펴본 베르누이 분포에서 MLE 를 통해 계수 추정을 해봅시다. 베르누이 분포의 pmf 는 아래와 같이 정의 됩니다. 

$$\theta^{x}(1-\theta)^{1-x}$$

베르누이 분포의 likelihood 는 아래와 같이 쓸 수 있습니다.

$$L(\theta) =  \theta^{\sum_{i=1}^{n}{x_i}} 1-\theta^{n- \sum_{i=1}^{n}{x_i} }$$   

그러면 log likelihood는 아래와 같습니다. 

$$l(\theta) = \log \theta \sum_{i=1}^{n}{x_i} + \log (1-\theta)(n-\sum_{i=1}^{n}{x_i} )$$

log likelihood 를 미분하면 아래와 같습니다. (이를 score 라고 부르기도 합니다.)

$$l'(\theta) = \frac{1}{\theta} \sum_{i=1}^{n}{x_i} - \frac{1}{1-\theta}( n-\sum_{i=1}^{n}{x_i} )$$

위 score를 0으로 만드는 theta (계수) 값은 표본 평균입니다. (이와 같이 많은 경우, MLE 는 표본 평균인 경우가 많습니다.)

$$\hat{\theta} = \bar{X}$$  

또한 이 경우 표본 평균의 기댓값이 모평균이기 때문에 unbiased estimation 이라고 볼 수 있습니다.

MLE 는 항상 unbiased 는 아니다. 

하지만 MLE 의 경우 모든 경우에 unbiased 인 것은 아닙니다. 예를 들어, uniform distribution(0,theta) 의 계수를 MLE 로 추정값은 관측한 데이터중 최댓값이 됩니다(관련된 자세한 설명은 생략하겠습니다). 그리고 이값의 기댓값은 아래와 같이 n에 의존하는 값이 되기 때문에 biased estimation 이라고 할 수 있습니다.(참고링크) 

MLE 의 분산과 Efficiency

그럼에도 불구하고 MLE 는 관측 데이터가 많으면 많을 수록 가능한 unbiased estimator 집합들이 가질 수 있는 분산의 최솟값으로 수렴한다는 큰 장점을 가지고 있습니다. 즉, MLE 는 이론적으로 나올 수 있는 가장 작은 분산을 가진다는 의미이며, 이를 "Asymptotically efficient 하다" 라고 부르기도 합니다. 

이 때, 이론적으로 나올 수 있는 분산의 최솟값을 Rao-Cramer Lower Bound 라고 부릅니다. 어떤 unbiased estimator Y가 있을 때, Y의 분산의 최솟값은 아래과 같습니다. 

$$var(Y) \ge \frac{1}{nI_1{\theta}}$$

여기서 I는 fisher information 이고 아래와 같습니다. I_1 은 데이터 하나로 구한 fisher information 입니다. 

$$I(\theta) = E(-l''(\theta))$$

$$I(\theta) = nI_1(\theta)$$

베르누이 분포에서 fisher information 과 rao-cramer lower bound 를 구해봅시다. 데이터 1개에서 log likelihood 를 구하면 fisher information 을 구하고 이를 통해 rao-cramer lower bound 를 구할 수 있습니다.

$$l(\theta) = xlog\theta + (1-x)log(1-\theta)$$

$$l'(\theta) = \frac{x}{\theta} - \frac{1-x}{1-\theta}$$

$$l''(\theta) = \frac{x}{\theta^2}-\frac{1-x}{(1-\theta)^2}$$

$$I_1(\theta) = \frac{1}{\theta(1-\theta)}$$

rao-cramer lower bound 는 아래와 같습니다. 

$$var(\hat{\theta}) \ge \frac{\theta(1-\theta)}{n}$$  

실제 계수 추정값의 분산을 구해봅시다. 

$$var(\hat{\theta}) = var(\bar{X}) = var(\frac{  \sum_{i=1}^{n}{x_i}  }{n}) = \frac{1}{n^2} \sum_{i=1}^{n}{var(x_i)} = \frac{\theta(1-\theta)}{n}$$

이 값은 rao-cramer lower bound 의 값과 같습니다. 따라서 베르누이 분포에서 MLE 는 최소 분산을 가지며, efficient estimator 라고 할 수 있습니다. 또한 unbiased 이기 때문에, minimum variance unbiased  estimator (MVUE) 라고 부르기도 합니다.