一般化線形混合モデル

これまで線形モデル(LM) ではモデルの作成に data.frameとformula、model が必要だということを学びました。さらに一般化線形モデル(GLM) では GLMを組むには「観測できる値が従う確率分布」をfamilyで指定すること、また結果を解釈をしやすくするためにはリンク関数の理解が必要であることを学びました。以下のコードの問、Q.1–3を見て復習しましょう。

# LM
# Q1. data.frame, formula, model はどれでしょうか。
f = Reaction ~ Days
m = lm(f, sleepstudy)

# GLM
# Q2. なぜ family=binomial と指定しているのでしょうか。
f = Correct ~ Days
m = glm(f, sleepstudy, family=binomial)

# Q3. Days と Estimate(-0.16017) の値は Correct とどういう関係でしょうか。
summary(m)
> Coefficients:
>              Estimate Std. Error z value Pr(>|z|)  
> (Intercept)  0.83782    0.29216   2.868  0.00413 **  
> Days        -0.16017    0.05442  -2.943  0.00325 **

Q1はOKとして、上で binomial を指定する理由は LM の Reaction と異なり、 Correct は0/1で二項分布に従うからでした。また結果の解釈には「リンク関数」が必要でした。 binomial の場合、logit関数がリンク関数だったことを思い出しましょう。解釈に不安が残る場合は前回の内容を再度実行して p/(1-p) = exp(ax)*exp(b) が何を示すのかを思い出しましょう。

前回までのところまででLMを一般化させる方向は達成したので、今回は「混合」の部分を足します。とはいえ前回の内容に比べればだいぶ簡単で、 Rにおける縦線 | (バーと読みます) の役割を確認するだけです¹。まずは実データに触れながら (i) 構造を表現する記号を確認します。そして(ii) 構造的なノイズを扱うメリットを抑え、最後に (iii) glmer を用いた一般化線形混合モデルを作成、分析していきます。

構造を表現する記号

まずRにおける | の役割を確認するために、とりあえず以下のコードを実行してプロットしてみましょう(Submitボタンを押す)。プロットの結果は右のペインの Plots を押すと別ウィンドウが開かれて、拡大されて Subject ごとに xyplot が作れていることが確認できるかと思います。どうしても拡大できない人は xyplot内の# を外してサブセットを見てみましょう。


    library(lattice); library(lme4)
    set.seed(43)
    # data
    rt = sleepstudy$Reaction
    logistic = function(x) {1 / (1 + exp(-x))}
    scale = function(x) (x-min(x))/(max(x)-min(x))
    ps = logistic(-5*scale(rt)+2)
    xs = as.integer(rbinom(ps, n=length(ps), size=1))
    sleepstudy$Correct = xs


    f = Reaction ~ Days | Subject
    xyplot(f, sleepstudy, type=c("p", "r"),
           # subset=Subject %in% c(370, 371)
           )


    f = Reaction ~ Days | Subject
    xyplot(f, sleepstudy, type=c("p", "r"),
           subset=Subject %in% c(370, 371))


    test_function("xyplot")
    success_msg("Great job!")

前回 xyplot を使った時の formula は f = Reaction ~ Days だったのに対し、今回のプロットは | Subject が足されており、 Subjectごとに描画している点に注意しましょう。これを足すと Subject ごとに図が書かれたことから、 | Subject は「Subject ごとに」という意味があることがわかります²。つまり、被験者ごとにDaysを条件付けてプロットしたいときの formula には | が使えるのです。

条件付けられた formula を受けとって xyplot が描画したのは被験者ごとの傾きや切片です。これをモデリングでもできれば、被験者ごとに傾きや切片をモデリングできそうなことは伝わるでしょうか。ただ図を見渡すと、傾きや切片の傾向に「全体的な効果」はあるものの、被験者ごとのデータには「構造的なノイズ」があってバラバラです。結論を言えば、データを「全体的な効果」と「構造的なノイズ」を明示的に混ぜてモデリングするのが混合モデルになります。

次にデータを「全体的な効果」と「構造的なノイズ」を明示的に混ぜることのメリットを知るため、これらを分けず単に個々の傾きや切片といったパラメータを推定、分析することの問題点を見てみます。

固定効果とランダム効果に分けるメリット

個人ごとの傾きや切片といったパラメータを推定、分析することは最終的な目的ではありません。被験者ごとのパラメータはバラバラなので、それらの報告はとりとめのないものになってしまいます。そうではなく、全体的な効果を「固定効果」、個人差などの構造的なノイズを「ランダム効果」として分離するのが一般化線形混合モデルになります。

まずは効果を分離せず、最初の「個人ごとの傾きや切片といったパラメータを推定、分析すること」がいかに悪い方法であるかを実際に試して確認しましょう。何が問題なのかすぐにわかります。 Rでは lmList という関数を使って Days | Subject つまり被験者ごとのDaysの効果を求められます。まずは先ほどと同様、とりあえず以下のコードを実行してみましょう。


    library(lattice); library(lme4)
    set.seed(43)
    # data
    rt = sleepstudy$Reaction
    logistic = function(x) {1 / (1 + exp(-x))}
    scale = function(x) (x-min(x))/(max(x)-min(x))
    ps = logistic(-5*scale(rt)+2)
    xs = as.integer(rbinom(ps, n=length(ps), size=1))
    sleepstudy$Correct = xs


    f = Reaction ~ Days | Subject # 前回の範囲
    models = lmList(f, sleepstudy)
    summary(models)


    f = Reaction ~ Days | Subject
    models = lmList(f, sleepstudy)
    summary(models)


    test_function("summary")
    success_msg("Great job!")

上を実行すると、以下のような被験者ごと(308, 309, …)の切片(Intercept)と Daysの効果が求められたと思います。

> Coefficients:
>    (Intercept) 
>     Estimate Std. Error  t value     Pr(>|t|)
> 308 244.1927   15.04169 16.23439 2.419368e-34
> 309 205.0549   15.04169 13.63244 1.067180e-27
> 310 203.4842   15.04169 13.52802 1.993900e-27
> :
>    Days 
>      Estimate Std. Error    t value     Pr(>|t|)
> 308 21.764702   2.817566  7.7246464 1.741840e-12
> 309  2.261785   2.817566  0.8027444 4.234454e-01
> 310  6.114899   2.817566  2.1702769 3.162541e-02
> :

確かにこれでReactionに対する「個人ごとの傾きや切片といったパラメータ」を推定できました。でも分析や報告するときは「大体の人によっては効果ありました。」のようなアバウトな物になってしまいそうです。これに対して、もし固定効果 (全体的な効果) とランダム効果 (個人差などの構造的なノイズ) を同時にモデリングできれば、「全体的に効果ありました。個人差は織り込み済みです。」と言えます。

これが「混合」、つまり固定効果とランダム効果を明示的に混ぜる動機になります。最後に (i) 構造的なノイズをモデリングする方法と (ii) Mixed の意味、 (iii) 一般化線形モデルの作り方の3点を次の節で確認しましょう。

glmer を用いた一般化線形混合モデル

まず、Rにおける線形混合モデルには lmerがあり³、これを一般化(generalized)させた一般化線形混合モデルが glmer です。チュートリアルの流れ的には線形モデル、一般化線形モデル、一般化線形混合モデルですが、ネーミング的には線形モデル、線形混合モデル、一般化線形混合モデルになっています。まずは lmerで構造的なノイズをモデリングする方法を確認しましょう。

構造的なノイズをモデリングする方法

線形モデルでは f = Reaction ~ Days で、 f = Reaction ~ Days | Subject だと効果を被験者ごとに推定してしまいます。そうではなく、(一般化)線形混合モデルでは f = Reaction ~ Days + (Days | Subject)のように、括弧でランダム効果を表現します。 Reactionが正規分布に従っていると前提を起き⁴、まずはmodel にlmerを使いましょう。

データは sleepstudyで formula はf = Reaction ~ Days + (Days | Subject) です。後で試す glmerの時はこのformulaを見ないで作れるようになりましょう。なお、この formula の意味は、固定効果 (全体的な効果) Days だけでなく「Subject ごとの Days」もランダム効果 (個人差などの構造的なノイズ) としてモデリングして、という意味です。


    library(lattice); library(lme4);
    set.seed(43)
    # data
    rt = sleepstudy$Reaction
    logistic = function(x) {1 / (1 + exp(-x))}
    scale = function(x) (x-min(x))/(max(x)-min(x))
    ps = logistic(-5*scale(rt)+2)
    xs = as.integer(rbinom(ps, n=length(ps), size=1))
    sleepstudy$Correct = xs


    f = Reaction ~ Days + (Days | Subject)
    m = lmer(f, sleepstudy)
    summary(m)


    f = Reaction ~ Days + (Days | Subject)
    m = lmer(f, sleepstudy)
    summary(m)


    test_function("summary")
    success_msg("Great job!")

さきほどの lmList では個人ごとにモデリングしたためバラバラでしたが、今回のモデルの summary では Fixed effects: と Random effects: に分かれており、全体的な効果は Fixed effects 側が持っています。 Fixed effects の Days が 10.467 であることから、 1日増えると反応時間が10.467増えることがわかります。ただ Random effects の方は本当に個人ごとの Daysの効果のバラつきを考慮できているのでしょうか。

lmList で確認した個人ごとのバラつき、つまり構造的なノイズは Random effects がモデリングしてます。込み入った話になるのですが、この「バラバラ感」は Random effects の Std.Dev で確認できます。 (Intercept) が平均 24.740 ばらついていて、 Days の効果が平均 5.922 ばらついていることからわかります。

つまり結果は以下のように解釈できます。まず、Fixed effects の Days が 10.467 であることから、 1日増えると反応時間が10.467増えることがわかります。ただし、その増え具合は被験者で平均 5.922 程度ばらつく、というものです。切片も同様で、基本は 251.405 だけれども被験者で平均 24.740 程度ばらつくよ、という解釈です。

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.09   24.740       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  251.405      6.825   36.84
Days          10.467      1.546    6.77

次にランダム効果を確認する方法を述べます。結果で報告する内容ではないので飛ばしてもらっても結構ですが、知っておいたほうがモデルの理解は深まると思います。

ランダム効果の実体

Random effects の内訳は作成したモデルを ranef(m) という名前そのまんまな関数に入れるとわかります。以下のコードを実行してみてください。


    library(lattice); library(lme4)
    set.seed(43)
    # data
    rt = sleepstudy$Reaction
    logistic = function(x) {1 / (1 + exp(-x))}
    scale = function(x) (x-min(x))/(max(x)-min(x))
    ps = logistic(-5*scale(rt)+2)
    xs = as.integer(rbinom(ps, n=length(ps), size=1))
    sleepstudy$Correct = xs


    f = Reaction ~ Days + (Days | Subject)
    m = lmer(f, sleepstudy)
    ranef(m)$Subject


    f = Reaction ~ Days + (Days | Subject)
    m = lmer(f, sleepstudy)
    summary(m)


    test_function("summary")
    success_msg("Great job!")

すると、以下のように被験者ごとのランダム効果が出たと思います。

$Subjects
    (Intercept)        Days
308   2.2585654   9.1989719
309 -40.3985770  -8.6197032
310 -38.9602459  -5.4488799
330  23.6904985  -4.8143313
:

例えば308番のデータを用いて式を書くと、 $a \cdot Days + b + (a_{308} \cdot Days + b_{308})$ となることになります。少し式を変えると、 $(a+a_{308}) \cdot Days + (b+b_{308})$ とも言えます。先ほどの固定効果を考えると、 $(10.467+9.1989719) \cdot Days + (251.405+2.2585654)$ が $Reaction_{308}$ の推定値になります。

さて、ここで確認問題です。 ranef(m) は被験者ごとのランダム効果を示し、固定効果(fixefとしましょう)は summaryのFixed effects にあります。モデルの係数はcoef(m) で確認でき、10.467+9.1989719 や 251.405+2.2585654、つまり固定効果とランダム効果を足したものが被験者ごとの係数になります。 fixef+ranef(m)=coef(m) となるのですが、coef(m)-ranef(m)はどうなるでしょうか。

上の答えに納得できれば、固定効果とランダム効果、係数の関係を理解していることになります。最後に lmer を一般化させた glmer でモデル作成、分析をしてみましょう。

一般化線形混合モデルの作り方

前節では Reaction が正規分布に従うという仮定を置きましたが、その仮定は Correct では成り立ちません。 Y が従う分布を一般化し、固定効果とランダム効果を混合してようやく「一般化線形混合モデル」になります。これで大体のデータは分析できるようになります。 Correct と Days の関係をランダム効果ありで作成し、 glmer 関数を使ってモデリングしましょう。なお、family は Correct が従う分布を選択してください。


    library(lattice); library(lme4)
    set.seed(43)
    # data
    rt = sleepstudy$Reaction
    logistic = function(x) {1 / (1 + exp(-x))}
    scale = function(x) (x-min(x))/(max(x)-min(x))
    ps = logistic(-5*scale(rt)+2)
    xs = as.integer(rbinom(ps, n=length(ps), size=1))
    sleepstudy$Correct = xs


    f = 
    m = glmer(f, sleepstudy, family=)
    summary(m)


    f = Correct ~ Days + (Days | Subject)
    m = glmer(f, sleepstudy, family=binomial)
    summary(m)


    test_function("summary")
    success_msg("Great job!")

これで一般化線形混合モデルを作れるようになりました。分析、報告の内容は前回と同じになります。

まとめ

これで統計のチュートリアルは終わりです。線形モデルで要因の効果を検証し、一般化して様々な分布に対応が可能であることを学び、構造的なノイズを組み込めるようになりました。今回は Subjects だけでしたが、利用したアイテムや環境などもランダム効果に組み込めます。

このチュートリアルを通じてみなさんが統計を勉強する時間を減らせたり、読める論文やできる実験が増えたりしたなら幸いです。ただ、今回のチュートリアルでは実験の組み方や要因の選び方に関しては全く触れていません。また、データ収集後のデータを変数の値に落とし込むのも一手間かかるかもしれません。ここは指導教官や先輩に事前に相談するとよいでしょう。

アンケートにもご協力お願いいたします。

back

脚注

バーとかパイプとか呼ばれていますが、Rのパイプは別のものを示す場合が多いので、バーとしましょう。 ↩
これは条件付確率 $P(X|Z)$ が「$Z$という文脈を与えた時の$X$の確率」を示すことと同じく、Days を被験者の文脈で条件付けていることと同じです。 ↩
lmerTestパッケージを使えばp値も出せます。 ↩
もちろん Gamma分布や対数正規分布過程してもいいが、その場合は glmer を使うことになるのと、 Gamma分布の際は解釈が難しくなる。対数正規分布は試したが(glmerでfamily=gaussian(link=log))、どうも切片のランダム効果が小さくなりすぎて警告を受けた。 ↩