Stataで潜在クラス分析ができるのでやってみた。
Stata15から標準機能で利用できるらしい。
Lanzaらのプラグインを使えばStata11から利用可能である。
Stataの標準機能だと機能が足りないように思うので、Lanzaらのパッケージを結局使うことになりそうだ。
LanzaらのパッケージはSASでも利用できる。SASといっても、無料で使えるSAS University Editionには対応していない。
ソフトウェアの価格の面から考えるとStataでの利用となりそうだ。
以下で実行しているデモはStata15から搭載されている標準機能を使ったものである。
サンプルデータ
Stataの用意しているサンプルデータを利用する。
http://www.stata-press.com/data/r16/gsem_lca1.dta
| 変数 | ラベル |
|---|---|
| accident | would testify against friend in accident case 事故の場合に友人に対して証言するだろう |
| play | would give negative review of friend's play 友人の遊びに否定的な評価をするだろう |
| insurance | would disclose health concerns to friend's insurance company 友人の健康上の問題を保険会社に開示するだろう |
| stock | would keep company secret from friend 友人から会社の秘密を守るだろう |
何のデータがいまいちわからないが、デモであれば特に問題ないだろう。
潜在クラス分析
Stataでは潜在クラスの式は最初に走らせる必要がある。とりあえず3クラスモデルで走らせている。
. use http://www.stata-press.com/data/r16/gsem_lca1.dta, clear . gsem (accident play insurance stock <- ), logit lclass(C 3)
コマンドはgsemである。一般化SEMのバリエーションの一つだ。フォーマットは、()の中に潜在クラスに利用する変数を入れて最後に<-をつける。次に、コンマを打ちlogit lclass(C 3)と潜在クラスのオプションをつける。(C 3)は３クラスモデルを意味している。
クラス構成割合
クラス構成割合のコマンドと表示は以下である。
. estat lcprob
Latent class marginal probabilities Number of obs = 216
--------------------------------------------------------------
| Delta-method
| Margin Std. Err. [95% Conf. Interval]
-------------+------------------------------------------------
C |
1 | .1560766 15.31178 2.1e-100 1
2 | .6346788 11.94001 2.55e-44 1
3 | .2092446 3.373463 1.17e-18 1
--------------------------------------------------------------
応答確率
条件付き応答確率は次のように出力する。
. estat lcmean
Latent class marginal means Number of obs = 216
--------------------------------------------------------------
| Delta-method
| Margin Std. Err. [95% Conf. Interval]
-------------+------------------------------------------------
1 |
accident | .3935409 25.44072 1.20e-91 1
play | .0466199 12.58716 4.3e-243 1
insurance | .2541585 3.758665 4.52e-18 1
stock | .0656737 2.330868 3.26e-34 1
-------------+------------------------------------------------
2 |
accident | .8220824 4.909919 1.23e-28 1
play | .455712 9.045093 7.63e-32 1
insurance | .4224278 5.13116 9.18e-19 1
stock | .1818264 3.94581 5.89e-24 1
-------------+------------------------------------------------
3 |
accident | .9963741 .3046207 4.65e-70 1
play | .9725123 1.448169 2.73e-45 1
insurance | .9850408 3.020941 2.1e-173 1
stock | .881904 6.083824 1.41e-49 1
--------------------------------------------------------------
フィッティング指標
潜在クラス分析ではクラス数を検討しなければならない。
estat lcgofコマンドで尤度比、尤度比検定、AIC、BICが出力できる。
. estat lcgof
----------------------------------------------------------------------------
Fit statistic | Value Description
---------------------+------------------------------------------------------
Likelihood ratio |
chi2_ms(1) | 0.387 model vs. saturated
p > chi2 | 0.534
---------------------+------------------------------------------------------
Information criteria |
AIC | 1034.602 Akaike's information criterion
BIC | 1081.856 Bayesian information criterion
----------------------------------------------------------------------------
フィッティング指標の比較
. gsem (accident play insurance stock <- ), logit lclass(C 2) estimates store twoclass . gsem (accident play insurance stock <- ), logit lclass(C 3) estimates store threeclass
それぞれの結果をtwoclassとthreeclassに格納する。この値は任意である。
この２つのフィッティング指標を比較することができる。
. estimates stats twoclass threeclass
Akaike's information criterion and Bayesian information criterion
-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
twoclass | 216 . -504.4677 9 1026.935 1057.313
threeclass | 216 . -503.3011 14 1034.602 1081.856
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] BIC note.
