Homework.1
ある1つのSNPと表現型の関連を調べたい。 - $y_k(k=1,\ldots,N)$: Sample $k$ の表現型（中心化: $\mathbb{E}\left[\mathbf{y}\right]=0$） - $x_k(k=1,\ldots,N)$: Sample $k$ の注目する遺伝型（中心化: $\mathbb{E}\left[\mathbf{x}\right]=0$） - $\varepsilon_k(k=1,\ldots,N)$: 残差項、$\varepsilon_k\sim_{\mathrm{i.i.d.}}\mathcal{N}\left(0, \sigma^2\right)$ とし、 $$y_k = \beta x_k + \varepsilon_k$$ のモデルを立てる。Q.1-1 （必須課題）
- $\beta$ の最小二乗（OLS; Ordinary least squares）推定量 $\hat{\beta}$ の値を $\varepsilon_k$ を⽤いて、$\hat{\beta}$ の分布を $\sigma^2$ を用いて求めてください。
- $\beta=0$ の帰無仮説の元で、$\hat{\beta}$ を線形変換して $z\sim\mathcal{N}(0,1)$ に従う確率変数 $z$ と $z^2$ が従う分布を求めてください。 （ちなみに、OLS推定量は最尤推定量と一致する。）
解答
この問題では、データの真の分布が上記の方程式で表され、$\varepsilon_k$ が誤差の実現値であるとみなす。したがって、求めるべきパラメータは $\beta$ であるが、観測されたデータは $(x,y)$ の一部であるため、観測されたデータから推定された $\hat{\beta}$ が真の値であるとは限らない。よって、$\hat{\beta}$ が残差項 $\varepsilon$ の確率分布の影響を受け、ある確率分布に従うと考える。
1
したがって、あるパラメータ $\beta$ の下でSample $k$ の表現型（$x_k$）から推定される遺伝型（$\hat{y_k}$）は、$\hat{y_k}=\beta x_k$ で表される。この推定値と観測値 $y_k$ の二乗誤差を最小にする $\beta$ が最小二乗推定量 $\hat{\beta}_{\text{OLS}}$ なので、二乗誤差
$$S = \sum_k\left(y_k-\hat{y_k}\right)^2 = \left(y_k-\beta x_k\right)^2$$を最小にすることを考えれば良い。よって、$S$ を $\beta$ について微分して、
$$ \begin{aligned} \left.\frac{\partial S}{\partial\beta}\right|_{\beta=\hat{\beta}_{\text{OLS}}} &= \sum_k\left(y_k-\beta x_k\right)\left(-x_k\right)\\ &= \sum_k\beta x_k^2 - \sum_kx_ky_k = 0\quad\cdots(\ast)\\ \therefore\hat{\beta}_{\text{OLS}} &=\frac{\sum_kx_ky_k}{\sum_k x_k^2}\\ \end{aligned} $$によって最小二乗推定量 $\hat{\beta}_{\text{OLS}}$ が求まった。ここで、先述の通り $\hat{\beta}_{\text{OLS}}$ は残差項 $\varepsilon$ の影響を受けてある確率分布に従っているので、
$$ \begin{aligned} \hat{\beta}_{\text{OLS}} &=\frac{\sum_kx_ky_k}{\sum_k x_k^2}\\ &=\frac{\sum_kx_k\left(\beta x_k + \varepsilon_k\right)}{\sum_k x_k^2}\\ &=\beta + \frac{\sum_kx_k\varepsilon_k}{\sum_k x_k^2} \end{aligned} $$と表せる。
ここで、$\hat{\beta}_{\text{OLS}}$ の確率分布を求めるために、統計量を調べる。
まず、$(\ast)$ の式から、
$$ \begin{cases} \begin{aligned} \sum_k\varepsilon_k &= y_k -\beta x_k = 0\\ \sum_kx_k\varepsilon_k &= \sum_kx_ky_k - \sum_k\beta x_k^2 = 0\\ \end{aligned} \end{cases} $$が成り立っていることがわかる。これは、残差の平均が $0$、残差と説明変数 $x_k$ が無相関となるように $\beta$ が推定されていることを表す。
これを用いて $\hat{\beta}_{\text{OLS}}$ の統計量を求めると、
$$ \begin{cases} \begin{aligned} \mathbb{E}\left[\hat{\beta}_{\text{OLS}}\right] &= \beta\\ \mathbb{V}\left[\hat{\beta}_{\text{OLS}}\right] &= \mathbb{E}\left[\left(\hat{\beta}_{\text{OLS}}-\beta\right)^2\right]\\ &= \mathbb{E}\left[\sum_k\varepsilon_k^2\right]\mathbb{E}\left[\frac{\sum_kx_k^2}{\left(\sum_kx_k^2\right)^2}\right]\\ &=\frac{\sigma^2}{\sum_kx_k^2} \end{aligned} \end{cases} $$したがって、$\varepsilon_k\sim_{\mathrm{i.i.d.}}\mathcal{N}\left(0, \sigma^2\right)$ より $\hat{\beta}_{\text{OLS}}$ が正規分布に従うことは明らかなので、上記の統計量より、$\hat{\beta}_{\text{OLS}}\sim\mathcal{N}\left(\beta, \sigma^2/\sum_kx_k^2\right)$
2
$\beta=0$ の帰無仮説の元で、$\mathcal{N}\left(0,1\right)$ に従うように $\hat{\beta}_{\text{OLS}}$ を線形変換するには、$\sqrt{\mathbb{V}\left[\hat{\beta}_{\text{OLS}}\right]}$ で $\hat{\beta}_{\text{OLS}}$ を割ってスケーリングすれば良いので、
$$ \begin{aligned} z &= \frac{\hat{\beta}_{\text{OLS}}}{\sqrt{\mathbb{V}\left[\hat{\beta}_{\text{OLS}}\right]}}\\ &= \frac{ \sum_kx_k\varepsilon_k / \sum_kx_k^2}{\sqrt{\sigma^2 / \sum_kx_k^2}}\\ &= \frac{\sum_kx_k\varepsilon_k}{\sigma\sqrt{\sum_kx_k^2}}\sim\mathcal{N}\left(0, 1\right) \end{aligned} $$したがって、$z$ は平均 $0$、分散 $1$ の正規分布に従う。
ここで、標準正規分布の確率密度関数は $p(x) = \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right)$ とかけるので、$y = f(x) = x^2$ と変数変換する。この時、逆関数が多価関数となるので、変換後の密度関数 $q(y)$ は、
$$q(y) = \sum_ip_i\left(f_i^{-1}(y)\right)\left|\frac{\partial f_i^{-1}(y)}{\partial y}\right|$$と全ての逆関数の和で表される。したがって、
$$ \begin{aligned} q(y)\Delta y &= p(x)\Delta x + p(-x)\Delta x\\ &= 2p(x)\Delta x\quad (y\geq0)\\ q(y) &= 2p\left(\sqrt{y}\right)\left|\frac{\partial \left(\sqrt{y}\right)}{\partial y}\right|\\ &= 2\cdot\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{y}{2}\right)\frac{1}{2\sqrt{y}}\\ &= \frac{1}{2^{\frac{1}{2}}\sqrt{\pi}}y^{-\frac{1}{2}}\exp\left(-\frac{y}{2}\right) \end{aligned} $$一方、自由度 $n$ のカイ二乗分布は、以下の確率密度関数 $(\ast)$ で表される分布であるので、$z^2$ が自由度 $1$ のカイ二乗分布に従うことがわかる。
$$ \begin{cases} \begin{aligned} f_n(x) &= \frac{1}{2^{\frac{n}{2}}\Gamma\left(\frac{n}{2}\right)}x^{\frac{n}{2}-1}\exp\left(-\frac{x}{2}\right)\quad\cdots\left(\ast\right)\\ \Gamma\left(\frac{1}{2}\right) &= \sqrt{\pi} \end{aligned} \end{cases} $$Q.1-2 （選択課題）
$\sigma^2$ を不偏推定量 $\hat{\sigma}^2=c\sum_k\left(y_k-\hat{\beta}x_k\right)^2$ （$c$ は定数）に置換するとき、$\hat{\sigma}^2$ の値を求めて下さい。また、$\beta=0$ の帰無仮説の元で $z_{\hat{\sigma}^2}$ が自由度 $N-1$ のt分布に従うことを示して下さい。
解答
全 $N$ サンプルの不偏分散なので、$c=N-1$ であり、
$$ \begin{aligned} \hat{\sigma}^2 &= \frac{1}{N-1}\sum_{k=1}^N\left(y_k - \hat{\beta}_{\text{OLS}}\right)^2\\ &= \frac{1}{N-1}\sum_{k=1}^N\left(y_k - \left(\beta + \frac{\sum_jx_j\varepsilon_j}{\sum_j x_j^2}\right)x_k\right)^2\\ &= \frac{1}{N-1}\sum_{k=1}^N\left(y_k - \beta x_k - \frac{\sum_jx_j\varepsilon_j}{\sum_j x_j^2}x_k\right)^2\\ &= \frac{1}{N-1}\sum_{k=1}^N\left(y_k-\beta x_k\right)^2\\ &= \frac{1}{N-1}\sum_{k=1}^N\varepsilon_k^2\\ &= \frac{\sigma^2}{N-1}\sum_{k=1}^N\left(\frac{\varepsilon_k}{\sigma}\right)^2 \end{aligned} $$よって、$\beta=0$ の帰無仮説の元で、$z$ は
$$ \begin{aligned} z_{\hat{\sigma}^2} &= \frac{\sum_kx_k\varepsilon_k}{\hat{\sigma}\sqrt{\sum_kx_k^2}}\\ &= \frac{\sum_kx_k\varepsilon_k}{\sqrt{\frac{\sigma^2}{N-1}\sum_{k=1}^N\left(\frac{\varepsilon_k}{\sigma}\right)^2}\sqrt{\sum_kx_k^2}}\\ &= \frac{\sum_kx_k\varepsilon_k}{\sigma\sqrt{\sum_kx_k^2}}/\sqrt{\sum_{k=1}^N\left(\frac{\varepsilon_k}{\sigma}\right)^2/\left(N-1\right)} \end{aligned} $$ここで、
$$ \begin{aligned} \frac{\sum_kx_k\varepsilon_k}{\sigma\sqrt{\sum_kx_k^2}} &\sim \mathcal{N}(0,1)\\ \sum_{k=1}^N\left(\frac{\varepsilon_k}{\sigma}\right)^2&\sim\mathcal{\chi}^2_N \end{aligned} $$なので、$z_{\hat{\sigma}^2}$ が自由度 $N-1$ のt分布に従うことは示せなかった。
どのようにして自由度 $N-1$ のカイ二乗分布に従う確率変数を導くのかわかりませんでした。
Q.1-3 （選択課題）
$\beta\neq0$ の仮定の元で、Q.1-1 の $z$ はどのような分布に従うかを $\beta,\sigma^2$ を用いて求めてください。また、$z^2$ はどのような分布に従うかを求めてください。
解答
$\beta\neq0$ の仮定のもとで Q.1-1 と同様の線形変換を行うと、
$$ \begin{aligned} z &= \frac{\hat{\beta}_{\text{OLS}}}{\sqrt{\mathbb{V}\left[\hat{\beta}_{\text{OLS}}\right]}}\\ &= \frac{\beta+\left( \sum_kx_k\varepsilon_k / \sum_kx_k^2\right)}{\sqrt{\sigma^2 / \sum_kx_k^2}}\\ &= \frac{\beta}{\sigma}\sqrt{\sum_kx_k^2}+\frac{\sum_kx_k\varepsilon_k}{\sigma\sqrt{\sum_kx_k^2}}\sim\mathcal{N}\left(\frac{\beta}{\sigma}\sqrt{\sum_kx_k^2}, 1\right) \end{aligned} $$したがって、$z$ は平均 $\frac{\beta}{\sigma}\sqrt{\sum_kx_k^2}$ 分散 $1$ の正規分布に従う。
ここで、上記の正規分布は、$\mu = \frac{\beta}{\sigma}\sqrt{\sum_kx_k^2}$ として $f(x)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^2}{2}\right)$ とかけるので、先ほどと同様に $y=x^2$ と変数変換すると、
$$ \begin{aligned} q(y) &= \frac{1}{2\sqrt{y}}\frac{1}{\sqrt{2\pi}}\left[\exp\left(-\frac{\left(y-2\sqrt{y}\mu+\mu^2\right)}{2}\right) + \exp\left(-\frac{\left(y+2\sqrt{y}\mu+\mu^2\right)}{2}\right)\right]\\ &= \frac{1}{2\sqrt{y}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{y + \mu^2}{2}\right)\left[\exp\left(\sqrt{y}\mu\right) + \exp\left(-\sqrt{y}\mu\right)\right]\\ &= \frac{1}{\sqrt{2\pi y}}\exp\left(-\frac{y + \lambda}{2}\right)\cosh\left(\sqrt{y\lambda}\right),\quad\left(\lambda = \mu^2\right)\\ &= \frac{1}{\sqrt{2\pi y}}\exp\left(-\frac{y + \lambda}{2}\right)\sum_{n=0}^{\infty}\frac{\left(\sqrt{y\lambda}\right)^{2n}}{2n!}\\ &= \frac{1}{\sqrt{2\pi y}}\exp\left(-\frac{y + \lambda}{2}\right)\sum_{n=0}^{\infty}\frac{\left(y\lambda\right)^{n}}{2n!} \end{aligned} $$一方、自由度 $k$、非心パラメータ $\lambda$ の $\chi^2$ 分布の確率密度関数が
$$f(x;k,\lambda) = \frac{1}{2}\exp\left(-\frac{x + \lambda}{2}\right)\left(\frac{x}{2}\right)^{\frac{k}{2}-1}\sum_{n=0}^{\infty}\frac{1}{n!\Gamma\left(\frac{k}{2} + n\right)}\left(\frac{\lambda z}{4}\right)^n$$と表される[Wiki]ことから、$z^2$ が自由度 $1$、非心パラメータ $\lambda=\mu^2 = \frac{\beta^2}{\sigma^2}\sum_kx_k^2$ の非心カイ二乗分布に従うことがわかる。
$k=1$ とした時、非心カイ二乗分布は、
$$ \begin{aligned} f(x;1,\lambda) &= \frac{1}{2}\exp\left(-\frac{x+\lambda}{2}\right)\left(\frac{x}{2}\right)^{-\frac{1}{2}}\sum_{n=0}^{\infty}\frac{1}{n!\Gamma\left(\frac{1}{2} + n\right)}\left(\frac{\lambda x}{4}\right)^n\\ &= \frac{1}{\sqrt{2x}}\exp\left(-\frac{x+\lambda}{2}\right)\sum_{n=0}^{\infty}\frac{1}{n!\Gamma\left(\frac{1}{2} + n\right)}\left(\frac{\lambda x}{4}\right)^n \end{aligned} $$と書けるので、先ほどの式と比べると、
$$\frac{1}{\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{\left(x\lambda\right)^{n}}{2n!} = \sum_{n=0}^{\infty}\frac{1}{n!\Gamma\left(\frac{1}{2} + n\right)}\left(\frac{\lambda x}{4}\right)^n$$が言えれば良い。ここで、
$$ \begin{aligned} \frac{1}{n!\Gamma\left(\frac{1}{2} + n\right)}\left(\frac{\lambda x}{4}\right)^n &= \frac{1}{n!\frac{\left(2n - 1\right)!!}{2^n}\sqrt{\pi}}\left(\frac{\lambda x}{4}\right)^n\\ &= \frac{1}{\sqrt{\pi}}\frac{1}{n!\left(2n-1\right)!!}\left(\frac{\lambda x}{2}\right)^n\\ &= \frac{1}{\sqrt{\pi}}\frac{1}{n!\frac{(2n)!}{2^nn!}}\left(\frac{\lambda x}{2}\right)^n\\ &= \frac{1}{\sqrt{\pi}}\frac{\left(\lambda x\right)^n}{2n!} \end{aligned} $$と変換できるので、正しさが確認された。
Q.1-4 （選択課題）
上記の解析（全tag SNPに対して独立に、表現型との間に線形な関係があると仮定して線形回帰を行う）を全SNPについて行うとする。$i$ 番目のSNPの効果量を $\beta_i$、$q_i:=\mathrm{Var}\left[\beta_ix_i\right]/\mathrm{Var}[y]$ としたとき、Q.1-3で求めた非心パラメータ $\lambda_i$ を $q_i^2$ を用いて表してください。
解答
$i$ 番目のSNPに関するモデルは以下のように表される。（$y_{ik}$ は、$i$ に依らず、Sample $k$ のみに依存）
$$y_{ik} = \beta_ix_{ik} + \varepsilon_{ik}$$したがって、誤差項が独立であることから分散の線型性を用いると、
$$ \begin{aligned} \mathbb{V}\left[\mathbf{y}_i\right] &= \mathbb{V}\left[\beta_i\mathbf{x}_{ik}\right] + \mathbb{V}\left[\boldsymbol{\varepsilon}_i\right]\\ &= \mathbb{V}\left[\beta_i\mathbf{x}_{ik}\right] + \sigma_i^2\\ \end{aligned} $$であるので、これを用いると、
$$ \begin{aligned} 1/q_i^2 &= \mathbb{V}\left[\mathbf{y}_i\right] / \mathbb{V}\left[\beta_i\mathbf{x}_{ik}\right]\\ &= 1 + \frac{\sigma_i^2}{\mathbb{V}\left[\beta_i\mathbf{x}_{ik}\right]}\\ &= 1 + \frac{\sigma_i^2}{\beta_i^2\sum_kx_{ik}^2} \end{aligned} $$という関係が導かれる。以上より、非心パラメータ $\lambda_i$ を $q_i^2$ で表すと
$$ \lambda_i= \frac{\beta_i^2}{\sigma_i^2}\sum_kx_{ik}^2 = \frac{q_i^2}{1-q_i^2} $$Q.1-5 （選択課題）
全 $M$ 個のSNPの $\chi^2$ 統計量(null仮説の時 $\chi^2$ 分布に従う統計量)のみからこの表現型の遺伝率 $h^2$ を推定するためにはどうすれば良いか、ただし、LD（連鎖不均衡）は存在しないものとし、causal variants は全て Genotyping されており、complex trait ($q_i\ll1$) と仮定して良い。
解答
遺伝率は、表現型の分散のうち遺伝型で説明される分散の割合を表す。つまり、
$$\underbrace{y_k}_{P\text{: Phenotype}} = \underbrace{\sum_i^M\beta_ix_{ik}}_{G\text{: Gene}} + \underbrace{\varepsilon_k}_{E\text{: Environment}}$$とおいた時に、$h^2:=\mathbb{V}\left[G\right]/\mathbb{V}\left[P\right]$ として定義される。
ここで、complex diseaseでは、GWAS解析において少数の有意なSNPと、多数の弱い関連SNPがある。そこで、$h_{\text{GWAS}}$（表現型の分散のうち、GWASで有意なSNPのみで説明できる分散の割合）を求める。これは、各SNPの効果量のGWAS推定値 $\hat{\beta_i}$ を重みとし、線形和として表現型の推定値 $\hat{\phi_i}$ を
$$\hat{\phi_i} = \sum_{\text{Significant SNP i}}\hat{\beta_i}x_{ik}$$で求め、$h_{\text{GWAS}}^2 = r^2\left(\hat{\boldsymbol{\phi}},\mathbf{y}\right)$ という関係を用いて求める。
Homework.2
与えられた身長のSimulation dataを見て、気が付いたことを答えてください。- PCAを行い、サンプル内にクラスタが存在するか確認せよ。（必須課題）
- QQplotを書いてみよ。（必須課題）
- 余裕があれば、後述のPrediction(Risk Score)も求めよ。
解答
| ファイル名 | 概要 | カラムの説明 |
|---|---|---|
simu.legend |
SNPに関する情報が書いてあります。1行が1つのSNPに対応しています。 |
|
simu.phe |
サンプルに関する情報が書いてあります。1行が1サンプルに対応しています。 |
|
simu.genot |
遺伝型情報が書いてあります。1行が1サンプルに、1列が1SNPに対応しています。順番はsimu.genot,simu.pheと同じです。 |
各要素はminor alleleの個数 |
import numpy as np
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt
df_legend = pd.read_csv("simu/simu.legend")
df_phe = pd.read_csv("simu/simu.phe")
df_genot = pd.read_csv("simu/simu.genot", header=None, sep=" ")
print(f"Total num of Samples (N): {len(df_phe)}")
print(f"Total num of SNPs (M): {len(df_legend)}")
print(f"(M,N) = {df_genot.shape}")
Y = df_phe.height
X = np.asarray(df_genot).T
# 標準化
X = (X - np.average(X, axis=0)) / np.std(X, axis=0)
Y = (Y - np.average(Y)) / np.std(Y)
PCAを行い、サンプル内にクラスタが存在するか確認せよ。¶
from sklearn.decomposition import PCA
print(f"(N,M) = {X.shape}")
pca = PCA(n_components=2)
pca.fit(X)
X_transformed = pca.transform(X)
X_pc1, X_pc2 = X_transformed.T
plt.figure(figsize=(6,6))
plt.scatter(X_pc1, X_pc2, s=10)
plt.title("PCA", fontsize=16), plt.xlabel("Principle Component 1"), plt.ylabel("Principle Component 2")
plt.show()
上のプロットより、サンプル内にクラスタが存在することが確認できた。そこで、試しにKMeans法を用いてクラスを分割する
from kerasy.utils import findLowerUpper
from sklearn.cluster import KMeans
model = KMeans(n_clusters=3, random_state=0)
model.fit(X_transformed)
Xlus = findLowerUpper(X_transformed, margin=5e-2, N=100)
X1lus,X2lus = np.meshgrid(*Xlus)
Xs = np.c_[X1lus.ravel(),X2lus.ravel()]
Zs = model.predict(Xs).reshape(X1lus.shape)
cls = model.predict(X_transformed)
plt.figure(figsize=(6,6))
plt.scatter(X_pc1, X_pc2, s=15, c=cls, cmap="jet")
plt.pcolor(X1lus,X2lus,Zs,alpha=0.2, cmap="jet")
plt.title("PCA", fontsize=16), plt.xlabel("Principle Component 1"), plt.ylabel("Principle Component 2")
plt.show()
QQplotを書いてみよ。¶
各SNPと表現型の間に以下のモデルを立てる。
$$y_k = \beta x_k + \varepsilon_k$$def calc_uniform_order_statistic_medians(sample_size):
v = np.zeros(sample_size, dtype=np.float64)
v[-1] = 0.5**(1.0 / sample_size)
v[0] = 1 - v[-1]
i = np.arange(2, sample_size)
v[1:-1] = (i - 0.3175) / (sample_size + 0.365)
return v
betas = np.asarray([np.sum(x*Y)/np.sum(x*x) for x in X.T])
sample_size = len(betas)
osm_uniform = calc_uniform_order_statistic_medians(sample_size)
osm = stats.norm.ppf(osm_uniform)
osr = np.sort(betas)
plt.scatter(osm, osr, color="b", s=10)
plt.title("QQ Plot"), plt.xlabel('Theoretical quantiles'), plt.ylabel("Ordered Values"), plt.grid()
plt.show()
Homework.3
Q.1では、各SNPごとに単回帰を行い、$\hat{\beta}_i$ を推定していたが、ここでは、全SNPによる効果と表現型の関連を同時に考えたい。 - $\mathbf{y} = \left(y_k\right)$: 従属変数ベクトル（表現型） - $\boldsymbol{\beta} = \left(\beta_i\right)$: 効果量ベクトル - $\mathbf{X} = \left(x_{k,i}\right)$: 独立変数行列 - $\boldsymbol{\varepsilon}=\left(\varepsilon_k\right)$: 残差項ベクトル $\left(\sim_{\text{i.i.d.}}\mathcal{N}\left(0,\sigma^2\mathbf{I}_n\right)\right)$ $$ \begin{aligned} y_{k} = \sum_{i=1}^M\beta_ix_{ik} + \varepsilon_{k}\quad \left(\underbrace{\mathbf{y}}_{(N)} = \underbrace{\mathbf{X}}_{(N\times M)}\underbrace{\boldsymbol{\beta}}_{(M)} + \underbrace{\boldsymbol{\varepsilon}}_{(N)}\right) \end{aligned} $$ の多重線形モデルを考える。また、この時 $x_{ik},y_{k}$ は標準化されているものとする。Q.3-1 （必須課題）
重回帰を用いて遺伝型 $\mathbf{X}$ と表現型 $\mathbf{y}$ から最小二乗推定値 $\tilde{\boldsymbol{\beta}}$ を求めてください。また、LDがないとき、$\tilde{\beta_i}$ の分布を $\beta_i$ を用いて求めてください。
解答
最小二乗推定量 $\hat{\boldsymbol{\beta}}_{\text{OLS}}$ は、残差平方和
$$s:= \left\|\mathbf{y}-\hat{\mathbf{y}}\right\|^2 = \left(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\right)^T\left(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\right)$$を最小化するので、
$$ \begin{aligned} \left.\frac{\partial s}{\partial\boldsymbol{\beta}}\right|_{\boldsymbol{\beta}=\hat{\boldsymbol{\beta}}_{\text{OLS}}} &= \frac{\partial}{\partial\boldsymbol{\beta}}\left(\mathbf{y}^T\mathbf{y} - 2\boldsymbol{\beta}^T\mathbf{X}^T\mathbf{y} + \boldsymbol{\beta}^T\mathbf{X}^T\mathbf{X}\boldsymbol{\beta}\right)\\ &= -2\mathbf{X}^T\mathbf{y} + 2\mathbf{X}^T\mathbf{X}\boldsymbol{\beta}\quad\cdots\left(\ast\right)\\ &= \mathbf{0}\\ \therefore\hat{\boldsymbol{\beta}}_{\text{OLS}} &=\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y} \end{aligned} $$ここで、$D:=\left(\mathbf{X}^T\mathbf{X}\right)$ の各要素は、
$$ D_{ij} = \mathbf{x}_i^T\mathbf{x}_j = \begin{cases} N\mathrm{Var}\left[\mathbf{x}_i\right] = N & (i=j)\\ N\mathrm{Cov}\left[\mathbf{x}_i,\mathbf{x}_j\right] & (i\neq j) \end{cases} $$LDが存在しない時、つまり、SNP間に相関がない時 $\left(\mathrm{Cov}\left[\mathbf{x}_i,\mathbf{x}_j\right] = 0\left(i\neq j\right)\right)$ 、全てのSNPが独立なので、$\hat{\boldsymbol{\beta}}_{\text{OLS}}=\mathbf{X}^T\mathbf{y}/N$ とかける。
これより、
$$ \begin{aligned} \hat{\boldsymbol{\beta}}_{\text{OLS}_i} &= \frac{1}{N}\left(\mathbf{X}^T\mathbf{y}\right)_i\\ &= \frac{1}{N}\left(\mathbf{X}^T\left(\mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\right)\right)_i\\ &= \left(\boldsymbol{\beta} + \frac{1}{N}\mathbf{X}^T\boldsymbol{\varepsilon}\right)_i \end{aligned} $$と分解できるので、
$$ \begin{cases} \begin{aligned} \mathbb{E}\left[\hat{\boldsymbol{\beta}}_{\text{OLS}_i}\right] &= \boldsymbol{\beta}_i\quad \left(\because \text{equation }(\ast)\right)\\ \mathbb{V}\left[\hat{\boldsymbol{\beta}}_{\text{OLS}_i}\right] &= \boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon} = \frac{1}{N}\mathbb{V}\left[\boldsymbol{\varepsilon}\right] \end{aligned} \end{cases} $$と $\boldsymbol{\beta}$ の統計量がわかる。この時、遺伝率が
$$\underbrace{y_k}_{P\text{: Phenotype}} = \underbrace{\sum_i^M\beta_ix_{ik}}_{G\text{: Gene}} + \underbrace{\varepsilon_k}_{E\text{: Environment}}$$において $h^2:=\mathbb{V}\left[G\right]/\mathbb{V}\left[P\right]$ と表されていたことを考えれば、$\mathbb{V}\left[\boldsymbol{\varepsilon}\right] = 1-h^2$ であることがわかる。以上より、
$$\hat{\boldsymbol{\beta}}_{\text{OLS}_i}\sim\mathcal{N}\left(\beta_i,\frac{1-h^2}{N}\right)$$Q.3-2 （選択課題）
LDがなく、Infinitesimal Modelの場合、$\beta_i|\tilde{\beta_i}$ が従う事後分布を求めてください。
解答
Q.3-2,Q.3-3 では、ベイズの定理を用いて $\beta_i$ の事後分布を求める。
$$p\left(\beta_i|\hat{\beta}_{\text{OLS}_i}\right)\propto p(\beta_i)p\left(\hat{\beta}_{\text{OLS}_i}|\beta_i\right)$$Infinitensimal Model では、全SNPが小さな効果量を持つと仮定し、全SNPに等しく分散を割り当てる。
$$\beta_i\sim\mathcal{N}\left(0,\frac{h^2}{M}\right)$$したがって、
$$ \begin{aligned} p\left(\beta_i|\hat{\beta}_{\text{OLS}_i}\right) &\propto p(\beta_i)p\left(\hat{\beta}_{\text{OLS}_i}|\beta_i\right)\\ &=\mathcal{N}\left(0,\frac{h^2}{M}\right)\mathcal{N}\left(\beta_i,\frac{1-h^2}{N}\right)\\ &\propto\exp\left(-\frac{1}{2\frac{h^2}{M}}\left(x-0\right)^2\right)\exp\left(-\frac{1}{2\frac{1-h^2}{N}}\left(x-\beta_i\right)^2\right)\\ &=\exp\left(-\frac{1}{2}\left(\frac{M}{h^2}x^2 + \frac{N}{1-h^2}\left(x^2-2\beta_ix + \beta_i^2\right)\right)\right)\\ &\propto\exp\left(-\frac{1}{2}\left(\frac{M(1-h)^2+Nh^2}{h^2(1-h)^2}\right)\left(x^2 - 2\frac{Nh^2\beta_i}{M(1-h)^2+Nh^2}\right)x\right)\\ &\propto\exp\left(-\frac{1}{2}\left(\frac{M(1-h)^2+Nh^2}{h^2(1-h)^2}\right)\left(x - \frac{Nh^2\beta_i}{M(1-h)^2+Nh^2}\right)^2\right)\\ &= \mathcal{N}\left(\frac{Nh^2\beta_i}{M(1-h)^2+Nh^2},\frac{h^2(1-h)^2}{M(1-h)^2+Nh^2}\right) \end{aligned} $$Q.3-3 （選択課題）
LDがなく、Non-infinitesimal Modelの場合、$\beta_i|\tilde{\beta_i}$ が従う事後分布を求めてください。
解答
Non-infinitensimal Model では、割合 $\rho$ のSNPのみ効果量を持つと仮定し、それらに等しく分散を割り当てる。
$$ \beta_i\sim \begin{cases}\begin{aligned} &\mathcal{N}\left(0,\frac{h^2}{M_{\rho}}\right) & \left(\text{with probability $\rho$}\right)\\ &0 & \left(\text{else}\right) \end{aligned}\end{cases} $$よって、Q.3-2と同様に考えて、
$$ p\left(\beta_i|\hat{\beta}_{\text{OLS}_i}\right)\sim \begin{cases}\begin{aligned} &\mathcal{N}\left(\frac{Nh^2\beta_i}{M_{\rho}(1-h)^2+Nh^2},\frac{h^2(1-h)^2}{M_{\rho}(1-h)^2+Nh^2}\right) & \left(\text{with probability $\rho$}\right)\\ &0 & \left(\text{else}\right) \end{aligned}\end{cases} $$