RのpoLCAパッケージで潜在クラス分析のエントロピーを計算してみたい。
相対的エントロピーと絶対的エントロピー
エントロピーという指標はMplusユーザーにとってはお馴染みのものである。
Mplusで出力されるのは正確に言うと、相対的エントロピー(Relative Entropy)と呼ばれるものである。RのpoLCAパッケージでもエントロピーが出力できるが、絶対的エントロピーであって、相対的エントロピーとは異なったものだ。
相対的エントロピーの文献はRamaswamyらの論文であろうか。
相対化エントロピーは最小値が0で1が最大値である。1に近い方がクラスの分類が正確だとされている。
相対化エントロピーの値は0.7以上で十分、0.8以上で良好な分類だと言われている。
Ramaswamyらの論文ではこの基準が書かれているかというと、似たような記述はあるが、決して基準の提唱をしているわけではない。Mplusの作者Muthen & MuthenがMplusのマニュアルで0.8以上の推奨をしていることが確認できた。
下記の本の304頁には0.7以上で十分、0.8以上で良好という記述がある。
- 作者: Sigrid Bloemeke,Feng-Jui Hsieh,Gabriele Kaiser,William H. Schmidt
- 出版社/メーカー: Springer
- 発売日: 2013/09/23
- メディア: ハードカバー
- この商品を含むブログを見る
該当の論文はResearch Gateにある。
本来は異なったクラス数のモデル間で、相対化していないエントロピーを比較するのが正しい使い方なのではないかと思うが、標準化してもこの機能は果たせるので、相対化エントロピーでも問題はない。ただ、0.7以上といった基準は個人的にはあまりどうでもいい気がする。この議論は相関係数がどの程度あれば強い相関なのか、といった議論と同種のものであるからだ。
絶対的エントロピーの計算
潜在クラス分析は先のエントリーで使用したものを利用する。先のエントリでは共変量を設定しているが下記のスクリプトでは共変量のないノーマルな潜在クラス分析を行っている。
library("poLCA") data("election",package="poLCA") # electionデータ読み込み election$GENDER <- factor(election$GENDER) # 因子型にする d1 <- election[,c(1:3,7:12,16:17)] # 道徳性・思いやり・知識が豊富のゴア、ブッシュの値、性別、政党を抜き出し d1 <- na.omit(d1) # 欠損値除去 fm <- cbind(MORALG,CARESG,KNOWG,MORALB,CARESB,KNOWB)~1 # 式formulaの作成 res3 <- poLCA(formula=fm,data=d1,nclass=3,nrep=10) # 潜在クラス分析の実行 poLCA.entropy(res3) # エントロピーの計算
エントロピーは6.62と計算される。
6.621753
この関数でエラーが出るという報告がある(こちら)。僕の環境ではエラーは出ていないので、特に問題ないのだが、エラーが出る人もいるのだろうか。
その代わりに次のような手計算の式が書かれてあるのだが、僕の場合はこちらの計算でエラーが出た。
poLCA.entropy.fix <- function (lc) { K.j <- sapply(lc$probs, ncol) fullcell <- expand.grid(sapply(K.j, seq, from = 1)) P.c <- poLCA.predcell(lc, fullcell) return(-sum(P.c * log(P.c), na.rm = TRUE)) } poLCA.entropy.fix(lca.polca) [1] 3.602397
相対的エントロピーの計算
RのpoLCAは相対的エントロピーの計算はしてくれないので、自力で行う必要がある。自力といっても、潜在クラス分析の結果はもちろん使って行うので、先に潜在クラス分析の計算はしなくてはならない。下のスクリプトこちらを参考にした。
n <- 1346 # サンプルサイズ K <- 3 #クラス数 nume.E <- -sum(res3$posterior * log(res3$posterior)) # 分子 deno.E <- n*log(K) #分母 Entro <- 1-(nume.E/deno.E) Entro
結果は以下のようになる。
0.7329062
他の計算法としては、Daniel Oberskiが提案しているものがある(PDF)。
entropy <- function (p) sum(-p*log(p)) error_prior <- entropy(res3$P) # クラス割合 error_post <- mean(apply(res3$posterior, 1, entropy)) R2_entropy <- (error_prior - error_post) / error_prior R2_entropy