LCA
Ryohei Yamamoto
2024/03/30
Latent class analysisの自己学習用に作成
潜在クラス分析
潜在クラス分析は、観察されたカテゴリー変数と連続変数の特定の組み合わせを使用して、人口統計学や臨床特性などの類似した特性を持つ人々のグループをクラスタリングする教師なし機械学習手法の一種である。「潜在的」という用語は、サブグループが潜在的に存在するが、直接観察することができないという概念に由来する。
この分析では、観察されたデータは、いくつかの潜在的なグループの混合分布から生成されると仮定される。各参加者の各サブグループに属する確率は、混合モデルの最尤推定に基づいて推定される。通常、グループの数の選択は、ベイズ情報量規準（BIC）に基づいて行われる。BICが最小のモデルが、最良適合モデルとして決定される。BICは、最尤推定での尤度関数の値を用いて計算される。（PMID:
20192788、29982203）。
このデモでは、リッカート尺度を用いた調査研究の例でコードを確認していく。これらの集団のサブグループをみつけたいというモチベーションの場合に、２つの方法がある。 １つは潜在クラス分析で、もう1つはクラスター分析である。 階層クラスタリングやk-meansクラスタリングなどのクラスター分析は、大規模なデータセットをグループ化するために選択される方法であある。クラスター分析は、ユークリッド距離のような、個体間の類似性／非類似性のグラフィカルな距離尺度を生成する。そして、個々の測定値のサブグループ内の変動（すなわち、最短距離）を最小にする一方で、グループ間の変動（すなわち、最大距離）を最大にするサブグループを作成しようとする（PMID 35636591）。
この2つの手法にはいくつかの重要な違いがある。第1に、k-meansクラスタ分析がオブザベーションを絶対的にグループに割り当てるのに対して、潜在クラス分析では確率に基づいてオブザベーションをグループに割り当てる。第2に、k-meansは、多くのソフトウェア・パッケージで簡単に利用できるが、それは連続データにのみ使用可能である。潜在クラスタ分析は、多くのソフトウェア・パッケージでそれほど広く利用可能ではないが、カテゴリ・データを取り扱うことができる。
潜在クラス分析ソフトウェア・パッケージは少ない。おそらく最もよく使われているのはLatent Goldであるが、ライセンスが高額である。これは、日常業務に潜在クラスモデリングが含まれていない場合に特に当てはまる。SASではproc lcaで可能である。オープンソースでは、RパッケージのpoLCAとMCLUSTなどがある。Latent Goldで利用可能な多くの機能を必要としない限り、データ分析にはこれらのRパッケージで一般的に十分と思われる。
LCAは、クラスター分析よりもいくつかの利点がある：
LCAは、測定タイプが混在する変数（2値、順序、区間尺度、および幅が変化する尺度）をよりよく管理できる
LCAは、欠損データをよりよく扱うことができる
LCAは、より高い分類精度のために使用できるモデルベースのパラメータを提供する
クラスター分析は、確率的アプローチ（例えば、2ステップ・クラスター分析）と組み合わせることができるが、LCAは、パラメータの推定において2ステップ・アプローチに従う必要がない
Sample size
- Nylundは情報量基準（IC）と尤度検定が500～1,000人の被験者に基づいて正しいモデルを正確に識別できるが200人以下では識別できない
- 別のレビューでは、参加者が70人未満のモデルは「実行不可能」としている
- 100人未満の被験者に基づくモデルは「注意して」行うべき
- ほとんどのLCAは300人以上のサンプルサイズを必要とする
(PMID 37286148)
変数選択
指標変数は多いほうが良い
いくつの変数を含めるべきかのコンセンサスはないが、４つから２０以上まで報告はまちまちである
含める変数にはモデルに含める理論的根拠があった方が良い。 例：タバコの喫煙パターンを見るのであればアルコール嗜好を含めるがアイスの嗜好はいらないかもしれない
R code
rm(list = ls())
set.seed(1234)
library(pacman)
p_load(e1071,poLCA,reshape2,tidyverse,knitr,ztable,
forcats)以下は、American National Election Study(ANES)のデータを分析する方法の例である。 データは2000年アメリカ国民選挙調査の調査データ。大統領候補のアル・ゴアとジョージ・W・ブッシュについて、様々な特徴（道徳的、思いやりがある、物知り、優れた指導者、不誠実、知的）がどの程度あてはまるか、回答者の意見を尋ねる6つの質問を2人の候補者それぞれで尋ねた。回答4段階のリッカート尺度で、(1) 非常に良い、(2) かなり良い、(3) あまり良くない、(4) 全く良くない、である。多くの回答者は、これらの変数にさまざまな欠損値を持つ。
VOTE3 は、(1)Gore、(2)Bush、(3)Other でコード化されている。
EDUCは、(1)8学年以下、(2)9～11学年、それ以上の学歴なし、(3)高卒または同等の学歴、(4)12年以上の学歴、それ以上の学歴なし、(5)短大またはコミュニティ・カレッジ・レベルの学歴、(6)学士レベルの学歴、上級学歴なし、(7)上級学歴、としてコード化される。
GENDERは、(1)男性、(2)女性。
PARTYは、(1) 強い民主党、(2) 弱い民主党、(3) 無所属-民主党、(4) 無所属-無所属、(5) 無所属-共和党、(6) 弱い共和党、(7) 強い共和党でダミー変数として格納されてる.
# Example dataset from the poLCA package
data(election)
str(election)## 'data.frame': 1785 obs. of 17 variables:
## $ MORALG : Factor w/ 4 levels "1 Extremely well",..: 3 4 1 2 2 2 2 2 2 1 ...
## $ CARESG : Factor w/ 4 levels "1 Extremely well",..: 1 3 2 2 4 3 NA 2 2 1 ...
## $ KNOWG : Factor w/ 4 levels "1 Extremely well",..: 2 4 2 2 2 3 2 2 2 1 ...
## $ LEADG : Factor w/ 4 levels "1 Extremely well",..: 2 3 1 2 3 2 2 2 3 1 ...
## $ DISHONG: Factor w/ 4 levels "1 Extremely well",..: 3 2 3 2 2 2 4 3 4 4 ...
## $ INTELG : Factor w/ 4 levels "1 Extremely well",..: 2 2 2 2 2 2 2 2 2 1 ...
## $ MORALB : Factor w/ 4 levels "1 Extremely well",..: 1 NA 2 2 1 2 NA 3 2 2 ...
## $ CARESB : Factor w/ 4 levels "1 Extremely well",..: 1 NA 2 3 1 NA 3 4 3 4 ...
## $ KNOWB : Factor w/ 4 levels "1 Extremely well",..: 2 2 2 2 2 3 2 4 2 2 ...
## $ LEADB : Factor w/ 4 levels "1 Extremely well",..: 2 3 3 2 2 2 2 4 2 3 ...
## $ DISHONB: Factor w/ 4 levels "1 Extremely well",..: 4 NA 3 3 3 2 4 3 3 3 ...
## $ INTELB : Factor w/ 4 levels "1 Extremely well",..: 2 1 2 1 2 2 NA 4 2 2 ...
## $ VOTE3 : num 2 NA 1 1 2 1 NA 1 1 1 ...
## $ AGE : num 49 35 57 63 40 77 43 47 26 48 ...
## $ EDUC : num 5 4 3 4 5 2 4 7 6 3 ...
## $ GENDER : num 1 2 2 1 2 1 1 2 2 2 ...
## $ PARTY : num 5 3 1 3 7 1 6 1 1 1 ...PoLCAを使ったLCAのコードは以下。
# build the model with PARTY as the covariate
# select variables
mydata <- election %>% dplyr::select(MORALG, CARESG, KNOWG, LEADG, DISHONG, INTELG, MORALB, CARESB, KNOWB, LEADB, DISHONB, INTELB)
# define function
f<-with(mydata, cbind(MORALG, CARESG, KNOWG, LEADG, DISHONG, INTELG, MORALB, CARESB, KNOWB, LEADB, DISHONB, INTELB)~1)“~1” はpoLCAに基本的な潜在クラスモデルを推定するように指示する。 “~1” の部分に共変量を設定すると、潜在クラス回帰分析となる。
クラスター数を指定して解析
クラスター数は３つ（賛成、中立、反対）のように事前の仮説がある場合はクラスタ数を事前に指定して解析する。
# Run LCA on the ANES 2000 dataset 3 classes
anes2000 <- poLCA(f,election,nclass=3,nrep=10)## Model 1: llik = -16714.66 ... best llik = -16714.66
## Model 2: llik = -16714.66 ... best llik = -16714.66
## Model 3: llik = -16714.66 ... best llik = -16714.66
## Model 4: llik = -16714.66 ... best llik = -16714.66
## Model 5: llik = -16714.66 ... best llik = -16714.66
## Model 6: llik = -16714.66 ... best llik = -16714.66
## Model 7: llik = -16714.66 ... best llik = -16714.66
## Model 8: llik = -16714.66 ... best llik = -16714.66
## Model 9: llik = -16714.66 ... best llik = -16714.66
## Model 10: llik = -16714.66 ... best llik = -16714.66
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $MORALG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.1013 0.5990 0.2552 0.0446
## class 2: 0.1851 0.3668 0.2350 0.2131
## class 3: 0.5224 0.4109 0.0390 0.0277
##
## $CARESG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0227 0.5090 0.3819 0.0864
## class 2: 0.0879 0.2546 0.3524 0.3051
## class 3: 0.3970 0.4605 0.0895 0.0529
##
## $KNOWG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0429 0.8058 0.1408 0.0105
## class 2: 0.2542 0.4253 0.2382 0.0822
## class 3: 0.5878 0.3543 0.0266 0.0312
##
## $LEADG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0207 0.4886 0.4268 0.0639
## class 2: 0.0747 0.2582 0.3856 0.2816
## class 3: 0.3747 0.4772 0.1085 0.0396
##
## $DISHONG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0519 0.2182 0.4998 0.2301
## class 2: 0.2007 0.3060 0.2774 0.2158
## class 3: 0.0426 0.0487 0.3371 0.5717
##
## $INTELG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0600 0.7953 0.1317 0.0130
## class 2: 0.2936 0.4407 0.1852 0.0805
## class 3: 0.5763 0.3698 0.0197 0.0341
##
## $MORALB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0642 0.7137 0.2164 0.0057
## class 2: 0.5797 0.3787 0.0248 0.0168
## class 3: 0.0966 0.3795 0.3638 0.1601
##
## $CARESB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0093 0.4847 0.4380 0.0679
## class 2: 0.3483 0.5362 0.0836 0.0319
## class 3: 0.0162 0.1125 0.4045 0.4668
##
## $KNOWB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0191 0.7583 0.2182 0.0045
## class 2: 0.4970 0.4899 0.0132 0.0000
## class 3: 0.0682 0.3115 0.3828 0.2374
##
## $LEADB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0537 0.7119 0.2217 0.0127
## class 2: 0.5106 0.4648 0.0187 0.0058
## class 3: 0.0236 0.2735 0.4524 0.2505
##
## $DISHONB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0134 0.1471 0.5683 0.2711
## class 2: 0.0365 0.0740 0.1993 0.6902
## class 3: 0.0862 0.3096 0.4111 0.1931
##
## $INTELB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0435 0.7741 0.1796 0.0028
## class 2: 0.5284 0.4656 0.0060 0.0000
## class 3: 0.1179 0.3391 0.3440 0.1990
##
## Estimated class population shares
## 0.4194 0.2608 0.3198
##
## Predicted class memberships (by modal posterior prob.)
## 0.4256 0.2578 0.3166
##
## =========================================================
## Fit for 3 latent classes:
## =========================================================
## number of observations: 1311
## number of estimated parameters: 110
## residual degrees of freedom: 1201
## maximum log-likelihood: -16714.66
##
## AIC(3): 33649.32
## BIC(3): 34218.96
## G^2(3): 15016.7 (Likelihood ratio/deviance statistic)
## X^2(3): 16748689622 (Chi-square goodness of fit)
## nclass: モデルで計算される潜在クラスの数を指定する。デフォルトは2潜在クラス。poLCAは、実行するたびに潜在クラスの1セットを仮定する。したがって、複数のクラスを考えている場合は、複数回仮定する潜在クラスの数を指定し、コマンドを何回も実行する。
maxiter: 収束のための最大反復数。この反復回数に達する前に収束しない場合は、エラー・メッセージが表示され、分析は終了される。
graphs: パラメータ推定値を示すグラフを作成するかどうかを指定。デフォルトはFALSE。大規模なデータセットで4つ以上の潜在クラスを持つ分析を実行するには、Rではかなり時間がかかる。したがって、グラフを作成せずに分析を実行し、必要であれば、結果のモデルが比較された後、最も適合したモデルについてのみグラフを作成する方が早い。
na.rm: これはpoLCAが欠損値を持つケースをどのように扱うかを指定する。TRUEを指定すると、これらのケースはモデル推定の前にリストワイズ削除によって取り除かれる。FALSEを指定すると、欠損値を持つケースは保持される。デフォルトは TRUE。poLCAは欠損値を持つケースを計算から除外するので、モデルを推定する前に欠損値を持つケースを削除する必要はない。
nrep: このオプションは、異なる開始値を用いてモデルを推定する回数を指定するために使用する。アルゴリズムが対数尤度関数の局所最大値ではなく、大域最大値を見つけることを保証するために、nrepを1より大きく設定することが望ましい。
verbose: モデルの結果を画面に出力するかどうかを指定。デフォルトはTRUE。
出力の最初には以下が表示される。
Model 1: llik = -16714.66 … best llik = -16714.66
Model 2: llik = -16714.66 … best llik = -16714.66
Model 3: llik = -16714.66 … best llik = -16714.66
Model 4: llik = -16714.66 … best llik = -16714.66
Model 5: llik = -16714.66 … best llik = -16714.66
Model 6: llik = -16714.66 … best llik = -16714.66
Model 7: llik = -16714.66 … best llik = -16714.66
Model 8: llik = -16714.66 … best llik = -16714.66
Model 9: llik = -16714.66 … best llik = -16714.66
Model 10: llik = -16714.66 … best llik = -16714.66
出力の最初の部分（Model 1 to Model 10 llik and best llik）は、潜在クラス・モデルが異なる開始値を用いて（nrep=10で指定）10回推定されたことを示している。割り当てられた結果は、対数尤度関数の最大値を持つモデルについて推定されたもので、モデルのフィッティングの最初の試みで、-16714.66というグローバルな最大対数尤度が見つかったことを意味する。
出力の次のセクションは、各クラスについて、結果変数による条件付き項目回答確率を提供している。この出力は、各潜在クラスの回答者が、指標変数の回答を提供する確率を示す。このケースでは、モデルは3つの潜在クラスを仮定したので、3つの行がある。列は、指標変数のカテゴリ（1 Extremely well 2 Quite well 3 Not too well 4 Not well at all）を示している。
$MORALGは、この分析で使用される道徳心（ゴア）の変数である。この分析の場合、各変数への回答は、それぞれ1、2、3、4とコード化されている。したがって、上記の出力では、最初の列は、Extermely wellとなる確率である。潜在クラス1の回答者がMORALGに「Extermely well」と回答する確率は10%で、Quite wellと回答する確立は60％である。
“Estimated class population shares”は、各潜在クラスに属するオブザベーションのシェアに対応する推定比率を提供している。したがって、3クラス・モデルの場合、オブザベーションのシェアは、潜在クラス1で41.9%、潜在クラス2で26.1%, 潜在クラス３で32%である。
出力のPredicted class membershipsは、潜在クラスのサイズを推定するもう1つの方法である。これは事後確率を用いて潜在クラスにオブザベーションを割り当てる。Estimated class population sharesとPredicted class membershipsの値が似ている場合、よいモデル適合の徴候を示す。しかし、この値の一致だけで、モデルの適合を評価すべきではない。他の指標は次の出力にある、AIC、BIC、G2、X2である。
オブザベーションの数 は，分析で使用された完全にオブザベーションされたケース（すなわち，欠損値のないケース）の数で、推定パラメータの数は、モデルで使用された自由度の数を示す。残差自由度”が負でないことを確認する。負の残差は、通常、エラー・メッセージになる。
AICは赤池情報量規準で、BICはベイズ情報量規準で、潜在クラスの数が異なるモデルが比較されるとき，通常，最小のAICとBICを持つモデルが，ベスト・フィット・モデルとして選ばれる。これは、情報量規準の値が低いほど、モデルの適合とパーシモンのバランスがよいことを示唆するからである (Lanza, S.T. & Rhoades, B.L., 2013)。時には、AICとBICが同じモデルを最適適合モデルとして示さないこともある。基本的な探索的LCAの場合、モデルが比較的単純であるため、通常はBICの方がより適切である（Lin and Dayton 1997; Forster 2000）。
尤度比統計量とカイ二乗統計量も、モデルの適合度を評価するために使用できる。AICやBICと同様に、尤度比統計量とカイ2乗統計量を最小化するモデルを選択する一方で、パラメータの数を少なく維持することが目的です。これらの検定は、カイ2乗分布の仮定に違反するため、スパース・データ（LCAでよくある問題）では役に立たないことに注意することが重要である (Linzer & Lewis, 2011; Nylund et. al, 2007)。Nylund et. al (2007)は、モデル適合度を評価するための様々な方法を調査した研究から、ブートストラップ尤度比検定がモデル適合度を評価するための最良の方法であり、次いでBICであると結論付けている。poLCAパッケージを使ってブートストラップ尤度比検定を計算することはできない。したがって、poLCAを使用する場合は、BICを使用して最も適合するモデルを選択することが推奨される。
上記を個別に見る場合は以下のコマンドで実行できる
各変数の値の条件付き確率
anes2000$probs # クラスタを付与した時の、クラスタリングに使った各変数の値の条件付き確率## $MORALG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.1012766 0.598972 0.25517643 0.04457489
## class 2: 0.1851216 0.366794 0.23498253 0.21310188
## class 3: 0.5223517 0.410870 0.03902982 0.02774846
##
## $CARESG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.02267347 0.5090242 0.38194661 0.08635567
## class 2: 0.08791577 0.2545822 0.35235266 0.30514933
## class 3: 0.39702602 0.4605172 0.08952707 0.05292976
##
## $KNOWG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.04293128 0.8057945 0.14075120 0.01052305
## class 2: 0.25419444 0.4253413 0.23823313 0.08223118
## class 3: 0.58784748 0.3542820 0.02662775 0.03124273
##
## $LEADG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.02067176 0.4886119 0.4267951 0.06392122
## class 2: 0.07465177 0.2581511 0.3856206 0.28157648
## class 3: 0.37469255 0.4772173 0.1085304 0.03955969
##
## $DISHONG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.05187972 0.21819340 0.4998461 0.2300807
## class 2: 0.20074126 0.30601567 0.2774392 0.2158039
## class 3: 0.04255801 0.04868314 0.3370558 0.5717031
##
## $INTELG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.05996179 0.7953262 0.13171894 0.01299303
## class 2: 0.29361616 0.4406839 0.18515365 0.08054628
## class 3: 0.57630157 0.3698074 0.01974336 0.03414764
##
## $MORALB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0641706 0.7137016 0.21641132 0.005716453
## class 2: 0.5797291 0.3786837 0.02480855 0.016778687
## class 3: 0.0966234 0.3795247 0.36377749 0.160074393
##
## $CARESB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.009327214 0.4847185 0.43802804 0.06792621
## class 2: 0.348319004 0.5362152 0.08357466 0.03189113
## class 3: 0.016181151 0.1125055 0.40447168 0.46684166
##
## $KNOWB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.01905866 0.7582515 0.21822758 4.462285e-03
## class 2: 0.49695736 0.4898728 0.01316988 2.102069e-230
## class 3: 0.06824680 0.3115366 0.38281017 2.374064e-01
##
## $LEADB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.05371343 0.7118851 0.22172329 0.012678182
## class 2: 0.51063341 0.4648412 0.01871262 0.005812796
## class 3: 0.02358063 0.2734986 0.45240814 0.250512667
##
## $DISHONB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.01341454 0.14711021 0.5683310 0.2711442
## class 2: 0.03651393 0.07397303 0.1993052 0.6902078
## class 3: 0.08619161 0.30962388 0.4111306 0.1930539
##
## $INTELB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.04352115 0.7741027 0.17955295 0.002823181
## class 2: 0.52835819 0.4656163 0.00602548 0.000000000
## class 3: 0.11788768 0.3391487 0.34395083 0.199012746クラスタごとの所属確率
anes2000$posterior # 各サンプルのクラスタごとの所属確率各サンプルの所属クラスタ
anes2000$predclass # 各サンプルの所属クラスタ（ハードクラスタリング的な解釈）クラスター数の検索
今度はクラスター数を指定せずに解析する。事前にクラスター数が存在しない場合である。２～１０クラスターを指定して、解析するコードを示す。
#------ run a sequence of models with 1-10 classes and print out the model with the lowest BIC
max_II <- -100000
min_bic <- 100000
for(i in 2:10){
lc <- poLCA(f, mydata, nclass=i, maxiter=3000,
tol=1e-5, na.rm=FALSE,
nrep=10, verbose=TRUE, calc.se=TRUE)
if(lc$bic < min_bic){
min_bic <- lc$bic
LCA_best_model<-lc
}
}poLCA-packageから最高のモデルの標準出力が得られる。 Negative alertが出る場合は自由度が負になっているので、適切ではないモデルを意味する。
LCA_best_model## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $MORALG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3122 0.5222 0.1363 0.0292
## class 2: 0.0554 0.1769 0.2709 0.4968
## class 3: 0.5736 0.3811 0.0234 0.0219
## class 4: 0.7614 0.2334 0.0052 0.0000
## class 5: 0.0969 0.3501 0.3428 0.2103
## class 6: 0.1687 0.7495 0.0818 0.0000
## class 7: 0.1806 0.7263 0.0852 0.0079
## class 8: 0.0054 0.3932 0.4574 0.1440
##
## $CARESG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.1512 0.4632 0.2737 0.1119
## class 2: 0.0098 0.0262 0.3103 0.6537
## class 3: 0.4752 0.4594 0.0405 0.0249
## class 4: 0.6658 0.2798 0.0381 0.0163
## class 5: 0.0157 0.1806 0.3451 0.4587
## class 6: 0.0788 0.7025 0.2100 0.0087
## class 7: 0.0455 0.7579 0.1650 0.0316
## class 8: 0.0072 0.1237 0.6743 0.1948
##
## $KNOWG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.4007 0.5017 0.0976 0.0000
## class 2: 0.0956 0.3485 0.3006 0.2554
## class 3: 0.6531 0.3469 0.0000 0.0000
## class 4: 0.8556 0.1395 0.0000 0.0049
## class 5: 0.1006 0.3789 0.3241 0.1964
## class 6: 0.1515 0.7922 0.0484 0.0079
## class 7: 0.0928 0.8682 0.0390 0.0000
## class 8: 0.0187 0.6711 0.2948 0.0155
##
## $LEADG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.1199 0.4730 0.3557 0.0514
## class 2: 0.0070 0.0304 0.2892 0.6734
## class 3: 0.4675 0.4698 0.0627 0.0000
## class 4: 0.6796 0.3081 0.0000 0.0123
## class 5: 0.0215 0.1883 0.4864 0.3038
## class 6: 0.0294 0.7405 0.2301 0.0000
## class 7: 0.0374 0.7505 0.2011 0.0111
## class 8: 0.0108 0.1028 0.7080 0.1785
##
## $DISHONG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0466 0.2329 0.4545 0.2661
## class 2: 0.4395 0.2532 0.1520 0.1553
## class 3: 0.0000 0.0653 0.2836 0.6511
## class 4: 0.0704 0.0275 0.2019 0.7003
## class 5: 0.1422 0.2374 0.3094 0.3109
## class 6: 0.0164 0.1041 0.4986 0.3809
## class 7: 0.0063 0.1184 0.5548 0.3205
## class 8: 0.1321 0.4331 0.3403 0.0944
##
## $INTELG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.4700 0.4968 0.0332 0.0000
## class 2: 0.1428 0.3056 0.3008 0.2508
## class 3: 0.7035 0.2796 0.0069 0.0099
## class 4: 0.8729 0.1120 0.0026 0.0125
## class 5: 0.0572 0.4526 0.2867 0.2035
## class 6: 0.1508 0.7800 0.0656 0.0036
## class 7: 0.0822 0.8905 0.0273 0.0000
## class 8: 0.0490 0.6794 0.2477 0.0239
##
## $MORALB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.6119 0.3629 0.0252 0.0000
## class 2: 0.7604 0.2098 0.0215 0.0084
## class 3: 0.0739 0.1701 0.4350 0.3209
## class 4: 0.1574 0.5398 0.2611 0.0417
## class 5: 0.0668 0.1304 0.4229 0.3800
## class 6: 0.0386 0.4738 0.4469 0.0407
## class 7: 0.0499 0.7963 0.1508 0.0030
## class 8: 0.1587 0.7616 0.0797 0.0000
##
## $CARESB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3512 0.5152 0.1053 0.0283
## class 2: 0.4976 0.4542 0.0483 0.0000
## class 3: 0.0078 0.0189 0.2411 0.7321
## class 4: 0.0136 0.1687 0.5252 0.2925
## class 5: 0.0265 0.0444 0.3684 0.5607
## class 6: 0.0075 0.1686 0.5340 0.2899
## class 7: 0.0016 0.4515 0.4343 0.1125
## class 8: 0.0380 0.7186 0.2369 0.0065
##
## $KNOWB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5557 0.4261 0.0182 0.0000
## class 2: 0.6735 0.3103 0.0162 0.0000
## class 3: 0.0125 0.0497 0.3686 0.5692
## class 4: 0.1587 0.5938 0.2475 0.0000
## class 5: 0.0974 0.4324 0.2215 0.2486
## class 6: 0.0036 0.1462 0.7617 0.0884
## class 7: 0.0267 0.9237 0.0367 0.0129
## class 8: 0.0440 0.8198 0.1362 0.0000
##
## $LEADB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5400 0.4213 0.0387 0.0000
## class 2: 0.6336 0.3664 0.0000 0.0000
## class 3: 0.0000 0.0133 0.4255 0.5612
## class 4: 0.0624 0.4193 0.4057 0.1125
## class 5: 0.0893 0.2268 0.3403 0.3436
## class 6: 0.0000 0.3063 0.5829 0.1107
## class 7: 0.0402 0.7672 0.1857 0.0070
## class 8: 0.1125 0.8258 0.0584 0.0033
##
## $DISHONB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0293 0.0926 0.2450 0.6331
## class 2: 0.0193 0.0000 0.1077 0.8730
## class 3: 0.1860 0.4076 0.2332 0.1732
## class 4: 0.0541 0.2533 0.4175 0.2750
## class 5: 0.2596 0.2828 0.2800 0.1776
## class 6: 0.0143 0.2579 0.5746 0.1532
## class 7: 0.0000 0.1537 0.5889 0.2574
## class 8: 0.0033 0.1194 0.4359 0.4414
##
## $INTELB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.6371 0.3536 0.0050 0.0043
## class 2: 0.6267 0.3733 0.0000 0.0000
## class 3: 0.0535 0.0821 0.4204 0.4440
## class 4: 0.2310 0.6180 0.1509 0.0000
## class 5: 0.1349 0.3792 0.2790 0.2069
## class 6: 0.0000 0.2804 0.6198 0.0999
## class 7: 0.0552 0.9270 0.0178 0.0000
## class 8: 0.0713 0.8422 0.0853 0.0012
##
## Estimated class population shares
## 0.1276 0.0677 0.0814 0.12 0.0648 0.1461 0.2058 0.1866
##
## Predicted class memberships (by modal posterior prob.)
## 0.1255 0.0672 0.0784 0.1261 0.0611 0.1445 0.2106 0.1866
##
## =========================================================
## Fit for 8 latent classes:
## =========================================================
## number of observations: 1785
## number of fully observed cases: 1311
## number of estimated parameters: 295
## residual degrees of freedom: 1490
## maximum log-likelihood: -19922.58
##
## AIC(8): 40435.16
## BIC(8): 42053.87
## G^2(8): 12818.57 (Likelihood ratio/deviance statistic)
## X^2(8): 22833756 (Chi-square goodness of fit)
## 8クラスのモデルがBIC最小になっているようだ。
個別に行ってみる。
## models with different number of groups without covariates:
set.seed(01012)
lc1<-poLCA(f, data=mydata, nclass=1, na.rm = FALSE, nrep=30, maxiter=3000) #Loglinear independence model.
lc2<-poLCA(f, data=mydata, nclass=2, na.rm = FALSE, nrep=30, maxiter=3000)
lc3<-poLCA(f, data=mydata, nclass=3, na.rm = FALSE, nrep=30, maxiter=3000)
lc4<-poLCA(f, data=mydata, nclass=4, na.rm = FALSE, nrep=30, maxiter=3000)
lc5<-poLCA(f, data=mydata, nclass=5, na.rm = FALSE, nrep=30, maxiter=3000)
lc6<-poLCA(f, data=mydata, nclass=6, na.rm = FALSE, nrep=30, maxiter=3000)上の結果をまとめる。 以下のコードは https://statistics.ohlsen-web.de/latent-class-analysis-polca/ から引用しているが、ABICの計算があってるかはわからなかった。
# generate dataframe with fit-values
results <- data.frame(Modell=c("Modell 1"),
log_likelihood=lc1$llik,
df = lc1$resid.df,
BIC=lc1$bic,
ABIC= (-2*lc1$llik) + ((log((lc1$N + 2)/24)) * lc1$npar),
CAIC = (-2*lc1$llik) + lc1$npar * (1 + log(lc1$N)),
likelihood_ratio=lc1$Gsq)
results$Modell<-as.integer(results$Modell)## Warning: 強制変換により NA が生成されましたresults[1,1]<-c("Modell 1")
results[2,1]<-c("Modell 2")
results[3,1]<-c("Modell 3")
results[4,1]<-c("Modell 4")
results[5,1]<-c("Modell 5")
results[6,1]<-c("Modell 6")
results[2,2]<-lc2$llik
results[3,2]<-lc3$llik
results[4,2]<-lc4$llik
results[5,2]<-lc5$llik
results[6,2]<-lc6$llik
results[2,3]<-lc2$resid.df
results[3,3]<-lc3$resid.df
results[4,3]<-lc4$resid.df
results[5,3]<-lc5$resid.df
results[6,3]<-lc6$resid.df
results[2,4]<-lc2$bic
results[3,4]<-lc3$bic
results[4,4]<-lc4$bic
results[5,4]<-lc5$bic
results[6,4]<-lc6$bic
results[2,5]<-(-2*lc2$llik) + ((log((lc2$N + 2)/24)) * lc2$npar) #abic
results[3,5]<-(-2*lc3$llik) + ((log((lc3$N + 2)/24)) * lc3$npar)
results[4,5]<-(-2*lc4$llik) + ((log((lc4$N + 2)/24)) * lc4$npar)
results[5,5]<-(-2*lc5$llik) + ((log((lc5$N + 2)/24)) * lc5$npar)
results[6,5]<-(-2*lc6$llik) + ((log((lc6$N + 2)/24)) * lc6$npar)
results[2,6]<- (-2*lc2$llik) + lc2$npar * (1 + log(lc2$N)) #caic
results[3,6]<- (-2*lc3$llik) + lc3$npar * (1 + log(lc3$N))
results[4,6]<- (-2*lc4$llik) + lc4$npar * (1 + log(lc4$N))
results[5,6]<- (-2*lc5$llik) + lc5$npar * (1 + log(lc5$N))
results[6,6]<- (-2*lc6$llik) + lc6$npar * (1 + log(lc6$N))
results[2,7]<-lc2$Gsq
results[3,7]<-lc3$Gsq
results[4,7]<-lc4$Gsq
results[5,7]<-lc5$Gsq
results[6,7]<-lc6$Gsq
entropy<-function (p) sum(-p*log(p))
results$R2_entropy## NULLresults[1,8]<-c("-")
error_prior<-entropy(lc2$P) # class proportions model 2
error_post<-mean(apply(lc2$posterior,1, entropy),na.rm = TRUE)
results[2,8]<-round(((error_prior-error_post) / error_prior),3)
error_prior<-entropy(lc3$P) # class proportions model 3
error_post<-mean(apply(lc3$posterior,1, entropy),na.rm = TRUE)
results[3,8]<-round(((error_prior-error_post) / error_prior),3)
error_prior<-entropy(lc4$P) # class proportions model 4
error_post<-mean(apply(lc4$posterior,1, entropy),na.rm = TRUE)
results[4,8]<-round(((error_prior-error_post) / error_prior),3)
error_prior<-entropy(lc5$P) # class proportions model 5
error_post<-mean(apply(lc5$posterior,1, entropy),na.rm = TRUE)
results[5,8]<-round(((error_prior-error_post) / error_prior),3)
error_prior<-entropy(lc6$P) # class proportions model 6
error_post<-mean(apply(lc6$posterior,1, entropy),na.rm = TRUE)
results[6,8]<-round(((error_prior-error_post) / error_prior),3)
# combining results to a dataframe
colnames(results)<-c("Model","log-likelihood","resid. df","BIC","aBIC","cAIC","likelihood-ratio","Entropy")
lca_results<-results
lca_results# Generate a HTML-TABLE and show it in the RSTUDIO-Viewer (for copy & paste)
view_kable <- function(x, ...){
tab <- paste(capture.output(kable(x, ...)), collapse = '\n')
tf <- tempfile(fileext = ".html")
writeLines(tab, tf)
rstudioapi::viewer(tf)
}
view_kable(lca_results, format = 'html', table.attr = "class=nofluid")
# Another possibility which is prettier and easier to do:
ztable::ztable(lca_results)BICが低いものを選択する。
# plot 1
# Order categories of results$model in order of appearance
results$Model <- factor(results$Model)
#convert to long format
results2<-tidyr::gather(results,Kriterium,Guete,4:7)
results2#plot
fit.plot<-ggplot(results2) +
geom_point(aes(x=Model,y=Guete),size=3) +
geom_line(aes(Model, Guete, group = 1)) +
theme_bw()+
labs(x = "", y="", title = "") +
facet_grid(Kriterium ~. ,scales = "free") +
theme_bw(base_size = 16, base_family = "") +
theme(panel.grid.major.x = element_blank() ,
panel.grid.major.y = element_line(colour="grey", size=0.5),
legend.title = element_text(size = 16, face = 'bold'),
axis.text = element_text(size = 16),
axis.title = element_text(size = 16),
legend.text= element_text(size=16),
axis.line = element_line(colour = "black")) # Achsen etwas dicker## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# save 650 x 800
fit.plot
BICをプロットして、変曲点またはプラトー点を求めることが有用である（エルボープロット）。この場合はModel2か３が選択されるだろう。
各クラスの人口シェアを調べる 各クラスの人口シェアに興味がある場合は、次のようにして取得できる：
round(colMeans(lc3$posterior)*100,2)## [1] 27.79 29.08 43.13あるいは、クラスメンバーズの推定人数を検査する：
table(lc3$predclass)##
## 1 2 3
## 486 507 792round(prop.table(table(lc3$predclass)),4)*100##
## 1 2 3
## 27.23 28.40 44.37潜在クラスの順序付け
潜在クラスは順序付けされていないので、どの潜在クラスが1番、2番、3番…になるかは任意である。潜在クラスは名目的であるはずなので、1つのクラスが最初になる理由はない。必要であれば、潜在クラスを並べ替える。潜在クラスを手動で並べ替える関数： poLCA.reorder()
まず、LCAモデルを実行し、開始値を抽出し、再度モデルを実行しますが、今度は手動で順序を設定する。
#extract starting values from our previous best model (with 3 classes)
probs.start<-lc3$probs.start
#re-run the model, this time with "graphs=TRUE"
lc<-poLCA(f, mydata, nclass=3, probs.start=probs.start,graphs=TRUE, na.rm=TRUE, maxiter=3000)
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $MORALG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5224 0.4109 0.0390 0.0277
## class 2: 0.1851 0.3668 0.2350 0.2131
## class 3: 0.1013 0.5990 0.2552 0.0446
##
## $CARESG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3970 0.4605 0.0895 0.0529
## class 2: 0.0879 0.2546 0.3524 0.3051
## class 3: 0.0227 0.5090 0.3819 0.0864
##
## $KNOWG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5878 0.3543 0.0266 0.0312
## class 2: 0.2542 0.4253 0.2382 0.0822
## class 3: 0.0429 0.8058 0.1408 0.0105
##
## $LEADG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3747 0.4772 0.1085 0.0396
## class 2: 0.0747 0.2582 0.3856 0.2816
## class 3: 0.0207 0.4886 0.4268 0.0639
##
## $DISHONG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0426 0.0487 0.3371 0.5717
## class 2: 0.2007 0.3060 0.2774 0.2158
## class 3: 0.0519 0.2182 0.4998 0.2301
##
## $INTELG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5763 0.3698 0.0197 0.0341
## class 2: 0.2936 0.4407 0.1852 0.0805
## class 3: 0.0600 0.7953 0.1317 0.0130
##
## $MORALB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0966 0.3795 0.3638 0.1601
## class 2: 0.5797 0.3787 0.0248 0.0168
## class 3: 0.0642 0.7137 0.2164 0.0057
##
## $CARESB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0162 0.1125 0.4045 0.4668
## class 2: 0.3483 0.5362 0.0836 0.0319
## class 3: 0.0093 0.4847 0.4380 0.0679
##
## $KNOWB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0682 0.3115 0.3828 0.2374
## class 2: 0.4970 0.4899 0.0132 0.0000
## class 3: 0.0191 0.7583 0.2182 0.0045
##
## $LEADB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0236 0.2735 0.4524 0.2505
## class 2: 0.5106 0.4648 0.0187 0.0058
## class 3: 0.0537 0.7119 0.2217 0.0127
##
## $DISHONB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0862 0.3096 0.4111 0.1931
## class 2: 0.0365 0.0740 0.1993 0.6902
## class 3: 0.0134 0.1471 0.5683 0.2711
##
## $INTELB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.1179 0.3391 0.3440 0.1990
## class 2: 0.5284 0.4656 0.0060 0.0000
## class 3: 0.0435 0.7741 0.1796 0.0028
##
## Estimated class population shares
## 0.3198 0.2608 0.4194
##
## Predicted class memberships (by modal posterior prob.)
## 0.3166 0.2578 0.4256
##
## =========================================================
## Fit for 3 latent classes:
## =========================================================
## number of observations: 1311
## number of estimated parameters: 110
## residual degrees of freedom: 1201
## maximum log-likelihood: -16714.66
##
## AIC(3): 33649.32
## BIC(3): 34218.96
## G^2(3): 15016.7 (Likelihood ratio/deviance statistic)
## X^2(3): 16748681883 (Chi-square goodness of fit)
## # If you don´t like the order, reorder them (here: Class 1 stays 1, Class 3 becomes 2, Class 2 becomes 1)
new.probs.start<-poLCA.reorder(probs.start, c(1,3,2))
#run polca with adjusted ordering
lc<-poLCA(f, mydata, nclass=3, probs.start=new.probs.start,graphs=TRUE, na.rm=TRUE)
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $MORALG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5224 0.4109 0.0390 0.0277
## class 2: 0.1013 0.5990 0.2552 0.0446
## class 3: 0.1851 0.3668 0.2350 0.2131
##
## $CARESG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3970 0.4605 0.0895 0.0529
## class 2: 0.0227 0.5090 0.3819 0.0864
## class 3: 0.0879 0.2546 0.3524 0.3051
##
## $KNOWG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5878 0.3543 0.0266 0.0312
## class 2: 0.0429 0.8058 0.1408 0.0105
## class 3: 0.2542 0.4253 0.2382 0.0822
##
## $LEADG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3747 0.4772 0.1085 0.0396
## class 2: 0.0207 0.4886 0.4268 0.0639
## class 3: 0.0747 0.2582 0.3856 0.2816
##
## $DISHONG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0426 0.0487 0.3371 0.5717
## class 2: 0.0519 0.2182 0.4998 0.2301
## class 3: 0.2007 0.3060 0.2774 0.2158
##
## $INTELG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5763 0.3698 0.0197 0.0341
## class 2: 0.0600 0.7953 0.1317 0.0130
## class 3: 0.2936 0.4407 0.1852 0.0805
##
## $MORALB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0966 0.3795 0.3638 0.1601
## class 2: 0.0642 0.7137 0.2164 0.0057
## class 3: 0.5797 0.3787 0.0248 0.0168
##
## $CARESB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0162 0.1125 0.4045 0.4668
## class 2: 0.0093 0.4847 0.4380 0.0679
## class 3: 0.3483 0.5362 0.0836 0.0319
##
## $KNOWB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0682 0.3115 0.3828 0.2374
## class 2: 0.0191 0.7583 0.2182 0.0045
## class 3: 0.4970 0.4899 0.0132 0.0000
##
## $LEADB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0236 0.2735 0.4524 0.2505
## class 2: 0.0537 0.7119 0.2217 0.0127
## class 3: 0.5106 0.4648 0.0187 0.0058
##
## $DISHONB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0862 0.3096 0.4111 0.1931
## class 2: 0.0134 0.1471 0.5683 0.2711
## class 3: 0.0365 0.0740 0.1993 0.6902
##
## $INTELB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.1179 0.3391 0.3440 0.1990
## class 2: 0.0435 0.7741 0.1796 0.0028
## class 3: 0.5284 0.4656 0.0060 0.0000
##
## Estimated class population shares
## 0.3198 0.4194 0.2608
##
## Predicted class memberships (by modal posterior prob.)
## 0.3166 0.4256 0.2578
##
## =========================================================
## Fit for 3 latent classes:
## =========================================================
## number of observations: 1311
## number of estimated parameters: 110
## residual degrees of freedom: 1201
## maximum log-likelihood: -16714.66
##
## AIC(3): 33649.32
## BIC(3): 34218.96
## G^2(3): 15016.7 (Likelihood ratio/deviance statistic)
## X^2(3): 16748681883 (Chi-square goodness of fit)
## lc## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $MORALG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5224 0.4109 0.0390 0.0277
## class 2: 0.1013 0.5990 0.2552 0.0446
## class 3: 0.1851 0.3668 0.2350 0.2131
##
## $CARESG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3970 0.4605 0.0895 0.0529
## class 2: 0.0227 0.5090 0.3819 0.0864
## class 3: 0.0879 0.2546 0.3524 0.3051
##
## $KNOWG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5878 0.3543 0.0266 0.0312
## class 2: 0.0429 0.8058 0.1408 0.0105
## class 3: 0.2542 0.4253 0.2382 0.0822
##
## $LEADG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.3747 0.4772 0.1085 0.0396
## class 2: 0.0207 0.4886 0.4268 0.0639
## class 3: 0.0747 0.2582 0.3856 0.2816
##
## $DISHONG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0426 0.0487 0.3371 0.5717
## class 2: 0.0519 0.2182 0.4998 0.2301
## class 3: 0.2007 0.3060 0.2774 0.2158
##
## $INTELG
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.5763 0.3698 0.0197 0.0341
## class 2: 0.0600 0.7953 0.1317 0.0130
## class 3: 0.2936 0.4407 0.1852 0.0805
##
## $MORALB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0966 0.3795 0.3638 0.1601
## class 2: 0.0642 0.7137 0.2164 0.0057
## class 3: 0.5797 0.3787 0.0248 0.0168
##
## $CARESB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0162 0.1125 0.4045 0.4668
## class 2: 0.0093 0.4847 0.4380 0.0679
## class 3: 0.3483 0.5362 0.0836 0.0319
##
## $KNOWB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0682 0.3115 0.3828 0.2374
## class 2: 0.0191 0.7583 0.2182 0.0045
## class 3: 0.4970 0.4899 0.0132 0.0000
##
## $LEADB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0236 0.2735 0.4524 0.2505
## class 2: 0.0537 0.7119 0.2217 0.0127
## class 3: 0.5106 0.4648 0.0187 0.0058
##
## $DISHONB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.0862 0.3096 0.4111 0.1931
## class 2: 0.0134 0.1471 0.5683 0.2711
## class 3: 0.0365 0.0740 0.1993 0.6902
##
## $INTELB
## 1 Extremely well 2 Quite well 3 Not too well 4 Not well at all
## class 1: 0.1179 0.3391 0.3440 0.1990
## class 2: 0.0435 0.7741 0.1796 0.0028
## class 3: 0.5284 0.4656 0.0060 0.0000
##
## Estimated class population shares
## 0.3198 0.4194 0.2608
##
## Predicted class memberships (by modal posterior prob.)
## 0.3166 0.4256 0.2578
##
## =========================================================
## Fit for 3 latent classes:
## =========================================================
## number of observations: 1311
## number of estimated parameters: 110
## residual degrees of freedom: 1201
## maximum log-likelihood: -16714.66
##
## AIC(3): 33649.32
## BIC(3): 34218.96
## G^2(3): 15016.7 (Likelihood ratio/deviance statistic)
## X^2(3): 16748681883 (Chi-square goodness of fit)
## 上の立体図はいまいちなので、別の図に作り変える。
lcmodel <- reshape2::melt(lc$probs, level=2)
zp1 <- ggplot(lcmodel,aes(x = L2, y = value, fill = Var2))+
geom_bar(stat = "identity", position = "stack")+
facet_grid(Var1 ~ .) +
scale_fill_brewer(type="seq", palette="Greys") +theme_bw()+
labs(x = "Item",y="Proportion", fill ="Answer")+
theme(axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
panel.grid.major.y=element_blank())+
guides(fill = guide_legend(reverse=TRUE))
print(zp1)
3クラスでいくので 以下でプロットしていく。
参考資料
全身性硬化症における肺動脈圧の潜在的軌跡モデリング：後方視的コホート研究 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9806097/
敗血症のLCA、K means https://pubmed.ncbi.nlm.nih.gov/31104070/
肥満手術の好みのLCA https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6439857/
方法 https://pubmed.ncbi.nlm.nih.gov/20192788/
Youtube https://cjvanlissa.github.io/visualizing_mixtures/ecdp_mixtures.html#1 https://youtu.be/6DStckg8NhE