潜在クラス分析[Stata]

Stataで潜在クラス分析ができるのでやってみた。
Stata15から標準機能で利用できるらしい。

Lanzaらのプラグインを使えばStata11から利用可能である。

www.methodology.psu.edu

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)は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コマンドで尤度比、尤度比検定、AICBICが出力できる。

. 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

それぞれの結果をtwoclassthreeclassに格納する。この値は任意である。
この2つのフィッティング指標を比較することができる。

. 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.

AICBICをみるといずれも2クラスモデルの方が値が小さいので、2クラスモデルの方がさそうである。