Finite Mixture Models の誤用に注意
太郎丸 博
2015年8月30日
library(knitr)
opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)Finite Mixture Models とか Latent Class Regression とか呼ばれる統計モデルを応用している研究を最近見かけたのだが、なんか使い方を間違っている気がしたので、少し考えたことをメモしておく。なお私はこのモデルについては、R の flexmix パッケージについている vignettes の文書 を読んだ程度の知識なので、以下の記述には間違いがあるかもしれないことに留意されたい。
1 Finite Mixture Models とは？
通常、回帰分析などで被説明変数の分布を予測する場合、すべてのケースが同じ回帰式にしたがって分布していると仮定する。例えば、 \[ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i \] という回帰式をたてる場合、すべてのケースがこの回帰式で示されるようなデータ生成プロセス (data generating process) に従うと考える。しかし、グループによって切片や傾きが異なる場合もある。その場合、そのグループを示すダミー変数と、\(X_i\)の交互作用効果を仮定したモデルをたてればよい。例えば、グループが2つあり、そのダミー変数を\(G_i\)とするならば、 \[ Y_i = \beta_0 + \beta_1 X_i + \beta_2 G_i + \beta_3 G_i X_i + \epsilon_i \] というモデルをたてればよい。ところが、\(G_i\)のような変数が測定されておらず、どのケースがどちらのグループに属しているのか、予めわからない場合もある。例えば、年齢とともに賃金が上がる労働者と、年齢が上がっても賃金は上がらない労働者がいるが、散布図を見ただけではそれらを明確に識別するのは、困難な場合もある。下の図1 は、そのような架空のデータである。
N <- 480
x0 <- rep(20:59, N/80)
set.seed(1986)
y1 <- 200 + rnorm(N/2, sd=50)
set.seed(1987)
y2 <- 9 * x0 + rnorm(N/2, sd=50)
y <- c(y1, y2)
x <- c(x0, x0)
plot(x, y)
図1 傾きと切片の異なる2つの潜在クラスからなるデータの散布図（架空）
これをみると、x が 40 をこえたあたりから、2つのグループに分岐していることがわかるが、40以下では両者が混じっていて、目測では恣意的にしか 2つの潜在クラスをわけられない。それゆえ、それぞれのグループの傾きと切片も正確な推定は困難である。
もちろん、男性、正規雇用の場合は年齢と賃金の関連は強く、女性、非正規雇用の場合は弱いといった予測は成り立つが、例外も多い。そのような場合、2つの観察されない潜在クラスが存在し、その潜在クラスによって、回帰係数が異なると仮定して、潜在クラスの比率とそれぞれの潜在クラスにおける回帰係数を予測するということが考えられる。このようなモデルを Finite Mixture Models とか Latent Class Regression と呼ぶようである。推定は最尤法だが初期値によっては local maxima の1つに到達してしまい、最尤値を得られない場合もあるので、初期値を変えて複数回推定することが推奨される。以下の例では各モデルそれぞれ１回しか推定していないが、それはモデルが単純で正解がわかっているからであって（モデルが単純だと最尤値を得やすい気がする）、実際の応用では複数回推定したほうが良かろう。
上の図1 のデータに関して潜在クラスの数を 2 と仮定して Finite Mixture Models を推定した結果が以下である。
library(flexmix)
m1 <- flexmix(y ~ x, k = 2)
summary(m1)##
## Call:
## flexmix(formula = y ~ x, k = 2)
##
## prior size post>0 ratio
## Comp.1 0.494 247 405 0.610
## Comp.2 0.506 233 414 0.563
##
## 'log Lik.' -2766.348 (df=7)
## AIC: 5546.697 BIC: 5575.913summary(refit(m1))## $Comp.1
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 190.00356 15.88686 11.9598 <2e-16 ***
## x 0.32251 0.35361 0.9121 0.3617
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.94712 16.41496 0.2405 0.81
## x 8.88910 0.37421 23.7542 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1上の出力の Comp.1, Comp.2 が２つの潜在クラスに付けられた名前であり（Finite Mixture Models では潜在クラスのことを Component と呼ぶようである）、最初の出力 (summary(m1)) の “prior” が各潜在クラスの推定比率、“size” はそれぞれの推定人数、“post>0” ある確率（デフォルトでは 0.0001 以上の確率で各潜在クラスに属するサンプルの数、“ratio” は “size” を “post>0” で割った値である。その下の summary(refit(m1)) の出力は２つの潜在クラスにおける切片と傾きの推定値とその検定結果が示されている。
このデータは私が乱数を回して生成した架空データなので、パラメータの「真の値」がわかっている。正解は、2つの潜在クラスの比率は 50% ずつで、一方の潜在クラスの切片は 200、傾きは 0、もう一方の潜在クラスは切片は 0 で傾きは 9 であるから、上の推定結果は、概ね正しいものになっていることがわかるだろう。
このように、複数の観察されないグループがあり、そのようなグループによって回帰式の係数が異なると考えられる場合、このモデルは役に立つと考えられる。
2 Finite Mixture Models の誤用
構造方程式モデリングや潜在クラス分析でも同じだが、潜在変数を使うモデルは柔軟性が高く、フィッティングが同程度の様々なモデルを考えることができるし、真のモデルとはかけ離れたモデルが高いフィッティングを示すこともある。それゆえ、理論や先行研究の知見、探索的分析の結果を十分に踏まえ、慎重にモデルを選択することが求められる。逆にいえば、安易に用いればとんでもない誤りを犯すことになる。いわゆる統計的捏造 (statistical artifacts) をしてしまうことになりかねない。
2.1 Nonlinearity: 非線形の関係なのに線形の関係を仮定してしまった場合
例えば、以下のような 2次曲線的な関係をしているデータに、機械的に線形モデルを仮定して Finite Mixture Model を当てはめると、とんでもない間違いを犯すことになる。
set.seed(1986)
y <- 700 -(x - 40)^2 + rnorm(N, sd=50)
plot(x, y)
図2 すべてのサンプルが２次曲線にしたがって分布しているデータ（架空）
m1 <- flexmix(y ~ x, k = 1)
m2 <- update(m1, k = 2)
m3 <- update(m1, k = 3)
m4 <- update(m1, k = 4)
m5 <- update(m1, k = 5)
BIC(m1, m2, m3, m4, m5)## df BIC
## m1 3 6018.105
## m2 7 5836.281
## m3 11 5780.477
## m4 15 5767.409
## m5 19 5772.767summary(m4)##
## Call:
## flexmix(formula = y ~ x, k = 4)
##
## prior size post>0 ratio
## Comp.1 0.309 154 307 0.502
## Comp.2 0.269 108 293 0.369
## Comp.3 0.221 124 341 0.364
## Comp.4 0.202 94 350 0.269
##
## 'log Lik.' -2837.401 (df=15)
## AIC: 5704.802 BIC: 5767.409summary(refit(m4))## $Comp.1
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1466.34525 48.60014 30.172 < 2.2e-16 ***
## x -18.29527 0.92683 -19.739 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -93.6974 31.7236 -2.9536 0.003141 **
## x 21.5870 1.1716 18.4258 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.3
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 781.58421 29.92911 26.1145 < 2.2e-16 ***
## x -3.33690 0.72389 -4.6097 4.033e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.4
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 317.74091 33.99267 9.3473 < 2.2e-16 ***
## x 8.24233 0.98103 8.4018 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1このデータの生成プロセスに潜在クラスなど存在しておらず、すべてのサンプルが、二次曲線にしたがって分布するように作られている。ところが、安易に線形関係を仮定して、Finite Mixture Models を推定すると、BIC を基準にしても 4つも潜在クラスがあるという推定結果になってしまう。ちなみに、「正解」は以下の通り。
m6 <- flexmix(y ~ x +I(x^2), k=1)
summary(m6)##
## Call:
## flexmix(formula = y ~ x + I(x^2), k = 1)
##
## prior size post>0 ratio
## Comp.1 1 480 480 1
##
## 'log Lik.' -2564.259 (df=4)
## AIC: 5136.518 BIC: 5153.213summary(refit(m6))## $Comp.1
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -848.662411 28.857730 -29.409 < 2.2e-16 ***
## x 77.199684 1.544275 49.991 < 2.2e-16 ***
## I(x^2) -0.962725 0.019383 -49.668 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1この１クラス二次曲線モデルは、上で最小の BIC を示していた 4クラス・モデルよりも 600以上小さな BIC を示しており、ずっと適合度が高いことがわかる。
この例から学ぶべき教訓は Finite Mixture Models を使う場合には、非線形関係について通常の回帰分析をする時以上に気をつけるべきだということであろう。
2.2 Distribution of Response: 被説明変数は正規分布していないのに正規分布していると仮定してしまった場合
また、被説明変数の残差が正規分布していない（つまり、被説明変数が条件付きで正規分布していない）のに、正規分布を仮定して Finite Mixture Models を推定すると、やはり存在してもいない潜在クラスを捏造してしまうことになる。
以下の図3 は、すべてのサンプルに関して、従属変数がポアソン分布している例であるが、これに正規分布を仮定して Finite Mixture Models を推定すると、やはり複数の潜在クラスを仮定したモデルのほうがフィッティングがいいことになってしまう。ただし、この例では3つ以上の潜在クラスを推定しようとするとエラーが出ることが多かった。以下の例でも潜在クラスの数を 4 と 5 に仮定したモデルは結果がおかしかった。
set.seed(1986)
y <- rpois(N, exp(scale(x)+0.5))
plot(x, jitter(y))
図3 被説明変数がポアソン分布している場合（架空データ）
m1 <- flexmix(y ~ x, k = 1)
m2 <- update(m1, k = 2)
m3 <- update(m1, k = 3)
m4 <- update(m1, k = 4)
BIC(m1, m2, m3, m4)## df BIC
## m1 3 2009.303
## m2 7 1829.599
## m3 11 1802.264
## m4 15 1804.468summary(m3)##
## Call:
## flexmix(formula = y ~ x, k = 3)
##
## prior size post>0 ratio
## Comp.1 0.233 77 480 0.160
## Comp.2 0.470 230 451 0.510
## Comp.3 0.297 173 350 0.494
##
## 'log Lik.' -867.1763 (df=11)
## AIC: 1756.353 BIC: 1802.264summary(refit(m3))## $Comp.1
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.79335 0.98646 -6.8866 5.715e-12 ***
## x 0.27667 0.02107 13.1308 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.8227921 0.3361448 -8.3975 < 2.2e-16 ***
## x 0.1332414 0.0092191 14.4528 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.3
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5653586 0.2914530 -5.3709 7.835e-08 ***
## x 0.0625399 0.0087994 7.1073 1.183e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1このデータに関するモデルの正解は、以下の通り。
m5 <- flexmix(y ~ x, k = 1, model = FLXMRglm(family="poisson"))
m6 <- update(m5, k= 2)
BIC(m5, m6)## df BIC
## m5 2 1553.563
## m6 5 1571.604summary(m5)##
## Call:
## flexmix(formula = y ~ x, k = 1, model = FLXMRglm(family = "poisson"))
##
## prior size post>0 ratio
## Comp.1 1 480 480 1
##
## 'log Lik.' -770.6079 (df=2)
## AIC: 1545.216 BIC: 1553.563ポアソン回帰を仮定すると、１グループで（つまり、潜在クラスを仮定しないモデルで）、 BIC が最小になっていることがわかる。ただし、正規分布を仮定したモデルとポアソン分布を仮定したモデルで、対数尤度を比較していいのか、ちょっとあやふやである（どちらも \(\mathrm{BIC} = -2 \mathrm{LogLik } + 2 \mathrm{NofParameters}\) で計算しているので比較していい気がするのだが）。
この例から学ぶべき教訓は、深く考えずに安易に正規分布を仮定していると、ありもしない潜在クラスを知らないうちに捏造してしまう事があるということである。Finite Mixture Models を使うときはいつも以上に残差の分布をまじめにチェックした方がいいだろう。
3 Discussion
このような Finite Mixture Models の誤用は、散布図を見たり、残差の分布をチェックするなどすれば避けられるはずである。というか、通常の回帰分析をするときにもとうぜんチェックすべきことなのだが、これをまじめにやっている社会学者がどれぐらいいるのだろうか。いずれにせよ、闇雲に Finite Mixture Models をあてはめて、潜在クラスの数しか検討しないのであれば失敗する可能性が高いので、モデルや残差の分布について慎重にチェックするべきだろう。
Finite Mixture Models はロジスティック回帰でもポアソン回帰でも可能なので、適当なリンク関数と分布の仮定をおければ有効活用はできるはずである。また、いろいろモデルの拡張もできるようなので、詳しくは flexmix パッケージについている解説論文を参照されたい。
https://cran.r-project.org/web/packages/flexmix/vignettes/mixture-regressions.pdf
https://cran.r-project.org/web/packages/flexmix/vignettes/regression-examples.pdf