Rで解析：ECMEアルゴリズムで線形混合モデルを高速に推定する「lmm」パッケージの紹介

医療や生物統計学の分野では、多くの因子が複雑に絡み合って生じるデータを分析することがあります。特に、個体間や測定時点間の変動を考慮したモデルを構築する必要があります。例えば、臨床試験での治療効果評価や遺伝子発現データの解析などです。このような作業では、線形混合モデルの適用が有効な手法とされることがあります。

「lmm」パッケージは、ECME（Expectation/Conditional Maximization Either）アルゴリズムによる最尤推定・制限付き最尤推定に加え、より高速に収束するアルゴリズム、そしてギブスサンプラーやMCMCによるベイズ推定を実装しています。特に、Schafer (1998) “Some improved procedures for linear mixed models” で示された手法を基盤としています。

同じ引数の枠組みの中で、最尤推定・制限付き最尤推定・ベイズ推定を目的に応じて切り替えられる点が、このパッケージの特長です。

パッケージバージョンは1.4。Windows 11 x64 (build 26200)のR version 4.6.0で確認しています。

＜おすすめのRに関する書籍です＞

遺伝統計学の基礎 Rによる遺伝因子解析・遺伝子機能解析

遺伝因子解析・遺伝子機能解析をRで行うための教科書。ゲノム解析に特有の統計的な考え方を学べる。

www.amazon.co.jp

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("lmm")
# パッケージの読み込み
library("lmm")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

ECME＝Expectation/Conditional Maximization Either。EMアルゴリズムを拡張した反復計算法で、Eステップの後に、パラメータの一部を固定した条件付き最大化（CM）ステップと通常の最大化（M）ステップを組み合わせて尤度関数を最大化します。ML＝最尤推定（maximum likelihood）。RML＝制限付き最尤推定（restricted maximum likelihood）。

marijuanaデータを用いて、最尤推定・制限付き最尤推定・ベイズ推定を一通り試す例です。

marijuanaデータは、9人の男性被験者に低用量・高用量・プラセボの3種類のマリファナ入りシガレットを投与した臨床試験のデータです（3×3のラテン方格デザイン）。投与15分後・90分後の心拍数の変化が記録されており、54個のデータのうち5個が欠測しています。

ecmeml()・fastml()・ecmerml()・fastrml()・mgibbs()・fastmcmc()に共通するオプション:

コマンド別に追加される主なオプション:

# marijuanaデータセットを読み込む
data(marijuana)

# 欠測値（NA）を除いた完全データのみを使用する
marijuana <- subset(marijuana,!is.na(y))

# 列名で変数を直接参照できるようにデータフレームをアタッチする
attach(marijuana)

# 固定効果・変量効果の説明変数となる計画行列を作成する
pred <- cbind(int,dummy1,dummy2,dummy3,dummy4,dummy5)

# Xi（固定効果）に使う列番号を指定する
xcol <- 1:6

# Zi（変量効果）に使う列番号を指定する（ここでは切片のみ）
zcol <- 1

# ECMEアルゴリズムで最尤推定を行う（212サイクルで収束）
ecmeml.result <- ecmeml(y,subj,pred,xcol,zcol)

# 高速アルゴリズムで最尤推定を行う（8サイクルで収束）
fastml.result <- fastml(y,subj,pred,xcol,zcol)

# 固定効果の推定値（切片＝最終時点の平均、残りは各時点との差）を取り出す
beta.hat <- fastml.result$beta

# 6時点分の推定平均値を求める
muhat <- c(beta.hat[2]+beta.hat[1], beta.hat[3]+beta.hat[1],
   beta.hat[4]+beta.hat[1], beta.hat[5]+beta.hat[1],
   beta.hat[6]+beta.hat[1], beta.hat[1]) 

# ECMEアルゴリズムで制限付き最尤推定を行う
ecmerml.result <- ecmerml(y,subj,pred,xcol,zcol)

# 高速アルゴリズムで制限付き最尤推定を行う
fastrml.result <- fastrml(y,subj,pred,xcol,zcol)

# 変量効果の経験ベイズ推定値を取り出す
b.hat <- as.vector(fastrml.result$b.hat)

# 分散パラメータの不確実性を考慮した新しい標準誤差を求める
se.new <- sqrt(as.vector(fastrml.result$cov.b.new))

# 従来の（分散パラメータを固定した）標準誤差を求める
se.old <- sqrt(as.vector(fastrml.result$cov.b))

# 従来法と改良法の95%区間を比較する表を作成する（Table 2を再現）
table2 <- cbind(round(b.hat,3),round(cbind(b.hat-2*se.old,b.hat+2*se.old,
      b.hat-2*se.new,b.hat+2*se.new),2),round(100*(se.new-se.old)/se.old))
dimnames(table2) <- list(paste("Subject",format(1:9)),
   c("Est.","Lower.old","Upper.old","Lower.new","Upper.new","Increase (%)"))
print(table2)

# 事前分布のハイパーパラメータを指定する
prior <- list(a=3*100,b=3,c=3,Dinv=3*5)

# 修正Gibbsサンプラーでベイズ推定を行う（5000サイクル）
gibbs.result <- mgibbs(y,subj,pred,xcol,zcol,prior=prior,seed=1234,iter=5000)

# 高速MCMCアルゴリズムでベイズ推定を行う（5000サイクル）
fmcmc.result <- fastmcmc(y,subj,pred,xcol,zcol,prior=prior,seed=2345,iter=5000)

# tsパッケージのacf()を使う場合はlibrary(ts)を読み込む（base Rのacf()でも動作する）
# 2つのアルゴリズムで得られたpsiの自己相関を比較するグラフを描画する
par(mfrow=c(2,1))
acf(log(gibbs.result$psi.series[1,1,]),lag.max=10, ylim=0:1)
acf(log(fmcmc.result$psi.series[1,1,]),lag.max=10, ylim=0:1)

# アタッチしたデータフレームをデタッチする
detach(marijuana)

実行結果：

Performing ECME...
Performing FAST-ML...
Performing ECME for RML estimation...
Performing FAST-RML...
            Est. Lower.old Upper.old Lower.new Upper.new Increase (%)
Subject 1  0.190     -2.39      2.77     -3.21      3.59           32
Subject 2 -0.168     -2.75      2.41     -3.43      3.09           26
Subject 3  0.011     -2.57      2.59     -2.59      2.61            1
Subject 4  0.215     -2.40      2.83     -3.46      3.89           40
Subject 5 -0.459     -3.06      2.14     -6.50      5.58          133
Subject 6 -0.288     -2.87      2.29     -4.54      3.96           65
Subject 7  0.668     -1.91      3.25     -7.52      8.86          218
Subject 8 -0.482     -3.06      2.10     -6.68      5.71          140
Subject 9  0.313     -2.31      2.93     -4.27      4.90           75
Performing modified Gibbs...
Performing FAST-MODE...
Performing FAST-MCMC...

＜おすすめのRに関する書籍です＞

RとPythonで学ぶ[実践的]データサイエンス&機械学習【増補改訂版】 | 有賀友紀, 大橋俊介

RとPythonを併走させながらデータサイエンスと機械学習を学ぶ。両方の言語を行き来する人に向く。

この記事が誰かの役に立ちますように。