2024-10-01

混合モデル方程式から回帰係数の期待値/分散を導出

はじめに

次のような混合モデルを考えます。

$y=X\beta+Zu+e$
where $u \sim MVN(0,K\sigma_{g}^2), e \sim MVN(0,I\sigma_{e}^2)$

この時、このモデルのβとuのBLUE及びBLUPは次の方程式を解くことで求めらることが分かっています(混合モデル方程式)。

$\begin{bmatrix} X^{T}X & X^{T}Z \\ Z^{T}X & Z^{T}Z+G^{-1} \\ \end{bmatrix}\begin{bmatrix} \beta \\ u \\ \end{bmatrix}=\begin{bmatrix} X^{T}y \\ Z^{T}y \\ \end{bmatrix}$

また、この時のパラメータの分散も次のように求まることが知られています。

$\begin{bmatrix} Var(\beta) & Cov(\beta, u)\\ Cor(\beta, u) & Var(u) \\ \end{bmatrix}=\begin{bmatrix} X^{T}X & X^{T}Z \\ Z^{T}X & Z^{T}Z+G^{-1} \\ \end{bmatrix}^{-1}\sigma_{e}^2$

非常に覚えやすい形をしていて個人的に気に入っているのですが、今回の私はこれに関する話題です。

一般に線形回帰 $y_{i}=\beta_{0}+\beta_{1}x_{i}+e_{i}$ の解は最小二乗法を解くことにより求まります。

そしてこの最小二乗解は
$\beta=(\beta_{0}, \beta_{1})^T , X= \begin{bmatrix}1&x_{1}\\1&x_{2}\\ \vdots & \vdots \\1&x_{n}\end{bmatrix}$ としたとき、
$X^{T}X \beta = X^{T}y$ を解くことにより求まるわけですがこれは混合モデル方程式のuが無いときに対応しています。

この解であるβの分散も知ることができれば傾きや切片自体に対して統計的検定(例えば傾きが有意に0と異なるか、など)を行えるのですが、この分散はやや分かりにくい形をしています。

$Var(\beta_{0})=\frac{\sigma_{e}^2\Sigma x_{i}^2}{n\Sigma(x_{i}-\bar{x})^2}$
$Var(\beta_{1})=\frac{\sigma_{e}^2}{\Sigma(x_{i}-\bar{x})^2}$

さて、最小二乗解が混合モデル方程式のuが無いバージョンに対応しているように、この分散も混合モデルの分散を求める式のuが無いところ $Var(\beta)=(X^{T}X)^{-1}\sigma_{2}^2$ に対応していると考えられます。
そこで今回は、この線形回帰の分散を混合モデルの分散を求める方程式から求めてみようと思います。

問題設定

$Var(\beta)=(X^{T}X)^{-1}\sigma_{2}^2$ 　から　 $Var(\beta_{0})=\frac{\sigma_{e}^2\Sigma x_{i}^2}{n\Sigma(x_{i}-\bar{x})^2}, Var(\beta_{1})=\frac{\sigma_{e}^2}{\Sigma(x_{i}-\bar{x})^2}$ 　を求める。

計算

まず、 $X= \begin{bmatrix}1&x_{1}\\1&x_{2}\\ \vdots & \vdots \\1&x_{n}\end{bmatrix}$ であることを利用して $X^{T} X$ を計算します

$X^{T}X=\begin{bmatrix} 1 & 1 & \dots & 1\\ x_{1} & x_{2} & \dots & x_{n}\end{bmatrix}\begin{bmatrix}1&x_{1}\\1&x_{2}\\ \vdots & \vdots \\1&x_{n}\end{bmatrix}=\begin{bmatrix} n & \Sigma x_{i} \\ \Sigma x_{i} & \Sigma x_{i}^2 \end{bmatrix}$

2×2の行列の逆行列はすでに公式が知られているのでこれに照らし合わせると
$(X^{T}X)^{-1}=\begin{bmatrix} n & \Sigma x_{i} \\ \Sigma x_{i} & \Sigma x_{i}^2 \end{bmatrix}^{-1}=\frac{1}{n\Sigma x_{i}^2-(\Sigma x_{i})^2}\begin{bmatrix} \Sigma x_{i}^2 & -\Sigma x_{i} \\ -\Sigma x_{i} & n \end{bmatrix}$

ここでよく知られた分散の公式から $n\Sigma x_{i}^2-(\Sigma x_{i})^2 = n^2(\Sigma x_{i}^2/n-(\Sigma x_{i}/n)^2=n\Sigma(x_{i}-\bar{x})$ であるので、

$(X^{T}X)^{-1}=\frac{1}{n\Sigma(x_{i}-\bar{x})}\begin{bmatrix} \Sigma x_{i}^2 & -\Sigma x_{i} \\ -\Sigma x_{i} & n \end{bmatrix}$

よって、

$Var(\beta)=\frac{\sigma_{e}^2}{n\Sigma(x_{i}-\bar{x})}\begin{bmatrix} \Sigma x_{i}^2 & -\Sigma x_{i} \\ -\Sigma x_{i} & n \end{bmatrix}$

ここで以下が成り立つ
$Var(\beta)=\begin{bmatrix} Var(\beta_{0}) & Cor(\beta_{0}, \beta_{1}) \\ Cor(\beta_{0}, \beta_{1}) & Var(\beta_{1}) \end{bmatrix}$

よって、
$Var(\beta_{0})=\frac{\sigma_{e}^2\Sigma x_{i}^2}{n\Sigma(x_{i}-\bar{x})^2}$
$Var(\beta_{1})=\frac{\sigma_{e}^2}{\Sigma(x_{i}-\bar{x})^2}$

まとめ

今回、マトリックスを分解することで、混合モデルの分散の式から線形回帰のパラメータの分散を求めることできるとわかりました。
混合モデルの分散の式を覚えるのはそんな難しくはないので、パラメータの分散の式を忘れたときは、今後はここに立ち戻って計算しようと思います

2024-07-14

GBLUPの変量効果からマーカー効果を推定②：マルチカーネル版

線形混合モデル

はじめに
ひとつ前の記事でGBLUPの変量効果からマーカー効果を推定する方法の背景にある計算を紹介しました。そこで紹介したGBLUPのモデルには変量効果が一つしか含まれませんでしたが、場合によってはこの変量効果が二つ以上になる場合がありえます。そのように変量効果が二つ以上入る混合モデルはマルチカーネルと呼ばれますが、今日はそのようなマルチカーネルになったGBLUPにおいて、マーカー効果を推定する方法を計算します。(先に結論を述べると先週紹介した計算式と全く同じになります)。なお、以下の計算は先日の方法をベースに拡張しています。まだ先日の記事を読んでいない方はそちらを読むことを強くお勧めします。

マルチカーネル
まずはマルチカーネルのモデル、及びその時の混合モデル方程式を紹介します。マルチカーネルのモデルは一般に以下のように定式化されます。
$y=X\beta+Zu+Wv+e$
$u\sim MVN(0,G_{u})$
$v\sim MVN(0,G_{v})$
$e\sim MVN(0,R)$
そしてその時の混合モデル方程式は次のように書かれます。
(この導出はいろいろなところで見ることができるので省略します)

$\begin{bmatrix} X^{T}R^{-1}X & X^{T}R^{-1}Z & X^{T}R^{-1}W \\ Z^{T}R^{-1}X & Z^{T}R^{-1}Z+G_{u}^{-1} & Z^{T}R^{-1}W\\ W^{T}R^{-1}X & W^{T}R^{-1}Z & W^{T}R^{-1}W+G_{v}^{-1}\\ \end{bmatrix}\begin{bmatrix} \beta \\ u \\ v\\ \end{bmatrix}=\begin{bmatrix} X^{T}R^{-1}y \\ Z^{T}R^{-1} \\ W^{T}R^{-1} \\ \end{bmatrix}$

これを行列表示せず、別々の式として表すと3つの式が出てきます。その中でも今回はuとvがかかわる次の2つに注目します。
$Z^{T}R^{-1}Z\beta+(Z^{T}R^{-1}Z+G_{u}^{-1})u+ Z^{T}R^{-1}Wv=Z^{T}R^{-1}$
$W^{T}R^{-1}X\beta+W^{T}R^{-1}Zu + (W^{T}R^{-1}W+G_{v}^{-1})v=W^{T}R^{-1}$

これをu,vそれぞれについて解くと次の式が出てきます。
$u=(Z^{T}R^{-1}Z+G_{u}^{-1}-Z^{T}R^{-1}W(W^{T}R^{-1}W+G_{v}^{-1})^{-1}(W^{T}R^{-1}Z))^{-1}Z^{T}R^{-1}(y-X\beta-W(W^{T}R^{-1}W+G_{v}^{-1})^{-1}W^{T}R^{-1}(y-X\beta))$
$v=(W^{T}R^{-1}W+G_{v}^{-1}-W^{T}R^{-1}Z(Z^{T}R^{-1}Z+G_{u}^{-1})^{-1}(Z^{T}R^{-1}W))^{-1}W^{T}R^{-1}(y-X\beta-Z(Z^{T}R^{-1}Z+G_{u}^{-1})^{-1}Z^{T}R^{-1}(y-X\beta))$
これが一般的なマルチカーネルの混合モデルのu,vに対する解となります。
以下ではより具体的なモデルを仮定し、その場合におけるuとvの具体形を求めていきます。

GBLUPの場合
GBLUPの場合、次のような式として記述することができます
$y=X\beta+u+v+e$
$u\sim MVN(0,K_{u}\sigma_{u}^2)$
$v\sim MVN(0,K_{v}\sigma_{v}^2)$
$e\sim MVN(0,R)$
先ほど提示した一般的な混合モデルにおいて、Z,Wを単位行列に、uとvが従う分散共分散行列を変形すればよいことが分かります。
従って、先ほどのuとvに対する2式は以下のように書き下せます。
$u=(R^{-1}+K_{u}^{-1}\sigma_{u}^{-2}-R^{-1}(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1})^{-1}R^{-1}(y-X\beta-(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1}(y-X\beta))$
$v=(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2}-R^{-1}(R^{-1}+K_{u}^{-1}\sigma_{u}^{-2})^{-1}R^{-1})^{-1}R^{-1}(y-X\beta-(R^{-1}+K_{u}^{-1}\sigma_{u}^{-2})^{-1}R^{-1}(y-X\beta))$

マーカー効果推定の場合
次にマーカー効果を推定する場合のモデルを考えます。それは以下のように書けるでしょう。
$y=X\beta+Z\mu+W\gamma+e$
$\mu\sim MVN(0,I\sigma_{\mu}^2))$
$\gamma\sim MVN(0,I\sigma_{\gamma}^2)$
$e\sim MVN(0,R)$
よって、
$\mu=(Z^{T}R^{-1}Z+I\sigma_{\mu}^{-2}-Z^{T}R^{-1}W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}(W^{T}R^{-1}Z))^{-1}Z^{T}R^{-1}(y-X\beta-W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}R^{-1}(y-X\beta))$
$\gamma=(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2}-W^{T}R^{-1}Z(Z^{T}R^{-1}Z+I\sigma_{\mu}^{-2})^{-1}(Z^{T}R^{-1}W))^{-1}W^{T}R^{-1}(y-X\beta-Z(Z^{T}R^{-1}Z+I\sigma_{\mu}^{-2})^{-1}Z^{T}R^{-1}(y-X\beta))$

両方のモデルの関係
今紹介した二つの式は
$u=Z\mu$
$v=W\gamma$
によって結ばれます。
また、これより自動的に以下の関係が導かれます。
$K_{u}\sigma_{u}^{2}=ZZ^{T}\sigma_{\mu}^{2}$
$K_{v}\sigma_{v}^{2}=WW^{T}\sigma_{\gamma}^{2}$

マーカー効果とGBLUPの変量効果の関係
さて、いまuとμの式に注目します。
$u=(R^{-1}+K_{u}^{-1}\sigma_{u}^{-2}-R^{-1}(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1})^{-1}R^{-1}(y-X\beta-(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1}(y-X\beta))$
$\mu=(Z^{T}R^{-1}Z+I\sigma_{\mu}^{-2}-Z^{T}R^{-1}W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}(W^{T}R^{-1}Z))^{-1}Z^{T}R^{-1}(y-X\beta-W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}R^{-1}(y-X\beta))$
両式の後半部分に注目すると次に示すよく似たパーツを持っていることが分かります。
$(y-X\beta-(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1}(y-X\beta))$
$(y-X\beta-W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}R^{-1}(y-X\beta))$
さて、今から
$W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}=(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}$
となることを示します。
$(AB)^{-1}=B^{-1}A^{-1}$ であることを用いて左辺の式変形を行うと
$(RHD)=W(W^{-T}W^{T}R^{-1}W+W^{-T}\sigma_{\gamma}^{-2})^{-1}$
すなわち
$(RHD)=W(R^{-1}W+W^{-T}\sigma_{\gamma}^{-2})^{-1}$
さらに変形すると
$(RHD)=(R^{-1}+W^{-T}W^{-1}\sigma_{\gamma}^{-2})^{-1}$
つまり
$(RHD)=(R^{-1}+(WW^{T})^{-1}\sigma_{\gamma}^{-2})^{-1}$
ここで、
$K_{v}\sigma_{v}^{2}=WW^{T}\sigma_{\gamma}^{2}$
を思い出して、代入すると、
$(RHD)=(R^{-1}+( K_{v}\sigma_{v}^{2}\sigma_{\gamma}^{-2})^{-1}\sigma_{\gamma}^{-2})^{-1}$
これを整理すると
$(RHD)=(R^{-1}+ K_{v}^{-1}\sigma_{v}^{-2})^{-1}$
これより以下の二式が同じであることが分かりました。
$(y-X\beta-(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1}(y-X\beta))$
$(y-X\beta-W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}R^{-1}(y-X\beta))$
次に前半部分
$(R^{-1}+K_{u}^{-1}\sigma_{u}^{-2}-R^{-1}(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1})^{-1}$
$(Z^{T}R^{-1}Z+I\sigma_{\mu}^{-2}-Z^{T}R^{-1}W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}(W^{T}R^{-1}Z))^{-1}Z^{T}$
極めて似ていることに着目し、両者の関係性を探っていきます。
ここでは2つ目の式
$(Z^{T}R^{-1}Z+I\sigma_{\mu}^{-2}-Z^{T}R^{-1}W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}(W^{T}R^{-1}Z))^{-1}Z^{T}$
を変形して1つ目の式に近づけていきます。
まず、末尾にくっついている $Z^{T}$ をカッコ内に吸収すると
$(R^{-1}Z+Z^{-T}\sigma_{\mu}^{-2}-R^{-1}W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}R^{-1}Z)^{-1}$
次に無理やりZをくくりだすと
$( (R^{-1}+Z^{-T}Z^{-1}\sigma_{\mu}^{-2}-R^{-1}W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}R^{-1})Z)^{-1}$
よって、
$Z^{-1} (R^{-1}+(ZZ^{T})^{-1}\sigma_{\mu}^{-2}-R^{-1}W(W^{T}R^{-1}W+I\sigma_{\gamma}^{-2})^{-1}W^{T}R^{-1})^{-1}$
先ほどと同様のカラクリを使えば
$Z^{-1} (R^{-1}+K_{u}^{-1}\sigma_{u}^{2}-R^{-1}(R^{-1}+K_{v}^{-1}\sigma_{v}^{-2})^{-1}R^{-1})^{-1}$
ここで $Z^{-1}$ 以降は1つ目の式と同じ形になっています。

以上をまとめるとuとμの間には次の関係があることが分かります。
$\mu=Z^{-1}u$
ここで、 $K_{u}\sigma_{u}^{2}=ZZ^{T}\sigma_{\mu}^{2}$ より
$Z^{-1}=Z^{T}K_{u}^{-1}\sigma_{u}^{-2}\sigma_{\mu}^{2}$
これを代入して
$\mu=(Z^{T}K_{u}^{-1}\sigma_{u}^{-2}\sigma_{\mu}^{2})u$
ここで
$\sigma_{u}^{2}=\sigma_{\mu}^{2}$ が成り立つとすると
$\mu=Z^{T}K_{u}^{-1}u$
これでマーカー効果とGBLUPの変量効果が結びつきました。
今の話はvとγにおいても同様に成り立ち、その対称性から
$\gamma=W^{T}K_{v}^{-1}v$
となることが直ちにわかります。
この結果は先週の記事の結論と同じであり、結局、マルチカーネルであろうとなかろうと対象の変量効果に注目して計算すればよいことが示されました。

2024-07-14

GBLUPの変量効果からマーカー効果を推定

線形混合モデル

はじめに

DNAマーカーの開発により遺伝子情報から表現型(作物の特性)を予想することが可能になりました。これは育種においては画期的な技術革新です。例えば果樹の品種改良を考えたとき、それまでは交配した種を実がなるまで育てる必要がありましたが、この技術を使えば小さい苗の段階でDNAを抽出、それをもとに果実形質を予測、果実形質が悪いと予測された個体を早期に淘汰することで育種にかかるコストや労力を大幅に減らす言ことが可能です。その予測に使われるのがいわゆる混合モデル(mixed model)と呼ばれるモデルです。ベースとなっているは回帰モデルで、その切片や傾きが確率分布にしたがっていると仮定します。このモデルは一般的に次のように記述されます。

$y=X\beta+Zu+e$

$u\sim MVN(0,G)$

$e\sim MVN(0,R)$

ここでyは応答変数に対応するベクトルでn×1の次元をもっています。β, uはそれぞれ固定効果、変量効果と呼ばれる変数で、前者が確率分布に従わないのに対し、後者は確率分布に従う確率変数になります。変量効果が従う確率変数は二つ目の式で規定されていて、Gは分散共分散行列と呼ばれます。それぞれr×1、m×1を持っており、yのn×1次元との調整をするために使われるのがXやZです。これらはデザイン行列と呼ばれる行列で、それぞれn×r次元、n×m次元を持っています。最後のeは誤差であり、これもまた確率分に従う確率変数と定義されます。Rはその確率分布がとる分散共分散行列ですが、一般的に単位行列の定数倍が使われます。

さて、混合モデルを使って遺伝子情報と表現型の関係をモデル化する際、uが遺伝的効果、yが表現型情報になるわけですが、このuをどう設定するかによって大きく二つに分かれます。一つ目はマーカーごとに効果を設定する方法、もう一つはマーカーの効果全体を丸っとまとめて扱う方法です。特に後者はGBLUPと呼ばれますが、このGBLUPで求まった個体が持つ遺伝的効果の総和からマーカー効果を求める方法を本日は紹介します。これ以降はそれぞれの方法の基本は知っている体で話を進めます。

GBLUPによるモデル

GBLUPによるモデルでは以下のようにモデル化されます。なお、個体反復はないものとします。

$y=X\beta+u+e$

$u\sim MVN(0,K\sigma^{2}_{u})$

$e\sim MVN(0,I\sigma^{2}_{e})$

ここでKはゲノム関係行列と呼ばれる行列でマーカー情報から計算されます。

マーカー効果を使ったモデル

マーカー効果をそれぞれモデルに組み込んだものは次のようにモデル化されます。先ほどと同様個体反復はないものとします。

$y=X\beta+Z\mu+e$

$\mu\sim MVN(0,I\sigma^{2}_{\mu})$

$e\sim MVN(0,I\sigma^{2}_{e})$

ここでμはマーカー分の長さを持つベクトルになります。

両者の関係

さて、いま二つのモデルを示しましたが、これらは

$u=Z\mu$

つまり、マーカー効果μの総和がuであると定義すると同値であることが分かります。今回やりたいことは推定されたuからμを求めることです。一見するとZの逆行列をuに書ければよいように思いますが、推定された値を使っているため実はそう簡単にはいきません。以下ではこの推定されたuからμを求める方法を混合モデル方程式と呼ばれる方法から導出します。

混合モデル方程式

さて、一番はじめに示した式に戻り一般的な解法を復習します。

$y=X\beta+Zu+e$

$u\sim MVN(0,G)$

$e\sim MVN(0,R)$

この式が成り立つとき、βやuは次の式を解くことで求められることが知られています(ここだけ黒抜きで表示されているかもしれませんが特に意味はありません)。

$\begin{bmatrix} X^{T}R^{-1}X & X^{T}R^{-1}Z \\ Z^{T}R^{-1}X & Z^{T}R^{-1}Z+G^{-1} \\ \end{bmatrix}\begin{bmatrix} \beta \\ u \\ \end{bmatrix}=\begin{bmatrix} X^{T}R^{-1}y \\ Z^{T}R^{-1} \\ \end{bmatrix}$

なお、この式は混合モデル方程式と一般的に呼ばれています。さて、この式を変形すると二つの式が得られますが、そのうち二つ目の式に注目します。

$Z^{T}R^{-1}X\beta+{Z^{T}R^{-1}Z+G^{-1}}u=Z^{T}R^{-1}y$

これをuについて解くと

$u=(Z^{T}R^{-1}Z+G^{-1})^{-1}Z^{T}R^{-1}(y-X\beta)$

ここで

$(AB)^{-1}=B^{-1}A^{-1}$

が一般的に成り立ちます。

$Z^{T}=((Z^{T})^{-1})^{-1}$

であることを用いると、

$u=(R^{-1}Z+(Z^{T})^{-1}G^{-1})^{-1}R^{-1}(y-X\beta)$

さらにRについても同様のことをすると

$u=(Z+R(Z^{T})^{-1}G^{-1})^{-1}(y-X\beta)$

さて、ここまでもとまったところで、次からは具体的にGBLUPの場合、各マーカーを計算する場合のモデルにこの式を適用していきます。

GBLUPの場合

GBLUPの場合、先ほどのモデルにおいて、

$G\to K\sigma^{2}_{u}$

$R\to I\sigma^{2}_{e}$

$Z\to I$

に変換すればいいことが分かります。従って先ほど求めた

$u=(Z+R(Z^{T})^{-1}G^{-1})^{-1}(y-X\beta)$

の式はGBLUPの場合次のように変形することができます

$u=(I+\sigma^{2}_{e}(K\sigma^{2}_{u})^{-1})^{-1}(y-X\beta)$

さらに変形すると

$u=(I+\sigma^{2}_{e}\sigma^{-2}_{u}K^{-1})^{-1}(y-X\beta)$

もう少し変形して

$u=K(K+I\sigma^{2}_{e}\sigma^{-2}_{u})^{-1}(y-X\beta)$

マーカー効果を使ったモデルの場合

マーカー効果を使った場合、先ほどのモデルにおいて、

$u \to \mu$

$G\to I\sigma^{2}_{\mu}$

$R \to I\sigma^{2}_{e}$

に変換すればよいことが分かります。よって

$\mu=(Z+\sigma^{2}_{e}(Z^{T})^{-1}\sigma^{-2}_{\mu})^{-1}(y-X\beta)$

さらにZ転置をくくりだして

$\mu=((Z Z^{T}+I\sigma^{2}_{e}\sigma^{-2}_{\mu})(Z^{T})^{-1})^{-1}(y-X\beta)$

くくりだしたものを外に出すと

$\mu=Z^{T}(Z Z^{T}+I\sigma^{2}_{e}\sigma^{-2}_{\mu})^{-1}(y-X\beta)$

さて、 $u=Z\mu$ より次が成り立ちます

$var(u)=var(\mu)=K\sigma_{u}^2=Z Z^{T}\sigma_{\mu}^2$

よって

$\mu=Z^{T}(K\sigma_{u}^2\sigma_{\mu}^{-2}+I\sigma^{2}_{e}\sigma^{-2}_{\mu})^{-1}(y-X\beta)$

もう少し整理して

$\mu=Z^{T}\sigma_{u}^2\sigma^{-2}_{\mu}(K+I\sigma^{2}_{e}\sigma^{-2}_{u})^{-1}(y-X\beta)$

最終変形

以上をまとめますと

GBLUPでは

$u=K(K+I\sigma^{2}_{e}\sigma^{-2}_{u})^{-1}(y-X\beta)$

各マーカー効果を使ったモデルでは

$\mu=Z^{T}\sigma_{u}^2\sigma^{-2}_{\mu}(K+I\sigma^{2}_{e}\sigma^{-2}_{u})^{-1}(y-X\beta)$

となることが分かりました

ここで両者に共通して以下の項が含まれています

$(K+I\sigma^{2}_{e}\sigma^{-2}_{u})^{-1}(y-X\beta)$ ]

これを消去すると

$\mu=\sigma_{u}^2\sigma^{-2}_{\mu}Z^{T}K^{-1}u$

多くの場合 $\sigma_{u}=\sigma_{\mu}$ が仮定できるので以下が導かれます

$\mu=Z^{T}K^{-1}u$

2023-10-30

非線形混合モデル

非線形混合モデル

量的遺伝学ではいわゆる線形混合モデルと呼ばれるモデルがよく使われていますが、成長曲線などをモデル化する場合は、この線形混合モデルを拡張した非線形混合モデルが必要になってきます。ちょうど研究の中で非線形混合モデルを使用しており、まとめて勉強する機会があったので、今日はひとまずそもそも非線形混合モデルがどう数式として定式化されるかを書こうと思います。

非線形であるとは

そもそも非線形とは何かをまずはじめにはっきりさせておきます。線形と聞くと"まっすぐ"、"直線"といったイメージを持ちがちなので、 $y=ax^2+bx+c$ でかかれる二次関数も非線形として捉えがちですが(そして文脈によってはそれは間違いではないのですが)、ここではパラメータに対して非線形であるモデルを非線形であると呼ぶことにします。 $y=ax^2+bx+c$ の例でいうと、この関数はパラメータも変数であると考えることによって、 $f(a,b,c,x)=ax^2+bx+c$ ととらえることができるわけですが、このa,b,cについては線形性が保たれているためこれは線形なモデルということになります。例えば $f(a,b,c,x)=a^2x^2+bx+c$ なんかは非線形なモデルというわけです。

植物分野でよく使われる非線形なモデル

非線形なモデルを調べるとよく薬学の話、特に薬の体内での濃度変化の話が取り上げられますが、植物分野では成長曲線がよく使われます。特に使われるのが、①logistic 曲線、②Gompertz 曲線、③von Bertalanffy 曲線の3つです。例えば前者二つはそれぞれ次のように数式化されます。

logistic 曲線

$f(a,b,c,t)=\frac{a}{1+b*exp(-ct)}$

Gompertz 曲線

$f(a,b,c,t)=a*b^{exp(-ct)}$

非線形なモデルの優位性

このようにxなどの変数に対しての非線形性は認められるため、線形なモデルでもかなり柔軟なモデルにすることができるわけですがそれでも非線形のモデルを使いたい理由はどこになるのでしょうか。よく取り上げられるのが次の3つです。①パラメータの解釈可能性。多くの非線形なモデルは微分方程式に基づいて導出されます(Gompertzモデルといった例外ももちろんありますが…)。微分方程式は物理量の関係性を表すモデルですので、そこに現れるパラメータには何らかの物理的解釈を加えることが可能です。線形なモデルではパラメータの解釈が難しく考察や問題解決がしにくいわけですが、非線形なモデルではそういうことが起こらないというメリットがあります。②倹約性。一般に同じ柔軟性を持つモデルを比較すると、非線形なモデルのパラメータ数は線形なもでるのそれよりも少ない傾向にあります。パラメータ数が増えると多重線形性といった問題が生じてくるのですが、非線形なモデルはそういった問題をあまり気にしなくて済みます。③外挿に対する安定性。機械学習では、学習データに含まれない領域のデータに対する予測精度が低いという問題がよく取り上げられます。非線形なモデルは先ほど書いたように微分方程式に基づき、それはどのデータ領域にも適応できるため、そのようなデータに対しても予測精度が下がりにくいというメリットがあります。

時系列データの非線形モデルによる定式化

さて、ここからは時系列データが非線形なモデルを使ってどのように定式化できるかを考えていきます。まず、品種がn品種あり、それぞれに対して時系列データがとられたという状況を想定します。この時、i番目の品種のj回目の測定データは $y_{ij}$ と表すことができます。ここで $i=1\sim n$ 、 $j=1\sim n_{i}$ であり、 $n_{i}$ はi番目の品種に対して測定された回数を表します。次に非線形なモデル $f(a,b,c,t)$ ですが、パラメータは品種によって違うと考えられるためi番目の品種に対する非線形なモデルは $f(a_{i},b_{i},c_{i},t)$ となります。さらに、 $\phi_{i}=(a_{i},b_{i},c_{i})$ とすると、 $f(\phi_{i},t)$ と書くことができます。測定された時刻を $t_{ij}$ と表すことにすると、 $y_{ij}$ は次のようになります：

$y_{ij}=f(\phi_{i},t_{ij})+e_{ij}$

ここで $e_{ij}$ はノイズで $e_{ij}\sim N(0,\sigma)$ に従います。これがまず第一段階の定式化です。しかし、これはあくまでも非線形モデルであって、非線形混合モデルではありません。次にこれを非線形混合モデルに拡張していきます。

非線形混合モデル

先ほどの式を改めて書くと次のようになります

$y_{ij}=f(\phi_{i},t_{ij})+e_{ij}$

where $e_{ij}\sim N(0,\sigma)$

これを混合モデル化するのは簡単で、パラメータ $\phi_{i}$ が母数効果と変量効果の線形和になっている、つまり次のように表せると仮定すればよいわけです。

$\phi_{i}=X_{i}\beta+Z_{i} u_{i}$

where $u_{i}\sim MVN(0,G)$

ここで、 $\beta$ は固定効果を、 $u_{i}$ は変量効果を表しています。イメージとしては前者が品種を平均した値で、後者が品種による平均からのずれといったところでしょうか。また変量効果が従う多変量正規分布の分散共分散行列 $G$ はパラメータが3つの場合は $(a_{i},b_{i},c_{i})$ 間の共分散を表しています。

よって、まとめると非線形混合モデルは次のように表されるモデルとなります。

$y_{ij}=f(\phi_{i},t_{ij})+e_{ij}$

where $\phi_{i}=X_{i}\beta+Z_{i} u_{i}$ 、 $u_{i}\sim MVN(0,G)$ 、 $e_{ij}\sim N(0,\sigma)$

実データにこれらのモデルをフィッテイングしていく際は、非線形性ならではの困難が伴うのですが、それはまた次回。やる気が出たときに書こうと思います。

2023-03-19

ハーディーワインベルグの法則

遺伝学の基礎

ハーディーワインベルグの法則は量的遺伝学において最も重要な法則の一つで、無作為に交配が起こっていればその集団の遺伝的構成はどの世代でも変わらないことを主張します。そのように構成が不変な状態は平衡状態とみなすことができるため、ハーディーワインベルグ平衡と言ったりもします。主張の詳しい説明は後にしますが、育種では無作為な交配と作為的な交配を繰り返すため、ハーディーワインベルグの法則が成り立つ状況なのかどうかを見極めることが非常に重要です。

さて、以下ではハーディーワインベルグの法則をマメの集団を例に説明していきます。このマメには丸型としわ型があり、それぞれ $A$ 、 $a$ という対立遺伝子に支配されているとします。また、 $A$ は $a$ に対して優性(顕性)であるとします。

言葉の定義

遺伝子型頻度(genomic frequency)

集団における $AA$ 、 $Aa$ 、 $aa$ を持つ個体の頻度。ここではそれぞれ $P$ 、 $Q$ 、 $R$ で表す。この時、 $P+Q+R=1$ を満たすことが直ちにわかる。

遺伝子頻度

集団における $A$ 、 $a$ の頻度。ここではそれぞれ $p$ 、 $q$ で表す。同様に、 $p+q=1$ を満たす。

遺伝子型頻度と遺伝子頻度の関係

ここで、遺伝子型頻度と遺伝子頻度の間には次の関係式が成り立つことが分かる。

$p=P+\frac{1}{2}Q$

$q=R+\frac{1}{2}Q$

例

ある1000個体のマメ集団において、 $AA$ 、 $Aa$ 、 $aa$ の遺伝子型を持つ個体がそれぞれ、500、100、400個体いた場合、遺伝子型頻度は次のように求まる。

$P=\frac{500}{1000}=0.50$

$Q=\frac{100}{1000}=0.10$

$R=\frac{400}{1000}=0.40$

また、遺伝子型頻度と遺伝子頻度の関係より、

$p=0.50+0.05=0.55$

$q=0.40+0.05=0.45$

さて、次はハーディーワインベルグの法則の主張を詳しく見ていきます。

主張

ハーディーワインベルグの法則はある条件を満たす集団における遺伝的構成に関する主張ですが、その"ある条件"が非常に大切なので先にこちらを見ていきます。

①無作為交配が起きていること

→完全にランダムに雄と雌が選ばれ交配されている状況です。自家受粉をする植物は、同じ個体に由来する雄(花粉)と雌(胚珠)の交配が優先されるため無作為交配にはなりません。

②考えている集団サイズが大きいこと

→AB型同士の親から連続してA型の子供が生まれるように、考えている集団サイズが小さい場合、たまたまある遺伝子型が偏って生じることが考えられます。このような効果を減らすため、集団サイズは大きい必要があります。

③他の集団からの流入、への流出がないこと

丸のマメしか存在しない集団から種が流入するようなことがあると必然的に丸のマメの頻度が多くなってしまいます。このような影響を考えないため集団は他の集団から孤立していることが求められます。

④突然変異が起こらないこと

丸としわのマメの集団の中に突然三角形のマメが生じるようなことがあると、集団の遺伝的構造が大きく変わるため、このような変異はないものと考えます。

⑤遺伝子型に有利不利がないこと

例えば考えている遺伝子がマメの丸・しわではなく、花の色などの場合どちらかが受粉に有利に働く場合があります。この場合、交配のランダム性がなくなってしまいます。

そして、この条件を満たす集団において次のことを主張します。

①遺伝子型頻度と遺伝子頻度は変化しない(ハーディーワインベルグ平衡)。

②ハーディーワインベルグ平衡には一世代で至ることができる。

ここで、①についてはそのままですが②は少しわかりにくいかもしれませんので少し説明を加えます。例えば、先の条件を満たしていない集団があったとします。この時、ハーディーワインベルグの法則は成り立たないため、遺伝子型頻度と遺伝子頻度は世代ごとに変化します。ところが、その集団が突然条件を満たすようになった場合、その次の世代以降はハーディーワインベルグ平衡に至り、遺伝子型頻度と遺伝子頻度が一定になる、という主張です。

証明

まず①の証明を行います。

遺伝子型頻度が $P_1$ 、 $Q_1$ 、 $R_1$ 、遺伝子頻度が $p_1$ 、 $q_1$ の集団を想定し、その次の世代の集団における遺伝子型頻度 $P_2$ 、 $Q_2$ 、 $R_2$ と遺伝子頻度 $p_2$ 、 $q_2$ を考えます。

遺伝子型頻度

例えば、 $Aa$ の個体は、オスから $A$ とメスから $a$ 、またはオスから $a$ とメスから $A$ をもらうことによって誕生しますが、無作為なが交配が起こっているため、オスから $A$ をもらう確率は遺伝子頻度 $p_1$ 、メスから $a$ をもらう確率は $q_1$ となります。逆も同様であり、したがって $Aa$ が誕生する確率は $2p_1q_1$ となるわけです。 $AA$ 、 $aa$ についても同様で $p_1^2$ 、 $q_1^2$ の確率で誕生します。集団に対する細かい条件が付いていた理由はまさにここにあり、 $A$ と $a$ が同様に確からしいことを保証していたわけです。

ここで、次の世代の集団が $N$ 個体から成り立っているとします。 $AA$ 、 $Aa$ 、 $aa$ は、それぞれ $p_1^2$ 、 $2p_1q_1$ 、 $q_1^2$ の確率で生じるので、次の世代の集団の分布は多項分布に従うことになります。つまり、全 $N$ 個体において $AA$ 、 $Aa$ 、 $aa$ がそれぞれ $n_{AA}$ 、 $n_{Aa}$ 、 $n_{aa}$ 個体生じる確率 $P(n_{AA},n_{Aa},n_{aa})$ は、

$P(n_{AA},n_{Aa},n_{aa})=\frac{N!}{n_{AA}!n_{Aa}!n_{aa}!}(p^2)^{n_{AA}}(2pq)^{n_{Aa}}(q^2)^{n_{aa}}$

となるわけです。ここで多項分布の理論より、 $AA$ が出る頻度 $frac{n_{AA}}{N}$ の期待値は $p_1^2$ 、 $Aa$ が出る頻度 $frac{n_{Aa}}{N}$ の期待値は $2p_1q_1$ 、 $aa$ が出る頻度 $frac{n_{aa}}{N}$ の期待値は $q_1^2$ となることがわかるので、

$P_2=p_1^2$

$Q_2=2p_1q_1$

$R_2=q_1^2$

が期待値の意味で成り立つことがわかります。

遺伝子頻度

次に遺伝子頻度ですが、 $p+q=1$ を使うと遺伝子型頻度と遺伝子頻度の関係より、

$p_2=p_1^2+2p_1q_1=p_1(p_1+q_1)=p_1$

$q_2=q_1^2+2p_1q_1=q_1(q_1+p_1)=q_1$

が成り立ちます。これは前の世代の集団と同じ値です。つまり、最初の集団の遺伝子頻度が $p_1$ 、 $q_1$ であれば、それ以降の世代の遺伝子型頻度は $P_n=p_1^2$ 、 $Q_n=2p_1q_1$ 、 $R_n=q_1^2$ 、遺伝子頻度は $p_n=p_1$ 、 $q_n=q_1$ となるわけです。これがハーディーワインベルグの法則の①の証明です。注目したいのは $P_1$ 、 $Q_1$ 、 $R_1$ は一切計算に出てこなかった点です。これが②の主張につながります。

②の証明

最初の世代を色々な場所からサンプリングして集めてきた集団だとしましょう。その集団を一か所に集めて無作為に交配を行ったとします。この時、先ほどの議論より2世代目以降の遺伝子型は $P_n=p_1^2$ 、 $Q_n=2p_1q_1$ 、 $R_n=q_1^2$ となりますが、 $P_1$ 、 $Q_1$ 、 $R_1$ はランダムに集められた集団のため、 $P_n=P_1$ 、 $Q_n=Q_1$ 、 $R_n=R_1$ は保証されません。逆にいえば、一世代分過ぎれば集団はハーディーワインベルグ平衡に至るということです。これは②の主張にほかなりません。

以上のことより、集団に対するある条件のもとハーディーワインベルグの法則が成り立つことがわかります。

応用

以下ではハーディーワインベルグの法則がどのようなところに応用されるかについて、ヘテロ接合個体の個体数推定を取り上げようと思います。

ヘテロ接合個体の個体数推定

遺伝子型頻度と遺伝子頻度の例では、 $AA$ 、 $Aa$ 、 $aa$ の遺伝子型を持つ個体をそれぞれ500、100、400個体としましたが、通常 $AA$ と $Aa$ は区別することができないため $Aa$ の個体数を知ることはできません。しかし、考えている集団がハーディーワインベルグ平衡にあるのであれば、そこから $Aa$ の個体数を導くことが可能です。

まず、全1000個体を調査したうち、40個体が $aa$ だったとします。この時 $R=\frac{40}{1000}=0.04$ となるわけですが、ハーディーワインベルグ平衡においては $R_n=q^2$ が保証されます。これより $q=0.2$ となりますが、 $p+q=1$ より、[p=0.8]がわかります。さらに $Q_n=2pq$ であるため、 $Q=0.16$ と $Aa$ の頻度がわかります。この頻度に1000をかけた160個体が $Aa$ の個体数となるわけです。このようにして、 $Aa$ の個体数を推測することが可能です。