Rでの実践編〜線形混合モデル(REML) - ryamadaの遺伝学・遺伝統計学メモ

linear mixed modelのためのRパッケージはいくつかあるようで、lme4パッケージなどが有名なようだが
linear mixed modelの頻繁な利用対象は、サンプルをいくつかの群にわけてそれぞれでことなるバラツキを想定したりするものらしく、すべてのサンプルに関する分散共分散行列を指定する方法は、「理論としてはカバー」していても、「実装が対応していない」ものがあり、lme4もその一つらしい
もともと、分散共分散行列で変量効果を扱うことを前提にしているパッケージがvarCompパッケージ
その名も、Variance Component法(この名前は、考え方によれば線形混合モデルであるし、被説明変数の分散を指定した構造の分散に分解する方法という意味で『分散分解法』と言ってもよいことからついている)
肝となる関数はvarComp()関数
これの作りは
- $y = X\beta e$ と線形回帰する
- ただし、 $e$ の項は、 $V$ なる分散共分散に従う多次元正規乱数項
- ここで $V=\sum_{j=1}^R \sigma_j^2 K_j + \sigma_e^2 W$ と分解される
- $K_j$ は半正定値行列、 $W$ は正の対角行列
- さらに $K_j = Z_j G_j Z_j^T$ と分解し(Zは対角行列、Gは分散共分散行列)
- varComp()関数では、ZとGの組(いくつも指定できる)とWの対角成分を指定して、変量効果の形を決められる。それらの対応する $\sigma_j,\sigma_e$ を推定対象とする
- そしてもちろん、 $\beta$ により固定項も推定する
実際、SNPによる相加モデルでの遺伝率推定にあたっては、Zは単位行列、Wも単位行列で、Gがゲノムワイドな個人の分散共分散行列、変量効果はこれだけでR=1。また、Xは既知リスクマーカージェノタイプ行列でその線形寄与係数が $\beta$
実際に簡単なデータを作って回してみる

library(varComp)
n <- 100
y <- rnorm(100)
n.marker <- 5
x <- matrix(rnorm(n*n.marker),ncol=n.marker)
x2 <- matrix(rnorm(n*n.marker),ncol=n.marker)
G <- cov(t(x2))
dim(G)
Z <- diag(rep(1,n))
W <- rep(1,n)
varComp.out <- varComp(fixed = y ~ x, random = ~ Z, varcov = G, weights = W)