2014-05-20

Latent class mixed effect model

R lcmm

参考はこちら

library('lcmm');

data(data_hlme)
head(data_hlme)

データは経時的に観察するYと観察時刻Timeが各行で異なり、それ以外の共変量は個人の番号IDと同様観測回数ごとに繰り返して入力されている形式。以下では、ID=(1,2)の２人分のデータが見えている

> head(data_hlme)
  ID      Y  Time X1 X2    X3
1  1 22.820 1.831  1  1 1.163
2  1 20.622 3.831  1  1 1.163
3  1 18.340 5.831  1  1 1.163
4  2 24.390 0.089  1  0 0.411
5  2 21.015 2.089  1  0 0.411
6  2 18.188 4.089  1  0 0.411

m<-lcmm(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,
subject='ID',ng=2,data=data_hlme,link="linear")
summary(m)

> summary(m)
# 説明文のタイトル
> summary(m)
General latent class mixed model 
     fitted by maximum likelihood method 
# 実行したコマンドの再掲 
lcmm(fixed = Y ~ Time * X1, mixture = ~Time, random = ~Time, 
    subject = "ID", classmb = ~X2 + X3, ng = 2, link = "linear", 
    data = data_hlme)
# モデルの説明 
Statistical Model: 
     Dataset: data_hlme # データセット名
     Number of subjects: 100 # IDで識別される個体数
     Number of observations: 326 # 観測数、データの行数
     Number of latent classes: 2 # 推定クラス数
     Number of parameters: 13  # モデルを構成する変数の数
     Link function: linear  # 時間変数・共変量を使って推定した値を従属変数と結びつけるために使う関数
# 推定がどんな具合だったか
Iteration process: 
     Convergence criteria satisfied # 収束したか
     Number of iterations:  17 # 何回、回して収束したか
     Convergence criteria: parameters= 1e-05 # 収束判定条件のうち、パラメタの変化幅の判定基準
                         : likelihood= 1.4e-05 # 尤度の変化幅の判定基準
                         : second derivatives= 7.7e-12 # 収束アルゴリズムで用いた二次微分の変化幅の判定基準
# 最終推定のあてはまりを数値で表すと…
Goodness-of-fit statistics: 
     maximum log-likelihood: -773.82  
     AIC: 1573.64  
     BIC: 1607.51  
 
# 最優推定の内訳 
Maximum Likelihood Estimates: 

# "classmb=~X2+X3"と指定したようにX2とX3はクラス分けに影響するパラメタとして指定されている
Fixed effects in the class-membership model:
(the class of reference is the last class) 

                     coef      Se   Wald p-value
intercept class1 -0.10518 0.33396 -0.315 0.75280
X2 class1         0.73939 0.48979  1.510 0.13114
X3 class1         0.06818 0.19284  0.354 0.72366

# 時間依存の要素(時間に沿って増える分)
# 前の要素にclass1の基礎値が入っているので、こちらではclass1は0
# クラス別の係数が出ているのはTimeとその切片。それは"mixture~Time"としているから
# 切片が入るのは、切片を入れるのは回帰式の書き方のルールがそうなっているから。切片を抜きたければ"Y~Time -1"のようにするのは、「回帰式」指定で共通するルール
# ただし、lcmm()では-1として取り除かないこと
# クラスは考慮しないが、"Y~Time*X1"と指定してあるので、それに応じた係数が推定されている切片(これはデフォルトで入ることになっている。
# "Time*X1"と回帰式で指定したら、切片+Time+X1+Time*X1のように、その配下の全組み合わせを盛り込んだことにするのも、「回帰式」の書き方のルール
# ただし、切片 + Time + X1 + Time*X1の切片は、クラス別の切片に吸収されている
Fixed effects in the longitudinal model:

                                     coef      Se    Wald p-value
intercept class1 (not estimated)  0.00000                        
intercept class2                 -8.97287 0.77084 -11.640 0.00000
Time class1                       0.13534 0.17844   0.758 0.44817
Time class2                      -1.03737 0.21216  -4.890 0.00000
X1                                2.53261 0.68591   3.692 0.00022
Time:X1                          -0.29719 0.20777  -1.430 0.15262

# "random ~ Time"と、指定したのでTimeにだけ、ばらつきに関する推定が返っている
Variance-covariance matrix of the random-effects:
          intercept      Time
intercept 4.8025533          
Time      0.3839559 0.8465551

# 固定のばらつき項
Residual standard error (not estimated) = 1

# 上記で出した値をさらにlink関数に渡さないといけない。そのlink関数はlinearなので、上記で推定した値xに対して、以下のLinear1 + x*Linear2とする…？？ちょっと不安
Parameters of the link function:

                         coef      Se   Wald p-value
Linear 1 (intercept) 29.39251 0.54904 53.534       0
Linear 2 (std err)    0.96808 0.06066 15.959       0