Rで潜在クラス分析をするパッケージとしてはpoLCAが有名だが、RandomLCAでも潜在クラス分析ができる。 poLCAではブートストラップ尤度比検定(Bootstrapped Likelihood Ratio Test: BLRT)ができないが、RandomLCAではできるらしい。
RandomLCAの特徴
RandomLCAは1) スタンダードな潜在クラス分析 2)シングルレベルランダム効果モデル 3)２レベルランダム効果モデル３つが使用できる。
The single level random effects model
Qu, Y., Tan, M. and Kutner, M.H. (1996) Random effects models in latent class analysis for evaluating accuracy of diagnostic tests. Biometrics, 52, 797–810. https://www.ncbi.nlm.nih.gov/pubmed/8805757The two level random effects model
Beath, Ken J., and Hellerm, Gillian Z. (2009) Latent trajectory modelling of multivariate binary data. Statistical Modelling, 9(3)), 199-213. https://journals.sagepub.com/doi/abs/10.1177/1471082X0800900302
ちなみにこのパッケージの作者はBeath自身である。 RandomLCAは他の統計パッケージではできない解析ができるようだ。
ランダム効果モデルの話はまたの機会にするとして、今回はとりあえず、スタンダードな潜在クラス分析を行ってみたい。
データの特徴
RandomLCAで利用するデータの形式は私たちが普段使っている形式とは異なっている。サンプルデータのHIV testing dataを例にとろう。
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より大きい割合である。１は修正値である。
print((sum(simlrt>=obslrt)+1)/(nsims+1)) # P値の計算
結果は下記のようになる。
0.001
P値が有意であった場合、２クラスモデルから３クラスモデルにしてモデルが改善した、つまり２ではなく3クラスモデルが支持されるという解釈を行う。
課題と考察
P値の計算式の分母に反復回数であるnsimsが入っているので、P値を細かく出そうとすると反復回数はそこそこ大きい数字の方が求められる。小数点3桁でP値を表示するには999+1回の反復が必要である。
この例ではシミュレートされたLRTが実際のLRTより大きかった回数(sum(simlrt>=obslrt)+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についてちゃんと勉強しないといけないようである。
