構造クラス分析のCAICとｐ値

poLCAでCAICとｐ値を算出する方法が分かったのでメモ。

library(poLCA)

data(gss82)	
x <- gss82
f <- cbind(PURPOSE,ACCURACY,UNDERSTA,COOPERAT)~1
minc <- 1 #最小クラス数
maxc <- 5 #最大クラス数

pol.out <- matrix(0, maxc, 5)

# 計算にかなり時間がかかるので注意
for (i in minc:maxc) {
	temp <- poLCA(f, x, nclass=i, maxiter=100000)
	caic <- (-2 * temp$llik) + ((log(temp$N) + 1) * temp$npar) #CAIC
	Cp <- round(pchisq(temp$Chisq, temp$resid.df, lower.tail=FALSE), 6) #ｐ値
	Gp <- round(pchisq(temp$Gsq, temp$resid.df, lower.tail=FALSE), 6) #ｐ値
	pol.out[i,] <- c(temp$aic, temp$bic, caic, Cp, Gp)
	colnames(pol.out) <- c("AIC","BIC","CAIC","Chi.p","G.p")
	rownames(pol.out) <- minc:maxc
	}
pol.out
#        AIC      BIC     CAIC    Chi.p      G.p
# 1 5756.459 5787.010 5793.010 0.000000 0.000000
# 2 5592.536 5658.729 5671.729 0.000000 0.000000
# 3 5551.478 5653.313 5673.313 0.057085 0.060503
# 4 5547.242 5684.719 5711.719 0.745440 0.642435
# 5 5558.009 5731.128 5765.128 0.081996 0.093657

matplot(minc:maxc, pol.out[,1:3], type="b", xla="classes", yla="bic", col=2:4,
	pch=c("A","B","C"),
	main="AIC、BIC、CAICの比較\nデータセット“gss82”の場合")
legend(2, max(pol.out), legend=c("AIC","BIC","CAIC"),col=2:4,lty=c(1:3))

うーん、AIC、BIC、CAICそれぞれの支持するクラス数が違う……
ｐ値から見て４クラスかなぁ。