読者です 読者をやめる 読者になる 読者になる

疫学と医療統計学と遺伝学と時々、大学院生活

疫学を専門とする大学院生の研究に関する備忘録的ページ。

R package "BLR"でベイジアンラッソ(Bayesian Lasso regression)

Genetic epidemiology GWAS R programming epidemiology statistics

あけましておめでとうございます(疫学会や予防早期医療創成センターのワークショップ等への参加もあり、一ヶ月ぶりの更新です)。
今回はRのパッケージ"BLR"を用いて、ベイジアンラッソを実行する。ちなみに参考にするのは、Cedric Gondroらによる『Genome Wide Association Study and Genomic Prediction』のP.299からである。BLR内のデータセットwheatを用いて、10-fold cross-validationして観測値と予測値の相関係数を求める。

実際に始める

wheatの中に入っているデータはこちらを参照。

library(BLR)
data(wheat)
set.seed(12345)
#表現型を選択する
y <- Y[,1]

#パラメータの設定
prior <- list(varE=list(df=5, S=1.5),
           lambda=list(type="random", value=20,
                       shape=2, rate=2e-4))

#クロスバリデーションの設定
cvCor <- numeric()
folds <- rep(1:10, 60)[order(runif(599))] #10回のCVを想定する

#実際に10-folds CVを実行する
for(i in 1:10)
{
  tst <- which(folds==i)
  yNA <- y; yNA[tst] <- NA
  fm <- BLR(y=yNA, XL=X, nIter=6000, burnIn=1000, prior=prior)
  cvCor[i] <- cor(y[tst], fm$yHat[tst])
}

#最後に10回分の実測値と予測値の相関を表示する
cvCor

結果(全体的に中程度かな、、、)

0.6299844 0.4754850 0.5333422 0.2822434 0.5821051 
0.5377218 0.5674720 0.6276386 0.5019760 0.5267373

振り返り

肝はCVみたいになっているコードですが、重要なのは、forのループ内にあるBLR関数です。BLR関数のインプットには、y(表現型)、XL(当てはめのための行列)、nIter(反復数)などが必要です。詳しくはBLRパッケージのpdfをご参照ください。

参考にした書籍

www.springer.com
f:id:ryosukefujii0320:20170202234346p:plain



20170202
RF