公式のマニュアルに記載された例を紹介する。
データ
データdata_SRHS_longはミシガン大学が実施したHealth and Retirement Studyに由来する自己申告健康状態に関するデータセットである。
t: 時点。
id ID。
gender: 男が1、女が2。
race: 「白人」が1、「黒人」が2、「その他」が3。
educational: 教育レベル。1は「高校」、2は「一般教育修了証」、3は「高校卒業」、4は「ある程度の大学」、5は「大学以上」。
age: 年齢。
srhs: 自己申告した健康状態。「優れている」を1、「非常に良い」を2、「良い」を3、「まあまあ」を4、「悪い」を5。
データ名にもあるが、このデータはロングデータと呼ばれている形である。
data("data_SRHS_long")
SRHS <- data_SRHS_long[1:2400,]
# Categories rescaled to vary from 0 (“poor”) to 4 (“excellent”)
SRHS$srhs <- 5 - SRHS$srhs
head(SRHS)
データの頭部分。
t id gender race education age srhs 1 1 1 1 1 3 56 1 2 2 1 1 1 3 58 1 3 3 1 1 1 3 60 2 4 4 1 1 1 3 62 2 5 5 1 1 1 3 64 1 6 6 1 1 1 3 66 2
コード
out <- lmest(responsesFormula = srhs ~ NULL,index = c("id","t"),
data = SRHS,k = 3,start = 1, modBasic = 1, seed = 123)
out
summary(out)
引数の解説。
index: 1 つ目はID、2つ目は時点
k: 潜在状態の数(潜在クラスの数)
modBasic=1: 均質な移行確率
seed: 乱数のシード
引数の最初にあるのはモデルの指定である。
responsesFormula・測定モデル
responsesFormula = y1 + y2 ~ NULL
共変量と2つのレスポンス(y1とy2)を含まないLMモデルが指定。responsesFormula = NULL
共変量なしのLMモデルを推定。id列とtime列を除くデータ中のすべての列がレスポンスとして用いられる。ResponsesFormula = y1 + y2 ~ x1 + x2
2つのレスポンス（y1とy2）と2つの共変量（x1とx2）を持つLMモデルが指定。
latentFormulal・潜在モデル
responsesFormula = y1 + y2 ~ NULL
latentFormula = ~ x1 + x2 | x3 + x4
2つのレスポンス(y1とy2)と初期確率に影響を与える2つの共変量(x1とx2)と移行確率に影響を与える残りの2つの共変量(x3とx4)を持つLMモデルが指定。responsesFormula = y1 + y2 ~ NULL
latentFormula = ~ 1 | x1 + x2
(orlatentFormula = ~ NULL | x1 + x2)
共変量は移行確率にのみ影響し、切片は初期確率に指定。responsesFormula = y1 + y2 ~ NULLlatentFormula = ~ x1 + x2
初期確率と移行確率の両方に影響を与える2つの共変量(x1とx2)を持つLMモデルを指定。responsesFormula = y1 + y2 ~ NULL
latentFormula = ~ NULL | NULL
(orlatentFormula = ~ 1 | 1)
初期確率と移行確率の切片のみを持つLMモデルを指定。
今回は上から2番目のid列とtime列を除くデータ中のすべての列を潜在変数に組み込むモデルになっている。
結果
Convergence info:
LogLik np k AIC BIC n TT
-2796.989 20 3 5633.979 5708.054 300 8
Coefficients:
Initial probabilities:
est_piv
[1,] 0.5608
[2,] 0.3119
[3,] 0.1273
Transition probabilities:
state
state 1 2 3
1 0.9318 0.0608 0.0074
2 0.0000 0.9464 0.0536
3 0.0040 0.0200 0.9760
Conditional response probabilities:
, , item = 1
state
category 1 2 3
0 0.0028 0.0049 0.2924
1 0.0000 0.1208 0.6178
2 0.0832 0.6093 0.0860
3 0.4979 0.2438 0.0038
4 0.4160 0.0212 0.0000
Transition probabilitiesをみると、ほとんどのクラスが維持されていることがわかる。
余談
連続変数のageを入れて潜在クラス分析を走らせるのは無理があるので、ageを抜いた形で潜在クラス分析をR上で行うと次のようになる。1波のみのデータで実施。
library(poLCA) d1 <- subset(SRHS, SRHS$t==1) # 1波のみ抽出 f <- cbind(gender, race, education, srhs)~1 d1$gender <- as.factor(d1$gender) d1$race <- as.factor(d1$race) d1$education <- as.factor(d1$education) d1$srhs <- as.factor(d1$srhs) lca.res <- poLCA(f,d1,nclass=3,maxiter=3000,nrep=100) lca.res
