はじめに

機械学習の勉強方法は色々あると思いますが、私はやっぱり勉強した内容をコードで実装するまで、ちゃんと理解した感じがしないので、本を読んだ後はなるべく実装するようにしています。機械学習におけるパラメータ推定で必須のMCMCについては、メトロポリスヘイスティングは実装したことがありましたが、ギブスサンプリングは今までやったことがありませんでした。
私が所属する研究室では、遺伝子情報を用いた解析をしいていますが、そこでは線形混合モデルと呼ばれるモデルがよく出てきます。そこで今回は、この線形混合モデルに出てくるパラメータの推定をギブスサンプリングを通して行っていきたいと思います。
本投稿は、大きく二つに分けて分かれています。一つ目が線形混合モデルのパラメータに対する完全事後分布を計算する理論編。もう一つが計算された完全事後分布を元にRで実装する実装編です。今回は後者の記事にあたります。前者の記事も合わせてお読みいただければ幸いです。
線形混合モデルのギブスサンプリングによるパラメータ推定① 理論編

目次

・データ生成
・ギブスサンプリングによるパラメータ推定

データ生成

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

library(gaston)
library(mvnfast)

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

# generate Genomic relationship matrix ( GRM )
Genome <- read.vcf( "Rice_Data/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 )

今回のデータの場合$\sigma_e$は次のような値になりました。

sigma <- sqrt( var )
print( sigma )
[1] 0.5306414

また、$\mathbf{u}$の真値と$\mathbf{y}$を対応付けてプロットすると以下のようになりました。遺伝率が高いので相関もかなり高めです。

plot( u, y , pch = 19, xlab = "genetic value", ylab = "phenotype" )

ギブスサンプリングによるパラメータ推定

上で生成されたデータを用いて各変数の推定を実際に行っていきます。とはいっても前の記事の最後の式をRで書き直しているだけなので、特に面白みはありません。サンプリングの部分はもう少しちゃんと考えるべきなんでしょうが、今回はburnoutだけを導入し、最初の100サンプリングは取り除くようにしました。

# initial values
init_b <- 0
init_u <- rep( 0, 100 )
init_prec_g <- 1
init_prec_e <- 1

# insert initial values to vectors or matrix to store samplings
b_vec <- matrix( init_b )
u_vec <- matrix( init_u )
prec_g_vec <- init_prec_g
prec_e_vec <- init_prec_e

# definition of N and n
N <- length( y )
n <- n

# scale and rate parameters for gamma distribution of precision g and e
a_g <- 1
b_g <- 1
a_e <- 1
b_e <- 1

# burnout and number of sampling
burnout <- 100
n_sample <- 10000

# difinition of X and Z
X <- matrix( rep( 1, 100 ) )
Z <- diag( 1, 100 )

for( i in 2:n_sample ){
  print( i )
  
  # sampling of b
  b_mean <- prec_e_vec[ i-1 ] * solve( t( X ) %*% X ) %*% t( X ) %*% ( y - Z %*% u_vec[ , i-1 ] )
  b_var <- solve( t( X ) %*% X * prec_e_vec[ i-1 ] )
  b_vec <- cbind( b_vec, mvnfast::rmvn( 1, b_mean, b_var ) )
  
  # sampling of u
  u_mean <- ( solve( ( t(Z) %*% Z * prec_e_vec[ i-1 ] + solve( GRM ) * prec_g_vec[ i-1 ] ) ) *
    prec_e_vec[ i-1 ] * t( Z ) ) %*% ( y - X %*% b_vec[ , i ] )
  u_var <- solve( ( t(Z) %*% Z * prec_e_vec[ i-1 ] + solve( GRM ) * prec_g_vec[ i-1 ] ) )
  u_vec <- cbind( u_vec, t( mvnfast::rmvn( 1, u_mean, u_var) ) )
  
  # sampling of precision of g
  prec_g_a <- N/2 + a_g
  prec_g_b <- 1/2 * t(y - Z %*% u_vec[ , i ] - X %*% b_vec[ , i ] ) %*% (y - Z %*% u_vec[ , i ] - X %*% b_vec[ , i ] )+ b_g
  prec_g_vec <- c( prec_g_vec, rgamma( 1, shape = prec_g_a, rate = prec_g_b) )
  
  # sampling of precision of e
  prec_e_a <- n/2 + a_e
  prec_e_b <- 1/2 * t( u_vec[ , i ] %*% solve( GRM )% * % u_vec[ , i ] ) + b_e
  prec_e_vec <- c( prec_e_vec, rgamma( 1, shape = prec_e_a, rate = prec_e_b ) )
  
}

b_estimated <- mean( b_vec[ 1, - ( 1:burnout ) ] )
u_estimated <- apply( u_vec[ ,- ( 1:burnout ) ], 1, mean )
prec_g_estimated <- mean( prec_g_vec[ - ( 1:burnout ) ] )
prec_e_estimated <- mean( prec_e_vec[ - ( 1:burnout) ] )

ということで各変数の推定が終わりました。最後に$\mathbf{u}$の真値と推定された$\mathbf{u}$の値を比較してみます。

plot( u, u_estimated, pch= 19 )

想像以上に良い推定をしてくれました($\mathbf{u}$を生成する際に平均値を加えているのでx軸の平均が0から大きく逸脱していますが、そんな大した問題ではありません)。
ギブスサンプリングは少し難しい印象を持っていましたが、実際やってみると思いのほかシンプルで、とっつきやすいことが分かりました。今後は自作のモデルなんかで推定をしてければと思います。