私たちが使い慣れているデータ(行にケース、列に変数という形式のもの)を使ってRでブートストラップ尤度比検定(Bootstrapped Likelihood Ratio Test)を計算する方法についてのエントリである。
poLCAには残念ながらBLRTの機能がないので、poLCAのpredcell関数でデータを加工してからRandomLCAパッケージでBLRTの計算をするという手順を踏む。
データ
poLCAのcarcinomaデータを利用する。
子宮頸部のがんのデータである。デモなのでデータの説明は省く。
library(poLCA) data(carcinoma) carcinoma[100:105,] #100~105行目を抽出
次のようなデータである。
A B C D E F G 100 2 2 2 2 2 1 2 101 2 2 2 2 2 1 2 102 2 2 2 2 2 1 2 103 2 2 2 2 2 2 2 104 2 2 2 2 2 2 2 105 2 2 2 2 2 2 2
predcell関数でデータ形式を変換
データ形式を下記の手順で変換する。
data(carcinoma) f <- cbind(A,B,C,D,E,F,G)~1 clc2 <- poLCA(f,carcinoma,nclass=2) carcinoma2<- clc2$predcell # スタック形式carcinoma2とする print(carcinoma2)
このようにするとRandomLCAで扱える形式のデータに変換できる。
A B C D F G observed expected 1 1 1 1 1 1 1 36 30.028 2 1 2 1 1 1 1 10 16.197 3 1 2 1 1 1 2 6 2.023 4 2 1 1 1 1 1 2 3.751 5 2 1 2 1 1 2 1 0.203 6 2 2 1 1 1 1 4 2.023 7 2 2 1 1 1 2 8 4.075 8 2 2 1 1 2 2 1 2.769 9 2 2 1 2 1 2 3 4.447 10 2 2 1 2 2 2 3 3.222 11 2 2 2 1 1 2 13 11.858 12 2 2 2 1 2 2 5 8.592 13 2 2 2 2 1 2 10 13.797 14 2 2 2 2 2 2 16 9.997
もっとスマートな方法は他にもありそうだが、とりあえずpredcell関数しか僕は発見できていないので、今回はこの方法を利用する。
randomLCAパッケージで潜在クラス分析
データcarcinoma2のA～Gはデータは1/2で構成されている。RandomLCAは0/1である必要があるため、memiscパッケージを利用して0/1にリコードする。
library(memisc) carcinoma2$A <- recode(carcinoma2$A , 0 <- 1, 1 <- 2) carcinoma2$B <- recode(carcinoma2$B , 0 <- 1, 1 <- 2) carcinoma2$C <- recode(carcinoma2$C , 0 <- 1, 1 <- 2) carcinoma2$D <- recode(carcinoma2$D , 0 <- 1, 1 <- 2) carcinoma2$E <- recode(carcinoma2$E , 0 <- 1, 1 <- 2) carcinoma2$F <- recode(carcinoma2$F , 0 <- 1, 1 <- 2) carcinoma2$G <- recode(carcinoma2$G , 0 <- 1, 1 <- 2)
データタイプは整数型でも実数型でも走るようだ。
今回は2クラスモデルと3クラスモデルを比較する。
比較する2クラスモデルと3クラスモデルを作成しておく。
library(randomLCA) carcinoma2.lca2 <- randomLCA(carcinoma2[,1:7],freq=carcinoma2$observed,nclass=2) carcinoma2.lca3 <- randomLCA(carcinoma2[,1:7],freq=carcinoma2$observed,nclass=3)
BLRT
BLRTのスクリプトに関する細かい説明は先のエントリを参照のこと。
まず、対数尤度の差の二乗。
## 3クラス、2クラスモデルの対数尤度の差の二乗 obslrt <- 2*(logLik(carcinoma2.lca3)-logLik(carcinoma2.lca2)) print(obslrt)
結果は以下。
'log Lik.' 47.09539 (df=23)
BLRT。
nsims <- 999 # ブートストラップの回数を指定
thesims <- simulate(carcinoma2.lca2, nsims) # シミュレーションを作成
simlrt <- rep(NA,nsims)
for (isim in 1:nsims) {
submodel <- refit(carcinoma2.lca2,newpatterns=thesims[[isim]]) # submodelは比較対象の2クラスモデル
fullmodel <- refit(carcinoma2.lca3,newpatterns=thesims[[isim]]) # fullmodelは3クラスモデル
simlrt[isim] <- 2*(logLik(fullmodel)-logLik(submodel)) # 2つのモデルの尤度の差を二乗
}
P値。
print((sum(simlrt>=obslrt)+1)/(nsims+1)) # P値の計算
結果は以下。
p-value=0.001
Mplusに比べて手順が煩雑なのが欠点であろうか。
