교차표에서 효과를 추정하는 방법은 아래 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) 라고 부르기도 합니다. 

LD score regression 

LD score regression 은 Genome-wide association study(GWAS) 에서 특정 trait 의 polegenicity 를 추정하기 위해 사용하는 방법이다. LD score regression 은 GWAS summary statatistics 를 기반으로 SNP-heritability 추정, SNP-heritability 기반 genetic correlation 의 추정 등 다양한 measure 들을 계산하는데에 활용되고 있다. 본 문서에서는 LD score regression 의 등장 배경과 의미에 대해 알아보고자 한다. 

polygenic trait 에 대해 GWAS 를 수행한 후에, 각 SNP 들의 p-value 의 분포를 살펴보면 null distribution 과 비교하여 값들이 낮게 나타나는 것을 확인할 수 있다. 이렇게 높게 나타나는 검정 통계량은 형질이 polygenic 함을 의미할 수도 있지만, confounding bias 나 population stratification 가 영향을 주었을 수도 있다. polygenicity 로부터 위와 같은 bias 를 분리해내는 방법이 LD score regression 이며, LD score regression 은 이 과정에서 Linkage Disequillibrium (LD) 와 검정 통계량 (test statistics) 의 관계를 이용한다. 

"Both polygenicity (i.e. many small genetic effects) and confounding biases, such as cryptic relatedness and population stratification, can yield inflated distributions of test statistics in genome-wide association studies (GWAS). "

LD score regression 의 아이디어

어떤 SNP j 에 대해서 이 LD 관계에 있는 SNP 들이 많을 수록, polygenic 한 trait 에 대해서는 test statistics 이 높게 나올 가능성이 높다.  LD 는 유전자 변이간의 연관성을 의미한다. 만약, 어떤 SNP 이 LD 관계에 있는 SNP 이 많다고 하면, causal variant 와 LD 관계일 가능성이 높고, causal variant 와 LD 관계라면, test statistics 가 높게 나온다. 따라서 LD 관계에 있는 SNP 이 많으면, test statistics 가 높게 나올 가능성이 높다. LD score regression 은 이렇게 LD 관계에 있는 SNP 이 많은 SNP 일 수록 test stat 이 높게 나올 가능성이 높다라는 관계를 이용하는 방법이다. 그리고, 이러한 경향성이 강한 trait 일 수록 polygenicity 가 강하다고 말할 수 있다. (즉, LD score 와 test stat 의 연관성이 강할 수록 polygenic effect 로 phenotype 을 설명할 수 있는 비중이 높다.)

GWAS test statistics

LD score regression 에서는 chi-square value 를 regression 의 종속변수로 선정한다. chi-square value 는 무엇일까? 일반적으로 GWAS 결과로 effect size (beta) 와 standard deviation (sd) 값이 나오게 된다. 이 때, beta/sd 를 z-value 라고 한다. (이는 beta = 0 이라는 귀무가설 하에 구한 z-score 이다.) z-value 관측된 beta 값이 0으로부터 몇 standard deviation 떨어져 있는지를 의미한다. 이 때 chi-square value 는 z value 의 제곱으로 계산된다. z-value 는 표준정규분포를 따르며, 표준 정규 분포의 확률변수 z 의 제곱은 자유도가 1인 chi-square 분포를 따르기 때문이다. 

LD 와 LD score 는 무엇일까?

먼저 LD 와 LD score 를 계산하는 방법을 간단히 알아보자. 일반적으로 두개의 SNP A,B 의 LD 와 관련된 지표 D 와 r^2 은 아래와 같이 계산된다. D 값이 높을 수록 A,B 변이는 함께 나타날 가능성이 높음을 의미 한다. 만약, P(A)=0.3, P(B)=0.4, P(A,B) = 0.15 라고 하면, D는 0.03 으로 계산되며, A,B 는 LD 관계가 아닐 것으로 판단된다. 

$$ D_{AB} = P(A \cap B) - P(A) P(B) $$ 

일반적으로 많이 사용되는 지표인 r^2 은 아래와 같이 계산된다. 

$$ r_{AB}^2 = \frac{D^2}{P(A)(1-P(A))P(B)(1-P(B))}  $$ 

특정 SNP j 에 대한 LD score 는 아래와 같이 계산된다. LD score 는 각각의 SNP 에 대해 다른 모든 SNP 들과의 LD 값 (r^2) 들을 더한 값으로 볼 수 있다. 

$$ l_j = 1 + \sum_{k \neq j}r_{jk}^2 $$ 

LD score 와 test statistics (chi square value) 의 관계를 아래와 같이 시각화해볼 수 있다. 아래 차트는 LD score 의 bin 과 평균 chi-square value 의 관계를 보여준다. 직선은 아래 점들을 대상으로 단순 선형 회귀 분석을 한 결과를 표현한다 (아래 차트에서 각 점들에 해당 하는 SNP 의 갯수에 가중치를 두어 regression 을 돌리면 결국 전체 snp 과 chi-square 를 대상으로 regression 을 돌린 것과 같은 값이 나오게 될 것이다). 

이를 LD score regression 이라고 하며, 이 선형 회귀 분석에서 기울기는 polygenicity 를 반영하고, 절편은 bias 를 반영한다. LD 와 chi-square 의 연관성 (기울기) 이 polygenicity 이며, 전반적으로 chi-square 가 inflation 이 된 정도 (절편) 가 bias 라는 것이다. 이러한 방법을 통해 polygenicity 와 bias 를 분해할 수 있게 된다. 

또한, LD score regression 에서 기울기는 heritability 를 반영한다. 아래 그림과 같이 heritability 가 높을 수록 LD score regression 의 기울기가 커지게 된다. 

또한, 기울기는 샘플 사이즈와 사용한 SNP의 전체 개수에도 영향을 받는다. 샘플 사이즈 N 이 커질 수록 chi-square 값이 커지고, 사용한 SNP 의 개수가 많아질 수록, LD score 의 값이 기본적으로 높아진다. 이를 고려하여, 특정 SNP j 의 test stat 을 설명하는 아래와 같은 regression 모델을 고려할 수 있다.

$$ E[\chi^2 | l_j] = Nh^2/M l_j + Na + 1 $$

기울기는 heritability 와 N, M 으로 분해하여 나타낸다. 또한 절편은 Na + 1 로 표현되는데, 이 때, a 가 population structure 또는 confounding bias 와 같은 요인으로 인해 test stat 이 inflation 된 정도를 의미한다. 1이 더해진 이유는, LD score 가 0 인 SNP (즉, 그 어떤 SNP 과도 LD 관계에 있지 않은 SNP) 의 경우, chi-square value 는 causal variant 가 아닌 이상 기댓값은 1일 것이다 (자유도가 1인 chi-square distribution 의 평균값). 따라서, 절편은 1에 가까울 것이며, 1에서 벗어난 만큼을 bias 로 판단하겠다는 의미를 가진다. 

Bivariate LD score regression

두개의 trait X, Y 에 대한 test statistics 를 이용해 LD score regression 을 하는 것을 Bivariate LD score regression 이라고 한다. 구체적으로, 두개의 trait X, Y 에 대한 각각의 z-value 의 곱에 대하여 LD score regression 을 한다. 앞선 LD score regression 에 대해서는 chi-square value 를 사용했는데, chi-square value 는 z-value 의 곱이다. 즉, z^2 대신에 z_x * z_y 를 넣어서 regression 을 한다는 것이다.

참고) X_n 이 표준정규분포로부터 추출된 random variable 일 때, X_n 의 제곱의 합은 자유도가 n인 chi-square 분포를 따른다. 

위에서는 자유도가 1인 chi-square value 이기 때문에 z^2 = chi-square 가 된다. 

$$ X^2_1 + X^2_2 + ... X^2_n \sim \chi^2(n) $$ 

이 때의 기울기는 무슨 의미를 가질까? 만약 두개의 z-value 간에 아무런 연관성이 없다면, 1을 중심으로 퍼져있는 분포를 나타내게 된다. (z^2 의 분포는 자유도가 1인 카이제곱 분포를 따르며, 자유도가 1인 카이제곱 분포의 기댓값은 1이기 때문) 하지만, z-value 간에 연관성이 있다면, 기울기를 갖게 되며, 이 때의 기울기는 두 trait 간의 유전적 연관성을 의미한다. 두 trait 간에 유전적 연관성을 나타내는 지표로 co-heritability 가 있다. 기울기는 co-heritability 를 반영한다고 볼 수 있다. 이를 모델링하면 아래와 같이 표현할 수 있다. 

$$ E[z_{xj}z_{yj}] = \frac{\sqrt{N_xN_y}h_{xy}^2}{M}l_j + \frac{\rho N_s}{\sqrt{N_xN_y}} $$ 

왼쪽 그림에서 검은색 선은 유전적 연관성이 없는 trait 에 대한 기울기를 보여주며, 오른쪽 그림에서의 검은색 선은 연관성이 있는 trait 에 대한 기울기를 보여준다. 기울기가 가파를 수록 두 trait 간에 유전적 연관성이 높다고 볼 수 있다. 

참고자료

- LD Score Regression Distinguishes Confounding from Polygenicity in Genome-Wide Association Studies (2015)

- https://cnsgenomics.com/data/teaching/GNGWS22/module4/Lecture11_from_pgc_stat_bulik_2015.pdf

- An atlas of genetic correlations across human diseases and traits (2015)

- https://annahutch.github.io/PhD/LD-score-regression.html

genetic score, heritability, co-heritability, genetic correlation 관련 개념 정리

phenotype Y 의 분산은 genetics 로 설명되는 분산과 environment 로 설명되는 분산으로 나누어진다.  유전율 (heritability) 는 phenotype(또는 trait) 의 분산에서 genetics 를 통해 설명되는 분산을 의미한다. heritability 는 0~1 사이의 값을 가진다. 

$$ Var(Y) = Var(G) + Var(E) $$

$$ h^2 = Var(G)/Var(Y) $$ 

genetics 로 설명되는 분산이란 무엇일까? phenotype 을 예측하기 위해 유전정보를 이용해 어떠한 score 를 만들고, 이를 genetic score 라고 하자. genetic score 는 phenotype 에 대한 예측 값이며, 이 값이 높을 수록 phenotype 의 값이 높을 가능성이 높음을 의미한다 (만약 질병과 같은 binary trait 인 경우, 질병의 걸릴 가능성이 높음을 의미한다.)

genetic score 는 유전체 정보를 이용해 구한 phenotype Y 에 대한 예측값이다. 따라서 아래와 같이 쓸 수 있는데 heritability 의 식이 결정계수의 식과 닮아 있음을 알 수 있다. 결정계수는 전체 분산중 어떠한 모델을 통하 예측값의 분산이 차지하는 비율이며, 이것이 곧, 모델을 통해 설명되는 분산을 의미한다. 

$$ h^2 = Var(\hat{Y}) / Var(Y) = r^2 $$ 

genetics 로 설명되는 분산은 genetic score 의 분산으로 정의할 수 있다. 만약 genetic score 를 구할 때, additive genetic effect 만 고려하여, additive genetic score 를 구해 heritability 를 구한 것을 narrow-sense heritability 라고 한다.  

$$ h_n^2 = Var(AG)/Var(Y) $$

만약 , Y 가 standardization 이 되어 있다고하면, Y의 평균은 0이고, Y의 분산은 1이다. 그러면, 간단히, additive genetic score 의 분산이 바로 narrow sense heritability 가 된다. 

"If the traits are standardized (that is, phenotypic variance = 1) and the genetic values consider only the additive genetic effects, then the genetic variances are narrow-sense heritabilities."

$$ h_n^2 = Var(AG) $$

두 가지 trait 의 유전적인 연관성을 정량적으로 표현하는 지표로 coheritability 라는 개념이 있다.

"Co-heritability is an important concept that characterizes the genetic associations within pairs of quantitative traits."

co-heritability 는 아래와 같이 정의되며, -1~1 사이의 값을 가진다. 

$$ h_{x,y} = \frac{Cov(g_x,g_y)}{Var(X)Var(Y)} $$

이 식의 의미를 살펴보면 분자의 covariance 에 Cov(X,Y) 가 오게 된다면, pearson 상관계수와 같음을 알 수 있다. 이 식은 Cov(X,Y) 대신에 X,Y 에 대한 genetic score 를 대입시킴으로써, 두 trait 의 유전적 상관성을 표현했다고 볼 수 있다. 여기서도 마찬가지로 trait X,Y 를 평균이 0이고 분산이 1인 표준화된 trait 을 사용했다면, Var(X) = Var(Y) = 1 이기 때문에 아래와 같다. 

$$ h_{x,y} = Cov(g_x,g_y) $$

두 가지 trait 의 유전적인 연관성을 정량적으로 표현하는 지표로 genetic correlation 이라는 개념도 있다. genetic correlation 은 아래와 같이 정의된다. 

"The genetic correlation is a quantitative genetic parameter that describes the genetic relationship between two traits"

$$ \rho_g = \frac{Cov(g_x,g_y)}{\sqrt{Var(g_x)Var(g_y)}} $$ 

위 식은 pearson 상관계수의 식과 같으며, genetic correlation 의 통계적인 의미는 X,Y 의 genetic score 상관성 (pearson 상관계수) 라고도 할 수 있다. 만약 두가지 trait, 예를 들어 키와 발가락 길이의 유전적 연관성이 높다라고 한다면, 유전자를 통해 예측한 키 (키에 대한 genetic score) 와, 예측된 발가락 길이 (발가락 길이의 genetic score) 의 연관성이 높을 것이다. 이를 수치화한 것이 genetic correlation 이라고 볼 수 있다. genetic correlation 도 마찬가지로 -1~1사이의 값을 가진다. 

genetic correlation 과 co-heritability 모두, 두가지 trait 의 유전적 연관성을 표현한다. 둘의 차이점은 무엇일까?  trait X,Y 가 표준화 되어있다고 하면 genetic correlation 은 아래와 같이 정의된다. 아래 식을 보면, genetic correlation 은 co-heritability 가 X,Y 각각의 trait 의 heritability 로 보정된 식임을 알 수 있다. 

$$ \rho_g = \frac{h_{x,y}}{\sqrt{h^2_x h^2_y}} $$ 

따라서, 두 trait 의 heritabilty 값이 작더라도, genetic correlation 은 높을 수 있다. 예를 들어, 발가락 길이와 키의 heritability 가 10% 라고 하자 (실제로는 더 높을 것이나 예시임). 즉, 전체 분산에서 genetic score 의 분산이 차지하는 부분이 10% 이다. 하지만, 두개의 genetic score 의 연관성이 높다라고 하면, genetic correlation 은 높게 추정될 수 있다. 따라서, genetic correlation 을 해석할 때, trait 을 genetics 가 설명하는 비중 (heritability) 도 함께 고려해야할 필요가 있다. 

참고자료

- Genetic correlations of polygenic disease traits: from theory to practice, Nature review genetics, 2020

- Optimal Estimation of Co-heritability in High-dimensional Linear Models

- Statistical methods for SNP heritability estimation and partition: A review

List 를 이용한 방법

# dictionary 생성
dict <- list(name = "John", age = 30, city = "New York")

# dictionary 사용
dict$name
# [1] "John"

dict$age
# [1] 30

dict$city
# [1] "New York"

vector 를 활용한 방법

ㄴ setNames 함수를 활용

# dictionary 생성
dict <- setNames(c("John", 30, "New York"), c("name", "age", "city"))

# dictionary 사용
dict$name
# [1] "John"

dict$age
# [1] 30

dict$city
# [1] "New York"

hash 함수를 활용한 방법

library(hash)
h <- hash() 
h[['a']] <- 'a'
h[['b']] <- 'b'
h[['c']] <- 'c'
h[['d']] <- 'd'

h[['a']]

R 에서 특정 변수를 카테고리화 하고 싶을 때가 많다. 

다양한 방법이 있지만, 

아래 cut 함수를 사용하는 코드로 0~5, 6~10, 11~15, ... >100 으로 카테고리화가 가능하다. 

cat <- seq(0,100,5)
df$cat <- cut(df$x, breaks = c(cat, Inf), labels = cat)
df$cat <- factor(df$cat, levels=cat)

-> breaks 의 element 보다 labels 의 elements 의 갯수가 1개 적다. 

좀 더 일반적으로는 다음과 같다.

# 예시 데이터 생성
set.seed(123)
data <- data.frame(id = 1:10, value = rnorm(10, mean = 50, sd = 10))

# 카테고리화
data$cat <- cut(data$value, breaks = c(0, 25, 50, 75, 100), labels = c("low", "medium-low", "medium-high", "high"))

Hmisc 의 Lag 변수를 통해 timeseries 데이터의 lag 변수를 만들 수 있다. 

만약, 그룹별 Lag 변수를 만들고 싶으면 dplyr group_by 를 통해 만들 수 있다. 

library(Hmisc)
data <- data %>% group_by(gender, age) %>% mutate(lag = Lag(variable, 1))