RandomLCAパッケージを用いた潜在クラス分析とBLRT[R]

Rで潜在クラス分析をするパッケージとしてはpoLCAが有名だが、RandomLCAでも潜在クラス分析ができる。 poLCAではブートストラップ尤度比検定(Bootstrapped Likelihood Ratio Test: BLRT)ができないが、RandomLCAではできるらしい。

cran.r-project.org

RandomLCAの特徴

RandomLCAは1) スタンダードな潜在クラス分析 2)シングルレベルランダム効果モデル 3)2レベルランダム効果モデル3つが使用できる。

ちなみにこのパッケージの作者はBeath自身である。 RandomLCAは他の統計パッケージではできない解析ができるようだ。

ランダム効果モデルの話はまたの機会にするとして、今回はとりあえず、スタンダードな潜在クラス分析を行ってみたい。

データの特徴

RandomLCAで利用するデータの形式は私たちが普段使っている形式とは異なっている。サンプルデータのHIV testing dataを例にとろう。

f:id:iDES:20191119174243p:plain

V1~V4までの0/1のパターンは2の4乗なので16通りあり、右側にその条件を満たす頻度が書かれている。もともとは潜在クラス分析の計算はこの形式にしてから計算をするので、慣れない感じはあるが、こちらのデータセット形式の方が正統的ではある。このブログでは一度も言及したことはずだが、統計パッケージのLEMもこのような形式のデータセットを用いる。

潜在クラス分析

サンプルデータとして利用するのはHandelman et al. (1986)の歯のX線データである。 1-5行目が変数であり、6行目はfreqという頻度のデータである。

library(randomLCA)
data(dentistry) 
dentistry.lca2 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=2)
summary(dentistry.lca2)

結果は次のようになる。

  Classes      AIC      BIC     AIC3    logLik penlogLik
        2 14952.77 15021.64 14963.77 -7465.385 -7465.446
Class probabilities 
Class  1 Class  2 
  0.1961   0.8039 
Outcome probabilities 
             V1     V2     V3     V4     V5
Class  1 0.4034 0.7129 0.5981 0.4888 0.9155
Class  2 0.0106 0.1020 0.0136 0.0316 0.3053

BLRT

BLRTは下記のように実施する。 今回は2クラスモデルと3クラスモデルを比較する。

dentistry.lca2 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=2) # 2クラスモデル
dentistry.lca3 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=3) # 3クラスモデル
obslrt <- 2*(logLik(dentistry.lca3)-logLik(dentistry.lca2)) # 3クラス、2クラスモデルの対数尤度の差の二乗
print(obslrt)

Printはこの出力は、観察された(ブートストラップでシミュレーションされたものではない)実際の対数尤度の差の二乗である。この値はP値を求める際に使用する。

'log Lik.' 108.3152 (df=17)

BLRTは下記の次のように行う。

nsims <- 999 # ブートストラップの回数を指定
thesims <- simulate(dentistry.lca2, nsims) # シミュレーションを作成
simlrt <- rep(NA,nsims) 
for (isim in 1:nsims) {
submodel <- refit(dentistry.lca2,newpatterns=thesims[[isim]]) # submodelは比較対象の2クラスモデル
fullmodel <- refit(dentistry.lca3,newpatterns=thesims[[isim]]) # fullmodelは3クラスモデル
simlrt[isim] <- 2*(logLik(fullmodel)-logLik(submodel)) # 2つのモデルの尤度の差を二乗
}

P値はシミュレートされたLRTが実際のLRTより大きい割合である。1は修正値である。

print((sum(simlrt>=obslrt)+1)/(nsims+1)) # P値の計算

結果は下記のようになる。

0.001

P値が有意であった場合、2クラスモデルから3クラスモデルにしてモデルが改善した、つまり2ではなく3クラスモデルが支持されるという解釈を行う。

課題と考察

P値の計算式の分母に反復回数であるnsimsが入っているので、P値を細かく出そうとすると反復回数はそこそこ大きい数字の方が求められる。小数点3桁でP値を表示するには999+1回の反復が必要である。

この例ではシミュレートされたLRTが実際のLRTより大きかった回数(sum(simlrt>=obslrt)+1)は1回もないので値は0である。Rでsum(simlrt>=obslrt)+1と打って返されるのは0であることで確認できる。そこに修正値の1を足しているので分子は1である。

反復回数が9回だった場合、このモデル比較だと、P値は0+1/(9+1)であり、0.1、つまり10%となる。頻度主義は5%以下をだいたいの場合は有意とするので、モデルの改善がないと判断される。つまり、BLRTでは反復回数は完全に自由に決めてよいというわけではないのだ。

修正値が1の場合は反復回数は少なくとも20回ないと5%有意にならない。20回というのは分子が0(クラス数の増加によつて明らかにモデル改善がある場合)だけであって、ぎりぎりの数値であり、実際には20回というのは少なすぎる。

修正値を0.1にすると、反復回数は100回で小数点3桁のP値が得られる。修正値は分子・分母(特に分母)を0にしないために入れられるものなので、1である必要はないのかもしれない。0.1であってもよさそうに思える。ブートストラップの反復回数と修正値は連動しているというのは少なくとも間違いないだろう。修正値として何が適切なのかを知るには、BLRTについてちゃんと勉強しないといけないようである。