Rで潜在クラス分析
以前書いたメモにトラックバックを頂戴しました。
id:kosugitti:20091203
本当にただのメモだったのでなんとも申し訳ない。
なのでお詫びに、lca() と poLCA() と randomLCA() の簡単な比較でもしてみましょう。
# 潜在クラス分析のパッケージいろいろ library(e1071) # lca() library(poLCA) # poLCA() library(randomLCA) # randomLCA() # 比較のために使うデータセット。poLCAパッケージに入っている。 data(carcinoma)
それぞれの関数で同じデータを分析してbicを比べてみます。
まずは分析の実行。
# e1071パッケージのlca()の例 xx <- countpattern(carcinoma - 1) # lca()で分析できる形に変換 lca.bic <- numeric(6) for (i in 1:6) {lca.bic[i] <- lca(xx, i+1, niter=5)$bic} # poLCAの例 f <- cbind(A,B,C,D,E,F,G)~1 pol.bic <- numeric(6) for (i in 1:6) {pol.bic[i] <- poLCA(f, carcinoma, nclass=i+1, maxiter=6000)$bic} # randomLCAの例 # randomLCAは度数表から分析するのですが、生のデータからrandomLCAで扱える度数表 # を作る方法が分からなかったので、e1071のcountpatternの出力を参考に直接入力。 x <- matrix(c( 0,0,0,0,0,0,0,34, 0,0,0,0,1,0,0,2, 0,1,0,0,0,0,0,6, 0,1,0,0,0,0,1,1, 0,1,0,0,1,0,0,4, 0,1,0,0,1,0,1,5, 1,0,0,0,0,0,0,2, 1,0,1,0,1,0,1,1, 1,1,0,0,0,0,0,2, 1,1,0,0,0,0,1,1, 1,1,0,0,1,0,0,2, 1,1,0,0,1,0,1,7, 1,1,0,0,1,1,1,1, 1,1,0,1,0,0,1,1, 1,1,0,1,1,0,1,2, 1,1,0,1,1,1,1,3, 1,1,1,0,1,0,1,13, 1,1,1,0,1,1,1,5, 1,1,1,1,1,0,1,10, 1,1,1,1,1,1,1,16),20,8,byrow=TRUE) colnames(x) <- c("A","B","C","D","E","F","G","freq") rl.bic <- numeric(6) for (i in 1:6) {rl.bic[i] <- randomLCA(x[,1:7], freq=x[,8] ,nclass=i+1)$bic[1]} # 警告が出るが一応答えは出てる?
それぞれのbicをまとめてグラフに描いてみます。
y <- cbind(lca.bic, pol.bic, rl.bic) matplot(2:7, y, type="b", xla="classes", yla="bic", col=2:4, pch=c("l","p","r"), main="lca、poLCA、randomLCAの比較\nデータセット“carcinoma”の場合") legend(2, max(y), legend=c("lca","poLCA","randomLCA"),col=2:4,lty=c(1:3))
うーん、poLCAとrundomLCAはほぼ一致して3クラスを支持していますが、lcaは2クラスを支持していますね。
どれが妥当なのか私にはちょっと判断が付きませんが……
私の場合は、生データを直接突っ込めるのとよく分からない警告が出ないのでpoLCAを使っています。
graphs=TRUE を付けるとなにやら楽しげなグラフがウニウニ動いたりしますし。
とはいえ、潜在クラス分析自体の使いどころがいまひとつ掴めていないので、出番が無いのですが。
クラスタ分析と言うことだともっぱらclusterSimパッケージを使っていたりします。