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パッケージを使っていたりします。