はじめに

社会学や生物学では線形混合モデルと呼ばれるモデルが広く使われています。この線形混合モデルとは簡単にいうと一般的な線形モデルの拡張モデルであり、変量効果と呼ばれるパラメータを持つことが特徴としてあげられます(対する通常のパラメータは固定効果と呼ばれる)。変量効果とは簡単にいえば確率変数に従い生成されるパラメータであり、本当の意味でのパラメータではありません。真の意味のパラメータは変量効果が従う確率変数のパラメータであり、パラメータを規定するパラメータという解層構造を取っています。このパラメータを推定する方法はいくつかありますが、このような階層構造を取るが故に、EMアルゴリズムと相性が良いことが知られています。そこで今回は混合モデルの中でも生物学(育種)でよく使われているGBLUPにターゲットを絞り、EMアルゴリズムでパラメータを推定してみようと思います。前回はアルゴリズムの具体形を計算しました(混合モデルをEMアルゴリズムで解く&Rで実装する①)。今回はその計算をもとにGBLUPのEMアルゴリズムをRで実装します。
なお、私の所属している研究室ではprmlの輪読および実装を行っています。Qiitaにもまとめていますのでもしよろしければ参考にしてください(PRML　演習問題　解答集　まとめ)

表現型生成
EMアルゴリズム

表現型生成

実データを用いても良いのですが、推定したパラメータが真値とどれだけ近いかを評価できないため、McCouchら（2016; Nature Communication7:10532）が用いたイネのゲノムデータを用いて表現型生成を行うことにしました。パラメータ$b=5, \sigma_g=2$とし、$\sigma$は遺伝率が0.8になるように決定しています。

n <- 100 # number of individuals
b <- 5 # global mean
sigma_g <- 2 # standard deviation for genome

# generate Genomic relationship matrix ( GRM )
Genome <- read.vcf( "HDRA-G6-4-RDP1-RDP2-NIAS.AGCT_MAF005MIS001NR_imputed.vcf.gz" )
G_mat <- as.matrix( Genome ) - 1
G_mat <- G_mat[ 1:n, ]
GRM <- G_mat %*% t( G_mat ) / ncol( G_mat )

# generate genetic value from GRM
u <- rmvn( 1, rep( b, n ), GRM * sigma_g )

# calculate the variance for residual so that it satisfies the heritability.
var_g <- var( as.vector( u ) )
var <- ( 1 - 0.8 ) / 0.8 * var_g

# generate residuals and phenotype
e <- rmvn( 1, rep(0, n ), diag( var, n ) )
y <- as.vector( u + e )

# demonstrate the standard deviation for residuals
sigma <- sqrt( var )
print( sigma )
plot( u, y , pch = 19, xlab = "genetic value", ylab = "phenotype" )

結果、誤差に対する標準偏差は0.484827となり、表現型と遺伝子型の関係は以下のようにプロットされました。

EMアルゴリズム

次にEMアルゴリズムを用いてパラメータ推定を行います。EMアルゴリズムは局所解に陥りやすいため、ここでは初期値としてかなり真値に近い値を設定しました。この計算で出てくる式の詳細は混合モデルをEMアルゴリズムで解く&Rで実装する①にまとめてありますのでどうぞご覧下さい。

# set initial values
b <- 3
var_g <- 3
var_e <- 1

# prepare a set of useful values
X <- matrix( rep( 1, n ) )
XTX <- t( X ) %*% X
XTX_inv <- solve( XTX )
GRM_inv <- solve( GRM )

i <- TRUE # Turns to FALSE if estimated parameters are converged
conv <- 0.001 # the criteria to decide whether values have converged

while( i ){
  
  # prepare a set of useful values
  var_g_inv <- as.numeric( var_g^( - 1 ) )
  var_e_inv <- as.numeric( var_e^( - 1 ) )
  
  # prepare a set of useful values
  A_inv <- ( diag( var_e_inv, n ) + GRM_inv * var_g_inv )
  A <- solve( A_inv )
  m <- var_e_inv * ( y - X %*% b )
  Am <- A %*% m
  
  # calculate updated value of the parameters
  b_new <- XTX_inv %*% t( X ) %*% ( y - Am ) 
  var_g_new <- 1 / n * ( sum ( diag( GRM_inv %*% A ) ) + t( Am ) %*% GRM_inv %*% Am )  
  var_e_new <- 1 / n * ( t( y - X %*% b_new ) %*% ( y - X %*% b_new ) + sum ( diag( GRM_inv %*% A ) ) + t( Am ) %*% GRM_inv %*% Am - 2 * t( y - X %*% b_new  ) %*% Am )
  
  # examine whether values have converged
  i <- ( abs( b - b_new ) > conv | abs( var_g - var_g_new ) > conv | abs( var_e - var_e_new ) > conv )
  
  # update the estimated parameters
  b <- b_new
  var_e <- var_e_new
  var_g <- var_g_new
  
}

最後に推定された値と真値との比較を行いました。

print( paste( "真値", c( 5, 4, sigma^2 ), "推定された値", c( b, var_g, var_e ) ) )

[1] "真値 5 推定された値 5.17462034596268"
[2] "真値 4 推定された値 0.825535441775579"
[3] "真値 0.212147208607806 推定された値 0.782133880676785"

bの値はよく推定されていますが他のパラメータに関しては壊滅的ですね。計算のどこかにミスがあったのか…今一度確認が必要ですね