poLCAでデータ加工をしてRandomLCAでBLRTの計算をする[R]

私たちが使い慣れているデータ(行にケース、列に変数という形式のもの)を使って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に比べて手順が煩雑なのが欠点であろうか。