はじめに

GBLUPは混合モデルの一種であるが、変量効果を含むモデルはそれ自体を潜在変数と扱うことができるためEMアルゴリズムと相性が良い。先日GBLUPにおけるEMアルゴリズムの導出を行ったのでここでその実装を行う。詳しい計算については混合モデルをEMアルゴリズムで解く&Rで実装する①を参考にせよ。

パッケージ読み込み

今回はゲノムデータの読み込みにgastonを多変量正規分布の生成にmvnfastを利用する。

library( gaston )

##  要求されたパッケージ Rcpp をロード中です

##  要求されたパッケージ RcppParallel をロード中です

## 
##  次のパッケージを付け加えます: 'RcppParallel'

##  以下のオブジェクトは 'package:Rcpp' からマスクされています:
## 
##     LdFlags

## Gaston set number of threads to 4. Use setThreadOptions() to modify this.

## 
##  次のパッケージを付け加えます: 'gaston'

##  以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     sigma

##  以下のオブジェクトは 'package:base' からマスクされています:
## 
##     cbind, rbind

library( mvnfast )

## Warning: パッケージ 'mvnfast' はバージョン 4.3.2 の R の下で造られました

表現型生成

実データを用いても良いが推定したパラメータが真値とどれだけ近いかを評価できないため、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" )

## ped stats and snps stats have been set. 
## 'p' has been set. 
## 'mu' and 'sigma' have been set.

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^2 )

# 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 )

## [1] 0.5680535

EMアルゴリズム

次にEMアルゴリズムを用いてパラメータ推定を行う。EMアルゴリズムは局所解に陥りやすいため、ここでは初期値としてかなり真値に近い値を設定した。詳しい計算は混合モデルをEMアルゴリズムで解く&Rで実装する①にまとめたのでここでの説明は省くがパラメータを更新しながら真のパラメータを発見するようなアルゴリズムになっている。

# set initial values
b <- 5
var_g <- 6
var_e <- 3

# vector to store the transition of estimated parameters
b_vec <- b
var_g_vec <- var_g
var_e_vec <- var_e

# 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 <- 1e-10 # the criteria to decide whether values have converged
rep <- 1

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
  
  b_vec <- c( b_vec, b )
  var_g_vec <- c( var_g_vec, var_g ) 
  var_e_vec <- c( var_e_vec, var_e ) 
  
  rep <- rep + 1
  
}

print(rep)

## [1] 1879

真値との比較

最後に真値との比較を行った。\(b\)の推定は比較的よくできている一方、分散についてはかなり推定精度が悪かった。EMアルゴリズムの計算に間違っている可能性も踏まえ再度検討が必要か。

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

## [1] "真値 5 推定された値 4.87176626470429"                 
## [2] "真値 4 推定された値 1.40651697220631"                 
## [3] "真値 0.322684829667875 推定された値 0.764691400759617"

値の遷移

plot( b_vec, ylab = "b" )

plot( var_g_vec, ylab = "var_g" )

plot( var_e_vec, ylab = "var_e" )

EMアルゴリズムによるGBLUPの実装

木村奏

2024-12-29