Monday, February 16, 2015

A Primer on Linear Models (chapter 1)

Chapter1 Examples of the General Linear Model

1.0 Contents

1.1 Introduction

1.2 One-Sample Problem

1.3 Simple Linear Regression

1.4 Multiple Regression

1.5 One-Way ANOVA

1.6 First Discussion

1.7 Two-Way Nested Model

1.8 Two-Way Crossed Model

1.9 Analysis of Covariance

1.10 Autoregression

1.11 Discussion


1.1 Introduction

\[y=Xb + e\]
\(y\) : N x 1 vector of observed reponses
\(X\): N x p matrix of fixed constants
\(b\): p x 1 vector of fixed, unknown parameters
\(e\): N x 1 vector of unobserved errors with zero mean

위 식이 general linear model 의 기본식. 1 장에서는 general linear model의 예를 알아본다.

1.2 One-Sample Problem

\[Xb = 1(\mu)\]
\(y_i\) 는 randomly sampled from a population with unknown mean \(\mu\) and variance \(\sigma^2\)

1.3 Simple Linear Regression

원래 식 \[y_i = \beta_0 +\beta_1x_i + e_i\]
general linear model 로 나타낸 식 \[y = \left[\begin{array}{c} y_1\\ y_2 \\ ... \\ y_N \end{array}\right], \quad Xb = \left[\begin{array}{cc} 1 & x_1 \\ 1 & x_2 \\ ... & ... \\ 1 & x_{N-1} \\ 1 & x_N \\ \end{array}     \right]\left[\begin{array}{c} \beta_0 \\ \beta_1 \end{array}\right] , \quad e = \left[\begin{array}{c} e_1\\ e_2 \\ ... \\ e_N \end{array}\right]  \]

1.4 Multiple Regression

원래 식 \[y_i = \beta_0 +\beta_1x_{i1} +\beta_2x_{i2} +\beta_kx_{ik} + e_i, \quad i = 1,...,N\]
general linear model 로 나타낸 식 \[y = \left[\begin{array}{c} y_1\\ y_2 \\ ... \\ y_N \end{array}\right], \quad Xb = \left[\begin{array}{cccc} 1 & x_{11} & x_{12} & ... & x_{1k} \\ 1 &  x_{21} & x_{22} & ... & x_{2k} \\ 1 & ... & ... & ... & ... \\  1 & x_{N1} & x_{N2} & ... & x_{Nk} \\ \end{array} \right] \left[\begin{array}{c} \beta_0 \\ \beta_1 \\ \beta_2 \\ ... \\ \beta_k \end{array}\right] , \quad e = \left[\begin{array}{c} e_1\\ e_2 \\ ... \\ e_N \end{array}\right]  \]

1.5 One-Way ANOVA

원래 식 \[y_{ij} = \mu + \alpha_i + e_{ij}\]
general linear model 로 나타낸 식 \[y = \left[\begin{array}{c} y_{11}\\ y_{12} \\ ... \\ y_{1n_1} \\ y_{21} \\ y_{22} \\ ... \\ y_{2n_2} \\ y_{a1} \\ y_{an_a}  \end{array}\right], \quad Xb = \left[\begin{array}{ccccc} 1_{n_1} & 1_{n_1} & 0 & ... & 0 \\ 1_{n_2} & 0 & 1_{n_2} & ... & 0 \\ ... & ... & ... & ... & ... \\  1_{n_a} & 0 & 0 & ... & 1_{n_a} \\ \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_1 \\ \alpha_2 \\ ... \\ \alpha_a \end{array}\right] , \quad e = \left[\begin{array}{c} e_{11}\\ e_{12} \\ ... \\ e_{1n_1} \\ e_{21} \\ e_{22} \\ ... \\ e_{2n_2} \\ ... \\ e_{a1} \\ ... \\ e_{an_a} \\ e_N \end{array}\right]  \]

\(i = 1,...,a\)
\(j = 1...,n_i\)

1.6 First Discussion

추후 정리 필요

1.7 Two-Way Nested Model

원래 식\[y_{ijk} = \mu + \alpha_i + \beta_{ij}+  e_{ijk}\]

\(i = 1,...,a\)
\(j = 1,...,b_i\)
\(k = 1,...,n_{ij}\)

\(\alpha \) 가 두 가지 F, P, 곧 \(\alpha_i = F, P\) 이고 각 \(\alpha\) 에 대해 \(\beta\) 가 \(D, M\) 과 \(A,B,E\),  곧 (\(\beta_F = {D, M}, \beta_P={A, B, E} \)) 일때


general linear model 로 나타낸 식  \[y = \left[\begin{array}{c} y_{F,D,1}\\ ... \\ y_{F,D,n_{11}} \\ y_{F,M,1} \\ ... \\ y_{F,M,n_{12}} \\ y_{P,A,1} /\\ ... \\ y_{P,A,n_{21}} \\ y_{P,B,1} \\ ... \\ y_{P,B,n_{22}} \\ y_{P,E,1} \\ ... \\ y_{P,E,n_{23}}  \end{array}\right], \quad Xb = \left[\begin{array}{cccccccc} 1_{n_{11}} & 1_{n_{11}} & 0 & 1_{n_{11}} & 0 & 0 & 0 & 0 \\ 1_{n_{12}} & 1_{n_{12}} & 0 & 0 & 1_{n_{12}} & 0 & 0 \\  1_{n_{21}} & 0 & 1_{n_{21}} & 0 & 0 &  1_{n_{21}} & 0 & 0 \\ 1_{n_{22}} & 0 & 1_{n_{22}} & 0 &  0 & 0 & 1_{n_{22}} & 0 \\ 1_{n_{23}} & 0 & 1_{n_{23}} & 0 & 0 & 0 & 0 & 1_{n_{23}}  \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_F \\ \alpha_P \\ \beta_{F, D} \\ \beta_{F, M} \\ \beta_{P, A} \\ \beta_{P, B} \\ \beta_{P, E}  \end{array}\right]  \]

1.8 Two-Way Crossed Model

원래 식 \[y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} + e_{ijk}\]
\(i = 1,...,a\)
\(j = 1,...,b\)
\(k = 1,...,n_{ij}\)

위 식에서 \(\beta\) 가 없다면 이는 two-way nested model 이 되고 \(\alpha\)와 \(\beta\) 사이에 interaction이 없다면 \(\gamma\) 는 0이 된다.

\(\alpha\)와 \(\beta\) 사이에 interaction이 없고 \(n_{ij}\) = 1, (곧 no replication), \(i = 1,...,a\), \(j = 1,...,b\) 이며 balanced case 라고 한다면 맨 처음 식은
\[y_{ij} = \mu + \alpha_i + \beta_j +  e_{ij}\]

general linear model 로 나타낸 식은 \[y = \left[\begin{array}{c} y_{11}\\ y_{12} \\ ... \\ y_{1b} \\ y_{21} \\ y_{22} \\ ... \\ y_{2b} \\ ... \\ y_{a1} \\ y_{ab}  \end{array}\right], \quad Xb = \left[\begin{array}{cccccc} 1_b & 1_b & 0 & ... & 0 & I_b \\ 1_b & 0 & 1_b & ... & 0 & I_b \\ 1_b & 0 & 0 & ... & 0 & I_b \\ ... & ... & ... & ... & ... & ... \\  1_b & 0 & 0 & ... & 1_b & I_b \\  \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_1 \\ \alpha_2 \\ ... \\ \alpha_a \\ \beta_1 \\ ... \\ \beta_b \end{array}\right] \]

1.9 Analysis of Covariance

원래 식 \[y_{ij}\ = \mu + \alpha_i + \beta x_{ij} + e_{ij}\]
\(i=1,...,a\)
\(j=1,...,n_i\)

general linear model 로 나타낸 식은 \[y = \left[\begin{array}{c} y_{11}\\ y_{12} \\ ... \\ y_{1n_1} \\ y_{21} \\ y_{22} \\ ... \\ y_{2n_2} \\ ... \\ y_{a1} \\ ... \\ y_{an_a}  \end{array}\right], \quad Xb = \left[\begin{array}{cccccc} 1_{n_1} & 1_{n_1} & 0 & ... & 0 & x_1 \\ 1_{n_2} & 0 & 1_{n_2} & ... & 0 & x_2 \\ 1_{n_3} & 0 & 0 & ... & 0 & x_3 \\ ... & ... & ... & ... & ... & ... \\  1_{n_a} & 0 & 0 & ... & 1_{n_a} & x_a \\  \end{array} \right] \left[\begin{array}{c} \mu \\ \alpha_1 \\ \alpha_2 \\ ... \\ \alpha_a \\ \beta \end{array}\right] \]

만약 treatment \(\alpha\) 에 따라 covariate 의 slope 이 다르다면 식은 \[y_{ij}\ = \mu + \alpha_i + \beta_i x_{ij} + e_{ij}\] 가 된다.

* analysis of covariance는 식이 two-way nested model과 비슷하나 error가 \(e_{ijk}\) 가 아닌 \(e_{ij}\) 이다. 이는 오히려 하나의 factor 에 의해 샘플이 구분되는  one-way ANOVA 와 비슷한데 여기서 같은 그룹내에 각 샘플이 갖는 개별적인 속성(covariate) 의 변수를 넣은 것이라고 이해하면 될 것이다.

1.10 Autoregression

simple linear trend model for data collected over time \[y_i = \beta_0 + \beta_1t + e_t, \quad t = 1,2,...,N\]
\(i\)대신 \(t\)를 사용한 이유는 time ordering을 강조하기 위함.
추후 정리 필요


1.11 Discussion

위에서와 같이 다양한 통계 모델은 general linear model의 special case이다.  \(e_i\) 가 uncorrelated random variables with zero mean and constant variance \(\sigma^2\) 인 경우가 있고 다른 경우는 correlated 한 경우가 있다.
추후 정리 필요

Wednesday, February 11, 2015

PostgreSQL homebrew로 9.3 에서 9.4 로 upgrade 하기

Mac home-brew update & upgrade를 하면 PostgreSQL 이 버젼업 되는데 아래와 같은 에러가 난다.

LOG:  skipping missing configuration file "/usr/local/var/postgres/postgresql.auto.conf"
FATAL:  database files are incompatible with server
DETAIL:  The data directory was initialized by PostgreSQL version 9.3, which is not compatible with this version 9.4.0.
이때 아래 사이트를 따라하면 해결됨.

http://rcmachado.github.io/postgres/postgis/2014/12/20/upgrading-postgresql-using-homebrew.html

Wednesday, February 4, 2015

Linear Algebra

rank가 \(r\) 인 \( A_{m\mbox{x}n}\)  matrix 의 4가지 vector space 의 property 정리.

  • C(A) : column space of A. Linear combination of columns of A
  • R(A) : row space of A. Linear combination of rows of A
  • N(A) : null space of A. \(N(A) = \{y: Ay = 0 \} \)
  • N(A') : null space of A' (다르게 표현하면 \(\textrm{C}(A)\) 의 orthogonal complement). \(\textrm{C}(A) = S\) 일때 \(S^{\perp} = \textrm{N}(A')\) 


vector space C(A) R(A) N(A) N(A')
vector dimension \(R^m\) \(R^n\) \(R^n\) \(R^m\)
space dimension \(r\) \(r\) \(n-r\) \(m-r\)
projection matrix \(AA^-\) \(A'A'^-\) \(I - A^-A\) \(I - A'A'^-\)
- - - - -
- - - - -

* \(A'\) : transpose of A
* \(A^-\) : generalized inverse (GI) of a matrix A
* projection matrix : a square matrix \(P\) is a projection matrix that projects onto the vector space \(S\), \(S \subseteq R^n\)

  1. \(P\) is idempotene (=> \(PP = P\))
  2. \(Px \in S \; \forall x \in R^n \)
  3. \( Pz = z \; \forall z \in S \)



01PreliminaryLinearAlgebra

Fact V1
  • \(x_1,x_2,...x_n\) 이 LD (linearly dependent)  \( \Leftrightarrow \) \( x_j \) 는 \( x_1,x_2,..x_{j-1}\) 의 LC (linear combination)
Fact V2
  • vector space \(S\) 의 basis의 갯수가 n개이고 \(S\) 의 LI (linearly independent) 한  k개의 vector가 있다면 k \(\leq\) n
Fact V3
  • \({a_1,...,a_n}\) 과 \({b_1,...,b_k}\) 이 각각 vector space S의 basis라고 한다면 \(n=k\).
Fact V4
  • \(a_1,...,a_n\) 이 vector space S, S \(\in dim(\mbox{S})=n\) 의 LI vector라면 이는 곧 S의 basis 가 된다.
Fact V5
  • \(a_1,...,a_k\) 이 vector space S, S \(\in dim(\mbox{S})=n\) 의 LI vector라면 \(a_1,...,a_k\) 를 포함한 basis가 존재한다.
Fact V6
  • \(a_1,...,a_k\)이 \(R^n\)의 LI vectors 라면 \(a_1,...,a_n\) 이 orthornormal LI 가 되는  \(a_{k+1},...,a_n\) 이 존재한다.
Result A.1
  • \(rank(AB) \le min(rank(A),rank(B))\)
Result A.2
  • \(A=BC\) 이면 \(\textrm{C}(A) \subseteq \textrm{C}(A)\)
  • \(\textrm{C}(A)\subseteq \textrm{C}(B)  \) 이면 \(A=BC\) 를 만족하는 \(C\)가 존재한다.
Result A.3
  • \(rank(A_{mxn})=n \Longleftrightarrow \textrm{N}(A)={0}\)
TheoremA.1
  • \(rank(A_{mxn})=r\) 이면 \(dim(\textrm{N}(A))=n-r\) 이고 다르게 표현하면 \(dim(\textrm{N}(A_{mxn}))+dim(\textrm{C}(A_{mxn}))=n\)

02OrthogonalComplements

Result A.4
  • vector space \(S \subseteq R^n \) 이 있을 때 모든 \(x \in R^n\) 는 \(s\in S, t \in S^{\perp}\) 에 의해 \(x = s+ t\)  로 나타내어진다. 그리고 이 \(s,t\) 의 decomposition은  unique하다.
Result A.5
  • \(A_{mxn}\) 이 있을 때 \(\textrm{C}(A)^{\perp} = \textrm{N}(A')\)
Result A.6
  • \(S_1,S_2\) 는 \(R^n\) 안의 vector space 이고 \(S_1 \subset S_2\) 이면 \(S_2^{\perp} \subseteq S_1^{\perp}\) 이다.

04ProjectionMatrices

Result P.1
  • \(P\) 가 idempotent matrix (=> \(PP=P\)) 이고 vector space \(S\) 의 projection matrix 라면 \(S = \textrm{C}(P)\)   
Result A.14 
  • \(AA^-\) 는 \(\textrm{C}(A)\) 의 projection matrix
Result A.15
  • \(I - A^-A\) 는 \(\textrm{N}(A)\) 의 projection matrix 
Result A.16
  • \(P\) 가 symmetric 이고 idempotent라면 \(P\)는 \(\textrm{C}(P)\) 의 unique symmetric projection matrix 이다. 곧 \(P\) 가 자신의 column space 의 projection matrix가 된다.
Corollary A.4
  • \(P\) 가 symmetric idempotent matrix 일때 \(I - P\) 는 \(\textrm{C}(P)^{\perp} = \textrm{N}(P) = \textrm{N}(P') \) 의 sysmmetric projection matrix 이다.

05IntroGLM

asf

07PxNormalEquations

Theorem 2.1

  • \(P_x\) 는 \(\textrm{C}(X)\) 의 orthogonal projection 일 때, \(P_x = X(X'X)^-X'\) 이고 이는 아래와 같은 성질을 갖는다.
    1. \([X(X'X)^-X'][X(X'X)^-X'] = [X(X'X)^-X'] \)
    2. \( X(X'X)^-X'y \in \textrm{C}(X) \; \forall y \in R^n \)
    3. \( X(X'X)^-X'x = x \;  \forall x \in \textrm{C}(X)\) 
    4. \(X(X'X)^-X' = [X(X'X)^-X']' \)
    5. \(X(X'X)^-X'\) 는 \(X'X\)의 \( \forall GI\) 와 같다
Lemma 2.1
  • \(\textrm{N}(X'X) = \textrm{N}(X)\)
Result 2.1
  • \(\textrm{C}(X'X)=textrm{C}(X')\)
Result 2.3
  • \(Q(b)\)를 최소화 하는 \(\hat\beta\) 값이 \(\textrm{NE}(X'Xb=X'y)\) 의 해이다. 
Result 2.4
  • \(X'XA=X'XB \Longleftrightarrow XA=XB\)
Result 2.5
  • \((X'X)^-\) 가 \(X'X\)의 GI 이면 \((X'X)^-X'\) 는 \(X\)의 GI 이다. 곧 \(X(X'X)^-X'X = X\)
Corollary to result 2.5
  • \((X'X)^-\) 가 \(X'X\)의 GI 이면 \(X(X'X)^-\) 는 \(X'\)의 GI 이다. 곧 \(X'X(X'X)^-X' = X'\)
Result 2.6
  • \(I-P_x = I-X(X'X)^-X'\) 는 \(\textrm{C}m(X)^{\perp}=\textrm{N}(X')\) 
Theorem 2.2
  • \(\textrm{C}(W) \subseteq \textrm{C}(X)\) 일 때, \(P_x-P_w\)는 \(\textrm{C}((I-P_w)X)\) 의 orthorgonal projection 이다.
Normal Equations (NE)
  • \(X'X \hat\beta = X'y\) 
  • NE의 해는 Result 2.3 과 같이 \(Q(b)\)를 최소화 하는 \(\hat\beta\) 값이 되고 이는 \(X(X'X)^-X'y\) 이다.
  • solution of NE \( \Longleftrightarrow \hat\beta \Longleftrightarrow \mbox{b minimized of }Q(b) \Longleftrightarrow X(X'X)^-X'y\)
Corollary 2.1
  • NE 는 consistent 하다. 곧 항상 해가 있다.
Corollary 2.2
  • \(rank(X'X)=rank(X)\)
Corollary 2.3
  • NE의 그 어떤 해 \(\hat\beta\) 의 \(X\hat\beta\) 의 값은 동일하다.

08Reparameterization

\(y = W\alpha+\epsilon\) 과 \(y = X\beta+\epsilon\) 는 동등하다 혹은 reparameterization iff \(\textrm{C}(X)=textrm{C}(W)\)
Result 2.8
  • \(\textrm{C}(X)=\textrm{C}(W)\) 이면 \(P_x=P_x\).
Corollary 2.4
  • \(\textrm{C}(X)=\textrm{C}(W)\) 이면 \(\hat y = P_{x}y = P_wy \) 이고 \(\hat \epsilon = (I-P_x)y = (I-P_w)y = y- \hat y\)
Result 2.9
  • \(\textrm{C}(W)=\textrm{C}(X)\) 이고 \(\hat \alpha\) 가 NE \(W'W\alpha=W'y\) 의 해일 때, \(\hat \beta = T \hat \alpha\) 는 NE \(X'X\alpha=X'y\) 의 해이고  \(T \ni W= XT\)




Tuesday, January 27, 2015

Estimated Breeding Values (EBV)


EBV (Estimated Breeding Value) 란 상대적 유전값 (relative genetic value)으로 육종(breeding), 즉 인간이 원하는 방향으로 생물을 개량하는 데 객관적인 유전적 지표로 사용되어 진다.  아래 전체적인 EBV를 구하는 방법과 그 적용 예를 들어본다.

Contents
  • Introduction
  • Theoretical Background
  • Statistical Model
  • Computing and Using EBV's
  • How to use ebb's to make selection decisions
  • Multiple-trait model




Introduction

EBV 란 statistical prediction of the relative genetic merit, 곧 개체의 유전적인 유리함이 어떤지에 대한 통계적 예측값이다. 이 값은 육종에 사용된다. 원하는 형질의 EBV가 높은 개체만이 생식하게 나고 그렇지 않은 개체의 생식은 막음으로서 population 이 특정 형질에 특화되는 방향으로 흘러가게 하는 것이다.


Theoretical Background

표현형(phenotype) 은 고기의 지방함량, 사료대비 성장률등의 관측가능한 육체적 형질을 의미한다. 이 phenotype은 아래식과 같이 크게 유전적 요인과 비유전적요인으로 구성되어진다. 유전적 요인은 자손에게 유전이 되며 비유전적 요인은 환경적요인과 그밖에 그 개체에 존재하는 특성들을 의미한다. 비유전적 요인은 크게 fixed effect (sex, year of birth, management condition) 과 그밖에 측정하기 힘든 effect로 나뉜다.

\[ \mbox{Phenotypic observation = Environmental effects + Genetic effects + Residual effects} \]
이를 통계적 모델로 표현하면 아래와 같다.
\[ y_{ij} = \mu_i + g_i + e_{ij} \quad  (1)\]

\( y_{ij}\) = i번째 개체의 j 번째 관측치
\( \mu_i \) = i번째 개체의 fixed  environmental effect
\( g_i \) = i번째 개체의 유전형의 addictive(\(g_a\)), dominance(\(g_d\)), epistatic(\(g_e\)) genectic value의 총 합
\( e_{ij} \) = i 번째 개체의 j번째 관측에서의 random environmental effect의 총합

여기서 \(g_i\)를 (\(g_a\)), (\(g_d\)), (\(g_e\))로 대체하게 되면
\[\begin{eqnarray} y_{ij} & = & \mu_i + (g_a)_i + (g_d)_i + (g_e)_i + e_{ij}  \\  & = & \mu_i + (g_a)_i + e_{ij}^*    & (2) \end{eqnarray} \]
\(e_{ij}^* = (g_d)_i + (g_e)_i + e_{ij}\)

위 2번 식이 \(\mu_i\) 를 알고 있다는 가정하 EBV를 계산하는 기본식이다. 그러나  \(\mu_i\) 를 정확하게 알고 있는 것은 거의 불가능하기 때문에 같은 데이터를 가지고  \(\mu_i\) 와  \((g_a)_i\)를 동시에 구해야 한다. 이 두가지 값을 동시에 구하는 방법이 BLUP (Best Linear Unbiased Prediction) 이다.

BLUP은 동시에  fixed effect에 의해 발생하는 차이 조절하고 집단에서의 가계도에 있는 유전적 관계를 설명하는 과정에서 EBV를 구한다.
BLUP을 이용하기 위해서는 형질에 대한 heritability를 알고 있어야 한다.
BLUP을 표현하는 식으로 mixed model equation 이 있다.
\[ Y = X\beta + Zu + e \quad (3)\]
\(Y \) = vector of observed value
\(\beta \) = vector of fixed effect
\(u \) = vector of random genetic effect (2번식에서 \((q_a)i\) 와 같은 의미)
\(e \) = vector of residual error
\(X, Z \) = design matrix

3번식은 여러 단계를 거쳐 아래와 같은 4번식으로 유도가 된다.
\[ \left[\begin{array}{cc} X'X & X'Z \\ Z'X & Z'Z + A^{-1}\alpha \end{array}\right] \left[\begin{array}{cc} \beta \\ u \end{array}\right] = \left[\begin{array}{cc} X'Y \\ Z'Y  \end{array}\right] \quad (4) \]
\(A^{-1}\) = \(A\)의 inverse matrix,  A는 가계도를 바탕으로 개체간의 유전적 관계를 수학적으로 나타낸 matrix
\( \alpha \) = (1 - \(h^2\)) / \(h^2\), \(h\) 는 heritability
개체간의 유전적 관계에서 breeding value을 추측할 수 있는 BLUP의 핵심은 \(A^{-1}\) 에 있다.


Animal ID Sire ID Dam ID Sex Contemporary Group physical trait Test Score 
1 Unknown Unknown M
2 Unknown Unknown F
3 1 Unknown M
4 1 Unknown F
5 Unknown Unknown M
6 3 2 M 1 2
7 3 2 F 2 4
8 5 4 M 1 3
9 5 4 F 2 5

table 1

9마리의 개체가 있고 개체간의 가계도는 표1에서 와 같다고 하자. \(A\)는 개체간의 유전적 관계를 나타내는 matrix 로 9X9 diagonal matrix 이다.  \(A\)를 구하는 방법은 아래와 같다.

먼저 위의 표와 같이 가족 관계를 표를 만드는데 부모가 되는 개체가 자식이 되는 개체보다 앞이 되도록 표를 만든다.

for i to n,
    for j to i,

  1. 두 부모(s, d)를 모두 아는 경우 \[ \left\{\begin{array}{r1} if \: j < i,  a_{ij} = a_{ji} = 0.5(a_{js} + a_{jd}) \\ if  \: j = i, a_{ii} = 1 + 0.5a_{sd}  \end{array}\right. \]
  2. 한쪽 부모(s)만 아는 경우  \[ \left\{\begin{array}{r1}if \: j < i,  a_{ij} = a_{ji} = 0.5(a_{js} \\ if  \: j = i, a_{ii} = 1  \end{array}\right. \]
  3. 양쪽 부모를 다 모르는 경우 \[ \left\{\begin{array}{r1}if \: j < i,  a_{ij} = a_{ji} = 0 \\ if  \:   = i, a_{ii} = 1 \end{array}\right. \]


위 방법을 이용하게 되면 \(A\)는 아래와 같이 구해진다.
table2

위의 값이 클수록 두 개체간의 유전적 관계가 밀접한 것을 의미한다. 개체수가 많아질 경우 A row, column의 수는 기하급수적으로 많아지고 \(A^{-1} \) 의 계산은 더욱 힘들어 진다. 1970 중반에 \(A\) 를 구하지 않고 \(A^{-1} \) 를 바로 구할 수 있는 방법이 고안되면서 수십만 개체의 genetic relationship matrix 의 계산이 가능해졌다.


Statistical Model

BLUP이 fixed effect와 genetic effect를 동시에 고려하기 때문에 타당한 fixed effect 의 선택이 되어야만 개체간의 올바른 유전적 비교가 가능하다. fixed effect가 제대로 선택되지 않으면 fixed effect에 대한 효과가 residual error에 포함되어지고 그만큼 모델의 설명력은 떨어지게 된다.  fixed effect에는 sex, 그리고 physical interpretation을 수행한 평가자, 나이, 평가 시기 등이 있다.


Computing and Using EBV's

heritability(h) 를 20%, fixed effect 를 contemporary group이라는 가정하에 예제 표1의 EBV를 계산하면 아래와 같다.
\[ Y = X\beta + Zu + e \]

\[ \left[\begin{array}{c} y_6 \\ y_7 \\ y_8 \\ y_9  \end{array}\right] = \left[\begin{array}{c} 2 \\ 4 \\ 3 \\ 5  \end{array}\right] = \left[\begin{array}{cc} 1 & 0 \\ 0 & 1 \\ 1 & 0 \\  0 & 1 \end{array}\right]  \left[\begin{array}{c} b_1 \\ b_2   \end{array}\right] +  \left[\begin{array}{r1} 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0  \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0  \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0  \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1   \end{array}\right] \left[\begin{array}{c} a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \\ a_7 \\ a_8 \\ a_9  \end{array}\right]  + \left[\begin{array}{c} e_6 \\ e_7 \\ e_8 \\ e_9  \end{array}\right]\]
가 된다.
\[ \left[\begin{array}{cc} X'X & X'Z \\ Z'X & Z'Z + A^{-1}\alpha \end{array}\right] \left[\begin{array}{cc} \beta \\ u \end{array}\right] = \left[\begin{array}{cc} X'Y \\ Z'Y  \end{array}\right] \quad (4) \]
을 이용하기 위해서 \(A^{-1}\)  와 \( \alpha \)를 구하게 되면

 \(A^{-1}\) 은 아래와 같고
table 3

heritability = 0.2 로 가정하였기 때문에 \( \alpha \) = (1 - 0.2) / 0.2 = 4 가 된다.

4번식에 \(X, Z, A^{-1}, \alpha \) 값을 넣고 계산하게 되면 아래와 같이 fixed effect의 값과 개체의 EBV를 구하게 된다.              
\[ \left[\begin{array}{c} \beta \\ u \end{array}\right] = \left[\begin{array}{c} \mbox{Contemporary group 1 average} \\ \mbox{Contemporary group 2 average} \\ EBV 1 \\ EBV 2 \\ EBV 3 \\ EBV 4 \\ EBV 5 \\ EBV 6 \\ EBV 7 \\ EBV 8 \\ EBV 9 \end{array}\right] = \left[\begin{array}{c} 2.5 \\ 4.5 \\ 0 \\ -0.09 \\ -0.07 \\ 0.07 \\ 0.09 \\ -0.13 \\ -0.13 \\ 0.13 \\ 0.13 \end{array}\right]  \]

How to use EBV's to Make Selection Decisions

EBV 는 그저 germplasm(= collection of genetic resources) 일 뿐 EBV가 높다고 항상 좋은 것은 아니다. phenotype을 disease susceptibility 로 두었다면 EBV 이 최소값을 같는 개체를 선발하는 것이 옳을 것이다.
위의 예의 physical trait를 사료 대비 성장률이라고 성장률을 높이길 원한다는 가정하면 먼저 성별로 개체를 나누고 EBV를 내림차순으로 정렬한다 (표 4).

Male ID Male EBV Female ID Female EBV
8 0.13 9 0.13
5 0.09 4 0.07
1 0.00 2 -0.09
3 -0.07 7
6 -0.13
table 4

trait 을 더욱 강화하고자 한다는 가정하에서는 male은 8번 개체를 female 개체는 9번 개체를 선택한다.



* 위 내응은 https://www.seeingeye.org/images/HowToUseEBV.pdf 를 바탕으로 함.

GWAS 2

일전에 GWAS (Genome-Wide Association Study)  에 대해 간략하게 포스팅을 한적이 있었는데 이번에 조금 더 자세하게 정리해 보려 한다.


Content
  • Important Questions in Human Genetics
  • Concepts Underlying the Study Design
  • Capturing Common Variation
  • Genotyping Technologies
  • Study Design
  • Association Test
  • Replication and Meta Analysis


Important Questions in Human Genetics

  • human genetics 의 핵심적인 목표는 common complex disease의 유전적 위험 요소를 찾는것. 여러 기술과 study design이 있는데 그중 하나가 GWAS. 
    • 유전적 위험 요소를 찾는 GWAS의 적용사례 중 하나가 age-related macular degeneration (AMD) 의 원인이 되는 Complement Factor H gene을 찾은 것
  • GWAS를 통한 personalized medicine 분야의 적용도 또다른 human genetics의 목표라 할 수 있음
    • GWAS를 통한 항응고제인 wafarin 복용량 결정이 한 예


Concetps Underlying the Study Design


  • 가장 일반적인 SNP을 이용함. SNP 은 일반적으로 MAF (minor allele frequency) 가 높은 single base substitution 을 뜻함
  • rare disorder 의 경우 특정 소수 혹은 단일의 유전자에 의한 것이 많기 때문에 몇몇의 가계도를 분석하는 linkage analysis를 통해 genetic risk factor 가 잘찾아짐. 하지만 common disease의 경우 linkage analysis 가 효과가 없었음
  • common disease는 다르게 접근해야 한다는 생각과 common disease의 원인 유전자가 몇몇 common variant (MAF가 높은 SNP) 이라는 발견과 맞물려 common disease의 원인은 common variant 라는  CD/CV (common disease / common variant) 가정이 성립됨. 이를 바탕으로 하는 것이 GWAS
    • CD/ CV 에 서 분파된 하나의 생각은 common variant가 disease에 영향을 주기 때문에 그 영향의 크기는 상대적으로(rare variant to rare disease에 비해) 적을 것이다(=low penetrance)
    • common variant가 low penentrance이지만 질병이 유전된기때문에 multiple SNP이 원인이 되는게 아닐까 생각하게 됨.
    • 위 두가지 생각으로 family-based에서 population-based study로 초점이 옮겨가게 됨
    • SNP의 allele frequency와 그 SNP에 의한 질병의 위험도(effect size)가 중요 요소가 genetic study design의 핵심. 위 그림에서 우측 하단이 GWAS 연구하기 좋은 케이스 
    • The National Human Genome Institute GWAS catalog 에 GWAS 를 통한 disease causal SNP 이 정리 되어 있음


Capturing Common Variation

CD / CV 가설을 제대로 테스트 하기 위해서는 몇가지 체계적 접근이 필요함

  • Common variant (SNP) 의 파악 필요. 
    • 연구해야 할 SNP들이 어떤것이 있는지 또한  genome 상의 위치 의 정보 필요 
  • population 간 SNP의 차이에 대한 이해 필요
    • population 간에 차이를 고려 하지 않으면 GWAS연구에서 target trait 관련 SNP이 아니라 population 간의 차이가 나는 SNP이 뽑힐 수 있음
  • SNP간의 correlation 정보 필요
    • 중복된 정보를 제거하기 위해 혹은 효율적인 연구를 위해
  • 위 3가지 목적에서 생긴 것이 The International HapMap Project
    • 1.6 M SNPs
    • 11 human population
    • LD (linkage disequilbrium) 측정
      • SNP 이 chromosome의 동일선상에 있을 경우 연관되어 있는 정도를 나타내는 것
      • african이 generation 이 다른 population에 비해 높기 때문에 LD 영역이 상대적으로 적음
      • recombination 이 일어나면 두 SNP은 LD가 깨짐, 곧 linkage equilibrium
      • LD의 지표로 \(D' \) 와 \( r^2\) 를 많이 사용
        • \(D'\)는 0 ~ 1 사이값을 갖음. 1 이 complete LD를 의미, 0 은 HWE 상태 하에 두 SNP이 독립일 경우.
        • \[ D' = \left\{\begin{array}{r1} \frac{\pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB}}{\mbox{min}(\pi_A\pi_b,\pi_a\pi_B)}   & if \pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB} > 0 \\ \frac{\pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB}}{\mbox{min}(\pi_A\pi_B,\pi_a\pi_b)}   & if \pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB} < 0   \end{array} \right\} \]
        • \( r^2\) 은 correlation 식
        • \[r^2 = \frac{(\pi_{AB}\pi_{ab} - \pi_{Ab}\pi_{aB})^2}{\pi_A\pi_b\pi_a\pi_B} \]
        • LD 분석을 통해 tagSNP 이 선택되어짐(LD block내 대표 SNP). population마다 LD가 다르기 때문에 tagSNP도 달라짐
        • 이 LD로 인해 GWAS로 뽑힌 결과가 꼭 질병과 direct하게 연관되어 있다는 보장은 없음. 뽑힌 SNP의 LD 관계에 있는 SNP이 질병의 직접적인 원인 SNP 일 수 있음.



Genotyping Technologies

Illumina 와 Affymerix 사의 SNP chip이 많이 사용됨. 요즘은 sequencing으로 많이 함.


Study Design

case control design 과 quantitative design

  • quantitative design
    • target trait (혹은 disease)의 지표가 quantitative 함. 곧 양적으로 측정됨
    • 통계적으로 유의한 genetic effect를 파악하기 쉬운 경향이 있음
    • heart disease의 지표가 될수 있는 HDL, LDL 등의 lipid level이 quantitative trait의 예
  • case control design
    • trait가 affected VS unaffected 식으로 이원화 된 경우의 study



Association Test 

사실상 GWAS분석이라는 것은 SNP 하나하나를 독립적으로 phenotype과 연관성이 있나 검사하는 것임

  • quantitative design
    • GLM(depedent variable를 함수로 변형시켜 linear 하게 만든 것) 또는 ANOVA로 분석
  • case control design
    • contingency table method(chi-square test, Fisher's exact test), logistic regression을 사용
    • logistic regression은 covariates 의 보정을 할수 있어 contingency table 보다  더 많이 쓰임

genotype data를 숫자로 변경하는 "data encoding"  하는 방식 선택도 테스트의 통계적 검정력에 영향이 있음(genotype-based group의 숫자에 따라 자유도 변화등).

  • allelic VS genotypic
    • allelic : SNP 의 하나의 allele과 phenotype의 연관성을 봄
    • genotypic : 한쌍의 allele 을 하나의 값으로 치환함. 
      • dominant, recessive, multiplicative, addictive model이 있음.
sex, age, study site, population substructure등의 covariate 보정필요

  • population substructure를 고려해야 하는 이유 : 예를 들어 폐암 환자그룹에 특정 population이 많이 있다면 폐암 관련 유전자가 아니라 population specific한 SNP이 뽑힐것임.
  • STRUCTURE나 ENGENSTART 방법으로 보정 필요. HapMap의 인종별 allele freuncy를 참고하거나 아니면 PCA를 이용.

SNP하나하나에 대한 통계적 분석이기 떄문에 multiple test correction 필요

  • Bonferroni correction : 모든  SNP이 independent 하다는 가정
  • FDR (false discover) 
  • permutation testing 
  • genome-wide significance : 특정 population 의 LD의 분포를 바탕으로 독립적인 "effective" genetic region 이 있음을 바탕으로  correction 을 해야하는 test 갯수를 기준으로 p-value threshold를 정함, 유럽인의 경우 7.2e-8


Replication and Meta Analysis


  • replication : GWAS 결과에서 얻어진 SNP을 검증하는 실험
  • meta analysis : 여러 study를 합쳐서 분석하는 방법



* 위 내용은 http://www.ncbi.nlm.nih.gov/pubmed/23300413 과 http://bioinformatics.org.au/ws09/presentations/ 의 day3 jstankovich의 pdf를 바탕으로 함

Wednesday, January 21, 2015

Statistical models in human genetics (lesson 1)

Contents 
  • genetic mapping
  • data for genetic study
  • 3 common question in genetic mapping
  • common designs for genetic studies



Genetic Mapping
형질의 유전패턴과 chromosomal region의 유전 패턴을 비교하여 "positional cloning" 하는 것
- positional cloning : 어떤 유전자인지를 아는 것과 상관없이 어디에 위치하는지를 찾는 것.


Data for a Genetic Study
  • pedigree : 관계가 알려진 개체 모임 (예: 가족)
  • genetic markers : 손쉽게 측정 가능한 genetic variant
  • phgenotype : mendelian phenotype 혹은 complex phenotype

3 Common Question in Genetic mapping
  • 어떠한 형질에 영향을 주는 유전자가 있느냐? Is a trait genetic? => Epidemiological study
      • 집단과 친척관계에서 형질의 분포를 조사함.
      • Crohn's disease 의 경우 전체 인구에서의 발병률과 twins of affected individuals 차이가 남 
    • 그 유전자는 어디에 있느냐? Where are those genes? => Linkage analysis
      • Crohn's 이 있는 가계도와 genetic marker를 비교해 보니 "D16S3136" marker가 연관됨을 확인. "D16S3136" marker 의 위치는 알려져 있음.
    • 그 유전자는 무엇이냐? What are those genes? => Association analysis
      • Crohn's 환자 중 NOD2의 유전자의 mutation 을 갖고 있는 환자가 많음.

    Common Designs for Genetic Studies
    • parametric linkage analysis
      • 특정 모델와 위치를 평가하는 방식(disease loci의 allele frequency, genotype 별 발병률)
      • Identify variants with relatively large contributions to disease risk 
    • allele-sharing analysis
    • association analysis
      • Identify genetic variants with relatively small individual contributions to disease risk 
      • 가장 단순한 방법으로는 case와 control의 allele frequency 비교
      • genome-wide하게 수백만 개의 marker를 가지고 살수도 있고 몇개의 candidate에 초점을 맞춰서 할수도 있음




    * 위 내용은 http://www.sph.umich.edu/csg/abecasis/class/666.01.pdf 를 바탕으로 함

    Tuesday, January 20, 2015

    ipython notebook to pdf

    1. ipython nbconvert sample.ipynb --to latex --post PDF
    trouble: 에러 발생 (원인 불명)

    2. ipython nbconvert sample.ipynb --to latex 후 pdflatex sample.tex
    trouble : 한글의 경우 utf-8 인대도 에러 발생

    3. brew install Caskroom/cask/wkhtmltopdf  깔고
    ipython nbconvert sample.ipynb --to html 후 wkhtmltopdf sample.html sample.pdf
    trouble :  plot이 page 넘어가는 부분에 있을 경우 plot이 쪼개짐